C语言处理rppg信号，计算心率

主要流程为：

1.对rppg信号进行滤波，用FIR滤波器进行滤波，对应函数为：bpfilter

2.对滤波后信号利用fftw库进行傅里叶变换并求功率谱密度

3.利用功率谱密度估计心率，呼吸率并且利用峰值检测法计算HRV（当前算法有问题，不准）

后续发现问题会继续更新

文件路径示意如下：

---rPPGc

        ----fftw

        ----tool

                ----physiology_postprocess.c

                ----physiology_postprocess.h

        ----main.c

        ----CMakeLists.txt

        

main.c

可以按照文件夹读取或者直接读取文件

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <fftw3.h>
#include <physiology_postprocess.h>
#include <dirent.h>

void calculateByFile(const char* filename,float* ans);
char** readLinesFromFile(const char* filename, int* numLines);
float CalculateMean(const float* array, int size);
float CalculateStandardDeviation(const float* array, int size);
void calculateByDir(const char* dirname);
void calculateWhole(float* finput,int numLines,float fs,int segment,float* ans);
void calculateSegment(float* finput,int numLines,float fs,int segment,float* ans);
void calculateAsPython(float* finput,int numLines,float fs,int segment,float* ans);

int main() {
    char* folderPath = "rppg_pure";
    calculateByDir(folderPath);
    
    // char* filepath = "0103ecg.txt";
    // float* ans = (float*)malloc(3*sizeof(float));
    // calculateByFile(filepath,ans);
    // printf("%s\n",filepath);
    return 0;
    
}

void calculateByDir(const char* dirname){
    char* folderPath = "/mnt/sn750_1/guoyunsheng/rppg/iai-PhysNet-rPPG/rppg_pure";
    // 打开目标文件夹
    DIR *dir;
    struct dirent *entry;
    // 打开文件夹
    dir = opendir(folderPath);

    if (dir == NULL) {
        perror("无法打开文件夹");
        return 1;
    }
    // 读取文件夹中的文件名称并组成绝对路径
    while ((entry = readdir(dir))) {
        if (entry->d_type == DT_REG) { // 仅处理普通文件
            char filepath[1000];
            snprintf(filepath, sizeof(filepath), "%s/%s", folderPath, entry->d_name);
            float* ans = (float*)malloc(3*sizeof(float));
            calculateByFile(filepath,ans);
            // printf("%s|%f|%f\n",entry->d_name,ans[0],ans[1]);
            // printf("%s\n",entry->d_name);
        }
    }
    // 关闭文件夹
    closedir(dir);
}

void calculateWhole(float* finput,int numLines,float fs,int segment,float* ans){
    // printf("Calculate By whole\n");
    //全体处理之后输入到算法中计算
    float* longinput = (float*)malloc(numLines*sizeof(float));
    float meanLong = CalculateMean(finput,numLines);
    float stdlong = CalculateStandardDeviation(finput,numLines);
    float minlong = 100,maxlong = -100;
    //归一化处理
    for(int i = 0;i < numLines;i++){
        longinput[i] = (finput[i] - meanLong)/stdlong;
        if(longinput[i] > maxlong) maxlong = longinput[i];
        if(longinput[i] < minlong) minlong = longinput[i];
    }
    //归一化处理
    for(int i = 0;i < numLines;i++){
        longinput[i] = (longinput[i] - minlong)/(maxlong-minlong);
    }
    //将计算的心率，呼吸率，心率变异性通过ans返回
    CalculatePhysiologyRate(longinput,numLines,segment,fs,ans);
    printf("HR:%f,HRV:%f,BR:%f\n",ans[0],ans[2],ans[1]);
    free(longinput);
}

