经验模态分解与边际谱：非平稳信号时频分析核心原理与工程实践

1. 项目概述

1.1 核心需求解析

先说结论：这个项目做的是"非平稳信号处理"，核心工具是经验模态分解（Empirical Mode Decomposition，简称EMD）和希尔伯特-黄变换（Hilbert-Huang Transform，简称HHT），最终目标是算出一张边际谱（Marginal Spectrum），用来分析信号在不同频率上的能量分布。

如果你手头有一堆振动数据、声音数据、心电信号，或者任何"看起来乱糟糟、频率随时间在变"的波形，传统的傅里叶变换往往只能告诉你"有哪些频率成分"，却说不出"这些频率是什么时候出现的"。EMD+HHT的组合拳解决的就是这个问题——先把复杂信号拆成若干个本征模态函数（IMF），再对每个IMF做希尔伯特变换得到瞬时频率和瞬时幅值，最后把所有瞬时信息汇总成边际谱，得到"哪个频率的能量贡献最大"。

我最初做这个项目，是因为当时的实验数据里有大量非平稳冲击成分，用短时傅里叶变换（STFT）处理后总感觉时频分辨率不够用，频率窗口怎么调都别扭。换成EMD之后，整个分析流程一下子顺了很多——虽然EMD本身也有不少坑，但总体而言，这套方法在机械故障诊断、地震信号分析、生物医学信号处理、海洋波浪能分析等领域都非常值得掌握。

这个项目适合谁？正在做信号处理相关课题的学生、需要分析振动或声音数据的工程师、想从时频分析角度换一种思路的科研人员。不需要你有多深的数学基础，但最好懂一点傅里叶变换的基本概念，知道什么是频率、什么是频谱，这样理解起来会快很多。

1.2 方案选型思路

其实在最开始，摆在我面前的不止EMD一个选择。当时我同时评估了小波变换（Wavelet Transform）、STFT和EMD三套方案，各自跑了一轮测试数据，最后才锁定EMD路线。

小波变换的问题是：母小波函数一旦选定了，相当于给信号强加了一个先验的"波形模板"。如果信号本身的形态和你选的母小波差了十万八千里，分析结果就会有偏差。而EMD不需要预先设定任何基函数，它是完全数据自适应的——信号长什么样，分解出来的IMF就长什么样。这一点在处理未知形态的实测信号时是巨大优势。

STFT的问题则在于时频分辨率的矛盾：窗口短，时间分辨率高但频率分辨率差；窗口长，频率分辨率高但时间分辨率差。海森堡测不准原理卡在那里，怎么调都有遗憾。EMD+HHT给出的瞬时频率是逐点计算的，从原理上绕开了这个矛盾。

当然，EMD也有自己的缺点，比如端点效应、模态混叠、对噪声敏感等，后面我会专门写一节怎么处理这些问题。但总体而言，对于非平稳、非线性信号，EMD+HHT确实是一个久经考验的经典框架，也是很多现代改进算法（如EEMD、CEEMDAN、VMD）的起点。

2. 核心原理拆解：从EMD到边际谱的逻辑链条

2.1 EMD到底在做什么

在讲EMD之前，先铺垫一个直觉：把一个复杂信号想象成"几层简单信号叠在一起"。比如一段音乐里，低音鼓的节奏是低频成分，人声的中频成分，镲片的沙沙声是高频成分，它们叠加在一起就是你能听到的完整音乐。EMD做的事，就是把这个混合物一层一层剥开，每一层都是一个"本征模态函数"（IMF）。

一个信号成为IMF需要满足两个条件：第一，在整个信号长度上，极值点（极大值和极小值）的数量与过零点的数量相差不超过一个；第二，在任意时间点上，由局部极大值定义的上包络和由局部极小值定义的下包络的均值必须为零。通俗说，IMF必须是"上下对称、局部波形相对单一"的分量。

