emd算法实现bcg信号心跳呼吸率提取

【实例简介】emd算法实现bcg信号心跳呼吸率提取

【实例截图】

【核心代码】

int imf(double x[N],double y[N],int b[N],double h[N]) //  输入：x，y，b； 输出：是否继续的条件；产生：imf分量；

{
// int go_on;  //继续分解成imf的标志位；
int i,j;
int n = 0;
int flag;
int count = 0;
double * temp_x_max;  
double * temp_y_max;
double * temp_x_min;
double * temp_y_min;
coefficient *coe_max;
coefficient *coe_min;
point p_max;
point p_min;
double m[N];

for(i=0; i<N; i )
{
h[i] = y[i];
}
/***************
**因为是采集一定的点之后再做EMD，而且每次做完EMD之后就要进行FFT来确定哪些分量是跟呼吸或者心跳有用的频率分量，所以采样的点数和频率要和FFT
**的相同。所以要采集512个点。但是此程序是测试程序，先将点的数目减少。
**第一步：模拟产生原始信号的横纵坐标；
**因为一次迭代不一定就求出了最高频，所以要反复迭代直到最后的平均包络为零；
****************/
do
{
flag = 0; //设置标志位，控制循环的条件；

// printf("\n第%d次结果如下 :\n",count);

//第二步：找极大值点，对极大值点进行包络，求出每两个极大值区间的包络函数；
n = extr_max(h,b);   //极大值点的个数；
if(n < 2)
return 1; //判断极值点是否大于等于二,如果不满足条件，就直接退出；

b[0] = 0;
b[n 1] = N-1; 

n = n 2; //这句话的意思是，要求整个函数的包络，要加上开头和结尾的点，故比极大值多两个点。
temp_x_max = (double *)malloc(n*sizeof(double));
temp_y_max = (double *)malloc(n*sizeof(double));
for(i=0; i<n; i )
{
temp_x_max[i] = x[b[i]];
temp_y_max[i] = h[b[i]];
}
// for(i=0; i<n; i )
// {
// printf("%f\t%f\n",temp_x_max[i],temp_y_max[i]);
// }
coe_max = (coefficient *)malloc((n - 1) * sizeof(coefficient));
p_max.xCoordinate = temp_x_max;
p_max.yCoordinate = temp_y_max;
p_max.f0 = 0;
p_max.fn = 0;
p_max.num = n;
p_max.con = NOTaKNOT;
p_max.coe = coe_max;
spline(&p_max);

//输出经过三次样条差值之后的每个区间的系数；
// for(i = 0; i < n - 1; i )
// printf("%f\t%f\t%f\t%f\n", coe_max[i].a3, coe_max[i].b2, coe_max[i].c1, coe_max[i].d0);

//求出极大值包络线的二分之一；
j = n - 2;
for(i=N-1; i>=0; i--)
{
if(x[i] >= x[b[j]])
{
m[i] = func(coe_max[j].a3,coe_max[j].b2,coe_max[j].c1,coe_max[j].d0,x[i],x[b[j]])/2;
}
else
{
j--;
m[i] = func(coe_max[j].a3,coe_max[j].b2,coe_max[j].c1,coe_max[j].d0,x[i],x[b[j]])/2;
}
}
/* printf("%极大值的平均包络值\n*****************************************\n"); //         
for(i=0; i<N; i )
{
printf("%f\n",m[i]); //——————————————————————————————每次迭代都输出最后的平均包络值；

}
*/
free(temp_x_max);
free(temp_y_max);
free(coe_max);

//第五步：找极小值点，对极小值点进行包络，求出每两个极小值区间的包络函数；
n = extr_min(h,b);   //极小值点的个数；

if(n < 2)
return 1; //判断极值点是否大于等于二,如果不满足条件，就直接退出；

b[0] = 0;
b[n 1] = N-1;
n = n 2; //这句话的意思是，要求整个函数的包络，要加上开头和结尾的点，故比极小值多两个点。
temp_x_min = (double *)malloc(n*sizeof(double));
temp_y_min = (double *)malloc(n*sizeof(double));
for(i=0; i<n; i )
{
temp_x_min[i] = x[b[i]];
temp_y_min[i] = h[b[i]];
}
/*
for(i=0; i<n; i )
{
printf("%f\t%f\n",temp_x_min[i],temp_y_min[i]);
}
*/

coe_min = (coefficient *)malloc((n - 1) * sizeof(coefficient));
p_min.xCoordinate = temp_x_min;
p_min.yCoordinate = temp_y_min;
p_min.f0 = 0;
p_min.fn = 0;
p_min.num = n;
p_min.con = NOTaKNOT;
p_min.coe = coe_min;
spline(&p_min);


// for(i = 0; i < n - 1; i )
// printf("%f\t%f\t%f\t%f\n", coe_min[i].a3, coe_min[i].b2, coe_min[i].c1, coe_min[i].d0);


//求出整体的二分之一，即包括极大值点。
j = n - 2; //j重新置数 这点很重要的！！
for(i=N-1; i>=0; i--)
{
if(x[i] >= x[b[j]])
m[i] = func(coe_min[j].a3,coe_min[j].b2,coe_min[j].c1,coe_min[j].d0,x[i],x[b[j]])/2;
else
{
j--;
m[i] = func(coe_min[j].a3,coe_min[j].b2,coe_min[j].c1,coe_min[j].d0,x[i],x[b[j]])/2;
}
}

free(temp_x_min);
free(temp_y_min);
free(coe_min);

// printf("%平均包络值\n*****************************************\n"); // 
for(i=0; i<N; i )
{
// printf("%f\n",m[i]); //——————————————————————————————每次迭代都输出最后的平均包络值；
if(m[i] >= 0.01) // 不知道判定条件是不是过于严谨，有说用系数的，看情况；
flag = 1;
}

//用原始的信号减去新产生的包络平均信号；
for(i = 0; i< N; i )
{
h[i] = h[i] - m[i]; 
}
//设置标志位，判断平均包络信号是否与x轴重合；
/* printf("%原始信号减去平均包络值\n*****************************************\n"); // 
for(i=0; i<N; i )
{
printf("%f\n",h[i]); //——————————————————————————————每次迭代都输出最后的平均包络值；
}
printf("flag = %d\n",flag);
*/
count ;
// if(count == 50)
// flag = 0;
}while(flag);  // 
// printf("count = %d\n ",count);//_________________________________________每迭代一次换行；
// getchar();
return 0;
}