void calculateSegment(float* finput,int numLines,float fs,int segment,float* ans){
    //input划分单个segment进行计算
    // printf("Calculate By segment\n");
    float* seginput = (float*)malloc(segment*sizeof(float));
    int segnum = numLines/segment;
    for(int i = 0;i < segnum;i++){
        for(int j = 0;j < segment;j++){
            seginput[j] = finput[j+i*segment];
        }
        float meanSeg = CalculateMean(seginput,segment);
        float stdSeg = CalculateStandardDeviation(seginput,segment);
        float segmin = 100.0;
        float segmax = -100.0;
        //归一化处理
        for(int j = 0;j < segment;j++){
            seginput[j] = (seginput[j]-meanSeg)/stdSeg;
            if(seginput[j] - segmax > 0){
                segmax = seginput[j];
            } 
            if(seginput[j] < segmin){
                segmin = seginput[j];
            } 
        }
        //printf("segmax:%f,segmin:%f\n",segmax,segmin);
        for(int j = 0;j < segment;j++){
            seginput[j] = (seginput[j]-segmin)/(segmax-segmin);
        }
        int window_size = 100;
        CalculatePhysiologyRate(seginput,segment,window_size,fs,ans);
        // printf("HR:%f\nHRV:%f\nBR:%f\n",ans[0],ans[2],ans[1]);
        printf("HR:%f,BR:%f\n",ans[0],ans[1]);
    }
    free(seginput);
}

void calculateAsPython(float* finput,int numLines,float fs,int segment,float* ans){
    //input划分单个segment进行计算
    // printf("Calculate By segment\n");
    int window_size = 100;
    float* seginput = (float*)malloc(segment*sizeof(float));
    float* window_input = (float*)malloc(window_size*sizeof(float));
    int segnum = numLines/segment;
    int window_len = numLines/window_size;
    int window_per_seg = segment / window_len;
    for(int i = 0;i < segnum;i++){
        int startindex = i * segment;
        float segmin = 100.0;
        float segmax = -100.0;
        for(int j = 0;j < window_per_seg;j++){
            //第j个窗口
            for(int k = 0;k < window_len;k++){
                window_input[k] = finput[startindex+j*window_len+k];
            }
            float mean_window = CalculateMean(window_input,window_len);
            float std_window = CalculateStandardDeviation(window_input,window_len);
            //均值方差归一化
            for(int k = 0;k < window_len;k++){
                window_input[k] = (window_input[k]-mean_window)/std_window;
                //归一化后的结果放到seginput中
                seginput[j*window_len+k] = window_input[k];
                if(window_input[k] - segmax > 0){
                    segmax = window_input[k];
                } 
                if(window_input[k] < segmin){
                    segmin = window_input[k];
                }
            }
        }
        //将所有的window组合起来形成一个segment计算最大最小值归一化
        for(int j = 0;j < segment;j++){
            seginput[j] = (seginput[j]-segmin)/(segmax-segmin);
        }
        CalculatePhysiologyRate(seginput,segment,window_size,fs,ans);
        // printf("HR:%f\nHRV:%f\nBR:%f\n",ans[0],ans[2],ans[1]);
        printf("%f|%f\n",ans[0],ans[1]);
    }
    free(seginput);
    free(window_input);
}

void calculateByFile(const char* filename,float* ans){
    int numLines;
    char** input = readLinesFromFile(filename, &numLines);
    if(numLines == 0){
        printf("file empty");
        return;
    }
    //printf("NumLines:%d\n",numLines);
    if (input == NULL) {
        printf("line null");
        return;
    }
    //转换成浮点型数据
    float* finput = (float*)malloc(numLines*sizeof(float));
    for(int i = 0;i < numLines;i++){
        float value;
        if (sscanf(input[i], "%f", &value) == 1) {
            finput[i] = value;
        } else {
            // 处理解析失败的情况，例如可以赋予默认值
            finput[i] = 0.0;
        }
    }
    float fs = 30;
    //分段计算心率
    int segment = 1000;
    // calculateWhole(finput,numLines,fs,segment,ans);
    // calculateSegment(finput,numLines,fs,segment,ans);
    calculateAsPython(finput,numLines,fs,segment,ans);
    free(finput);
    free(input);
}

//计算均值
float CalculateMean(const float* array, int size) {
    float sum = 0.0;
    for (int i = 0; i < size; i++) {
        sum += array[i];
    }
    return sum / size;
}
//计算方差
float CalculateStandardDeviation(const float* array, int size) {
    float mean = CalculateMean(array, size);
    float sumOfSquares = 0.0;

    for (int i = 0; i < size; i++) {
        float diff = array[i] - mean;
        sumOfSquares += diff * diff;
    }

    float variance = sumOfSquares / size;
    float stdDeviation = sqrt(variance);

    return stdDeviation;
}


