下载
加入VIP
  • 专属下载特权
  • 现金文档折扣购买
  • VIP免费专区
  • 千万文档免费下载

上传资料

关闭

关闭

关闭

封号提示

内容

首页 matlab实现HHT程序

matlab实现HHT程序.doc

matlab实现HHT程序

马君琦
2019-06-25 0人阅读 举报 0 0 暂无简介

简介:本文档为《matlab实现HHT程序doc》,可适用于IT/计算机领域

clearallx=load('TXT')fs=N=length(x)t=:fs:(N)fsz=xc=emd(z)计算每个IMF分量及最后一个剩余分量residual与原始信号的相关性m,n=size(c)fori=:ma=corrcoef(c(i,:),z)xg(i)=a(,)endxgfori=:m计算各IMF的方差贡献率定义:方差为平方的均值减去均值的平方均值的平方imfp=mean(c(i,:),)^平方的均值imfp=mean(c(i,:)^,)各个IMF的方差mse(i)=mean(c(i,:)^,)mean(c(i,:),)^endmmse=sum(mse)fori=:mmse(i)=mean(c(i,:)^,)mean(c(i,:),)^方差百分比,也就是方差贡献率mseb(i)=mse(i)mmse*显示各个IMF的方差和贡献率end画出每个IMF分量及最后一个剩余分量residual的图形figure()fori=:mdisp('imf',intstr(i))disp(mse(i)mseb(i))endsubplot(m,,)plot(t,z)set(gca,'fontname','timesNewRoman')set(gca,'fontsize',)ylabel('signal','Amplitude')fori=:msubplot(m,,i)set(gcf,'color','w')plot(t,c(i,:),'k')set(gca,'fontname','timesNewRoman')set(gca,'fontsize',)ylabel('imf',intstr(i))endsubplot(m,,m)set(gcf,'color','w')plot(t,c(m,:),'k')set(gca,'fontname','timesNewRoman')set(gca,'fontsize',)ylabel('r',intstr(m))画出每个IMF分量及剩余分量residual的幅频曲线figure()subplot(m,,)set(gcf,'color','w')f,z=fft(t,z)plot(f,z,'k')set(gca,'fontname','timesNewRoman')set(gca,'fontsize',)ylabel('initialsignal',intstr(m),'Amplitude')fori=:msubplot(m,,i)set(gcf,'color','w')f,z=fft(t,c(i,:))plot(f,z,'k')set(gca,'fontname','timesNewRoman')set(gca,'fontsize',)ylabel('imf',intstr(i),'Amplitude')endsubplot(m,,m)set(gcf,'color','w')f,z=fft(t,c(m,:))plot(f,z,'k')set(gca,'fontname','timesNewRoman')set(gca,'fontsize',)ylabel('r',intstr(m),'Amplitude')hx=hilbert(z)xr=real(hx)xi=imag(hx)计算瞬时振幅sz=sqrt(xr^xi^)计算瞬时相位sx=angle(hx)计算瞬时频率dt=diff(t)dx=diff(sx)sp=dxdtfigure()plot(t(:N),sp)title('瞬时频率')计算HHT时频谱和边际谱A,fa,tt=hhspectrum(c)E,tt=toimage(A,fa,tt,length(tt))figure()disphhs(E,tt)二维图显示HHT时频谱,E是求得的HHT谱pausefigure()fori=:size(c,)faa=fa(i,:)FA,TT=meshgrid(faa,tt)三维图显示HHT时频图surf(FA,TT,E)title('HHT时频谱三维显示')holdonendholdoffE=flipud(E)fork=:size(E,)bjp(k)=sum(E(k,:))*fsendf=(:N)N*(fs)figure()plot(f,bjp)xlabel('频率Hz')ylabel('信号幅值')title('信号边际谱')要求边际谱必须先对信号进行EMD分解functionA,f,tt=hhspectrum(x,t,l,aff)error(nargchk(,,nargin))ifnargin<t=:size(x,)endifnargin<l=endifnargin<aff=endifmin(size(x))==ifsize(x,)==ifnargin<t=:size(x,)endendNmodes=elseNmodes=size(x,)endlt=length(t)tt=t((l):(ltl))fori=:Nmodesan(i,:)=hilbert(x(i,:)')'f(i,:)=instfreq(an(i,:)',tt,l)'A=abs(an(:,l:endl))ifaffdisprog(i,Nmodes,max(Nmodes,))endendfunctiondisphhs(im,t,inf)DISPHHS(im,t,inf)displaysinanewfigurethespectrumcontainedinmatrix"im"(amplitudesinlog)inputs:im:imagematrix(eg,outputof"toimage")t(optional):timeinstants(eg,outputof"toimage")inf(optional):dynamicrangeindB(wrtmax)default:inf=utilisation:disphhs(im)disphhs(im,t)disphhs(im,inf)disphhs(im,t,inf)figurecolormap(bone)colormap(colormap)ifnargin==inf=t=:size(im,)endifnargin==iflength(t)==inf=tt=:size(im,)elseendendifinf>=error('infdoitetre<')endM=max(max(im))im=log(imMe)inf=infimagesc(t,fliplr((:size(im,))(*size(im,))),im,inf,)set(gca,'YDir','normal')xlabel('time')ylabel('normalizedfrequency')title('HilbertHuangspectrum')functionf,z=fftfenxi(t,y)L=length(t)N=^nextpow(L)fft默认计算的信号是从开始的t=linspace(t(),t(L),N)deta=t()t()m=:Nf=(N*deta)*m下面计算的Y就是x(t)的傅里叶变换数值Y=exp(i**pi*f)*fft(y)将计算出来的频谱乘以exp(i**pi*f)得到频移后,之间的频谱值Y=fft(y)z=sqrt(Y*conj(Y))

用户评价(0)

关闭

新课改视野下建构高中语文教学实验成果报告(32KB)

抱歉,积分不足下载失败,请稍后再试!

提示

试读已结束,如需要继续阅读或者下载,敬请购买!

文档小程序码

使用微信“扫一扫”扫码寻找文档

1

打开微信

2

扫描小程序码

3

发布寻找信息

4

等待寻找结果

我知道了
评分:

/12

matlab实现HHT程序

VIP

在线
客服

免费
邮箱

爱问共享资料服务号

扫描关注领取更多福利