希尔伯特黄变换HHT原理与Python实战:从EMD分解到瞬时频率分析
2026/9/7 7:46:39 网站建设 项目流程

简介:希尔伯特黄变换(HHT)是分析非线性非平稳信号的重要方法,这份代码包面向信号处理工程师、科研人员与学生,基于MATLAB与C语言实现经验模态分解、希尔伯特变换与瞬时频率提取,可用于生物医学信号、地震数据、机械故障诊断等场景。包内共72个文件、约94KB,包含40个MATLAB函数脚本、17个C源文件及11个头文件,另有shell配置脚本与mat数据文件,兼顾易用性与计算效率,附示例目录便于快速对照运行。已有5071人学习下载。除完整的EMD分解、希尔伯特谱计算与可视化工具外,还提供安装/卸载脚本与修复脚本,降低配置门槛;C版本可处理较大数据量或性能敏感任务,通过阅读源码可理解筛分停止条件、包络构造、模态混叠处理等关键细节,适合算法研究与工程落地参考。

1. 从一次振动信号处理说起:HHT解决的正是FFT最头疼的问题

我第一次用希尔伯特黄变换(HHT)处理旋转机械振动信号时,最大的震撼不是它算得多快,而是它把一段频谱图里看不出先后关系的频率成分,清清楚楚地展现在了时间轴上。当时手头的数据是一段包含冲击和变频成分的轴承振动信号,用FFT做出来只有几根峰,峰虽然好看,但我根本说不清这些峰对应的故障特征出现在哪个时刻。直到换成HHT,EMD把信号一层层拆开,瞬时频率随时间的变化像看心电图一样直观,很多原本要凭经验猜的信息,终于有了可解释的依据。

简单说,HHT是一种针对非线性、非平稳信号的自适应时频分析方法,由黄锷在1998年前后正式提出。它不依赖预先选好的基函数,而是通过经验模态分解(EMD)把复杂信号拆成若干个本征模态函数(IMF),再对每个IMF做希尔伯特变换,求出带有物理意义的瞬时频率。机械故障诊断、地震波分析、生物医学信号处理、金融时序分析,凡是信号特征会随时间和工况变化的领域,几乎都能看到它的身影。这篇我会从原理讲到代码,再把实战里最常踩的坑单独拎出来说,希望能帮你在自己的数据上少走弯路。

1.1 非平稳信号的“时变”本质

所谓非平稳,最直白的理解就是:信号的统计特性随时间变化。你去看一段平稳信号,比如50Hz的正弦波,它的均值、方差、频率基本恒定;但真实工程信号很少这么乖,轴承磨损到一定程度会出现周期性冲击,转速变化会导致频率持续爬升,噪声背景的强度也可能忽大忽小。

这类信号最大的麻烦在于,频率本身成了一个“变量”。上一秒振动的主要频率可能是100Hz,下一秒就可能跳到180Hz。我们需要的不是“这段信号里有哪些频率”,而是“每个频率在什么时间出现、强度怎么变化”。HHT的设计目标,就是回答后面这个问题。

1.2 FFT和STFT为什么总是差点意思

傅里叶变换是最经典的频谱分析工具,但它把整段信号展开到一组固定频率的正弦波上,得到的是一个全局平均结果。一个频率成分只要在信号里出现过哪怕很短一段时间,FFT都会“记上一笔”,但它不会告诉你这笔账记在哪个时间节点。

短时傅里叶变换(STFT)试图补救,做法是把信号切成小段再做FFT。问题在于窗长是固定的:窗太短,频率分辨率差,频率相近的成分分不开;窗太长,时间定位就模糊,突变细节被抹掉。这就是海森堡测不准原则在时频分析里的体现,你不可能同时获得无限精确的时间分辨率和频率分辨率。小波变换做了尺度自适应,但需要提前选择母小波,换成不同小波,结果往往差不少。本质上,它还是“拿模板去套数据”。

1.3 HHT的自适应思路:先分解再算瞬时频率