//读取文件中的内容
/***
 * 将filename中的数据写入到numLines中
 * 参数说明
 * filename：读取的文件名称
 * numLine：读取的行数
**/
char** readLinesFromFile(const char* filename, int* numLines) {
    FILE* file = fopen(filename, "r");
    if (file == NULL) {
        printf("无法打开文件\n");
        return NULL;
    }
    
    char** lines = NULL;
    int capacity = 0;
    int size = 0;
    
    char buffer[256];
    while (fgets(buffer, sizeof(buffer), file) != NULL) {
        char* line = malloc(strlen(buffer) + 1);
        strcpy(line, buffer);
        
        if (size >= capacity) {
            capacity += 10;
            lines = realloc(lines, capacity * sizeof(char*));
        }
        lines[size] = line;
        size++;
    }
    
    fclose(file);
    
    *numLines = size;
    return lines;
}

calclulate_HR接受一个文件名作为输入，按行读取该文件中的rppg信号值，返回计算好的心率。

physiology_postprocess.h头文件如下：

#ifndef PHYSIOLOGYPOSTPROCESS
#define PHYSIOLOGYPOSTPROCESS

/***
 * 功能：将input划分为多个segment大小的段分别计算心率和呼吸率，根据多段的心率计算心率变异性并返回
 * input:输入原始的rppg信号
 * numLine：rppg信号的长度
 * segment:需要分割的每个段的大小
 * fs:采样频率
 * ans:用来返回计算的结果,ans[0]：心率,ans[1]呼吸率，ans[2]心率变异性
*/
void CalculatePhysiologyRate(const float* input,const int numLines,const int segment,const float fs,float* ans);

#endif

tool.c文件如下：

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

//写入文件
void writeFile(const char* filename,float* numbers,int length){

    FILE* file = fopen(filename, "w");
    if (file == NULL) {
        printf("无法打开文件。\n");
    }

    for (int i = 0; i < length; i++) {
        float x = numbers[i];
        char a[50];
        sprintf(a,"%f\n", x);
        fprintf(file, "%s", a);
    }

    fclose(file);
    //printf("write to %s already done.\n",filename);

}

static float bessel0(float x)
{
    int i;
    float d;
    float y;
    float d2;
    float sum;
    y = x / 2.0;
    d = 1.0;
    sum = 1.0;
    for (i=1; i<=25; i++)
    {
        d = d * y / i;
        d2 = d * d;
        sum = sum + d2;
        if (d2 < sum*(1.0e-8))
        {
            break;
        }
    }
    return(sum);
}

static float kaiser(int i, int n, float beta)
{
    float a;
    float w;
    float a2;
    float b1;
    float b2;
    float beta1;

    b1 = bessel0(beta);
    a = 2.0 * i / (float)(n - 1) - 1.0;
    a2 = a * a;
    beta1 = beta * sqrt(1.0 - a2);
    b2 = bessel0(beta1);
    w = b2 / b1;
    return(w);
}

//n：窗口长度 type：选择窗函数的类型 beta：生成凯塞窗的系数
static float window(int type, int n, int i, float beta)
{
    int k;
    float pi;
    float w;
    pi = 4.0 * atan(1.0);  //pi=PI;
    w = 1.0;

    switch (type)
    {
    case 1:
    {
        w = 1.0;  //矩形窗
        break;
    }
    case 2:
    {
        k = (n - 2) / 10;
        if (i <= k)
        {
            w = 0.5 * (1.0 - cos(i * pi / (k + 1)));  //图基窗
        }
        if (i > n-k-2)
        {
            w = 0.5 * (1.0 - cos((n - i - 1) * pi / (k + 1)));
        }
        break;
    }
    case 3:
    {
        w = 1.0 - fabs(1.0 - 2 * i / (n - 1.0));//三角窗
        break;
    }
    case 4:
    {
        w = 0.5 * (1.0 - cos( 2 * i * pi / (n - 1)));//汉宁窗
        break;
    }
    case 5:
    {
        w = 0.54 - 0.46 * cos(2 * i * pi / (n - 1));//海明窗
        break;
    }
    case 6:
    {
        w = 0.42 - 0.5 * cos(2 * i * pi / (n - 1)) + 0.08 * cos(4 * i * pi / (n - 1));//布莱克曼窗
        break;
    }
    case 7:
    {
        w = kaiser(i, n, beta);//凯塞窗
        break;
    }
    }
    return(w);
}

