Matlab实现HHT程序(源码-非常珍贵)

Mat lab实现HHT程序（源码-非常珍贵）

clear all;

x=load (")；

fs=1000000;

N=length(x);

t=O:l/fs:(N-l)/fs;

z=x;

c=emd(z);

%计•算每个IMF分量及最后一个剩余分量residual与原始信号的相关性[m,n]=size(c);

for i=l:m;

a=corrcoef(c(i,:),z);

xg(i)=a(l,2);

end

xg；

for i=l:m-l

%计算各IMF的方差贡献率

%定义：方差为平方的均值减去均值的平方

%均值的平方

%imfp2=mean(c(i,:),2).A2

%平方的均值

%imf2p=mean(c(i,:).A2/2)

%各个IMF的方差

mse(i)=mean(c(i/:).A2/2)-mean(c(i,:)z2).A2;

end;

mmse=sum(mse);

for i=l:m-l

mse(i)=mea n(c(i/:).A2/2)-mean(c(i/:)z2).A2;

%方差百分比，也就是方差贡献率

mseb(i)=mse(i)/mmse*100;

%显示各个IMF的方差和贡献率

end;

%画出每个IMF分量及最后一个剩余分量residual的图形figure(l)

for i=l:m-l

disp(['imf'/int2str(i)]) ;disp([mse(i) mseb(i)]);

end;

subplot(m+lU)

plot(t’z)

set(gca/fontname7times New RomarV)

set(gca/f on tsize1,

ylabelfl'signar/Amplitude'])

for i=l:m-l

subplot(m+l z l,i+l);

setfgcf/colorVw1)

plot(t,c(i,:),k)

set(gca/fontname1,'times New Roman‘)

set(gca/f on tsize1,

ylabel(['imf'/int2str(i)])

end

subplot(m+l z l,m+l);

setlgcf/color'/w')

plot(tjC(m,:),'k‘)

set(gca/fontname1,'times New Roman‘)

set(gca/f on tsize1,

ylabel(['r\int2str(m-l)])

%画出每个IMF分量及剩余分量residual的幅频曲线figure(2)

subplot(m+blj)

setfgcf/color'/w1)

[f,z]=fft(t,z);

plot(f,z；k')

set(gca/fontname1,'times New Roman‘)

set(gca/f on tsize1,

ylabel『initial signal',int2str(m-l),'Amplitude'])

for i=l:m-l

subplot(m+l z l,i+l);

settgcf/color'/w1)

[f,z]=fft(t/c(i/:))；

plot(f,z,'k')

set(gca/fontname7times New RomarV)

set(gca/f on tsize1,

ylabel(['imf'/int2str(i)/'Amplitude,])

end

subplot(m+l,l,m+l);

settgcf/color'/w1)

[f/z]=fft(t,c(m/:));

plot(f,z；k')

set(gca/fontname7times New Roman‘)

set(gca/f on tsize1,

ylabel([,r,/int2str(m-l)；Amplitude,])

hx=hilbert(z);

xr=real(hx);xi=imag(hx);

%计算瞬时振幅

sz=sqrt(xr.A2+xi.A2);

%计算瞬时相位

sx=a ng le(hx);

%计算瞬时频率

dt=diff(t);

dx=diff(sx);

sp=dx./dt;

figure(6)

plot(t(l:N-l)/sp)

title/瞬时频率J

%讣算HHT时频谱和边际谱

[A,fa,tt]=hhspectrum(c);

[E/ttl]=toimage(A,fa,tt,length(tt));

figure(3)

disp_hhs(E,ttl) %二维图显示HHT时频谱,E是求得的HHT谱pause

figure(4)

for i=l:size(c,l)

faa=fa(i,:);

[FA/TTl]=meshgrid(faa,ttl);%三维图显示HHT 时频图

surf(FA,TTl,E)

title('HHT时频谱三维显示]

hold on

end

hold off

E=flipud(E);

for k=l:size(E,l)

bjp(k)=sum(E(k/))*l/fs;