HHT和上面这些方法的根本区别,在于它不预设任何基函数。EMD会从信号自身的极值分布出发,按时间尺度从小到大逐层剥离出不同频率成分,得到一个IMF序列和一个趋势残差。每个IMF理论上都是单分量信号,可以直接做希尔伯特变换,从而得到瞬时幅值和瞬时频率。

这个“先分解、再求瞬时频率”的流程,等于把“频率随时间变化”这件事从隐式变成了显式。时频平面的分辨率不再受窗函数约束,而是由数据本身的局部变化决定。这也是为什么HHT特别适合处理突变、跳频和非线性效应明显的信号。不过,自适应是它的优势,也是它的软肋,分解过程带有较强的经验性,用不好容易出问题。后面我会重点讲这些坑。

2. EMD分解:把信号拆成一堆IMF的完整逻辑

EMD是整个HHT的基石。它做的事情可以概括成一句话:把一个复杂信号表示成若干个本征模态函数IMF加上一个残差。每个IMF对应一种“固有振动模式”,分解过程完全由信号本身的极值点时间尺度驱动,不需要任何先验知识。

2.1 IMF的两个硬性条件

一个函数要想被称为IMF,必须同时满足两个条件:

  • 在整个数据段内,极值点的个数与过零点的个数必须相等,或最多相差一个。这保证IMF基本是围绕零轴上下振荡的“窄带”分量。
  • 在任意时刻,由局部极大值确定的上包络与由局部极小值确定的下包络的均值为零。换句话说,IMF要局部对称,不能有明显的“偏置”。

这两条本质上是在保证信号“局部单分量”。只有满足这两个条件,后续做希尔伯特变换得到的瞬时频率才有物理意义。如果不满足,算出来的瞬时频率会出现负值或者剧烈跳动,根本没法解释。

2.2 筛分流程到底在做什么

EMD的筛分过程,我用大白话串一遍:

  1. 找出信号的全部局部极大值点和局部极小值点。
  2. 用三次样条插值分别拟合出上包络线和下包络线。
  3. 计算上下包络的平均值m1,相当于信号当前的中心漂移。
  4. 原始信号x(t)减去m1,得到候选分量h1。
  5. 检查h1是否满足IMF的两个条件。如果不满足,就把h1当作新的原始信号,重复1到4步。
  6. 一旦h1满足条件,就记为第一个IMF,记为c1。
  7. 原始信号减去c1,得到残差r1;把r1当成新信号继续筛分,直到残差变成单调曲线或幅值小于设定阈值。

为什么这里用三次样条而不是线性插值?因为瞬时频率需要对相位求导,如果包络线不光滑,求导结果会充满毛刺。三次样条能保证包络和信号曲面足够平滑,这是工程实现里很关键的一处细节。

2.3 为什么说EMD是“数据驱动”的分解

FFT把信号硬塞给一组三角函数,小波需要人工挑小波基,而EMD是从数据里“长”出分解结果来。极值点之间的时间间隔,决定了这一层信号尺度;信号本身有多快,分解就从哪个尺度开始。转速高的成分先被剥离,然后一层层往下筛,最后剩下的残差就是趋势项。

这种自适应能力让EMD在处理非线性调频、间歇性冲击和随机噪声混合的信号时,常常比传统方法更贴近物理过程。但代价是计算量明显更大,而且分解结果对极值点分布、筛选次数和端点处理方式都很敏感。所以EMD不是“输入信号、输出真理”的傻瓜工具,而是一个需要使用者理解参数含义的分析框架。

3. 希尔伯特谱与边际谱:从IMF到瞬时频率的临门一脚

EMD只完成了信号分解,真正让HHT“封神”的是后面这一步:对每个IMF做希尔伯特变换,得到随时间变化的瞬时频率和瞬时幅值。这两者组合在一起,就构成了希尔伯特谱。

3.1 解析信号与瞬时频率的定义

对任意一个IMF,记作c(t),先做希尔伯特变换:

