52matlab技术网站,matlab教程,matlab安装教程,matlab下载

 找回密码
 立即注册
搜索
热搜: 活动 交友 discuz
查看: 10580|回复: 1
打印 上一主题 下一主题

采用hilbert变换法寻找R波

[复制链接]

6

主题

8

帖子

63

积分

版主

Rank: 7Rank: 7Rank: 7

积分
63
跳转到指定楼层
楼主
发表于 2015-5-19 17:12:52 | 只看该作者 回帖奖励 |倒序浏览 |阅读模式
根据版主的提示,对心电信号的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
回复

使用道具 举报

3

主题

7

帖子

53

积分

版主

Rank: 7Rank: 7Rank: 7

积分
53
沙发
发表于 2015-5-19 17:30:57 | 只看该作者
相当实用和完整的帖子{:soso_e179:}
回复 支持 反对

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册

本版积分规则

Archiver|手机版|小黑屋|52matlab技术网站 ( 粤ICP备14005920号-5 )

GMT+8, 2024-11-28 06:23 , Processed in 0.081635 second(s), 19 queries .

Powered by Discuz! X3.2 Licensed

© 2001-2013 Comsenz Inc.

快速回复 返回顶部 返回列表