/***
 *  参数说明
    n：滤波器的阶数
    band：滤波器的类型,取值1-4,分别为低通、带通、带阻、高通滤波器
    wn：窗函数的类型，取值1-7，分别对应矩形窗、图基窗、三角窗、汉宁窗、海明窗、布拉克曼窗和凯塞窗
    fs：采样频率
    h：存放滤波器的系数
    kaiser：beta值
    fln：带通下边界频率
    fhn：带通上边界频率 
***/

void firwin(int n, int band, int wn, float fs, float *h, float kaiser, float fln, float fhn)
{
    int i;
    int n2;
    int mid;
    float s;
    float pi;
    float wc1;
    float wc2;
    float beta;
    float delay;
    beta = kaiser;
    pi = 4.0 * atan(1.0);   //pi=PI;
    
    if ((n % 2) == 0)/*如果阶数n是偶数*/
    {
        n2 = (n / 2) - 1;
        mid = 1;
    }
    else
    {
        n2 = n / 2;//n是奇数,则窗口长度为偶数
        mid = 0;
    }
    
    delay = n / 2.0;
    wc1 = 2 * pi * fln;
    wc2 = 2 * pi * fhn;
    
    switch (band)
    {
    case 1:
    {
        for (i=0; i<=n2; ++i)
        {
            s = i - delay;
            h[i] = (sin(wc1 * s / fs) / (pi * s)) * window(wn, n+1, i, beta);//低通,窗口长度=阶数+1，故为n+1
            h[n - i] = h[i];
        }
        if (mid == 1)
        {
            h[n / 2] = wc1 / pi;//n为偶数时，修正中间值系数
        }
        break;
    }
    case 2:
    {
        for (i=0; i<=n2; i++)
        {
            s = i - delay;
            h[i] = (sin(wc2 * s / fs) - sin(wc1 * s / fs)) / (pi * s);//带通
            h[i] = h[i] * window(wn, n+1, i, beta);
            h[n-i] = h[i];
        }
        if (mid == 1)
        {
            h[n / 2] = (wc2 - wc1) / pi;
        }
        break;
    }
    case 3:
    {
        for (i=0; i<=n2; ++i)
        {
            s = i - delay;
            h[i] = (sin(wc1 * s / fs) + sin(pi * s) - sin(wc2 * s / fs)) / (pi * s);//带阻
            h[i] = h[i] * window(wn, n+1, i, beta);
            h[n - i] = h[i];
        }
        if (mid == 1)
        {
            h[n / 2] = (wc1 + pi - wc2) / pi;
        }
        break;
    }
    case 4:
    {
        for (i=0; i<=n2; i++)
        {
            s = i - delay;
            h[i] = (sin(pi * s) - sin(wc1 * s / fs)) / (pi * s);//高通
            h[i] = h[i] * window(wn, n+1, i, beta);
            h[n-i] = h[i];
        }
        if (mid == 1)
        {
            h[n / 2] = 1.0 - wc1 / pi;
        }
        break;
    }
    }
}



//bp滤波器
void BpFilter(const float* input, float* filt,const int num,const float fs){
    //滤波器的开始频率
    float minfq = (0.8 * 2) / fs;
    //滤波器的截止频率
    float maxfq = (3 * 2) / fs;
    //n为阶数
    int n = num / 10;
    //动态分配内存,h为最后生成的滤波器系数
    float* h = (float *)malloc(sizeof(float) * (n+1));
    
    if(h == NULL){
        perror("bpfilter realloc memory error");
    }
    firwin(n, 2, 5, fs, h, 0,  minfq,  maxfq);
    //FIR滤波器处理过程
    for(int i = 0; i < num; i++){
        filt[0] = 0;
        for(int k = 0;k < n+1; k++){
            if(i - k >= 0)
                filt[i] += input[i-k]*h[k];
        }
    }
    // printf("-----------parameter-----------------\n");
    // for(int i = 0;i < n+1;i++){
    //     printf("%f\n",h[i]);
    // }
    free(h);
}

