Matlab实现HHT程序源码非常珍贵Word文档下载推荐.docx
《Matlab实现HHT程序源码非常珍贵Word文档下载推荐.docx》由会员分享，可在线阅读，更多相关《Matlab实现HHT程序源码非常珍贵Word文档下载推荐.docx（5页珍藏版）》请在冰豆网上搜索。
end
xg;
m-1
%--------------------------------------------------------------------
%计算各IMF的方差贡献率
%定义：
方差为平方的均值减去均值的平方
%均值的平方
%imfp2=mean（c（i,:
）,2）.^2
ylabel（['
signal'
'
Amplitude'
]）
subplot（m+1,1,i+1）;
set（gcf,'
color'
w'
）
plot（t,c（i,:
）,'
k'
set（gca,'
fontname'
timesNewRoman'
fontsize'
14.0）
imf'
int2str（i）]）
subplot（m+1,1,m+1）;
plot（t,c（m,:
r'
int2str（m-1）]）
%画出每个IMF分量及剩余分量residual的幅频曲线
figure
（2）
subplot（m+1,1,1）
[f,z]=fft（t,z）;
plot（f,z,'
initialsignal'
int2str（m-1）,'
[f,z]=fft（t,c（i,:
））;
int2str（i）,'
[f,z]=fft（t,c（m,:
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）
size（c,1）
faa=fa（i,:
[FA,TT1]=meshgrid（faa,tt1）;
%三维图显示HHT时频图
surf（FA,TT1,E）
HHT时频谱三维显示'
holdon
holdoff
E=flipud（E）;
fork=1:
size（E,1）
bjp（k）=sum（E（k,:
））*1/fs;
f=（1:
N-2）/N*（fs/2）;
figure（5）
plot（f,bjp）;
xlabel（'
频率/Hz'
ylabel（'
信号幅值'
信号边际谱'
）%要求边际谱必须先对信号进行EMD分解
function[A,f,tt]=hhspectrum（x,t,l,aff）
error（nargchk（1,4,nargin））;
ifnargin<
2
t=1:
size（x,2）;
3
l=1;
4
aff=0;
ifmin（size（x））==1
ifsize（x,2）==1
x=x'
;
t=1:
Nmodes=1;
else
Nmodes=size（x,1）;
lt=length（t）;
tt=t（（l+1）:
（lt-l））;
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））
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"
）
%-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;
size（im,2）;
ifnargin==2
iflength（t）==1
inf=t;
inf=-20;
ifinf>
=0
error（'
infdoitetre<
0'
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]）;
YDir'
normal'
xlabel（['
time'
normalizedfrequency'
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））;
升级会员 