EMD通过一种叫"筛选"（sifting）的迭代过程来提取IMF。具体做法是：先用三次样条插值把信号的局部极大值点连成上包络，把局部极小值点连成下包络，取两者的均值作为一条"中心线"，从原始信号里减掉这条均值线，得到第一次筛选结果；然后对这个结果反复重复同样操作，直到满足IMF条件为止，这样就得到了第一个IMF。接着用原始信号减去第一个IMF得到残余分量，对残余分量继续做同样的筛选，得到第二个IMF……一直筛到残余分量变成一个单调函数或足够小的值，分解结束。

这个"包络均值"的思想是整个EMD的灵魂。它不依赖任何全局基函数，而是根据信号的局部极值自适应地构造分解方向，所以EMD特别擅长处理那些"频率成分随着时间在漂移"的信号。

2.2 希尔伯特变换与瞬时频率

拿到IMF之后，下一步是计算每个IMF的瞬时频率。这里用到希尔伯特变换：对任意实信号x(t)，构造一个复信号z(t) = x(t) + j·H[x(t)]，其中H[x(t)]是x(t)的希尔伯特变换。z(t)的幅值就是IMF的瞬时幅值，z(t)相位的导数就是瞬时频率。

关键点是：希尔伯特变换只有在信号是"窄带"信号时，给出的瞬时频率才有物理意义。所谓窄带，就是信号的频谱只集中在一个较窄的频带内。IMF恰好大致满足这个条件——这也是为什么要先做EMD而不是直接对原始信号做希尔伯特变换的原因。如果直接对原始复杂信号做希尔伯特变换，算出来的瞬时频率会因为多个频率成分纠缠在一起而变得毫无意义。

2.3 边际谱的定义与物理含义

每个IMF经过希尔伯特变换后，都能得到一条瞬时频率随时间变化的曲线和一条瞬时幅值随时间变化的曲线。把所有这些IMF的信息都算出来之后，就可以构造希尔伯特谱H(ω, t)——它是频率、时间、幅值的三维谱图，描述的是"在时间t、频率ω处，信号的瞬时幅值有多大"。

边际谱h(ω)则是对希尔伯特谱在时间维度上的积分，数学表达式为：

h(ω) = ∫ H(ω, t) dt

它在工程上非常好解释：在观察的整个时间段内，频率ω对应的累计幅值能量是多少。它很像傅里叶频谱，但是两者有本质区别：傅里叶频谱里的"频率"是指一个在整个时间轴上恒定存在的正弦分量，而边际谱里的"频率"可以是随时变化的瞬时频率，它只是把所有时刻的瞬时能量按频率归位了而已。因此，边际谱在分析频率随时间漂移的信号时，常常能给出比傅里叶谱更锐利、更清晰的频率峰。

这里我举一个直观的例子：假设有一个信号，前半段时间是10Hz的正弦波，后半段是20Hz的正弦波。傅里叶谱会分别显示出10Hz和20Hz两个峰，但它看不出"10Hz只存在前半段"这个信息；边际谱同样会给出10Hz和20Hz两个能量峰，但除此之外，结合希尔伯特谱的时频图，你能清楚看到这两个频率成分各自出现的时间区间。这就是边际谱加希尔伯特谱组合分析的价值所在。

3. 实操过程：从零开始算边际谱

3.1 准备数据与环境

做这个项目我用的语言是Python，环境是基于Anaconda搭建的虚拟环境，核心依赖库如下：

如果你是第一次接触，建议先装PyEMD，用现成的库跑通整个流程，之后再考虑自己实现EMD来深入理解原理。我用的是PyEMD库里的EMD和Visualisation模块，计算核心是C语言加速的，速度还可以，几万点的数据几秒钟就能出结果。

测试数据用了一段模拟的非平稳信号，生成方式如下：

import numpy as np fs = 1000 # 采样率 1000Hz t = np.linspace(0, 1, fs, endpoint=False) # 信号1：10Hz正弦波，只存在于0到0.5秒 x1 = np.sin(2 * np.pi * 10 * t) * (t < 0.5) # 信号2：25Hz正弦波，只存在于0.5到1秒 x2 = 0.8 * np.sin(2 * np.pi * 25 * t) * (t >= 0.5) # 信号3：高斯白噪声，幅值较小 noise = 0.1 * np.random.randn(len(t)) x = x1 + x2 + noise