//计算功率谱密度
void CalSpectrum(float *filt,float *powerSpectrum,int numLines, int fs,int N){
    //filt = readLinesFromFile("s2f.txt", &numLines);
    //printf("-----------fft start-----------------\n");
    N = numLines*10;
    //傅里叶变换
    fftw_complex *in,*out;
    out = (fftw_complex*) fftw_malloc(sizeof(fftw_complex)*N);
    in = (fftw_complex*) fftw_malloc(sizeof(fftw_complex)*N);
    fftw_plan my_plan;
    //确认傅里叶变换输入
    for (int i = 0;i < numLines;i++){
        in[i][0] = filt[i];
        in[i][1] = 0.0;
        //printf("%f\n",filt[i]);
    }
    //补充不足的部分为0
    for(int i = numLines;i < N;i++){
        in[i][0] = 0.0;
        in[i][1] = 0.0;
    }
    my_plan = fftw_plan_dft_1d(N, in, out, FFTW_FORWARD, FFTW_ESTIMATE);
    fftw_execute(my_plan);
    //printf("-----------fft end-----------------\n");
    //printf("-----------calculate spectrum start-----------------\n");
    //傅里叶变换后，取幅值平方除以N，以此作为对真实功率谱的估计
    //计算功率谱密度
    
    float fremax = 0.0;
    for (int i = 0; i < N; i++) {
        float real = out[i][0];
        float imag = out[i][1];
        //printf("real:%f,imag:%f",real,imag);
        //printf("%d\n",i);
        powerSpectrum[i] = (real * real + imag * imag) / (N*fs);
        //printf("%f\n",powerSpectrum[i]);
        if(fremax < powerSpectrum[i]){
            fremax = powerSpectrum[i];
        }
        //printf("%f\n",powerSpectrum[i]);
    }
    //printf("-----------calculate spectrum end-----------------\n");
    //释放资源
    fftw_destroy_plan(my_plan);
    fftw_free(in); 
    fftw_free(out);
}
/***
 * welch法计算功率谱密度
 * signal为信号序列
 * signal_length为需要计算的信号总长度
 * segment_length为子段长度
 * len为信号原本长度
 * powerSpectrum保留记录出的计算功率谱密度
 * overlap重叠比例
*/
float* WelchPsd(float* signal, int signal_length, int segment_length,int len,int overlap) {
    int num_segments = signal_length / (segment_length * (1 - overlap));
    int i, j;
    // 分配FFTW输入和输出数组
    fftw_complex* in = (fftw_complex*)fftw_malloc(sizeof(fftw_complex) * segment_length);
    fftw_complex* out = (fftw_complex*)fftw_malloc(sizeof(fftw_complex) * segment_length);

    // 创建FFTW计划
    fftw_plan plan;

    // 初始化FFTW计划
    plan = fftw_plan_dft_1d(segment_length, in, out, FFTW_FORWARD, FFTW_ESTIMATE);
    //将signal补充到segment_length的长度
    float* temp = (float*)realloc(signal, signal_length * sizeof(float));
    //将powerSpectrum补充到足够的长度
    float* powerSpectrum = (float*)malloc((segment_length /2) * sizeof(float));

    for(int i = len;i < signal_length;i++){
        temp[i] = 0.0;
    }
    int k = 0;
    printf("num_segments:%d\n",num_segments);
    // 计算每个子段的功率谱密度
    for (i = 0; i < num_segments; i++) {
        int start = i * segment_length * (1 - overlap);

        // 将信号加载到输入数组
        for (j = 0; j < segment_length; j++) {
            in[j][0] = temp[start + j];
            in[j][1] = 0.0;
        }

        // 执行FFT
        fftw_execute(plan);

        // 计算功率谱密度
        float* psd = (float*) malloc(sizeof(float) * (segment_length / 2));
        for (j = 0; j < segment_length / 2; j++) {
            float real = out[j][0];
            float imag = out[j][1];
            psd[j] = (real * real + imag * imag) / segment_length;
            //记录在powerSpectrum中
            powerSpectrum[k] = psd[j];
            k++;
        }
        powerSpectrum = realloc(powerSpectrum, ((i+2) * segment_length / 2) * sizeof(float));
        // 打印功率谱密度
        printf("Segment %d PSD:\n", i + 1);
        for (j = 0; j < segment_length / 2; j++) {
            printf("%f\n", psd[j]);
        }
        printf("\n");
    }
    printf("powerSpectrum length:%d\n",k);

    // 释放内存和FFTW计划
    fftw_free(in);
    fftw_free(out);
    fftw_destroy_plan(plan);
    return powerSpectrum;
}

