脉搏电压测量程序改进版--实时反馈

cal-heart.c

/* [2013-07-06 六 12:18]  IF.Q
 * 脉搏检测
 
 * */

#include <stdio.h>
#include <stdlib.h>

typedef int beatype;

/* 临时这个么搞，读数据的缓冲区大小 */
#define MAX_LEN  999999
/* 每次送入检测的数据量 10秒数据, main程序中的使用的 */
#define INDATA_LEN    (SCAN_FREQ * 3)

/* 每秒采样次数,根据实际系统的采样率调整 */
#define SCAN_FREQ  78
/* 快速射血期 的扫描次数 */
#define QUICK_BLOOD_INJECT  ((int)(SCAN_FREQ * 0.11))
/* 检测到一次之后的跳跃点数 */
#define JUMP_INTERVAL   ((int)(SCAN_FREQ * 0.382))
/* 过滤用的门限值 */
#define THRESH_RATE   (float)(1 - 0.718)


/* 读取文件的buff */
beatype alldata[MAX_LEN];

/* 算法全局变量 */
beatype g_dif_max=0;            /* 导数最大值 */
beatype g_thr_max;              /* 导数阀值 */

/* is_diff_peak
 *    beatype *diff : 导数缓冲区数组
 *    int pos       : 缓冲区的偏移量
 *    beatype thr_max： 过滤用的阀值
 *
 * 判断导数值波峰的函数，就是找尖尖。
 * 尖尖附近的几个点的平均值应该大于阀值。通过这个来找。
 * 因此在不同的场合可能需要调整 THRESH_RATE 的值。
 *              .                              .
 *             ...              .            ......
 *           ... ..            ..         ....    ..
 *          ..    ..          ....      ...        ...
 *        ...      ...      ..   ....  ..            ..   .....
 *  .......          .......        ....              .....   ....
 *
 * 传入的pos 之后还要有QUICK_BLOOD_INJECT / 2个数据空间，然后会越界
 * */
int inline is_diff_peak(beatype *diff, int pos, beatype thr_max)
{
    beatype tmp;
    beatype *iter;
    int i;

    if (*(diff + pos) < thr_max) {
        return 0;
    }

    iter = diff + pos - (int)(QUICK_BLOOD_INJECT / 2);
    tmp = 0;
    for (i = 0; i < QUICK_BLOOD_INJECT; ++i, ++iter) {
        tmp += *iter;
    }

    tmp /= QUICK_BLOOD_INJECT;
    if (tmp < thr_max) {
        return 0;
    }

    return 1;
}

/************************************************************
 * calculate_beat()   计算心跳的主函数
 *    beatype *fifo: 传感器电压raw数据的数组
 *    int len      : 数组长度
 * 可对任意长度的数组进行计算，返回数组区域内的心跳数
 *            ..          ....
 *          ..  ....... ...  ....
 *          .         ...       .     ....      ......
 *          ..                 ..   ..   .......      .
 *           ...             ...    .       ...     ...
 *              . .        ...      ..         ......
 *                ....  ....         ...........
 *                   ....
 *  注意： 每次调用会更新全局变量 g_thr_max, 如果换人或拆装传感器
 *         可能导致电压幅度变化，应该清空两个全局变量。
 *  ************************************************************/
int calculate_beat(beatype *fifo , int len)
{

    int i;

    beatype *diff;               /* 导数缓冲区 */
    int  beat = 0;               /* 计算出的心跳数 */

    beatype dif_max = g_dif_max; /* 导数最大值 */
    beatype thr_max = g_thr_max; /* 导数阀值 */

    //printf ("len = %d\n", len);
    diff = malloc(sizeof(beatype) * (len+2));
    if (!diff) {
        perror("malloc fail\n");
        return (-1);
    }

    *diff=0;                    /* 头两个字节置0 */
    *(diff+1)=0;

    /* 按DIFF(i)=f(i+1)-f(i-1)+2*f(i+2)-2*f(i-2)公式计算导数 */
    for (i=2;i<len-2;i++)
    {
        *(diff+i) = (*(fifo+i+1)) + 2 * (*(fifo+i+2))
            - (*(fifo+i-1)) - 2 * (*(fifo+i-2));

        if ((dif_max) < *(diff+i)) {
            dif_max= *(diff+i);
            thr_max= dif_max * THRESH_RATE;       //设定阈值
        }
        //printf ("%d\n", diff[i]);
    }

    *(diff+len-1)=0;            /* 置0 */
    *(diff+len)=0;

    //printf ("max dif=%d, thr=%d\n", dif_max, thr_max);

    /* 检测 */
    for (i = QUICK_BLOOD_INJECT / 2; i < len - 2; i++)
    {
        if (is_diff_peak(diff, i, thr_max)) {
            beat++;
            i += JUMP_INTERVAL;
        }
    }

    free(diff);

    g_dif_max = dif_max;            /* 更新全局值 */
    g_thr_max = thr_max;

    return beat;
}


int main(int argc, char **argv)
{
    beatype tmp;
    int reallen = 0;
    int totlebeat;
    float time;
    int i;

    if (argc < 3) {
        printf ("argc not enough!\n");
        exit(-1);
    }

    time = atof(argv[2]);

    FILE *fp = fopen(argv[1],"r");

    /* 根据数据的类型input */
    while (fscanf(fp, "%d", &tmp) != EOF) {
        alldata[reallen++] = tmp > 255 ? 255 : tmp;
    }

    fclose(fp);
#if 0
    for (i = 0; i < reallen; ++i) {
        printf ("%f\n", alldata[i]);
    }
#endif

#if 1
    totlebeat = 0;
    /* 测量过程中，每次凑够 INDATA_LEN 的数据后就调用函数计算短时心率 */
    for (i = reallen; i > 0; i -= INDATA_LEN) {
        totlebeat += calculate_beat(alldata + reallen - i,
                                 i > INDATA_LEN ? INDATA_LEN : i);
        /* 返回值怎么用，根据程序情况自定 */
        printf ("--- allbeat = %d, rate = %f\n", totlebeat,
                totlebeat/(((float)(reallen-i)/reallen) * time));
    }
#endif

    /* 测量结束，从头到尾算一次 */
    totlebeat = calculate_beat(alldata, reallen);
#ifdef DBG_SW
    printf("dbg diff=%d, range=%d\n",cnt_diff_under, cnt_range_under);
#endif

    printf("total beat = %d, rate = %f\n", totlebeat, totlebeat/time);

    return 0;
}