用希尔伯特黄变换HHT求时频谱和边际谱.docx

1、用希尔伯特黄变换HHT求时频谱和边际谱【原创】用希尔伯特黄变换（HHT）求时频谱和边际谱寒假将至，精心将自己最近做的东西总结了一下，能跟大家分享讨论是我的荣幸。源代码也贴出来了，希望大家能提出宝贵意见顺祝大家寒假快乐，新年快乐1.什么是HHT？HHT就是先将信号进行经验模态分解（EMD分解），然后将分解后的每个IMF分量进行Hilbert变换，得到信号的时频属性的一种时频分析方法。2.EMD分解的步骤。截图07.png (51.93 KB2010-2-5 18:30截图08.png (31.36 KB2010-2-5 18:30截图09.png (30.93 KB2010-2-5 18:30E

2、MD分解的流程图如下：101.PNG (79.44 KB2010-2-5 18:493.实例演示。给定频率分别为10Hz和35Hz的两个正弦信号相叠加的复合信号，采样频率fs=2048Hz的信号，表达式如下：y=5sin(2*pi*10t+5*sin(2*pi*35t(1为了对比，先用fft对求上述信号的幅频和相频曲线。 复制内容到剪贴板 代码:function fftfenxiclear;clc;N=2048;%fft默认计算的信号是从0开始的t=linspace(1,2,N;deta=t(2-t(1;1/detax=5*sin(2*pi*10*t+5*sin(2*pi*35*t;% N1=

3、256;N2=512;w1=0.2*2*pi;w2=0.3*2*pi;w3=0.4*2*pi;% x=(t=-200&t-200+N1*deta&t-200+N2*deta&t=200.*sin(w3*t;y = x;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;plot(f(1:100,z(1:100;title(幅频曲线xiangwei=an

4、gle(Y;figure(2plot(f,xiangweititle(相频曲线figure(3plot(t,y,r%axis(-2,2,0,1.2title(原始信号102.PNG (57.26 KB2010-2-5 18:42103.PNG (24.85 KB2010-2-5 18:42104.PNG (25.65 KB2010-2-5 18:42（2）用Hilbert变换直接求该信号的瞬时频率 复制内容到剪贴板 代码:clear;clc;clf;%假设待分析的函数是z=t3N=2048;%fft默认计算的信号是从0开始的t=linspace(1,2,N;deta=t(2-t(1;fs=1/

5、deta;x=5*sin(2*pi*10*t+5*sin(2*pi*35*t;z=x;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;plot(t(1:N-1,sptitle(瞬时频率106.PNG (35.92 KB2010-2-5 18:42小结：傅里叶变换不能得到瞬时频率，即不能得到某个时刻的频率值。Hilbert变换是求取瞬时频率的方法，但如果只用Hilbert变换求出来的瞬时频率也不准确。（出现负频，

6、实际上负频没有意义！）（3）用HHT求取信号的时频谱与边际谱 复制内容到剪贴板 代码:function HHTclear;clc;clf;N=2048;%fft默认计算的信号是从0开始的t=linspace(1,2,N;deta=t(2-t(1;fs=1/deta;x=5*sin(2*pi*10*t+5*sin(2*pi*35*t;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;endxg;for i=1:m-1%-%计算各IMF的方差贡献率%定

7、义：方差为平方的均值减去均值的平方%均值的平方%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-1mse(i=mean(c(i,:.2,2-mean(c(i,:,2.2; %方差百分比，也就是方差贡献率mseb(i=mse(i/mmse*100;%显示各个IMF的方差和贡献率end;%画出每个IMF分量及最后一个剩余分量residual的图形figure(1for i=1:m-1disp(imf,int2s

8、tr(i ;disp(mse(i mseb(i;end;subplot(m+1,1,1plot(t,zset(gca,fontname,times New Romanset(gca,fontsize,14.0ylabel(signal,Amplitudefor i=1:m-1subplot(m+1,1,i+1;set(gcf,color,wplot(t,c(i,:,kset(gca,fontname,times New Romanset(gca,fontsize,14.0ylabel(imf,int2str(iendsubplot(m+1,1,m+1;set(gcf,color,wplot(t

9、,c(m,:,kset(gca,fontname,times New Romanset(gca,fontsize,14.0ylabel(r,int2str(m-1%画出每个IMF分量及剩余分量residual的幅频曲线figure(2subplot(m+1,1,1set(gcf,color,wf,z=fftfenxi(t,z;plot(f,z,kset(gca,fontname,times New Romanset(gca,fontsize,14.0ylabel(initial signal,int2str(m-1,Amplitudefor i=1:m-1subplot(m+1,1,i+1;s

10、et(gcf,color,wf,z=fftfenxi(t,c(i,:;plot(f,z,kset(gca,fontname,times New Romanset(gca,fontsize,14.0ylabel(imf,int2str(i,Amplitudeendsubplot(m+1,1,m+1;set(gcf,color,wf,z=fftfenxi(t,c(m,:;plot(f,z,kset(gca,fontname,times New Romanset(gca,fontsize,14.0ylabel(r,int2str(m-1,Amplitudehx=hilbert(z;xr=real(h

11、x;xi=imag(hx;%计算瞬时振幅sz=sqrt(xr.2+xi.2;%计算瞬时相位sx=angle(hx;%计算瞬时频率dt=diff(t;dx=diff(sx;sp=dx./dt;figure(6plot(t(1:N-1,sptitle(瞬时频率%计算HHT时频谱和边际谱A,fa,tt=hhspectrum(c;E,tt1=toimage(A,fa,tt,length(tt;figure(3disp_hhs(E,tt1 %二维图显示HHT时频谱，E是求得的HHT谱pausefigure(4for i=1:size(c,1faa=fa(i,:;FA,TT1=meshgrid(faa,t

12、t1;%三维图显示HHT时频图surf(FA,TT1,Etitle(HHT时频谱三维显示hold onendhold offE=flipud(E;for k=1:size(E,1bjp(k=sum(E(k,:*1/fs; endf=(1:N-2/N*(fs/2;figure(5plot(f,bjp;xlabel(频率 / Hz;ylabel(信号幅值;title(信号边际谱%要求边际谱必须先对信号进行EMD分解function A,f,tt = hhspectrum(x,t,l,afferror(nargchk(1,4,nargin;if nargin 2t=1:size(x,2;endif

13、nargin 3l=1;endif nargin 4aff = 0;endif min(size(x = 1if size(x,2 = 1x = x;if nargin = 0error(inf doit etre 0endM=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,normalxlabel(timeylabel(normalized frequencytitle(Hilbert-Huang spectrumfunct

14、ion f,z=fftfenxi(t,yL=length(t;N=2nextpow2(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;107.PNG (88.8 KB2010-2-5 18:42108.PNG (93.91 KB2010-2-5 18:42109.PNG (55.65

15、KB2010-2-5 18:42110.PNG (108.1 KB2010-2-5 18:42111.PNG (37.52 KB2010-2-5 18:424.总结。（1）边际谱与傅里叶谱的比较： 意义不同：边际谱从统计意义上表征了整组数据每个频率点的累积幅值分布，而傅里叶频谱的某一点频率上的幅值表示在整个信号里有一个含有此频率的三角函数组分。 作用不同：边际谱可以处理非平稳信号，如果信号中存在某一频率的能量出现，就表示一定有该频率的振动波出现，也就是说，边际谱能比较准确地反映信号的实际频率成分。而傅里叶变换只能处理平稳信号。（2） HHT与Hilbert变换的比较： Hilbert变换只是单纯地求信号的瞬时振幅，频率和相位，有可能出现没有意义的负频率；HHT变换先将信号进行EMD分解，得到的是各个不同尺度的分量，对每一个分量进行Hilbert变换后得到的是有实际意义的瞬时频率。PS:运行上面的程序需要装时频工具箱，我仅将用到的emd分解的程序贴到下面。