[ \hat{c}(t) = \frac{1}{\pi} \text{PV} \int_{-\infty}^{\infty} \frac{c(\tau)}{t-\tau} d\tau ]

然后构造解析信号:

[ z(t) = c(t) + j\hat{c}(t) = a(t)e^{j\theta(t)} ]

其中瞬时幅值 (a(t)=\sqrt{c^2+\hat{c}^2}),瞬时相位 (\theta(t)=\arctan(\hat{c}/c))。瞬时频率定义为:

[ f(t) = \frac{1}{2\pi}\frac{d\theta(t)}{dt} ]

这个定义看起来简单,但前提是信号必须是单分量窄带信号。这也正是EMD必须先存在的理由:只有IMF满足单分量条件,相位求导出来的频率才不会出现无意义的大幅波动。如果拿原始复杂信号直接做希尔伯特变换,瞬时频率往往会乱成一团,甚至出现负数。

3.2 希尔伯特谱和边际谱怎么看

把所有IMF的瞬时频率、瞬时幅值画在同一张时间-频率图上,横轴是时间,纵轴是频率,颜色或亮度代表瞬时幅值,就得到了希尔伯特谱 (H(\omega,t))。它展示的是信号能量在时频平面上的分布,能非常直观地看出频率随时间移动的轨迹。

如果把希尔伯特谱沿时间轴积分,就得到边际谱:

[ h(\omega) = \int H(\omega,t)dt ]

边际谱的物理意义是:信号中每个频率成分在整个观测时间段里累计贡献了多少幅值或能量。和FFT谱相比,边际谱不需要信号满足平稳性假设,对频率成分的定位通常更尖锐,也更接近真实的物理成分。很多论文里用边际谱做故障特征提取,就是因为同一个故障频率在边际谱上往往比在FFT谱里更好区分。

3.3 一个标准算例:啁啾信号在两种频谱下的差距

我用一个线性调频信号举个例子。信号频率从10Hz线性爬到50Hz,持续1秒。用FFT做频谱分析,会看到一个从10到50Hz的“平台”,你可能知道信号包含了这段频率,但看不出频率是爬升的、下降的还是随机跳变的。而在HHT谱上,你看到的是一条从10Hz斜着升到50Hz的曲线,时间关系一目了然。

这就是HHT在非平稳信号上的核心价值:比FFT多了时间维度,比STFT少了窗长的纠结。边际谱则可以进一步告诉你,哪个频率在整个观察期内贡献的能量最大,这对定位设备故障特征频率特别有用。

4. 用Python完整跑一遍HHT:代码、参数和结果解读

理论说得再漂亮,不上手试一次都是空的。下面我用一段模拟信号演示完整流程。模拟信号包含10Hz低频正弦、一个频率从30Hz扫到80Hz的线性调频分量,以及0.5秒后出现的100Hz间歇振动,再加少量白噪声。

4.1 环境准备与模拟信号构造

Python环境里最常用的库是PyEMD,安装命令如下:

pip install EMD-signal

这里有个容易踩的小坑:安装包名字叫EMD-signal,但import的时候用的是PyEMD。新版PyEMD的接口和旧版略有差异,建议安装完成后先跑一下版本和基本调用确认无误。

构造模拟信号:

import numpy as np from PyEMD import EMD import matplotlib.pyplot as plt fs = 1000 t = np.linspace(0, 1, fs) f0 = 10 S = (np.sin(2 * np.pi * f0 * t) + np.sin(2 * np.pi * (30 * t + 25 * t**2)) + np.where(t >= 0.5, 0.5 * np.sin(2 * np.pi * 100 * t), 0) + 0.05 * np.random.randn(len(t)))

第二条分量是线性调频,相位是 (2\pi(30t+25t^2)),瞬时频率就是对内层求导得到的 (30+50t),正好从30Hz升到80Hz。最后加了一点白噪声,模拟真实采集环境。

4.2 EMD分解的代码实现与IMF解释

调用EMD分解:

