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

1、精品文档用希尔伯特黄变换(hht求时频谱和边际谱1 .什么是hht？hht1是先将信号进行经验模态分解(emd解)，然后将分解后的每个imf分量进行hilbert变换，得到信号的时频属性的一种时频分析方法。2 .emd分解的步骤。“防解的流程图如下：3 . 实例演示。给定频率分别为 10hz 和 35hz 的两个正弦信号相叠加的复合信号， 采样频率 fs=2048hz 的信号，表达式如下： y=5sin(2*pi*10t)+5*sin(2*pi*35t)(1) 为了对比，先用 fft 对求上述信号的幅频和相频曲线。function fftfenxiclear;clc;n=2048;%fft 默

2、认计算的信号是从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=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)%

3、将计算出来的频谱乘以 exp(i*4*pi*f) 得到频移后-2,2 之间的频谱值y=fft(y);z=sqrt(y.*conj(y);plot(f(1:100),z(1:100);title( 幅频曲线 )xiangwei=angle(y);figure(2)plot(f,xiangwei)title( 相频曲线 )figure(3)plot(t,y,r)%axis(-2,2,0,1.2)title( 原始信号 ) 复制代码( 2 )用hilbert 变换直接求该信号的瞬时频率clear;clc;clf;设待分析的函数是z=ta3n=2048;%fft 默认计算的信号是从0 开始的t=lin

4、space(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;hx=hilbert(z);xr=real(hx);xi=imag(hx);%计算瞬时振幅sz=sqrt(xr.a2+xi.a2);%计算瞬时相位sx=angle(hx);%计算瞬时频率dt=diff(t);dx=diff(sx);sp=dx./dt;plot(t(1:n-1),sp)title( 瞬时频率 )复制代码hilbert 变换是求取（出现负频，实际小结：傅里叶变换不能得到瞬时频率，即不能得到某个时刻的频率值。瞬时频率的方法，但如

5、果只用 hilbert 变换求出来的瞬时频率也不准确。上负频没有意义！ )(3)用hht求取信号的时频谱与边际谱function hht clear;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); endx

6、g;for i=1:m-1%计算各imf 的方差贡献率%定义：方差为平方的均值减去均值的平方%均值的平方%imfp2=mean(c(i,:),2).a2%平方的均值%imf2p=mean(c(i,:).a2,2)%各个imf 的方差mse(i)=mean(c(i,:).a2,2)-mean(c(i,:),2).a2;end;mmse=sum(mse);for i=1:m-1mse(i)=mean(c(i,:).a2,2)-mean(c(i,:),2).a2;%方差百分比，也就是方差贡献率mseb(i)=mse(i)/mmse*100;%显示各个imf 的方差和贡献率end;%画出每个imf 分

7、量及最后一个剩余分量residual 的图形figure(1)for i=1:m-1disp(imf,int2str(i) ;disp(mse(i) mseb(i);end;subplot(m+1,1,1)plot(t,z)set(gca,fontname,times new roman)set(gca,fontsize,14.0)ylabel(signal,amplitude)for i=1:m-1subplot(m+1,1,i+1);set(gcf,color,w)plot(t,c(i,:),k)set(gca,fontname,times new roman)set(gca,fontsi

8、ze,14.0)ylabel(imf,int2str(i)endsubplot(m+1,1,m+1);set(gcf,color,w)plot(t,c(m,:),k)set(gca,fontname,times new roman)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=fftfenxi(t,z);plot(f,z,k)set(gca,fontname,times new roman)set(gca,

9、fontsize,14.0)ylabel(initial signal,int2str(m-1),amplitude)for i=1:m-1subplot(m+1,1,i+1);set(gcf,color,w)f,z=fftfenxi(t,c(i,:);plot(f,z,k)set(gca,fontname,times new roman)set(gca,fontsize,14.0)ylabel(imf,int2str(i),amplitude)endsubplot(m+1,1,m+1);set(gcf,color,w)f,z=fftfenxi(t,c(m,:);plot(f,z,k)set(

10、gca,fontname,times new roman)set(gca,fontsize,14.0)ylabel(r,int2str(m-1),amplitude)hx=hilbert(z);xr=real(hx);xi=imag(hx);%计算瞬时振幅sz=sqrt(xr.a2+xi.a2);%计算瞬时相位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,f

11、a,tt,length(tt);figure(3)disp_hhs(e,tt1) % 二维图显示hht时频谱，e是求得的hht谱 pausefigure(4)for i=1:size(c,1)faa=fa(i,:);fa,tt1=meshgrid(faa,tt1);%三维图显示 hht寸频图surf(fa,tt1,e)title(hht 时频谱三维显示)hold onendhold offe=flipud(e);for k=1:size(e,1)bjp(k)=sum(e(k,:)*1/fs;endf=(1:n-2)/n*(fs/2);figure(5)plot(f,bjp);xlabel( 频

12、率 / hz);ylabel( 信号幅值 );title( 信号边际谱)要求边际谱必须先对信号进行em防解function a,f,tt = hhspectrum(x,t,l,aff)error(nargchk(1,4,nargin);if nargin 2t=1:size(x,2);endif nargin 3l=1;endif nargin 4aff = 0;endif min(size(x) = 1if size(x,2) = 1x = x;if nargin = 0error(inf doit etre 0)endm=max(max(im);im = log10(im/m+1e-300

13、);inf=inf/10;imagesc(t,fliplr(1:size(im,1)/(2*size(im,1),im,inf,0);set(gca,ydir,normal)xlabel(time)ylabel(normalized frequency)title(hilbert-huang spectrum)function f,z=fftfenxi(t,y)l=length(t);n=2anextpow2(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

14、就是 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);复制代码4 . 总结。( 1 )边际谱与傅里叶谱的比较：意义不同：边际谱从统计意义上表征了整组数据每个频率点的累积幅值分布，而傅里叶频谱的某一点频率上的幅值表示在整个信号里有一个含有此频率的三角函数组分。作用不同：边际谱可以处理非平稳信号，如果信号中存在某一频率的能量出现，就表示一定有该频率的振动波出现， 也就是说， 边际谱能比较准确地反映信号的实际频率成分。而傅里叶变换只能处理平稳信号。(2) hht 与 hilbert