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

HHT

(51.93KB

2010-2-518:30

(31.36KB

2010-2-518:30

(30.93KB

2010-2-518:30

(79.44KB

2010-2-518:49

实例演示。

给定频率分别为10Hz和35Hz的两个正弦信号相叠加的复合信号，采样频率

fs=2048Hz的信号，表达式如下：y=5sin(2*pi*10t+5*sin(2*pi*35t

(1为了对比，先用fft对求上述信号的幅频和相频曲线。

复制内容到剪贴板

代码:

functionfftfenxiclear;clc;N=2048;

%fft默认计算的信号是从开始的

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=-200t=-200+N1*deta.*sin(w1*t+(t-200+N1*detat=-200+N2*deta.*sin(w2*t+(t-200+N2*detat=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=angle(Y;figure(2plot(f,xiangweititle(相频曲线figure(3

plot(t,y,r

%axis([-2,2,0,1.2]

title(原始信号

用Hilbert变换直接求该信号的瞬时频率

复制内容到剪贴板

代码:

clear;clc;clf;

%假设待分析的函数是z=t^3N=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;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(瞬时频率

小结：傅里叶变换不能得到瞬时频率，即不能得到某个时刻的频率值。

Hilbert变换是求取瞬时频率的方法，但如果只用Hilbert变换求出来的

瞬时频率也不准确。（出现负频，实际上负频没有意义！）

用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的方差贡献率

%定义：方差为平方的均值减去均值的平方

%均值的平方

%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-1

mse(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-1

disp([imf,int2str(i] ;disp([mse(i m