float CalculateMean(const float* array, int size) {
    float sum = 0.0;
    for (int i = 0; i < size; i++) {
        sum += array[i];
    }
    return sum / size;
}

float CalculateStandardDeviation(const float* array, int size) {
    float mean = CalculateMean(array, size);
    float sumOfSquares = 0.0;

    for (int i = 0; i < size; i++) {
        float diff = array[i] - mean;
        sumOfSquares += diff * diff;
    }

    float variance = sumOfSquares / size;
    float stdDeviation = sqrt(variance);

    return stdDeviation;
}

float CalculateMeanSquaredDeviation(const float* array, int size) {
    float mean = CalculateMean(array, size);
    float sumOfSquares = 0.0;

    for (int i = 0; i < size; i++) {
        float diff = array[i] - mean;
        sumOfSquares += diff * diff;
    }

    float msd = sumOfSquares / size;

    return msd;
}
/***
 * 将功率谱滤波计算心率
 * powerSpectrum:是功率谱密度
 * signal_length:功率谱密度的长度
 * fs:采样频率
 * lowfre-highfre:频率范围
*/
float FiltSpectrum(const float *powerSpectrum,float signal_length,int fs,float lowfre,float highfre){
    //有效功率谱范围
    int realfre_count = (highfre-lowfre)/(fs/signal_length)+3;
    //printf("realfre_count:%d\n",realfre_count);
    float *real_power = (float*)malloc(sizeof(float)*realfre_count);
    //初始化
    for(int i = 0;i < realfre_count;i++){
        real_power[i] = 0.0;
    }
    float *real_frequency = (float*)malloc(sizeof(float)*realfre_count);
    for(int i = 0;i < realfre_count;i++){
        real_frequency[i] = 0.0;
    }
    //计算频率轴的坐标
    float *frequecy = (float*)malloc(sizeof(float)*signal_length/2);
    int j = 0;
    //记录功率谱的最大值，以及最大值对应的index
    float powermax = 0.0;
    int maxindex = 0;
    for(int i = 0;i < signal_length/2;i++){
        frequecy[i] = i * fs / signal_length;
        if(frequecy[i] >= lowfre && frequecy[i] <= highfre){
            //printf("frequecy%d:%f\n",i,frequecy[i]);
            //对功率谱进行分析，首先排除不在指定范围内的数据，正常心率范围为lowfre-highfre,按照频率取功率谱
            //real_power[j] = powerSpectrum[i];
            if(powerSpectrum[i] > powermax){
                powermax = powerSpectrum[i];
                maxindex = j; 
            }
            real_frequency[j] = frequecy[i];
            //printf("real_frequency%d:%f\n",j,real_frequency[j]);
            j++;
        }
    }
    float hr = real_frequency[maxindex];
    //printf("max index:%d\nmax signal power:%f\nfrequency:%f\n",maxindex,real_power[maxindex],real_frequency[maxindex]);
    //printf("Heart rate:%f\n",hr);
    // char *frequencyfile = "frequency.txt";
    // writeFile(frequencyfile,frequecy,signal_length/2);
    // char *real_frequency_file = "real_frequency.txt";
    // writeFile(real_frequency_file,real_frequency,realfre_count);
    free(real_frequency);
    //free(real_power);
    return hr;
}
/***
 * input:输入的rppg信号
 * fs:采集的频率
 * numLines:rppg信号的长度
 * rate:记录心率和呼吸率
*/
void CalcHrBr(const float* input,const float fs,const float numLines,float* rate){
    //filt中存储经过滤波的数据
    float* filt = (float*)malloc(sizeof(float) * numLines);
    //初始化
    for(int i = 0;i < numLines;i++){
        filt[i] = 0.0;
    }
    BpFilter(input,filt,numLines, fs);
    //writeFile("filt.txt",filt,numLines);
    float signal_length = numLines*10;
    int N = numLines;
    float* powerSpectrum = (float*) malloc(sizeof(float) * numLines *10);
    for(int i = 0;i < numLines*10;i++){
        powerSpectrum[i] = 0.0;
    }
    CalSpectrum(filt, powerSpectrum, numLines,fs,&N);
    //计算心率
    float hr_lowfre = 1,hr_highfre = 2.5;
    float HR = FiltSpectrum(powerSpectrum,signal_length,fs,hr_lowfre,hr_highfre) * 60;
    //计算呼吸率
    float br_lowfre = 0.2,br_highfre = 0.3;
    float BR = FiltSpectrum(powerSpectrum,signal_length,fs,br_lowfre,br_highfre) * 60;
    // 释放内存时需要释放每行字符串的内存
    free(filt);
    free(powerSpectrum);
    rate[0] = HR;
    rate[1] = BR;
}

