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

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代码测试 参考 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工具箱 参考 这个真的相当强大！ 原始信号如下，有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

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

4. 图片背景的改善

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

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)

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