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: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);---------------------------------------------------------------------------------------------------------------------------------------------------------------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;----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------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(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.','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