// 寻找数组最小值的下标
int argmin(int* index, int index_len)
{
	int min_index = 0;
	int min = index[0];
	for (int i = 1; i < index_len; i++)
	{
		if (index[i] < min)
		{
			min = index[i];
			min_index = i;
		}
	}
	return min_index;
}

//寻找极值点函数
// data是存放数据的数组
//index是存放峰值点下标的数组
//len_index是峰值个数，即index数组长度
void AMPD(float* data,int* index,int *len_index,int size)
{
	int* p_data = (int*)malloc(sizeof(int) * size); //size可以最大为数组长度
	int* arr_rowsum = (int*)malloc(sizeof(int) * size);
	int min_index, max_window_length;
	for (int i = 0; i < size; i++)
	{
		p_data[i] = 0;
	}
	for (int k = 1; k <= size / 2 + 1; k++)
	{
		int row_sum = 0;
		for (int i = k; i <= size - k; i++)
		{
			if ((data[i] > data[i - k]) && (data[i] > data[i + k]))
				row_sum -= 1;
		}
		*(arr_rowsum + k - 1) = row_sum;
	}
	/*for (int i = 0; i < size/2; i++)
	{
		printf("%d\n", arr_rowsum[i]);
	}*/
	min_index = argmin(arr_rowsum, size/2); //此处为最大的窗口
	//printf("%d\n", min_index);
	max_window_length = min_index;
	for (int k = 1; k < max_window_length + 1;k++)
	{
		for (int i = 1; i < size - k; i++)
		{
			if ((data[i] > data[i - k]) && (data[i] > data[i + k]))
				p_data[i] += 1;
		}
	}
	for (int i_find = 0; i_find < size; i_find++)
	{
		if (p_data[i_find] == max_window_length)
		{
			index[*len_index] = i_find;
			(*len_index) += 1;
		}
	}
	free(p_data);
	free(arr_rowsum);
}


/**
 * 函数返回HRV的值，同时将峰值的时间点记录在output指针中
 * 通过计算rppg信号中的峰值计算HRV
 * input为输入的rppg信号
 * numLines为rppg信号长度
 * fps为采样频率
*/
float FindPeaksCalcHRV(const float* input,int numLines, int fps) {
    int* result = (int*)malloc(numLines * sizeof(int));
    int resultCount = 0;
    float* filt = (float*)malloc(sizeof(float) * numLines);
    //初始化
    for(int i = 0;i < numLines;i++){
        filt[i] = 0.0;
    }
    BpFilter(input,filt,numLines, fps);
    //writeFile("hr62br14hrv54.txt",input,numLines);
    AMPD(filt,result,&resultCount,numLines);
    printf("resultCount:%d\n",resultCount);
    //计算每个波峰之间的差值
    float* de_peak = (float*)malloc(numLines * sizeof(float));
    for(int i = 0;i < resultCount-1;i++){
        de_peak[i] = (result[i+1] - result[i])/fps;
    }
    //计算差值的均值
    float peakmean = CalculateMean(de_peak, resultCount-1);
    float hrv = CalculateStandardDeviation(de_peak,resultCount-1);
    float hr = 60 / (peakmean);
    printf("Peak Calculate HR:%f\n",hr);
    printf("Peak Calculate HRV:%f\n",hrv);
    return hr;
}

