MATLAB希尔伯特黄变换HHT

原创已于 2025-06-19 16:52:26 修改 · 2.8w 阅读

311 ·

CC 4.0 BY-SA版权

文章标签：

#EMD #HHT

于 2019-07-23 20:19:16 首次发布

这两天在学习希尔伯特黄变换，也就是HHT，趁着学习的劲赶紧整理整理，用的是MATLAB进行编程，所用到的工具箱便是EMD工具箱，链接如下，请自行下载。

希尔伯特黄变换HHT_HHT-电信代码类资源-优快云下载

EMD工具箱安装过程：MATLAB中添加emd文件夹的路径 → 进入emd文件夹 → 命令窗口输入 install_emd 并执行

如果在仿真的过程中，出现“错误使用 instfreq”警告，下载下方函数，并保存至emd文件中即可。
https://download.youkuaiyun.com/download/Arry_W/91070508 用于解决下述报错：

错误使用 instfreq
需要的 TFD 应为 非负。

首先先介绍下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中，有这么一个关系， $x=\sum imf+rn$ ,而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谱') ;