EMD分解

clcclearallcloseall%[x,Fs]=wavread('Hum.wav');%Ts=1/Fs;%x=x(1:6000);Ts=0.001;Fs=1/Ts;t=0:Ts:1;x=sin(2*pi*10*t)+sin(2*pi*50*t)+sin(2*pi*100*t)+0.1*randn(1,length(t));imf=emd(x);plot_hht(x,imf,1/Fs);k=4;y=imf{k};N=length(y);t=0:Ts:Ts*(N-1);[yenvelope,yfreq,yh,yangle]=HilbertAnalysis(y,1/Fs);yModulate=y./yenvelope;[YMf,f]=FFTAnalysis(yModulate,Ts);Yf=FFTAnalysis(y,Ts);figuresubplot(321)plot(t,y)title(sprintf('IMF%d',k))xlabel('Time/s')ylabel(sprintf('IMF%d',k));subplot(322)plot(f,Yf)title(sprintf('IMF%d的频谱',k))xlabel('f/Hz')ylabel('|IMF(f)|');subplot(323)plot(t,yenvelope)title(sprintf('IMF%d的包络',k))xlabel('Time/s')ylabel('envelope');subplot(324)plot(t(1:end-1),yfreq)title(sprintf('IMF%d的瞬时频率',k))xlabel('Time/s')ylabel('Frequency/Hz');subplot(325)plot(t,yModulate)title(sprintf('IMF%d的调制信号',k))xlabel('Time/s')ylabel('modulation');subplot(326)plot(f,YMf)title(sprintf('IMF%d调制信号的频谱',k))xlabel('f/Hz')ylabel('|YMf(f)|');findpeaks.m文件functionn=findpeaks(x)%Findpeaks.找极大值点，返回对应极大值点的坐标n=find(diff(diff(x)0)0);%相当于找二阶导小于0的点u=find(x(n+1)x(n));n(u)=n(u)+1;%加1才真正对应极大值点%图形解释上述过程%figure%subplot(611)%x=x(1:100);%plot(x,'-o')%gridon%%subplot(612)%plot(1.5:length(x),diff(x)0,'-o')%gridon%axis([1,length(x),-0.5,1.5])%%subplot(613)%plot(2:length(x)-1,diff(diff(x)0),'-o')%gridon%axis([1,length(x),-1.5,1.5])%%subplot(614)%plot(2:length(x)-1,diff(diff(x)0)0,'-o')%gridon%axis([1,length(x),-1.5,1.5])%%n=find(diff(diff(x)0)0);%subplot(615)%plot(n,ones(size(n)),'o')%gridon%axis([1,length(x),0,2])%%u=find(x(n+1)x(n));%n(u)=n(u)+1;%subplot(616)%plot(n,ones(size(n)),'o')%gridon%axis([1,length(x),0,2])plot_hht.m文件functionplot_hht(x,imf,Ts)%PlottheHHT.%::Syntax%ThearrayxistheinputsignalandTsisthesamplingperiod.%Exampleonuse:[x,Fs]=wavread('Hum.wav');%plot_hht(x(1:6000),1/Fs);%Func:emd%imf=emd(x);fork=1:length(imf)b(k)=sum(imf{k}.*imf{k});th=unwrap(angle(hilbert(imf{k})));%相位d{k}=diff(th)/Ts/(2*pi);%瞬时频率end[u,v]=sort(-b);b=1-b/max(b);%后面绘图的亮度控制%Hilbert瞬时频率图N=length(x);c=linspace(0,(N-2)*Ts,N-1);%0:Ts:Ts*(N-2)fork=v(1:2)%显示能量最大的两个IMF的瞬时频率figureplot(c,d{k});xlim([0c(end)]);ylim([01/2/Ts]);xlabel('Time/s')ylabel('Frequency/Hz');title(sprintf('IMF%d',k))end%显示各IMFM=length(imf);N=length(x);c=linspace(0,(N-1)*Ts,N);%0:Ts:Ts*(N-1)fork1=0:4:M-1figurefork2=1:min(4,M-k1)subplot(4,2,2*k2-1)plot(c,imf{k1+k2})set(gca,'FontSize',8,'XLim',[0c(end)]);title(sprintf('第%d个IMF',k1+k2))xlabel('Time/s')ylabel(sprintf('IMF%d',k1+k2));subplot(4,2,2*k2)[yf,f]=FFTAnalysis(imf{k1+k2},Ts);plot(f,yf)title(sprintf('第%d个IMF的频谱',k1+k2))xlabel('f/Hz')ylabel('|IMF(f)|');endendfiguresubplot(211)plot(c,x)set(gca,'FontSize',8,'XLim',[0c(end)]);title('原始信号')xlabel('Time/s')ylabel('Origin');subplot(212)[Yf,f]=FFTAnalysis(x,Ts);plot(f,Yf)title('原始信号的频谱')xlabel('f/Hz')ylabel('|Y(f)|');emd.m文件functionimf=emd(x)%EmpiricialModeDecomposition(Hilbert-HuangTransform)%EMD分解或HHT变换%返回值为cell类型，依次为一次IMF、二次IMF、...、最后残差x=transpose(x(:));imf=[];while~ismonotonic(x)x1=x;sd=Inf;while(sd0.1)||~isimf(x1)s1=getspline(x1);%极大值点样条曲线s2=-getspline(-x1);%极小值点样条曲线x2=x1-(s1+s2)/2;sd=sum((x1-x2).^2)/sum(x1.^2);x1=x2;endimf{end+1}=x1;x=x-x1;endimf{end+1}=x;%是否单调functionu=ismonotonic(x)u1=length(findpeaks(x))*length(findpeaks(-x));ifu10u=0;elseu=1;end%是否IMF分量functionu=isimf(x)N=length(x);u1=sum(x(1:N-1).*x(2:N)0);%过零点的个数u2=length(findpeaks(x))+length(findpeaks(-x));%极值点的个数ifabs(u1-u2)1u=0;elseu=1;end%据极大值点构造样条曲线functions=getspline(x)N=length(x);p=findpeaks(x);s=spline([0pN+1],[0x(p)0],1:N);FFTAnalysis.m文件%频谱分析function[Y,f]=FFTAnalysis(y,Ts)Fs=1/Ts;L=length(y);NFFT=2^nextpow2(L);y=y-mean(y);Y=fft(y,NFFT)/L;Y=2*abs(Y(1:NFFT/2+1));f=Fs/2*linspace(0,1,NFFT/2+1);endHilbertAnalysis.m文件%Hilbert分析function[yenvelope,yf,yh,yangle]=HilbertAnalysis(y,Ts)yh=hilbert(y);yenvelope=abs(yh);%包络yangle=unwrap(angle(yh));%相位yf=diff(yangle)/2/pi/Ts;%瞬时频率end