hilbert提取心音包络matlab程序

216 篇文章 52 订阅
35 篇文章 23 订阅
该博客详细介绍了心音信号的处理流程,包括重采样、低通滤波、经验模式分解(EMD)等步骤。通过MATLAB代码展示了如何从原始心音信号中提取固有模态函数(IMFs),并进行多次黄变换以分离信号成分。最终,进行了心音(S1和S2)的定位分析,为心率计算和心音特征提取奠定了基础。
摘要由CSDN通过智能技术生成

clear
close all;
%fid=fopen('out.txt','w');
%fclose(fid);
format long
[y,f,bit] = wavread('[Input the name of original signal here]');  % 其中y中是文件的数据,值在[-1,1],f表示的是采样频率,bit表示的是信号的位数
y=y(1:400000);
fs=f
figure(1)
subplot(3,1,1)
plot(y)
title('原心音信号');
n=length(y);               %获取数据的总长度
m=1;
% 开始重采样:
for i=1:15:n               %用for循环对数据重新采样,每15点取一个数据,这样的话频率就变为11025/15=735Hz,zhuanY()中的数据就是通过降低采样频率后的数据
    zhuanY(m)=y(i);
    m=m+1;
end
subplot(3,1,2)
plot(zhuanY)  %zhuanY中的数据就是经过采样后的数据。在保证不混叠的情况下,降低了数据量,有利于后续的处理
%开始低通滤波 
 
title('重采样后的信号');
[b,a]=cheby1(5,0.5,150/500);      %采用了cheby1型低通滤波器,截至频率设定在150Hz,经过观察后感觉视觉效果明显改善
x=filter(b,a,zhuanY);             %  数组x中就是通过低通滤波后的数据     x还是心音图       数组X中的数据是经过预处理后的数据
subplot(3,1,3)   
plot(x)                           %画出滤波过后的波形
axis([0 7000 -0.3 0.3]);

title('滤波后的信号');
Nxunhuan=0;
%% 下面开始经验模式分解  x是要处理的数据
[h,NmaxX,NminX,Nzero,mean_baoluo]=EMD(x);
C1=h;
while(~(((NmaxX+NminX==Nzero)|(abs(((NmaxX+NminX)-Nzero)-1)))&(max(abs(mean_baoluo))<0.01)))
    [h,NmaxX,NminX,Nzero,up_baoluo,mean_baoluo]=EMD(C1);
     C1=h;
    Nxunhuan=Nxunhuan+1;
end   
 figure(2)  
 subplot(7,2,1)
 plot(C1)                       %第一阶固有模态函数的波形图
 axis([0 4000 -0.3 0.3]);
title('第一阶固有模态函数的波形图');
LenC=length(C1);
%%下面开始第二次黄变换 即上面一步的操作得到了C1,下面一部要在C1的基础上取得C2
residue=x-h;
subplot(7,2,2)
plot(residue)             %残差函数的波形图
title('残差函数的波形图');
[h,NmaxX,NminX,Nzero,up_baoluo,mean_baoluo]=EMD(residue);
C2=h;
Nxunhuanr=0;                %用来记录循环的次数

while(~(((NmaxX+NminX==Nzero)|(abs(((NmaxX+NminX)-Nzero)-1)))&(max(abs(mean_baoluo))<0.01)))
    
 [h,NmaxX,NminX,Nzero,up_baoluo,mean_baoluo]=EMD(C2);
    C2=h;
    Nxunhuanr=Nxunhuanr+1;
end   
subplot(7,2,3)
plot(C2)                      %第二阶固有模态函数的波形图
axis([0 4000 -0.3 0.3]);
title('第二阶固有模态函数的波形图');

residue2=residue-C2;
subplot(7,2,4)
plot(residue2)            %第二次黄变换的残差
title('第二次黄变换的残差');
%%
[h,NmaxX,NminX,Nzero,up_baoluo,mean_baoluo]=EMD(residue2);
C3=h;
Nxunhuanr=0;                %用来记录循环的次数

while(~(((NmaxX+NminX==Nzero)|(abs(((NmaxX+NminX)-Nzero)-1)))&(max(abs(mean_baoluo))<0.01)))
    
 [h,NmaxX,NminX,Nzero,up_baoluo,mean_baoluo]=EMD(C3);
    C3=h;
    Nxunhuanr=Nxunhuanr+1;
end  
subplot(7,2,5)
plot(C3)                      %第三阶固有模态函数的波形图
axis([0 4000 -0.3 0.3]);
title('第三阶固有模态函数的波形图');

