資源描述：
《emd 希爾伯特黃變換程序》由會(huì)員上傳分享，免費(fèi)在線閱讀，更多相關(guān)內(nèi)容在行業(yè)資料-天天文庫。
1、（一）簡單的EMD程序functionimf=emd(x)%EmpiricialModeDecomposition(Hilbert-HuangTransform)%imf=emd(x)%Func:findpeaksx??=transpose(x(:));%轉(zhuǎn)置imf=[];while~ismonotonic(x)%當(dāng)x不是單調(diào)函數(shù)，分解終止條件??x1=x;??sd=Inf;%均值%直到x1滿足IMF條件，得c1??while(sd>0.1)
2、~isimf(x1)%當(dāng)標(biāo)準(zhǔn)偏差系數(shù)sd大于0.1或x1不是固有模態(tài)函數(shù)時(shí)，分量終止條件????s1=gets
3、pline(x1);%上包絡(luò)線????s2=-getspline(-x1);%下包絡(luò)線????x2=x1-(s1+s2)/2;%此處的x2為文章中的h????????sd=sum((x1-x2).^2)/sum(x1.^2);????x1=x2;??end????imf{end+1}=x1;??x??????=x-x1;endimf{end+1}=x;%FUNCTIONSfunctionu=ismonotonic(x)%u=0表示x不是單調(diào)函數(shù)，u=1表示x為單調(diào)的u1=length(findpeaks(x))*length(findpeaks(-x)
4、);ifu1>0,u=0;else,????u=1;endfunctionu=isimf(x)%u=0表示x不是固有模式函數(shù)，u=1表示x是固有模式函數(shù)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)%三次樣條函數(shù)擬合成元數(shù)據(jù)包絡(luò)線N=length(x);p=findpeaks(x);s=spline([
5、0pN+1],[0x(p)0],1:N);--------------------------------------------------------------------------------functionn=findpeaks(x)%Findpeaks.找到極值%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)%??雙邊帶調(diào)幅信號(hào)的EMD分解%PlottheHHT.%plot_hht(x,Ts)%%::Syntax%??ThearrayxistheinputsignalandTsisthesamplingperiod.%??Exampleonuse:[x,Fs]=wavread('Hum
7、.wav');%??????????????plot_hht(x(1:6000),1/Fs);%Func:emd%GetHHT.clearall;closeall;Ts=0.0005;t=0:Ts:1;??????????????%采樣率2000HZ%調(diào)幅信號(hào)x=sin(2*pi*t).*sin(40*pi*t);s1=getspline(x);%上包絡(luò)線s2=-getspline(-x);%上包絡(luò)線x1=(s1+s2)/2;%此處的x2為文章中的hfigure;plot(t,x);xlabel('Time'),ylabel('Amplitude');
8、title('雙邊帶調(diào)幅信號(hào)');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
9、);%figure;fork=v(1:2)??plot(c,d{k},'k.','Color'
