【转】Matlab Hilbert-Huang 变换分析总结

最新推荐文章于 2024-07-29 10:13:29 发布

信号处理学渣

最新推荐文章于 2024-07-29 10:13:29 发布

阅读量1.2w

点赞数 13

MATLAB C/C++ 专栏收录该内容

15 篇文章 3 订阅

订阅专栏

http://blog.sina.com.cn/s/blog_84024a4a01019pfw.html

标签： emdhilbert时频分析imf希尔伯特-黄分类：信号处理
这是对几篇参考文章的整理和总结，参考文章在后面会给出链接。
一、Hilbert变换测试
参考http://blog.sina.com.cn/s/blog_84024a4a01019opg.html
1. hilbert函数
matlab中，由hilbert函数得到的信号是合成的复信号，其虚部才是书上定义的Hilbert变换。这一点在基本概念一文中有说明。在上面的参考博文中，有这样几句代码：
y = sin(2*pi*f*t);
yh = hilbert(y); % matlab函数得到信号是合成的复信号
yi = imag(yh); % 虚部为书上定义的Hilbert变换

2. hilbert函数可以提取包络和调制信号频率

这里主要看第五张图和第七张图
调制信号的产生如下
ts = 0.001;
fs = 1/ts;
N = 200;
k = 0:N-1;
t = k*ts;
% 原始信号
f1 = 10;
f2 = 30;

f3 = 200;
a = 2 + cos(2*pi*f1*t);
m = sin(2*pi*f2*t); % 调制信号

x3 = sin(2*pi*f3*t);
y = a.*m; % 信号调制

（1）包络的提取
第五张图中蓝色的线就是通过hilbert变换提取的包络，关键语句如下（y代表的是调制结果信号，也就是我们一般情况下所测得的信号）：
yh = hilbert(y);
aabs = abs(yh); % 包络的绝对值
plot(t, aabs)

（2）高频信号的提取
这一点主要基于两个信号乘积的hilbert变换取决于高频信号的hilbert变换。
看第七张图，就是调制信号的瞬时频率，主要代码如下：
yh = hilbert(y);
aangle = unwrap(angle(yh)); % 包络的相位
af = diff(aangle)/2/pi; % 包络的瞬时频率，差分代替微分计算
plot(t(1:end-1), af*fs)

这部分代码主要有两点需要注意
（i）unwrap
网上有如下解释：
对于复数函数，写为指数形式A(t)*exp(j*B(t))，程序中处理的是些离散点，由于用angle求得的相角不是B(t)，而是C=rem(B(t),2*Pi)（限制在-Pi到Pi之间），unwrap就是根据C求B（会差一个常数），依据是数据足够密，所以相邻两点之间的相位变化不超过Pi。
（ii）af是归一化的结果，转换到实际频率要乘以fs
从结果上看，主要的高频信号为70Hz，与原信号相符
（3）有三个信号的情况
这个是我后来又增加的200Hz的调制信号，我们只看一下瞬时频率和包络，发现与前面的结论一致。
ts = 0.001;
fs = 1/ts;
N = 200;
k = 0:N-1;
t = k*ts;
% 原始信号
f1 = 10;
f2 = 30;

f3 = 200;
a = 2 + cos(2*pi*f1*t); % 包络2
m = sin(2*pi*f2*t); % 调制信号
x3 = sin(2*pi*f3*t);
y = a.*m.*x3; % 信号调制

二、Matlab 官网提供的EMD代码测试
参考http://blog.sina.com.cn/s/blog_84024a4a01019opp.html
1. IMF代表固有模态函数
2. 第四次分解信号的Hilbert分析

（1）IMF4是分解出来的第四个固有模式函数
（2）IMF4的包络
yh = hilbert(y);
yenvelope = abs(yh); % 包络

就是前面说的，做hilbert变换后求幅值。
（3）瞬时频率
yf = diff(yangle)/2/pi/Ts; % 瞬时频率

这个在前面也说过
（4）调制信号
yModulate = y./yenvelope;

这个也很好理解，信号结果 = 包络 * 调制信号，变换一下就是上面的公式。
三、G-Rilling EMD工具箱
参考http://blog.sina.com.cn/s/blog_84024a4a01019oqe.html
这个真的相当强大！
原始信号如下，有2个频率，10Hz和100Hz
fs = 1000;
ts = 1 / fs;
t = 0: ts: 0.3;
z = 2*sin(2*pi*10*t) + 5.*sin(2*pi*100*t);

1. EMD分解得到固有模式函数
imf = emd(z);

2. 绘制EMD分解图
emd_visu(z, t, imf)

z是原始信号，t是时间序列，imf就是emd函数的结果
解释一下产生图片中c2f和f2c，c2f代表由coarse到fine，就是由最后分解的结果不断向上叠加生成原始信号，f2c正好相反，是由第一分解的结果向下叠加。
3. 绘制时频图
[A, f, tt] = hhspectrum(imf);
[im, tt] = toimage(A, f);
disp_hhs(im)

其结果的坐标都是序列或归一化后的，想要得到实际时间和频率，进行简单变换即可

t = 横坐标 * ts

f = 纵坐标 * fs