这个信号的频谱结构很简单清晰：前半段只有10Hz，后半段只有25Hz，加上一点点噪声。用这个信号来检验EMD和边际谱算法，任何偏差都能一眼看出来。

3.2 EMD分解步骤详解

在实际计算时，我用的封装函数是PyEMD里的标准EMD。关键代码如下：

from PyEMD import EMD emd = EMD() IMFs = emd(x)

如果你对内部原理感兴趣，其实核心步骤完全可以手动复现，流程是：

1. 找出信号的全部局部极大值点和局部极小值点。
2. 用三次样条函数分别对极大值点、极小值点插值，得到上包络线和下包络线。
3. 计算上下包络的均值m(t)。
4. 用原始信号x(t)减去均值m(t)，得到候选IMF h(t)。
5. 检查h(t)是否满足IMF条件。如果不满足，把h(t)当作新信号，回到步骤1继续筛选。
6. 满足条件后，把h(t)作为第i个IMF，用原始信号减去它得到残余分量，对残余分量再次执行步骤1到5。
7. 重复直到残余分量单调或足够小，分解结束。

这里有两个细节值得注意。

第一个是停止条件的设置。如果筛选次数太多，IMF会退化为等幅调频信号，所有物理信息被抹平；如果筛选次数太少，IMF又不够规范。我在实现时采用了较经典的做法：设定筛选次数上限为100次，同时检查候选信号与原始信号的方差比，当变化足够小时提前终止。

第二个是三次样条插值的端点处理。默认情况下，样条插值在信号的两端容易出现大幅度的摆动，这会引起著名的"端点效应"——分解出来的IMF首尾两端会发生畸变。常用的缓解方法是信号延拓，比如在两端分别镜像一段数据，或者用特征波延拓、多项式外推等。我在项目中采用的是镜像延拓，效果中规中矩，对于一般分析场景够用。

分解完成后，检查一下得到的IMF个数和残余分量。以我的模拟信号来说，正常情况下应该分解出两到三个有效IMF，再加一条残余趋势项。其中第一个IMF主要包含25Hz成分和一部分噪声，第二个IMF主要包含10Hz成分，残余分量接近零线。

3.3 希尔伯特谱与边际谱的计算实现

EMD分解完成后，对每个IMF做希尔伯特变换。这里需要注意：如果你的某个IMF不是严格的窄带信号，直接做希尔伯特变换得到的瞬时频率会有较大的抖动，甚至出现负频率——这在物理上不合理。常见的处理方式有两种：一是先对IMF做平滑滤波再把瞬时频率算出来；二是直接使用PyEMD的希尔伯特谱接口，它内部做了瞬时频率提取和去野值处理。

实际的边际谱计算代码大致如下：

from PyEMD import Hilbert # 创建希尔伯特-黄变换对象 hht = Hilbert() # 计算瞬时频率、瞬时幅值 inst_freq, inst_amp = hht.transform(IMFs, fs) # 构建频率轴，映射瞬时频率和幅值到边际谱 from scipy.signal import welch freq_resolution = 0.5 # 频率分辨率，单位Hz max_freq = fs / 2 freq_bins = np.arange(0, max_freq, freq_resolution) marginal_spectrum = np.zeros(len(freq_bins)) for i, (freq, amp) in enumerate(zip(inst_freq, inst_amp)): for j in range(len(freq)): if not np.isnan(freq[j]) and amp[j] > 0: # 找最近的频点 idx = int(freq[j] / freq_resolution) if 0 <= idx < len(freq_bins): marginal_spectrum[idx] += amp[j]

这段代码本质上就是把每个时间点的瞬时幅值按瞬时频率累加到对应的频率仓上，积分变成了累加。如果你的数据量比较大，这个逐样本循环会慢一些，可以用numpy的histogram函数一次性完成累加：

