Hilbert EMD 分解总结（转载）_丶Young_

http://blog.sina.com.cn/u/3039219953

首页博文目录关于我

个人资料

微博

加好友发纸条

写留言加关注

博客等级：
博客积分：

博客访问：
关注人气：
获赠金笔：0支
赠出金笔：0支
荣誉徽章：

正文字体大小：大中小

Hilbert EMD 分解总结（转载）

(2017-03-18 11:05:16)

标签：

matlab

hilbert

emd

分类：算法

这是对几篇参考文章的整理和总结，参考文章在后面会给出链接。

一、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函数可以提取包络和调制信号频率

http://s3/middle/84024a4a4d95cc65bb2b2&690Hilbert-Huang 变换分析总结" TITLE="[转载]Matlab Hilbert-Huang 变换分析总结" ACTION-DATA="http://s3/middle/84024a4a4d95cc65bb2b2&690" ACTION-TYPE="show-slide" STYLE="margin: 0px; padding: 0px; list-style: none;" />

这里主要看第五张图和第七张图

调制信号的产生如下

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; % 信号调制

http://s9/middle/84024a4a4d95cc6632858&690Hilbert-Huang 变换分析总结" TITLE="[转载]Matlab Hilbert-Huang 变换分析总结" STYLE="margin: 0px; padding: 0px; list-style: none;" />

http://s10/middle/84024a4a4d95cc66aaf19&690Hilbert-Huang 变换分析总结" TITLE="[转载]Matlab Hilbert-Huang 变换分析总结" STYLE="margin: 0px; padding: 0px; list-style: none;" />

二、Matlab 官网提供的EMD代码测试

参考http://blog.sina.com.cn/s/blog_84024a4a01019opp.html

1. IMF代表固有模态函数

2. 第四次分解信号的Hilbert分析

http://s13/middle/84024a4a4d95cc672f5cc&690Hilbert-Huang 变换分析总结" TITLE="[转载]Matlab Hilbert-Huang 变换分析总结" STYLE="margin: 0px; padding: 0px; list-style: none;" />

（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

http://s8/middle/84024a4a4d95cc66e2797&690Hilbert-Huang 变换分析总结" TITLE="[转载]Matlab Hilbert-Huang 变换分析总结" ACTION-DATA="http://s8/middle/84024a4a4d95cc66e2797&690" ACTION-TYPE="show-slide" STYLE="margin: 0px; padding: 0px; list-style: none;" />

这个时频图有如下含义，横轴是时间，纵轴是频率，颜色代表幅值，可以看出，在0.1（即100Hz）和0.01（10Hz）和0.005（5Hz）有明显的谱线。其中5Hz为最后的残差频率。查看im矩阵可知幅值，主要有两个，一个是4.5左右，另一个是 0.8左右，幅值可以通过边缘谱查看。

4. 图片背景的改善

（这是另一个信号的分析，只有一个频率10Hz）

colormap(flipud(gray))

http://s2/middle/84024a4a4d95cc6701e01&690Hilbert-Huang 变换分析总结" TITLE="[转载]Matlab Hilbert-Huang 变换分析总结" ACTION-DATA="http://s2/middle/84024a4a4d95cc6701e01&690" ACTION-TYPE="show-slide" STYLE="margin: 0px; padding: 0px; list-style: none;" />

5. 自己绘制时频图

imagesc(tt*ts, [0, fs/2], im)

ylabel('frequency/Hz' )

set(gca, 'YDir' , 'normal' );

xlabel('time/s' )

title('Hilbert-Huang spectrum' )

6. 边缘谱的绘制

边缘谱其实就是对幅值进行时间上的积分。

edge_trum = mean(im, 2);

ff = (1: size(im, 1)) / size(im, 1) * fs/2;

figure, plot(ff, edge_trum)

http://s15/middle/84024a4a4d95cc673258e&690Hilbert-Huang 变换分析总结" TITLE="[转载]Matlab Hilbert-Huang 变换分析总结" ACTION-DATA="http://s15/middle/84024a4a4d95cc673258e&690" ACTION-TYPE="show-slide" STYLE="margin: 0px; padding: 0px; list-style: none;" />

可以看到跟傅里叶频谱很类似。有两个频率，一个是10Hz，一个100Hz。

阅读┊ 收藏 ┊ 喜欢 ▼ ┊打印┊举报/Report

前一篇：彻头彻尾理解傅里叶变换算法（很优秀的文章，转载）

后一篇：word 分两栏后如何加脚注

新浪BLOG意见反馈留言板　欢迎批评指正