emd = EMD() IMF = emd.emd(S, t) plt.figure(figsize=(10, 8)) for i, imf in enumerate(IMF): plt.subplot(len(IMF), 1, i + 1) plt.plot(t, imf) plt.ylabel(f'IMF{i + 1}') plt.xlabel('时间 (s)') plt.tight_layout() plt.show()

跑出来的IMF序列通常第一个是高频间歇成分,对应100Hz那一段;中间某个IMF对应临时频率从30到80Hz的扫频分量;后面某个IMF对应10Hz低频;最后一条接近趋势或极低频残差。

这里不推荐机械地把每个IMF和物理信号一一对应,因为当存在模态混叠时,一个物理成分可能被拆到多个IMF里。但模拟信号里各分量频率差异足够大,分解结果会比较干净,正好适合用来理解每个IMF的物理含义。

4.3 绘制希尔伯特谱的代码与读图方法

自定义一个快速求希尔伯特谱的片段:

from scipy.signal import hilbert plt.figure(figsize=(10, 6)) for imf in IMF: analytic = hilbert(imf) amp = np.abs(analytic) phase = np.unwrap(np.angle(analytic)) freq = np.diff(phase) / (2 * np.pi * np.diff(t)) times = (t[:-1] + t[1:]) / 2 valid = freq > 0 plt.scatter(times[valid], freq[valid], c=amp[:-1][valid], s=1, cmap='jet', vmin=0, vmax=np.max(amp)) plt.xlabel('时间 (s)') plt.ylabel('瞬时频率 (Hz)') plt.ylim(0, 150) plt.colorbar(label='瞬时幅值') plt.show()

读图时重点看三条“亮线”:一条水平在10Hz附近,一条从30Hz爬升到80Hz,一条在0.5秒后出现在100Hz附近。颜色越亮,表示该时刻该频率的能量越强。如果改用FFT看这幅图对应的原始信号,很难把这些特征分开,但HHT谱可以直接讲出信号变化的“故事”。

5. 实战中最容易翻车的三个细节:端点效应、模态混叠和停止准则

HHT用起来最痛苦的从来不是原理,而是那些藏在细节里的坑。我实际用过的项目里,翻车最频繁的就是下面三个问题。

5.1 端点效应:筛分时信号两端在“乱飘”

三次样条包络需要极值点作为支撑,但信号端点处哪边都没有额外极值点,包络会不受约束地向外发散。这会导致一个结果:同一个信号,只截取不同长度的段落,IMF在端点附近往往会变形,瞬时频率也会在首尾处出现明显的大幅摆动。

我常用的缓解办法有三个:

  • 端点镜像延拓:把端点附近的数据镜像翻转一段,人为补几个极值点,让包络在端点处更自然。
  • 极值延拓:根据端部附近极值点的时间间隔和幅值趋势,向外预测延拓。
  • 舍弃影响区:如果只关心中间段,就直接丢弃两端各一小段,尤其是做瞬时频率分析时。

注意,延拓方式本身就会影响分解结果,所以不要迷信某一种延拓是“绝对正确”的。可靠做法是人工检查IMF在端部的表现,如果端部畸变明显,要意识到这部分结果是不可信的。

5.2 模态混叠:物理单一频率被拆到多个IMF里

模态混叠是EMD最被诟病的问题。最常见的情况是:一个本来应该集中在单个IMF里的振动模式,被分散到了多个IMF中;或者一个IMF里同时混入了时间尺度差异很大的成分。典型诱因包括间歇性高频信号、脉冲冲击和强噪声。

缓解方案里最成熟的是集合经验模态分解(EEMD)和它的改进版CEEMDAN。思路是给原始信号加入有限幅值的白噪声,多次分解后对结果取平均,利用白噪声的统计均匀性把混叠的模态“拉”开来。噪声幅值一般取信号标准差的0.1到0.4倍,集合次数视计算资源而定,几百次是常见选择。

