Matlab实现HHT程序源码,非常珍贵).docVIP

PAGE PAGE 1 clearall; x=load(06514135360001170106.TXT); fs=1000000; N=length(x); t=0:1/fs:(N-1)/fs; z=x; c=emd(z); %计算每个IMF分量及最后一个剩余分量residual与原始信号的相关性 [m,n]=size(c); fori=1:m; a=corrcoef(c(i,:),z); xg(i)=a(1,2); end xg; fori=1:m-1 % %计算各IMF的方差贡献率 %定义：方差为平方的均值减去均值的平方 %均值的平方 %imfp2=mean(c(i,:),2).^2 %平方的均值 %imf2p=mean(c(i,:).^2,2) %各个IMF的方差 mse(i)=mean(c(i,:).^2,2)-mean(c(i,:),2).^2; end; mmse=sum(mse); fori=1:m-1 mse(i)=mean(c(i,:).^2,2)-mean(c(i,:),2).^2; %方差百分比，也就是方差贡献率 mseb(i)=mse(i)/mmse*100; %显示各个IMF的方差和贡献率 end; %画出每个IMF分量及最后一个剩余分量residual的图形 figure(1) fori=1:m-1 disp([imf,int2str(i)]);disp([mse(i)mseb(i)]); end; subplot(m+1,1,1) plot(t,z) set(gca,fontname,timesNewRoman) set(gca,fontsize,14.0) ylabel([signal,Amplitude]) fori=1:m-1 subplot(m+1,1,i+1); set(gcf,color,w) plot(t,c(i,:),k) set(gca,fontname,timesNewRoman) set(gca,fontsize,14.0) ylabel([imf,int2str(i)]) end subplot(m+1,1,m+1); set(gcf,color,w) plot(t,c(m,:),k) set(gca,fontname,timesNewRoman) set(gca,fontsize,14.0) ylabel([r,int2str(m-1)]) %画出每个IMF分量及剩余分量residual的幅频曲线 figure(2) subplot(m+1,1,1) set(gcf,color,w) [f,z]=fft(t,z); plot(f,z,k) set(gca,fontname,timesNewRoman) set(gca,fontsize,14.0) ylabel([initialsignal,int2str(m-1),Amplitude]) fori=1:m-1 subplot(m+1,1,i+1); set(gcf,color,w) [f,z]=fft(t,c(i,:)); plot(f,z,k) set(gca,fontname,timesNewRoman) set(gca,fontsize,14.0) ylabel([imf,int2str(i),Amplitude]) end subplot(m+1,1,m+1); set(gcf,color,w) [f,z]=fft(t,c(m,:)); plot(f,z,k) set(gca,fontname,timesNewRoman) set(gca,fontsize,14.0) ylabel([r,int2str(m-1),Amplitude]) hx=hilbert(z); xr=real(hx);xi=imag(hx); %计算瞬时振幅 sz=sqrt(xr.^2+xi.^2); %计算瞬时相位 sx=angle(hx); %计算瞬时频率 dt=diff(t); dx=diff(sx); sp=dx./dt; figure(6) plot(t(1:N-1),sp) title(瞬时频率) %计算HHT时频谱和边际谱 [A,fa,tt]=hhspectrum(c); [E,tt1]=toimage(A,fa,tt,length(tt)); figure(3) disp_hhs(E,tt1)%二维图显示HHT时频谱，E是求得的HHT谱 pause figure(4) fori=1:size(c,1) faa=fa(i,:); [FA,TT1]=meshgrid(faa,tt1);%三维图显示HHT时频图 surf(FA,TT1,E) title(HHT时频谱三维显示) holdon end holdoff E=