资源描述：
《emd 希尔伯特黄变换程序》由会员上传分享，免费在线阅读，更多相关内容在行业资料-天天文库。
1、（一）简单的EMD程序functionimf=emd(x)%EmpiricialModeDecomposition(Hilbert-HuangTransform)%imf=emd(x)%Func:findpeaksx =transpose(x(:));%转置imf=[];while~ismonotonic(x)%当x不是单调函数，分解终止条件 x1=x; sd=Inf;%均值%直到x1满足IMF条件，得c1 while(sd>0.1)
2、~isimf(x1)%当标准偏差系数sd大于0.1或x1不是固有模态函数时，分量终止条件 s1=gets
3、pline(x1);%上包络线 s2=-getspline(-x1);%下包络线 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不是单调函数，u=1表示x为单调的u1=length(findpeaks(x))*length(findpeaks(-x)
4、);ifu1>0,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([
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)% 双边带调幅信号的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%调幅信号x=sin(2*pi*t).*sin(40*pi*t);s1=getspline(x);%上包络线s2=-getspline(-x);%上包络线x1=(s1+s2)/2;%此处的x2为文章中的hfigure;plot(t,x);xlabel('Time'),ylabel('Amplitude');
8、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
9、);%figure;fork=v(1:2) plot(c,d{k},'k.','Color'