/***
 * 功能：将input划分为多个segment大小的段分别计算心率和呼吸率，根据多段的心率计算心率变异性并返回
 * input:输入原始的rppg信号
 * numLine：rppg信号的长度
 * segment:需要分割的每个段的大小
 * fs:采样频率
 * ans:用来返回计算的结果,ans[0]：心率,ans[1]呼吸率，ans[2]心率变异性
*/
void CalculatePhysiologyRate(const float* input,const int numLines,const int segment,const float fs,float* ans){
    int numseg = numLines/segment;
    float* HRList = (float*) malloc(numseg*sizeof(float));
    float* BRList = (float*) malloc(numseg*sizeof(float));
    float* HR_BR = (float*) malloc(2*sizeof(float));
    float* newarray = (float*)malloc(segment * sizeof(float));
    for(int i = 0;i < numseg;i++){
        for(int j = 0;j < segment;j++){
            int index = i*segment + j;
            newarray[j] = input[index];
        }
        //对newarray中的数据求心率和呼吸率
        // if(i == 0){
        //     writeFile("orig_whole_0.txt",newarray,segment);
        //     CalcHrBr(newarray,fs,segment,HR_BR);
        // }
        CalcHrBr(newarray,fs,segment,HR_BR);
        //printf("%f\n",HR_BR[0]);
        HRList[i] = HR_BR[0];
        BRList[i] = HR_BR[1];
    }
    float meanHR = CalculateMean(HRList,numseg);
    float meanBR = CalculateMean(BRList,numseg);
    //float std = CalculateMeanSquaredDeviation(HRList,numseg);
    float std = CalculateStandardDeviation(HRList,numseg);
    ans[0] = meanHR;
    ans[1] = meanBR;
    ans[2] = std;
    // printf("---------------------------\n");
    // float HRV = FindPeaksCalcHRV(input,numLines,fs);
    // printf("---------------------------\n");
    free(newarray);
    free(HR_BR);
    free(HRList);
    free(BRList);
}

CMakeList如下：

cmake_minimum_required(VERSION 3.16)
project(rPPG)

set(COMPILE_CLANG OFF) 

if(COMPILE_FILE)
    message("generate by clang")
    set (CMAKE_C_COMPILER "/mnt/wdssd0/xuao/android-ndk-r25c/toolchains/llvm/prebuilt/linux-x86_64/bin/armv7a-linux-androideabi22-clang")
    set (CMAKE_CXX_COMPILER "/mnt/wdssd0/xuao/android-ndk-r25c/toolchains/llvm/prebuilt/linux-x86_64/bin/armv7a-linux-androideabi22-clang++")
endif()

SET(CMAKE_BUILD_TYPE "Debug")  # 定义编译类型
SET(CMAKE_CXX_FLAGS_DEBUG "$ENV{CXXFLAGS} -O0 -Wall -g2 -ggdb -fpermissive -pthread") # 定义Debug编译参数
SET(CMAKE_CXX_FLAGS_RELEASE "$ENV{CXXFLAGS} -O3 -Wall") # 定义Release编译参数

set(FFTW_PATH ${PROJECT_SOURCE_DIR}/fftw-3.3.10/build)
set(FFTWINCLUDE ${FFTW_PATH}/include)

# 这一块添加fftw的库
set(FFTW_LIB ${FFTW_PATH}/lib/libfftw3.a)

#添加静态库
add_library(PhysiologyPostprocess SHARED ${PROJECT_SOURCE_DIR}/tool/physiology_postprocess.c)

# 添加包含头文件的目录
target_include_directories(PhysiologyPostprocess PUBLIC ${PROJECT_SOURCE_DIR}/tool)

add_executable(rPPG rPPG.c)

target_include_directories(rPPG PRIVATE ${FFTWINCLUDE})
target_link_libraries(rPPG PhysiologyPostprocess m ${FFTW_LIB} )