/**********************************************EMD*******************************************/

void emd(double *x, double *y, double fc)//fc为对给进的初始信号进行fft之后求得的频谱最大值对应的频率
{
int go_on;
double **im; 
int count = -1;
int count1 = 0;
int count2 = 0;
int i = 0;
int j = 0;
double total = 0;
double e_asp = 0;
double e_hea = 0;
// double asp[N],hea[N]; //  存储的呼吸和心跳重构信号；
double h[N];  //存储y的值，保证y的值不被破坏；
int b[N];
double general[N];
comp data[N]; //申请FFT变换的静态内存；
FILE * fp;
for(i=0; i<N; i )
{
h[i] = y[i];
}

do   //找到需要分解成几个IMF分量，以确定要申请的内存空间；
{
go_on = imf(x,h,b,general);
for(i=0; i<N; i )
{
h[i] = h[i] - general[i];
}
count ;
}while(!go_on);

// printf("###################\nThe number of imf is %d\n",count);
// 动态申请二维数组，存储IMF分量；
im = (double **)malloc(count * sizeof(double*));

for(i=0; i<count; i )
{
im[i] = (double *)malloc(N * sizeof(double));
}
//再来一次，所以h要重新赋值；
for(i=0; i<N; i )
{
h[i] = y[i];
}
//存储IMF分量值；
for(i=0; i<count; i )
{
go_on = imf(x,h,b,im[i]);
for(j=0; j<N; j )
{
h[j] = h[j] - im[i][j];
}
}