Matlab实现HHT程序(源码,非常珍贵)[参考].pdfVIP

精品文档 ·可编辑版 clear all; 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); for i=1:m; a=corrcoef(c(i,:),z); xg(i)=a(1,2); end xg; for i=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); for i=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) for i=1:m-1 disp([imf,int2str(i)]) ;disp([mse(i) mseb(i)]); end; subplot(m+1,1,1) plot(t,z) set(gca,fontname,times New Roman) set(gca,fontsize,14.0) ylabel([signal,Amplitude]) for i=1:m-1 subplot(m+1,1,i+1); 1 / 5 精品文档 ·可编辑版 set(gcf,color,w) plot(t,c(i,:),k) set(gca,fontname,times New Roman) 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,times New Roman) 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,times New Roman) set(gca,fontsize,14.0) ylabel([initial signal,int2str(m-1),Amplitude]) for i=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,times New Roman) 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,times New Roman) 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=angl