采用hilbert变换法寻找R波

yangshangbin

根据版主的提示，对心电信号的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;
}

风子师兄

相当实用和完整的帖子{:soso_e179:}