import numpy as np hist, _ = np.histogram(np.concatenate(inst_freq), bins=freq_bins, weights=np.concatenate(inst_amp)) marginal_spectrum = hist

两种写法结果等价，后者在数据量大时快很多。

做边际谱的时候，还有一个参数值得反复调：频率分辨率freq_resolution。分辨率设得越细，谱峰越尖锐，但噪声也容易混进来。通常以信号总时长的倒数作为基本分辨率下限，比如1秒的信号，分辨率为1Hz；如果你想看得更细，可以放宽到0.5Hz，但低于0.1Hz意义就不大了，因为底层算法从时间积分来看也很难区分更细的频率差异。

3.4 参数设置与效果观察

跑完上面流程后，我对比了傅里叶频谱和边际谱。直观感受是：在10Hz和25Hz附近，傅里叶谱和边际谱都有清晰的能量峰；但从峰形的锐度来看，边际谱的峰要明显窄一些、背景噪声更低。原因前面说过——边际谱的能量归位是按瞬时频率逐点统计的，噪声在频率轴上被均匀摊开了，而真实成分的瞬时频率集中在目标频率附近，所以峰更突出。

如果改用STFT来做同样的分析，窗口长度必须提前拍板。窗口太短，25Hz和10Hz在频率轴上都糊成一团；窗口太长，时间分辨力降低了，分不清前半段和后半段的区别。EMD+HHT把"时间和频率分辨率不可兼得"的约束松开了，这种体验在对比过后会非常明显。

不过我还是要泼一盆冷水：边际谱虽然好用，但它的幅值不是严格的"能量"概念，它本质上是瞬时幅值的累积，单位取决于你原始信号和瞬时幅值定义，并不适合直接和傅里叶频谱的幅值做绝对数值比较。做工程分析时，大家通常只看它的相对大小和峰位信息。

4. 常见问题与排查技巧实录

4.1 端点效应导致首尾畸变

做完第一次EMD分解后，我画出来的IMF首尾两端经常出现明显的"翘尾巴"现象——尤其是第一个IMF，两端会出现比较大的摆动，甚至让整个波形看起来不合理。

排查思路分两步。第一步，确认是不是原始数据本身首尾就存在异常——如果信号本身在两端有突变，EMD会把这种突变当成一个极值点来处理，从而产生虚假的IMF。第二步才是确认延拓方式是否有效，我推荐的延拓方式是"镜像延拓"加"特征波延拓"结合。前者实现简单，把端点附近的波形对称翻转出去，适合平稳段和缓变信号；后者思路更复杂一点，把端点附近的波形模式往外延拓，适合有明显节奏性的信号。

如果说延拓都做了，两端还有畸变，那还有一个非常常见的解决办法：分解前先把信号两端各裁剪掉一部分，分析完再把对应区间的时间轴平移回来。代价是丢失一点点端点附近的信息，但换来的是更可靠的IMF分解结果。实测中我用这个方法处理过一段采集时有直流偏置漂移的数据，效果立竿见影。

4.2 模态混叠问题

模态混叠是EMD最出名的一个坑。简单说，就是不同频率成分没有被干净地分离开，一个IMF里混着多个频率分量，或者在同一个IMF序列里突然跳变到另一个频率的被污染状态。

典型场景是：信号里同时存在一个持续的小幅高频成分和一段短暂的大幅低频冲击，EMD很容易把冲击时段的高频能量全部分配到低频IMF里，导致高频IMF在冲击点附近剧烈变形。

解决模态混叠的办法，最经典的是EEMD（Ensemble EMD，集合经验模态分解）。思路非常直观：给原始信号反复加上不同的白噪声，每次做一次EMD分解，最后把所有分解结果按IMF序号平均。加噪声的目的是让不同尺度的信号被"均匀标记"到对应的IMF中去，平均后噪声就相互抵消了，只剩下信号的真实成分。

EEMD代码用PyEMD里的EEMD模块即可：

from PyEMD import EEMD eemd = EEMD() IMFs = eemd(x)

