|
根据版主的提示,对心电信号的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[0] = in;
din[1] = 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[0] = 1;
h[len/2] = 1;
for (i = 1; i <= len/2-1; i++)
{
h = 2;
}
}
else
{
h[0] = 1;
for (i = 1; i < (len + 1) / 2 ; i++)
{
h = 2;
}
}
for (i = 0; i < len; i++)
{
dout[0] *= h;
dout[1] *= 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[0] / len, din[1] / len);
}
for (i = 0; i < len; i++)
{
out = sqrt((din[0] * din[0] + din[1] * din[1]) / (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[10] = {0,1,2,3,4,5,6,7,8,9};
double dst[10];
int i;
// test();
hilbert(buf, 10, dst);
for (i = 0; i < 10; i++)
{
printf("%f\r\n", dst);
}
printf("\r\n");
getchar();
return 0;
}
|
本帖子中包含更多资源
您需要 登录 才可以下载或查看,没有帐号?立即注册
x
|