但EEMD不是万能的。白噪声加太多会淹没真实弱信号,加太少又起不到抑制混叠的作用;集合次数太少则平均效果差、结果不稳定。用之前最好先拿模拟信号标定一组参数,再应用到真实数据上。

5.3 筛分停止准则:多次迭代不是越少越好

EMD筛分是迭代过程,迭代到什么程度停止,直接影响IMF质量。如果筛分次数太少,IMF可能不满足条件,瞬时频率会出现负值;如果筛分次数太多,IMF会被过度平滑,变成近似常幅值的纯调频信号,丢失真实的幅值调制信息。

经典做法是通过相邻两次筛分结果的标准化差SD来判断停止,比如:

[ SD = \frac{\sum_{k=0}^{T} |h_{i-1}(t) - h_i(t)|^2}{\sum_{k=0}^{T} h_{i-1}^2(t)} ]

当SD落在0.2到0.3之间时停止,是很多论文沿用的一种经验准则。实际使用中我还会额外设置最大迭代次数,防止程序在某些病态数据上无限循环。重要的是,不要盲目照抄参数,要多观察IMF的波形和瞬时频率曲线,确认结果在物理上说得通。

6. HHT适合做什么、不适合做什么:选型思路与能力边界

每次有人问我HHT好不好用,我都会说:它是一个很有价值的时频分析工具,但不是万能的。理解了它的边界,你才知道什么时候该用它。

6.1 典型适用场景

HHT在机械故障诊断里尤其好用。滚动轴承、齿轮箱出现局部故障时,会产生周期性冲击和频率调制成分,这些信号天然是非线性的,FFT经常被背景噪声和频率混叠搞得很头疼,HHT反而能比较清晰地把瞬时频率变化刻画出来。

地震信号处理也是HHT的重要阵地,地震波到达前后频率成分剧烈变化,需要精确到时频定位。生物医学领域里,脑电、心电、肌电信号普遍非平稳,经常用HHT做特征提取。金融时间序列分析也偶有使用,不过金融数据的信噪比太低,用得更多是EMD做预处理,而不是直接上希尔伯特谱。

6.2 与FFT、小波变换的选型对比

方法时间分辨率基函数依赖适合信号主要短板
FFT固定三角函数平稳信号没有时间定位能力
STFT受窗长限制固定窗函数缓变信号时频分辨率无法兼得
小波变换较好需要选择母小波部分非平稳信号结果依赖基函数选择
HHT完全自适应非线性、非平稳信号经验性强、有端点与混叠问题

选型逻辑很简单:信号平稳就用FFT,效率高、解释直接;信号缓变且需要一定时间信息就考虑STFT;小波适合已有明确先验并愿意调试基函数的场景;当你面对的是复杂变频、冲击、调制信号,同时又不想预设基函数时,HHT的优势就体现出来了。

6.3 HHT的局限到底在哪里

HHT最大的局限是“经验性”。EMD缺少像傅里叶变换那样的严格数学基础,很多步骤依赖启发式准则,不同参数设置可能得到不同结果。端点效应和模态混叠很难彻底消除,只能缓解。计算开销也明显大于FFT,处理长时间高频采样数据时会比较慢。

所以我个人建议,HHT更适合作为“分析手段”而不是“全自动黑盒”。用它之前,先构造一段和真实信号特征相近的模拟信号,把分解质量、瞬时频率曲线、谱图形态都确认没大问题,再上真实数据。按照这个流程做下来,HHT通常不会让你失望。

最后再分享一个我自己的习惯:每次拿到新数据,我会先用HHT和FFT各跑一遍,对照着看。FFT给出全局频率骨架,HHT补上时间演变细节,两张图互相印证,远比单靠一种方法就能把信号“看透”要可靠。你如果刚开始接触HHT,也建议从模拟信号开始,把IMF和瞬时频率的对应关系摸熟了,再挑战真实工程数据。

本文还有配套的精品资源,点击获取

需要专业的网站建设服务?

联系我们获取免费的网站建设咨询和方案报价,让我们帮助您实现业务目标

立即咨询