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); %调用了求心率的函数