EMD算法的matlab程序介绍解析

%此版本为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(sd0.1|~isimf(x1%当标准偏差系数sd大于0.1或x1不是固有模态函数时,分量终止条件s1=getspline(x1;%上包络线s2=-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;ifu10,u=0;else,u=1;endfunctionu=isimf(x%u=0表示x不是固有模式函数,u=1表示x是固有模式函数N=length(x;u1=sum(x(1:N-1.*x(2:N0;u2=length(findpeaks(x+length(findpeaks(-x;ifabs(u1-u21,u=0;else,u=1;endfunctions=getspline(x%三次样条函数拟合成元数据包络线N=length(x;p=findpeaks(x;s=spline([0pN+1],[0x(p0],1:N;---------------------------------------------------------------------------------------------------------------------------------------------------------------functionn=findpeaks(x%Findpeaks.找到极值,n为极值点所在位置%n=findpeaks(xn=find(diff(diff(x00;u=find(x(n+1x(n;n(u=n(u+1;----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------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%GetHHT.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';imf=emd(x;fork=1:length(imfb(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:2plot(c,d{k},'k.','Color',b([kkk],'MarkerSize',3;holdon;set(gca,'FontSize',8,'XLim',[0c(end],'YLim',[050];%设置x、y轴句柄xlabel('Time',ylabel('Frequency';title('原信号时频图';end%SetIMFplots.M=length(imf;N=length(x;c=linspace(0,(N-1*Ts,N;fork1=0:4:M-1figurefork2=1:min(4,M-k1,subplot(4,1,k2,plot(c,imf{k1+k2};set(gca,'FontSize',8,'XLim',[0c(end];title('EMD分解结果';endxlabel('Time';end