不过EEMD也有代价，一个是计算量成倍增加，另一个是需要调两个参数：集合次数和噪声幅值。我的经验是集合次数取100到200次，噪声幅值取原始信号标准差的0.1到0.3倍，基本能得到稳定结果。如果噪声幅值太小，模态混叠改善不明显；太大又会让分解结果过平滑，丢失真实细节。

4.3 瞬时频率出现负值或剧烈跳变

瞬时频率出现负值在数学上是可以的，但物理上通常意味着提取失败。这类问题大多源于IMF不满足窄带条件，尤其在信号频率成分靠得很近时。

我遇到过一次特别离谱的情况：模拟信号里有两个频率只差1Hz的成分，分解出的某个IMF瞬时频率在目标频率附近来回跳变，幅度足有5Hz以上。后来我把眼光放回IMF本身，发现该IMF在波形上有一个明显的包络调制——幅值忽大忽小。对这种信号，直接希尔伯特变换做瞬时频率必然出问题。

改进办法有两个：第一，对IMF做包络解调，先从IMF提取包络信号，再用包络信号去归一化原始IMF，得到"调频+调幅"解耦后的纯调频信号，再算瞬时频率；第二，用平滑窗去瞬时频率曲线上做中值滤波或移动平均，滤掉野值，保留趋势项。前者物理意义更严谨，后者实现更快，看你的需求。

4.4 边际谱高频段出现"假峰"

有段时间我的边际谱在很高频段（接近奈奎斯特频率）总是冒出一个莫名其妙的峰，源信号里根本没有这个频率成分。排查了半天，最后锁定在噪声上。

EMD分解时，噪声成分会优先进入前几个IMF，高频噪声的瞬时幅值可能特别大，累加到边际谱高频段就形成了假峰。解决办法有两种：一是分解前对原始信号做低通滤波，把明显高于目标频段的噪声干掉；二是对算出来的瞬时幅值做阈值过滤，幅值特别小的样本直接忽略。

实际经验是两者结合效果好：分解前滤波可以显著抑制高频假峰，阈值过滤则可以压掉残余的微弱噪声贡献。阈值取原始信号幅值标准差的0.5倍比较合适，太小没啥用，太大又会把弱小的真实成分误杀。

5. 工具与库选型心得

5.1 现有开源库横向对比

如果你不想从零造轮子，工程复现时直接选用开源库是最快路径。我试过几个主流库，各自的脾气摸得比较清：

用PyEMD入门非常快，三行代码就能跑通分解和边际谱计算。但是要注意版本兼容问题，我之前在某Python版本上升级依赖时碰到过接口不一致的情况，建议把环境固定下来，避免"昨天还能跑，今天报错"的尴尬。

5.2 从库走向自实现的时机

开源库虽然快，但如果项目要求你改算法的核心机制，比如修改筛选过程、定制包络构造，那就必须自己实现EMD了。我自己花了大概一周时间写了一个基于numpy和scipy的简化EMD，核心代码不到200行，但写完之后对EMD的理解绝对上了一个台阶。

自实现时有一个性能陷阱必须提：三次样条插值如果对每一个筛选取包络均值都全局重建一次，耗时是很大的。我的优化方式是预先算出所有极值点和包络采样点，然后在每次筛选中只对局部区间做更新，避免重复全局插值，速度大概能快三到四倍。如果你的数据是百万点级的，这种优化几乎是必须的。

6. 项目扩展与应用场景思考

6.1 机械故障诊断中的边际谱应用

这个项目做完之后，我第一时间想到的落地场景是机械故障诊断。滚动轴承出现局部损伤时，振动信号里会产生周期性冲击成分，这些冲击在时间上很窄、在频域上覆盖很宽。传统的傅里叶频谱只告诉你有高频能量和低频能量，但是冲击的周期性特征不明显。EMD+HHT可以逐圈追踪瞬时频率的变化，边际谱上对应故障特征频率处的峰会格外突出，故障识别准确率明显高于只用频谱做阈值判断的方案。

