采用hilbert变换法寻找R波
根据版主的提示，对心电信号的2阶差分再做一次hilbert变化后R波就变得很明显了，效果如图所示。查看matab的hilbert.m源码发现希尔伯特变化只需要通过FFT变化和IFF变化和一些简单的运算即可。因此解决此问题的首要问题就是编写一个混合基的FFT算法。目前比较有名气的开源FFT算法是FFTW（www.fftw.org)。当然KISS FFT也很不错，但是输入必须是实数。在此我们选择采用FFTW，FFTW安装方法请查阅附件。希尔伯特变化的代码如下。#include "stdafx.h"
#include <malloc.h>
#include "fftw3.h"
#include <stdio.h>
#include <string.h>
#include <math.h>
void hilbert(double *in, unsigned short len, double *out)
{
int i;
fftw_complex *din,*dout;
fftw_plan p;
double *h;
din= (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * len);
dout = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * len);
h = (double*)malloc(len * sizeof(double));
for (i = 0; i < len; i++)
{
din = in;
din = 0;
}
p = fftw_plan_dft_1d(len, din, dout, FFTW_FORWARD,FFTW_ESTIMATE);
fftw_execute(p); /* repeat as needed */
fftw_destroy_plan(p);
fftw_cleanup();
memset(h, 0, len * sizeof(double));
if (2*floor(len/2.0) == len)
{
h = 1;
h = 1;
for (i = 1; i <= len/2-1; i++)
{
h = 2;
}
}
else
{
h = 1;
for (i = 1; i < (len + 1) / 2 ; i++)
{
h = 2;
}
}
for (i = 0; i < len; i++)
{
dout *= h;
dout *= h;
}
p = fftw_plan_dft_1d(len, dout, din, FFTW_BACKWARD,FFTW_ESTIMATE);
fftw_execute(p); /* repeat as needed */
fftw_destroy_plan(p);
fftw_cleanup();
for(i=0;i<len;i++)/*OUTPUT*/
{
printf("%f,%fi\n",din / len, din / len);
}
for (i = 0; i < len; i++)
{
out = sqrt((din * din + din * din) / (len * len));
}
if (h != NULL) free(h);
if(din!=NULL) fftw_free(din);
if(dout!=NULL) fftw_free(dout);
}
int _tmain(int argc, _TCHAR* argv[])
{
double buf = {0,1,2,3,4,5,6,7,8,9};
double dst;
int i;
// test();
hilbert(buf, 10, dst);
for (i = 0; i < 10; i++)
{
printf("%f\r\n", dst);
}
printf("\r\n");
getchar();
return 0;
}
相当实用和完整的帖子{:soso_e179:}
页:
[1]