希尔伯特谱、边际谱、包络谱、瞬时频率/幅值/相位—

上一篇文章对希尔伯特-黄变换（HHT）的前世今生进行了介绍。

不过在研究中通常并不是到希尔伯特-黄变换就停止了。

而是要用到诸如希尔伯特谱、包络谱、边际谱、瞬时频率/幅值/相位等方法进一步分析。

这些方法究竟是什么含义，以及怎样使用、怎样实现呢？

一、希尔伯特谱（Hilbert Spectrum）

希尔伯特谱是希尔伯特-黄变换得到的最直观结果，（上一篇文章中倒数第二张图——经过HHT的“疯狂的M”），其反映的是信号时间、瞬时频率和幅值之间的关系。该图谱可以用于分析包含混合分量信号中各分量随时间变化的规律，以识别局部特征。需要注意的是，希尔伯特谱有时是将emd分解后所有imf分量作为分析对象的，而有些会有针对性地挑选出某个或某几个imf分量分析，具体怎样操作是要结合具体研究内容有针对性地选择的。

下边使用一个模拟的轴承故障信号，使用希尔伯特谱寻找缺陷[1]。

模拟在轴承外圈引入一个缺陷，该缺陷导致一系列冲击，导致磨损逐渐加剧。这些冲击将在轴承的外环故障频率(BPFO)处重复出现。

模拟出的轴承故障振动特征信号如下图：

在4~6秒之间的尖峰是由于轴承振动激发的共振引起的[1]，而从5~10秒之间，信号中叠加了幅值逐渐增加的冲击信号。

对该信号进行EMD分解，其中的部分分量如下图。

对IMF1做希尔伯特谱，得到下图结果。随着轴承磨损的逐渐加剧，冲击能量增加，在2500~3500Hz范围内希尔伯特谱幅值随之增大。

对IMF3做希尔伯特谱，得到下图结果。这里将轴承振动激发的共振信号分离出来，其频率范围在0~100Hz之内。

需要再次强调的是，希尔伯特谱是一种时频谱，可以与之类比的是连续小波变换、短时傅里叶变换这种同样是时频分析的方法的谱图。这种谱反应的是信号频率成分随着时间的变化，是做非平稳信号（例如例子中的故障信号）的重要手段。使用这种类型的分析方法强调的就是“变化”，即特征在时间尺度上的改变——因为如果信号没有随时间发生变化，使用频域分析手段就够了。正因如此，HHT的方法在在生物医学（如血压变化）、地球物理（如地震、海浪分析）、工程领域（故障诊断等）的非平稳信号为主要研究对象的领域广泛应用。

二、边际谱（Marginal Spectrum）

边际谱是建立在已经算出希尔伯特谱的基础上的，计算方法是将希尔伯特谱在时间轴上进行积分，使之从幅值-时间-频率三者间的关系转变为幅值-频率两者间的关系，描述的是幅值（或能量）在频率轴上的分布。公式为：

现在问题来了：同样是幅频谱，傅里叶谱和边际谱有什么区别呢？

在傅里叶谱中，在某一频率上存在着能量意味着具有该频率的正弦或余弦波存在于信号的整个持续时间内；
而在边际谱中，在某一频率上存在着能量意味着具有该频率的波在信号的整个持续时间内某一时刻出现的可能性较高。
因此在一定程度上，Hilbert边际谱具有一定的概率意义，Hilbert边际谱可以看作是一种加权的联合幅值-频率-时间分布，而赋予每个时间-频率单元的权重即为局部幅值，从而在Hilbert边际谱中，在某一频率上存在着能量就意味着具有该频率的振动存在的可能性，而该振动出现的具体时刻在Hilbert谱中给出。[2]

简单来说，傅里叶谱和边际谱有一定的相关性，但是在处理非平稳信号时，更适合使用边际谱，因为“傅里叶变换为了在数学上拟合原始数据的非平稳波形，不得不引入大量高频的'伪'谐波分量，这会导致傅里叶谱对低频能量的低估[2]”。

下图为上边故障轴承信号的傅里叶谱和边际谱：

三、包络谱（Envelope Spectrum）

包络谱最常见的应用场景就是机械产品故障诊断（尤其是轴承）。不过需要区别于边际谱的是，包络谱不是基于希尔伯特谱。包络谱的求法是：目标信号→希尔伯特变换→得到解析信号→求解析信号的模→得到包络信号→傅里叶变换→得到Hilbert包络谱。

包络谱是一种解调方法，在某些轴承故障中（例如表面损伤）会在轴承运行中激发出一些列周期性冲击信号，这些信号会与高频固有振动发生调制。包络谱分析能够有效地将这种低频冲击信号进行解调提取。

不过包络分析结果往往会收到低频噪声的影响，在分析前需要进行带通滤波以消除噪声干扰，而带通滤波器的参数常常难以选择，因为事先并不知道共振频带的范围。此时就可以结合EMD的方法，将原始信号进行EMD分解，选出包含共振频带的前几个IMF分量重构信号，再进行包络分析[3]。

