这两天在学习希尔伯特黄变换,也就是HHT,趁着学习的劲赶紧整理整理,用的是MATLAB进行编程,所用到的工具箱便是EMD工具箱,链接如下,请自行下载。
首先先介绍下matlab中经验模态分解所用到的函数emd。
imf = emd(x)
imf = emd(x,...,'option_name',option_value,...)
imf = emd(x,opts)
[imf,ort,nb_iterations] = emd(...)
做仿真的时候只是简单的对信号做一个模态分解,因此只用到了第一种方法,imf=emd(x),这其中若x为实向量,则计算经验模态,若x为复向量,则计算二维经验模式。对于x为实向量的情况下,得到的是1-n个imf(内部模态函数),在加上一个剩余项。
在得到IMF后,可以直接作图,也可以利用工具箱中的函数emd_visu作图,可以得到模态分解结果以及两种方向相反的重构信号过程,该函数的用法为
emd_visu(x,t,imf)
其中,x为原信号,t为时间序列,imf为上一步得到的结果,该函数可得三个图,第一个为emd图,后两个为信号重构过程,分别为f2c、c2f。
clear;close all;clc
N=1024;
fs=1024;
t=(0:N-1)/fs;
f1=70;f2=120;
x=cos(2*pi*f1*t)+2*cos(2*pi*f2*t);
imf=emd(x); %经验模态分解
emd_visu(x,t,imf) %作图
所得结果如下
第一个图所表示便是emd的过程了,7个imf加一个剩余项rn,后两张图就是原信号的重构过程了,我们知道在HHT中,有这么一个关系,,而c2f图和f2c图的区别就是,c2f从imf7往imf1方向重构原信号,而f2c则是从imf1往imf7方向进行的。
在emd结束后,便可以开始画HHT谱图了,此处我结合这两天查看的资料,给出两个所用到的函数用法的简单表述,如下:
[A,f,t] = hhspectrum(imf)
用于计算HHT谱,imf为前面所得结果,A为每个imf的瞬时振幅,f、tt分别为对应的瞬时频率和时间序列。
[im,tt,ff] = toimage(A,f,t,splx,sply)
其中A,f,t与hhspecturm函数的相关定义一致,splx、sply则分别对应输出im的列数和行数,splx一般与t的长度一致;
im为频谱的二维图像,tt为时间序列号(此处划重点,是“序列号”,而不是“序列”,也就是需要除以总长度得到时间序列),ff的注释解释为 “centers of the frequency bins”(频率门的中心值?)。
下面直接给出程序和结果。
[A,f,tt] = hhspectrum(imf(1:end-1,:));
[im,tt,ff] = toimage(A,f,tt,length(tt));
figure;
set(gcf,'Color','w');
imagesc(tt/N,[0,0.5*fs],im);
set(gca,'YDir','normal')
colormap('jet');colorbar;
xlabel('时间(s)');
ylabel('频率(Hz)');
axis([0 1 0 200]);title('HHT谱') ;