前言
在希尔伯特变换中我们已经知道，直接对"解析信号"做3瞬属性的分析时，"瞬时频率"会出现很多负频率，是无意义的！因此为解决这个问题并对非稳态信号做时频分析，才有了"希尔伯特-黄变换"！
思路：在对原始信号进行希尔伯特变换之前，要先把信号分解成瞬时频率具有意义的各个分量/函数！把信号进行分解的方法，在希尔伯特-黄变换理论中称为"经验模态分解"；分解出来的各个分量，其实是一类函数，称为"固有模态函数"，利用这类函数的局部特性，可以使得函数在任意一点的瞬时频率都有意义。
求瞬时频率(无意义)
经验模态分解，得到一类固有模态函数
对固有模态函数做希尔伯特变换，求瞬时频率(有意义)
所以，对希尔伯特-黄变换做简要总结(其实一点都不复杂)：
目标：求原信号的"有物理意义"的瞬时频率；
步骤：1. 经验模态分解；2. 希尔伯特变换/谱分析。
下面对2大步进行详细的说明，并给出对应的matlab手写程序。
经验模态分解
-
为什么希尔伯特-黄变换那么好？
个人认为它的关键就在于"经验模态分解"这一步。因为这种信号的分解，完全是"方法适应于数据"的，即不同的数据，方法会根据数据的情况做变换和调整，而不是始终用固定的那几种基函数(小波变换、短时傅里叶变换)。
这种自适应思想下的工具，一定是非常优越的！ -
什么是固有模态函数？
经验模态分解的目标就是将原始信号分解成一类固有模态函数，
固有模态函数其实只是满足下面两个条件的函数而已，没啥复杂的：
- 函数极值点的数目与零点的数目必须相等，或最多相差1；
- 极大值和极小值形成的包络，包络的均值为0；
其中均值为0的条件对应离散数据不好实现，一般平均值趋近于0（一般和0做差<0.1）即可。
-
经验模态分解何时结束？
经验模态分解就是为了分解出一系列的固有模态函数，所以当剩余数据已经无法再分出来固有模态函数时，经验模态分解这一步就彻底结束了。
如何判断剩余数据无法再分解了？
剩余数据极大值或极小值点数目有一个为0时(均值条件没法满足了)，或剩余数据已经是单调时(没法做包络线了)，就可以结束了。
分解的流程
设原始信号为，分解步骤如下：
- 步1：提取输入信号的所有极大值和极小值点，并用三次样条曲线连接出上、下包络线，设包络线均值为
，用输入信号减去均值可得
：
-
步2：判断步1得到的
-
步3：设得到的第k个固有模态函数为
，将其赋值给
：
将从剩余原始序列中分离出来，得到新的剩余项：
- 步4：判断新的剩余项是否满足经验模态分解的结束条件，不满足将剩余项带回到步1；满足则结束经验模态分解。
分解完毕后，原始信号可以用n个固有模态函数 + 1个剩余项的形式表示：
下面给出各步的matlab程序：
功能函数：
% 非主函数: 求极值点个数
% 注意: 若要极大值点数直接输入x，若要极小值则输入-x
function n = peaks(x)
[value,local] = findpeaks(x);
n = length(local); % 只要极值点个数
% 非主函数: 获得数据的极大值、极小值包络线(三次样条插值后的数据)
function s = getspline(x,t)
N = length(x);
[value,local] = findpeaks(x);
% 原始插值点: 极值点
t1 = t(local);
x1 = x(local);
s = spline(t1,x1,t);
% 非主函数: 判断当前函数是否为imf —— 零点个数与极值点个数对比
function u = isimf(x)
N = length(x);
u1 = sum( x(1:N-1).*x(2:N) < 0 ); % 零点个数
u2 = peaks(x) + peaks(-x); % 极大值 + 极小值总点数
if abs(u2-u1) > 1
u = 0; % 不是imf函数
else
u = 1; % 是imf函数
end
% 非主函数: 剩余信号的单调性判断(是否整个程序结束)
function u = ismono(x)
u1 = peaks(x)*peaks(-x);
if u1 ~= 0
u = 0; % 非单调: 不存在极大值点或极小值点
else
u = 1; % 单调：有一个极值不是0，就不是单调的
end
主调用函数：固有模态函数用二维矩阵记录
% 主函数: 经验模态分解
clear; clc;
x = xlsread('shuju.xlsx');
x = x(1001:1001+1023)'; % 有效数据长必须是2^n，所以我取1024，最后10几个点是0
N = length(x);
fs = 100; % 采样频率 = 1/采样间隔
t = (0:N-1)/fs; % 时间刻度
imf = [];
num = 1; % 计数器
xx = x; % 代替原始数据被操作而已
while ~ismono(xx)
x1 = xx;
sd = inf; % 用sd代替包络均值为0(<0.01)的条件
% 当sd不满足<0.01，或x1不是imf函数时，进入while内层循环
while (sd > 0.01) || ~isimf(x1)
s1 = getspline(x1,t); % 极大值包络插值结果
s2 = -getspline(-x1,t); % 极小值包络插值结果
x2 = x1 - (s1+s2)/2;
sd = sum( (x1-x2).^2 )/sum(x1.^2);
x1 = x2;
end
imf(num,:) = x1;
xx = xx - x1;
num = num + 1;
end
imf(num,:) = xx; % 最后一个余项也放进去
% 经验模态分解出的函数(最后一个是余项)
figure(1);
merge = 0;
[row,col] = size(imf);
for n = 1:row
% 把经验分解的内容再加起来
merge = merge + imf(n,:);
% 每一个画图
subplot(row,1,n)
plot(t,imf(n,:));
axis([min(t) max(t) -inf inf]);
end
suptitle('经验模态分解函数');
经验模态分解效果：
图1：经验模态分解得到的固有模态函数
希尔伯特变换/谱分析
对经验模态分解得到的固有模态函数依次进行希尔伯特变换，并求它们的3瞬属性即可。注意：一般这里就不要第一步保留的剩余项了。
每一个固有模态函数可以表示为：
我们把采样点时间、瞬时频率、瞬时振幅放在3维空间里：
注意：最后面一项是根据欧拉公式得到的。
可以看出：时间、瞬时频率、瞬时振幅3者之间是有联系的！
此时我们就可以做时频谱分析了！也就是得到"希尔伯特谱"。
为了方便起见，做谱分析这一步就用世界通用的两个函数来直接完成！一个是hhtspectrum.m和toimage.m函数。
emd经验模态分解效果：
图2：emd分解的希尔伯特-换变换时频谱
ceemdan分解效果：
图2：ceemdan分解的希尔伯特-黄变换时频谱
说明：希尔伯特-黄变换最初的理论是采用emd的经验模态分解，目前已经改进到采用ceemdan的模态分解方式，可以看到改进的方法精度进一步提高了！
后记
参考文献如下：
- 屈梁生, 张西宁, 沈玉娣. 机械故障诊断理论与方法[M]. 西安交通大学出版社, 2009.
- 薛雅娟. 地震信号时频分析及其在储层含气性检测中的应用研究[D]. 成都理工大学, 2014.
参考博客：希尔伯特-黄变换理论介绍。
本文所有程序下载地址：希尔伯特-黄变换。
网友评论