residue3=residue2-C3;
subplot(7,2,6)
plot(residue3)            %第三次黄变换的残差
title('第三次黄变换的残差');
%%
[h,NmaxX,NminX,Nzero,up_baoluo,mean_baoluo]=EMD(residue3);
C4=h;
Nxunhuanr=0;                %用来记录循环的次数

while(~(((NmaxX+NminX==Nzero)|(abs(((NmaxX+NminX)-Nzero)-1)))&(max(abs(mean_baoluo))<0.01)))
    
 [h,NmaxX,NminX,Nzero,up_baoluo,mean_baoluo]=EMD(C4);
    C4=h;
    Nxunhuanr=Nxunhuanr+1;
end  
subplot(7,2,7)
plot(C4)                      %第三阶固有模态函数的波形图
axis([0 4000 -0.3 0.3]);
title('第四阶固有模态函数的波形图');

residue4=residue3-C4;
subplot(7,2,8)
plot(residue4)            %第三次黄变换的残差
title('第四次黄变换的残差');
%%
[h,NmaxX,NminX,Nzero,up_baoluo,mean_baoluo]=EMD(residue4);
C5=h;
Nxunhuanr=0;                %用来记录循环的次数

while(~(((NmaxX+NminX==Nzero)|(abs(((NmaxX+NminX)-Nzero)-1)))&(max(abs(mean_baoluo))<0.01)))
    
 [h,NmaxX,NminX,Nzero,up_baoluo,mean_baoluo]=EMD(C5);
    C5=h;
    Nxunhuanr=Nxunhuanr+1;
end  
subplot(7,2,9)
plot(C5)                      %第三阶固有模态函数的波形图
axis([0 4000 -0.3 0.3]);
title('第五阶固有模态函数的波形图');

residue5=residue4-C5;
subplot(7,2,10)
plot(residue5)            %第三次黄变换的残差
title('第五次黄变换的残差');
%%
[h,NmaxX,NminX,Nzero,up_baoluo,mean_baoluo]=EMD(residue5);
C6=h;
Nxunhuanr=0;                %用来记录循环的次数

while(~(((NmaxX+NminX==Nzero)|(abs(((NmaxX+NminX)-Nzero)-1)))&(max(abs(mean_baoluo))<0.01)))
    
 [h,NmaxX,NminX,Nzero,up_baoluo,mean_baoluo]=EMD(C6);
    C6=h;
    Nxunhuanr=Nxunhuanr+1;
end  
subplot(7,2,11)
plot(C6)                      %第三阶固有模态函数的波形图
axis([0 4000 -0.3 0.3]);
title('第六阶固有模态函数的波形图');

residue6=residue5-C6;
subplot(7,2,12)
plot(residue6)            %第三次黄变换的残差
title('第六次黄变换的残差');
% %%
% [h,NmaxX,NminX,Nzero,up_baoluo,mean_baoluo]=EMD(residue6);
% C7=h;
% Nxunhuanr=0;                %用来记录循环的次数

% while(~(((NmaxX+NminX==Nzero)|(abs(((NmaxX+NminX)-Nzero)-1)))&(max(abs(mean_baoluo))<0.01)))
%     
%  [h,NmaxX,NminX,Nzero,up_baoluo,mean_baoluo]=EMD(C7);
%     C7=h;
%     Nxunhuanr=Nxunhuanr+1;
% end  
% subplot(7,2,13)
% plot(C7)                      %第三阶固有模态函数的波形图
% axis([0 4000 -0.3 0.3]);
% title('第七阶固有模态函数的波形图');

% residue7=residue6-C7;
% subplot(7,2,14)
% plot(residue7)            %第三次黄变换的残差
% title('第七次黄变换的残差');

pw=zeros(1,7);
for i=1:LenC
    pw(1)=pw(1)+(-C1(i)*C1(i)*log(C1(i)*C1(i)));%香农能量
    pw(2)=pw(2)+(-C2(i)*C2(i)*log(C2(i)*C2(i)));
    pw(3)=pw(3)+(-C3(i)*C3(i)*log(C3(i)*C3(i)));
    pw(4)=pw(4)+(-C4(i)*C4(i)*log(C4(i)*C4(i)));
    pw(5)=pw(5)+(-C5(i)*C5(i)*log(C5(i)*C5(i)));
    pw(6)=pw(6)+(-C6(i)*C6(i)*log(C6(i)*C6(i)));
