HHT仿真实例

本人针对www.ilovematlab.cn论坛中楼主cwjy的原创“【原创】用希尔伯特黄变换（HHT）求时频谱和边际谱”（见http://www.ilovematlab.cn/thread-65516-1-1.html），进行了验证，结果如下。

注：原始程序有所改动！

给定频率分别为10Hz和35Hz的两个正弦信号相叠加的复合信号，采样频率fs=2048Hz的信号，表达式如下：y=5sin(2*pi*10*t)+5*sin(2*pi*35*t)。

具体Matlab程序如下：

% 用希尔伯特黄变换（HHT）求时频谱和边际谱 clear;clc; N=2048;

% fft默认计算的信号是从0开始的 t=linspace(1,2,N); delta=t(2)-t(1); % fs=1/deta;

x=5*sin(2*pi*10*t)+5*sin(2*pi*35*t); %%

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

% 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).*sin(w1*t)+(t>-200+N1*deta&t<=-200+N2*deta).* % sin(w2*t)+(t>-200+N2*deta&t<=200).*sin(w3*t); m=0:N-1;

f=1./(N*delta)*m;% 可以查看课本就是这样定义横坐标频率范围的 % 下面计算的Y就是x(t)的傅里叶变换数值 % Y=exp(i*4*pi*f).*fft(y)

% 将计算出来的频谱乘以exp(i*4*pi*f)得到频移后[-2,2]之间的频谱值 y=fft(x);

z=sqrt(y.*conj(y)); figure;

plot(f(1:100),z(1:100));grid; title('图2 幅频曲线');% xiangwei=angle(y); figure;

plot(f,xiangwei);grid; title('图3 相频曲线'); figure;

plot(t,x);grid; title('图1 原始信号'); % axis([-2,2,0,1.2]) %%

1

%---------------------------------------------- % （2）用Hilbert变换直接求该信号的瞬时频率 % clear;clc;clf;

% % 假设待分析的函数是z=t^3 % 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); hx=hilbert(x);

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;

plot(t(1:N-1),sp);grid; %

title('图4 瞬时频率'); %%

%----------------------------------------------

% （3）用HHT求取信号的时频谱与边际谱

c=emd(x); % 计算每个IMF分量及最后一个剩余分量residual与原始信号的相关性 [m,n]=size(c); for i=1:m;

a=corrcoef(c(i,:),x); % 各个IMF与原始信号的相关性 xg(i)=a(1,2); end xg;

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

2

mse(i)=mean(c(i,:).^2,2)-mean(c(i,:),2).^2; % 方差百分比，也就是方差贡献率 mseb(i)=mse(i)/mmse*100; % 显示各个IMF的方差和贡献率 end

% 画出每个IMF分量及最后一个剩余分量residual的图形

%--------------------------------------------------------- % figure; %

% for i=1:m-1

% disp(['imf',int2str(i)]);disp([mse(i) mseb(i)]); % end

% subplot(m+1,1,1);plot(t,x);ylabel(['signal','Amplitude']); %

% for i=1:m-1

% subplot(m+1,1,i+1);plot(t,c(i,:));ylabel(['imf',int2str(i)]); % % 线条设置为黑色plot(t,c(i,:),'k') % end

% subplot(m+1,1,m+1);plot(t,c(m,:));ylabel(['r',int2str(m-1)]); %----------------------------------------------------------- figure;

for i=1:m-1

disp(['imf',int2str(i)]); disp([mse(i) mseb(i)]);

subplot(m+1,1,i+1);plot(t,c(i,:));ylabel(['imf',int2str(i)]); end

subplot(m+1,1,1);plot(t,x);ylabel('原始信号','fontsize',9); % ylabel(['原始','信号']);

title('图5 EMD得到的各分量');

subplot(m+1,1,m+1);plot(t,c(m,:));ylabel(['r',int2str(m-1)]); %------------------------------------------------------------ % 画出每个IMF分量及剩余分量residual的幅频曲线 % for i=1:m-1 % end

y1(1,:)=fft(c(1,:));

z1(1,:)=sqrt(y1(1,:).*conj(y1(1,:))); figure;

plot(f,z1);grid;

title('图6 imf1幅频曲线'); figure;

plot(f(1:100),z1(1:100));grid; title('图7 imf1幅频曲线(1:100)');

% 以上两条语句截取坐标为1-100;相当于放大.

3