需要注意的是，包络谱与频谱结果差异较大，包络谱更适用于做故障特征提取。

还是用轴承故障来举例说明，仿真生成一段轴承外环故障频率（83.33Hz）的轴承外圈故障信号[7]：

对这段信号直接做频谱图，只能观察到25Hz的转频频率分量。

使用包络谱则可以提取到该故障信号：

四、瞬时频率/幅值/相位

这里说的瞬时频率、瞬时幅值和瞬时相位这三个概念都是基于希尔伯特变换构造，相关概念不清楚的可以再看一下这里。

这些方法在地质[4]、电气[5]、生物医学[6]等领域的应用，不过这些不是本人的研究方向，且概念相对直观不难理解，就不做展开了。

另外瞬时频率的定义和估计方法不止这一种，想要了解更多的同学可以看这里：信号瞬时频率的估计方法及其应用。

五、MATLAB代码实现

上文提到的希尔伯特谱、边际谱、包络谱以及各种瞬时量的MATLAB编程实现方法在网上可以搜到一些，不过大多比较杂乱，注释不清，有些还存在错误（甚至官方库中都有错）。同时还存在MATLAB自带函数与第三方库共存，难以选择取舍等情况，对于新手使用不太友好。

博主重新将上述分析方法进行了编程实现，使之兼具了MATLAB自带函数和第三方库函数，兼容了各个版本的MATLAB软件，可以快速配置，快速上手。获取连接如下：

希尔伯特谱、边际谱、包络谱、瞬时频率/幅值/相位画图 | 工具箱文档

上述功能已全部封好，直接导入数据就能使用。

希尔伯特谱：

function hhtSpec(data,fs,type)
% 绘制信号希尔伯特谱
% 输入：
% data： 待分析信号
% fs：   采样频率，当fs未知时，可以将fs设置为空（[]），并将type设置为2。此时hht图纵轴将经过标准化
% type： 当type的值为1时，则优先采用MATLAB自带库中的hht函数进行画图，采用的绘图风格也与自带hht函数保持一致，此时hht函数的用法与MATLAB自带函数一致；
%        当type的值为2时，则强制采用第三方库文件中的希尔伯特谱函数进行画图，不过可能会存在内存不够无法顺利画图的可能
%        如果type没有传入参数，则该函数内将type设置为1

边际谱：

function [mgS,f] = marginalSpec(imf,fs,type)
% 求边际谱并画图
% 输入：
% imf：  imf分量，注意方向：imf是每行一个信号分量的矩阵
% fs：   采样频率
% type： 当type的值为1时，则优先采用MATLAB自带库中的函数进行运算；
%        当type的值为2时，则强制采用第三方库文件中的函数进行运算，不过可能会存在内存不够无法顺利画图的可能（数组超过预设的最大数组大小）
%        如果type没有传入参数，则该函数内将type设置为1
%        经测试，使用两种库函数做出的边际谱存在幅值差异，频率点基本一致
% 输出：
% mgS：  边际谱幅值。
% f：    边际谱的频率轴。
%
% 典型用例：
% imf = kEMD(y);
% marginalSpec(imf,fs);

包络谱：

function [envS,f,xEnv] = envSpec(data,fs)
% 求信号包络谱并画图
% 输入：
% data： 待分析信号
% fs：   采样频率
% 输出：
% envS： 包络谱数值
% f：    包络谱频率轴
% xEnv： 包络线

瞬时频率、瞬时相位和瞬时幅值：

[insF,insP,insA] = InsFPA(yBad,fs);
% 求信号的瞬时频率、瞬时相位和瞬时幅值
% 输入：
% data： 目标信号
% fs：   目标信号的采样频率
% 输出：
% insF： 瞬时频率
% insP： 瞬时相位
% insA： 瞬时幅值

参考资料：

[1] https://ww2.mathworks.cn/help/signal/ref/hht.html#mw_025b9602-e0a8-4b16-9224-21b3b95eb156

[2] 张郁山. 希尔伯特—黄变换(HHT)与地震动时程的希尔伯特谱——方法与应用研究[D]. 中国地震局地球物理研究所, 2003.

[3] 朱可恒. 滚动轴承振动信号特征提取及诊断方法研究[J]. 2013.

[4] 张恩厚, 陈前新. 瞬时振幅,瞬时相位,瞬时频率在地质雷达工程中的应用[J]. 安徽地质, 2004, 014(001):58-60.

[5]韩晓慧, 杜松怀, 苏娟, et al. 基于局部均值分解的触电故障信号瞬时参数提取[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2015, 31(17):221-227.

[6] 韩晓慧, 杜松怀, 苏娟, et al. 基于局部均值分解的触电故障信号瞬时参数提取[J]. 农业工程学报, 2015, 31(017):221-227.

[7] Envelope spectrum for machinery diagnosis