资源描述：
《EMD算法的matlab程序介绍.docx》由会员上传分享，免费在线阅读，更多相关内容在行业资料-天天文库。
1、%此版本为ALAN版本的整合注释版functionimf=emd(x)%EmpiricialModeDecomposition(Hilbert-HuangTransform)%imf=emd(x)%Func:findpeaksx=transpose(x(:));%转置为行矩阵imf=[];while~ismonotonic(x)%当x不是单调函数，分解终止条件x1=x;sd=Inf;%均值%直到x1满足IMF条件，得c1while(sd>0.1)
2、~isimf(x1)%当标准偏差系数sd大于0.1或x1不是固有模态函数时，分量终止条件s1=getspline(x1);%上包络线s2=
3、-getspline(-x1);%下包络线x2=x1-(s1+s2)/2;%此处的x2为文章中的hsd=sum((x1-x2).^2)/sum(x1.^2);x1=x2;endimf{end+1}=x1;x=x-x1;endimf{end+1}=x;%FUNCTIONSfunctionu=ismonotonic(x)%u=0表示x不是单调函数，u=1表示x为单调的u1=length(findpeaks(x))*length(findpeaks(-x));ifu1>0,u=0;else,u=1;endfunctionu=isimf(x)%u=0表示x不是固有模式函数，u=1表示x是固有
4、模式函数N=length(x);u1=sum(x(1:N-1).*x(2:N)<0);u2=length(findpeaks(x))+length(findpeaks(-x));ifabs(u1-u2)>1,u=0;else,u=1;endfunctions=getspline(x)%三次样条函数拟合成元数据包络线N=length(x);p=findpeaks(x);s=spline([0pN+1],[0x(p)0],1:N);-------------------------------------------------------------------------------
5、--------------------------------------------------------------------------------functionn=findpeaks(x)%Findpeaks.找到极值,n为极值点所在位置%n=findpeaks(x)n=find(diff(diff(x)>0)<0);u=find(x(n+1)>x(n));n(u)=n(u)+1;-----------------------------------------------------------------------------------------------
6、-----------------------------------------------------------------------------------functionplot_hht00(x,Ts)%双边带调幅信号的EMD分解%PlottheHHT.%plot_hht(x,Ts)%%::Syntax%Thearray（列）xistheinputsignalandTsisthesamplingperiod（取样周期）.%Exampleonuse:[x,Fs]=wavread('Hum.wav');%plot_hht(x(1:6000),1/Fs);%Func:emd%G
7、etHHT.clearall;closeall;Ts=0.0005;t=0:Ts:10;%采样率2000HZ%调幅信号%x=sin(2*pi*t).*sin(40*pi*t);x=sin(2*pi*t);s1=getspline(x);%上包络线s2=-getspline(-x);%上包络线x1=(s1+s2)/2;%此处的x2为文章中的hfigure;plot(t,x);xlabel('Time'),ylabel('Amplitude');title('双边带调幅信号');holdon;plot(t,s1,'-r');plot(t,s2,'-r');plot(t,x1,'g');i
8、mf=emd(x);fork=1:length(imf)b(k)=sum(imf{k}.*imf{k});th=angle(hilbert(imf{k}));d{k}=diff(th)/Ts/(2*pi);end[u,v]=sort(-b);b=1-b/max(b);%Settime-frequencyplots.N=length(x);c=linspace(0,(N-2)*Ts,N-1);%figure;fork=v(1:2)plot(c,d{k},'k