实际工程中还有一个额外优势：机械设备的转速通常不是恒定不变的，升降速过程中的振动信号是典型的非平稳信号，用STFT或传统频谱分析很难处理，EMD的自适应性就能派上用场。如果你想做设备健康监测，EMD+HHT几乎是可以直接套用的框架。

6.2 生物医学信号分析

心电信号（ECG）里的P波、QRS波群、T波各自占据不同的时频区域，RR间期变化携带重要的生理信息。EMD可以把ECG分解成若干IMF，其中特定IMF的瞬时频率曲线能反映心率变异性。脑电信号（EEG）处理也是一个热门方向，不同频段的节律波在EMD分解后更容易分离。

不过生物信号有个大前提：噪声和伪迹特别多，基线漂移、肌肉电噪声、电极接触不良产生的尖峰，都会让EMD产出一大堆没意义的IMF。我的实践经验是：先用带通滤波器把目标频段之外的成分滤掉，再做EMD分解，能大幅减少模态混叠和虚假IMF。别偷懒，预处理做得好，后面的一切都会轻松很多。

6.3 数据驱动的自动化分析

进一步扩展，可以把EMD+HHT包装成一个自动化分析管线。比如输入一段振动信号，自动完成去噪、EMD分解、瞬时频率提取、边际谱计算，然后根据边际谱上特征峰的位置和幅度自动判断设备状态。这里需要处理一个细节：不同贵细数据之间的IMF数量可能不一样，自动对齐时要设一个能量占比阈值，只保留累计能量达到85%以上的前几个IMF，丢弃后续能量贡献极低的尾巴。

另外一个比较好的思路是用EEMD替代EMD作为管线的基础模块，虽然计算量大了，但模态混叠问题能缓解不少，自动化场景下不用人工干预去检查每个IMF的分离质量。

6.4 小波包与VMD等替代方法的对比

现在VMD（变分模态分解）已经被越来越多的人使用，它本质上把模态分解问题转化为变分优化问题，在数学上更漂亮，模态混叠的抑制能力也更强，但代价是必须提前设定模态个数。模态个数设错了，结果会非常糟糕。EMD不需要设定模态个数，用起来更"懒人友好"。

我个人的体会是：如果信号很干净、频率成分分离明显，直接用EMD就很省心；如果信号复杂、频率靠得近，可以先试EEMD；如果你的场景允许人工调参且对精度要求苛刻，再考虑VMD。没有万能的方法，关键是理解自己信号的特点。

7. 实操总结与个人经验

最后说几点我在整个项目过程中积累的实在经验。

第一，数据预处理比算法本身更重要。我反复验证过，一个干净的带通滤波，效果远好于任何花哨的算法优化。滤波虽然会引入滤波器的相位响应，但对于边际谱分析来说，影响通常可以接受，因为边际谱关注的是能量分布，不是精确的波形重建。

第二，分解结果的检查别只看IMF曲线形态，一定要看瞬时频率曲线。瞬时频率是判断IMF分离质量的硬指标，如果瞬时频率在目标频率附近波动超过百分之十几，说明这个IMF不干净，要么换EEMD，要么重新调整预处理参数。

第三，边际谱的能量峰位置对频率分辨率参数并不敏感，但峰幅值对噪声很敏感。做定量比较时，尽量在同一套参数下跑所有数据，不然峰高差异可能来自参数而不是信号本身。

第四，工程上时间紧迫时，别沉迷于自实现EMD。先用PyEMD跑通整个流程，确认边际谱对信号特征确实有效，再决定要不要深入研究算法内部。项目周期里最贵的其实是调试时间，不是那几行代码。

这个项目做完后，我确实对"频率到底是什么"这个问题有了更深的理解。傅里叶说频率是永恒的，但现实世界里几乎没有永恒的信号。EMD和边际谱提供了一种更尊重信号真实性的分析语言——允许频率随时间变化，然后在这个基础上重新定义频谱。用起来会有点挑剔，但一旦摸清它的脾气，它会成为你手里的利器。