%     pw(7)=pw(7)+(-C7(i)*C7(i)*log(C7(i)*C7(i)));
end
figure(3)
plot(pw)
% fid=fopen('out.txt','a+');
% fprintf(fid,'%6.4f %6.4f %6.4f %6.4f %6.4f %6.4f %6.4f',pw(1),pw(2),pw(3),pw(4),pw(5),pw(6),pw(7));
% fprintf(fid,'\n');
% fclose(fid);
%%
%[h,NmaxX,NminX,Nzero,up_baoluo,mean_baoluo]=EMD(residue7);
%C8=h;
%Nxunhuanr=0;                %用来记录循环的次数

%while(~(((NmaxX+NminX==Nzero)|(abs(((NmaxX+NminX)-Nzero)-1)))&(max(abs(mean_baoluo))<0.01)))
    
% [h,NmaxX,NminX,Nzero,up_baoluo,mean_baoluo]=EMD(C8);
%    C8=h;
%    Nxunhuanr=Nxunhuanr+1;
%end  
%subplot(8,2,15)
%plot(C8)                      %第三阶固有模态函数的波形图
%axis([0 4000 -0.3 0.3]);
%title('第八阶固有模态函数的波形图');

%residue8=residue7-C8;
%subplot(8,2,16)
%plot(residue8)            %第三次黄变换的残差
%title('第八次黄变换的残差');
%%
%经过两次黄变换后,残余的部分对原信号的意义不大,到此结束
xh=C1+C2;
figure(4)
plot(xh)    %画出第一阶和第二阶固有模态函数之和的波形图
title('第一阶和第二阶固有模态函数之和的波形图');
HX=hilbert(xh);
HC1=hilbert(C1);
HC2=hilbert(C2);
%baoluoC1=sqrt(real(HC1).*real(HC1)+imag(HC1).*imag(HC1)); %用两阶IMF之和进行运算
baoluoC1=abs(C1+HC1*j);
baoluoC2=sqrt(real(HC2).*real(HC2)+imag(HC2).*imag(HC2));
%baoluoX=sqrt(real(HX).*real(HX)+imag(HX).*imag(HX));
 figure(5)
plot(baoluoC1)

%title('标注S1和S2');

%figure(10)
%plot(baoluoC2)
%figure(11)
%plot(baoluoX)

%%%下面对baoluoC1进行3次样条插值运算
[h,NmaxX,NminX,Nzero,up_baoluo,mean_baoluo]=EMD(baoluoC1);
C1_baoluo=up_baoluo;

% figure(12)


%%%下面开始对C1_baoluo进行心音的定位分析。 %建立一个函数来进行定位分析,函数名称为DW。

L=length(C1_baoluo);
[NT1,T1]=DW(C1_baoluo);  %通过DW()函数实现了对S1和S2的初步定位   找出了大于阈值1的点所造的位置


%%%%下面开始对S1和S2开始精确定位,也就是要找到S1和S2的起始为止和终止位置
[NT2,T2,NT3,T3]=DW2(T1,C1_baoluo);%通过这个函数的调用得到了S1和S2的起始和终止位置,起始位置保存在T2中,终止位置保存在T3中

%如果T2,T3中有相同的坐标,应该剔除,并对数组进行重新赋值

plot(C1_baoluo)
hold on
[T2,T3]=jianyan(T2,T3)    %T2表示的是起始位置,T3表示的是终止位置
l=length(T2);
m=length(T3);
for i=1:l
   plot(T2(i),C1_baoluo(T2(i)),'r*')      %标注起始点的位置
    hold on
end
for i=1:m
    plot(T3(i),C1_baoluo(T3(i)),'go')      %标注终止点的位置 
    hold on
end
axis([0 4000 0 1]);
title('标注S1和S2');

%得到起始和终止位置后,首先就可以获取心率这个参数,分别从T2和T3计算心率,得到值后求平均,这样可以保证正确率
[heart_rate]=xinlv(T2,T3);     %调用了求心率的函数
 

  • 3
    点赞
  • 16
    收藏
    觉得还不错? 一键收藏
  • 打赏
    打赏
  • 0
    评论
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包

打赏作者

fpga和matlab

你的鼓励将是我创作的最大动力

¥1 ¥2 ¥4 ¥6 ¥10 ¥20
扫码支付:¥1
获取中
扫码支付

您的余额不足,请更换扫码支付或充值

打赏作者

实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值