Matlab实现HHT程序源码非常珍贵.docx
《Matlab实现HHT程序源码非常珍贵.docx》由会员分享，可在线阅读，更多相关《Matlab实现HHT程序源码非常珍贵.docx（5页珍藏版）》请在冰豆网上搜索。
Matlab实现HHT程序源码非常珍贵
Matlab实现HHT程序（源码-非常珍贵）
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
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=flipud（E）;
fork=1:
size（E,1）
bjp（k）=sum（E（k,:
））*1/fs;
end
f=（1:
N-2）/N*（fs/2）;
figure（5）
plot（f,bjp）;
xlabel（'频率/Hz'）;
ylabel（'信号幅值'）;
title（'信号边际谱'）%要求边际谱必须先对信号进行EMD分解
function[A,f,tt]=hhspectrum（x,t,l,aff）
error（nargchk（1,4,nargin））;
ifnargin<2
t=1:
size（x,2）;
end
ifnargin<3
l=1;
end
ifnargin<4
aff=0;
end
ifmin（size（x））==1
ifsize（x,2）==1
x=x';
ifnargin<2
t=1:
size（x,2）;
end
end
Nmodes=1;
else
Nmodes=size（x,1）;
end
lt=length（t）;
tt=t（（l+1）:
（lt-l））;
fori=1:
Nmodes
an（i,:
）=hilbert（x（i,:
）'）';
f（i,:
）=instfreq（an（i,:
）',tt,l）';
A=abs（an（:
l+1:
end-l））;
ifaff
disprog（i,Nmodes,max（Nmodes,100））
end
end
functiondisp_hhs（im,t,inf）
%DISP_HHS（im,t,inf）
%displaysinanewfigurethespectrumcontainedinmatrix"im"
%（amplitudesinlog）.
%
%inputs:
-im:
imagematrix（e.g.,outputof"toimage"）
%-t（optional）:
timeinstants（e.g.,outputof"toimage"）
%-inf（optional）:
-dynamicrangeindB（wrtmax）
%default:
inf=-20
%
%utilisation:
disp_hhs（im）;disp_hhs（im,t）;disp_hhs（im,inf）
%disp_hhs（im,t,inf）
figure
colormap（bone）
colormap（1-colormap）;
ifnargin==1
inf=-20;
t=1:
size（im,2）;
end
ifnargin==2
iflength（t）==1
inf=t;
t=1:
size（im,2）;
else
inf=-20;
end
end
ifinf>=0
error（'infdoitetre<0'）
end
M=max（max（im））;
im=log10（im/M+1e-300）;
inf=inf/10;
imagesc（t,fliplr（（1:
size（im,1））/（2*size（im,1）））,im,[inf,0]）;
set（gca,'YDir','normal'）
xlabel（['time']）
ylabel（['normalizedfrequency']）
title（'Hilbert-Huangspectrum'）
function[f,z]=fftfenxi（t,y）
L=length（t）;N=2^nextpow2（L）;
%fft默认计算的信号是从0开始的
t=linspace（t
（1）,t（L）,N）;deta=t
（2）-t
（1）;
m=0:
N-1;
f=1./（N*deta）*m;
%下面计算的Y就是x（t）的傅里叶变换数值
%Y=exp（i*4*pi*f）.*fft（y）%将计算出来的频谱乘以exp（i*4*pi*f）得到频移后[-2,2]之间的频谱值
Y=fft（y）;
z=sqrt（Y.*conj（Y））;
升级会员 