经验模态分解(EMD)原理与MATLAB实战：从IMF到希尔伯特谱

搞信号分析的人，总绕不开一个名字：经验模态分解（EMD）。我最早接触它是为了处理一段振动信号，FFT做完之后只能在频谱上看到几个模糊的峰，完全看不出故障冲击发生的时刻，那种束手无策的感觉到现在还记得。后来在MATLAB里把EMD跑通，看到一条条固有模态函数（IMF）从复杂波形里被拆出来，才真正体会到什么叫“自适应时频分析”。

这篇文章我不打算把EMD写成数学论文，而是站在一个实际用过它的人的角度，把原理、MATLAB代码、参数调节、踩坑经验一次讲透。如果你是做机械故障诊断、生物医学信号分析、地震数据处理、气象时间序列预测，或者只是想把非平稳信号里的成分看清楚，这篇文章都值得花时间读完。

1. EMD是什么，为什么要折腾它

1.1 传统时频分析工具的短板

过去分析信号，用得最多的是傅里叶变换。它把一个时间信号分解成许多固定频率的正弦波之和，听起来很美，但实际问题往往不配合。真实世界的信号几乎都是非平稳的：说话时声带振动频率会漂移，机械轴承损坏时冲击响应是瞬时的，脑电波的节律是动态变化的。傅里叶变换只能告诉你“整体上包含哪些频率”，却无法告诉你“某个频率在什么时间出现”。于是人们又造出了短时傅里叶变换和小波变换，但它们说到底都依赖一个预先选好的窗函数或基函数。

窗选窄了，频率分辨率差；窗选宽了，时间分辨率差。这就是著名的海森堡测不准原理在信号分析中的体现。更麻烦的是，小波变换需要你提前选定小波基，而真实信号里的模态往往是某个频率周围动态变化的模式，很难用一个固定的基函数去完美匹配。所以当黄锷在1998年提出经验模态分解时，很多做信号分析的人像找到了新大陆：它不需要选基函数，不需要设窗长，完全根据信号自身的局部特征来分解。

1.2 EMD拆信号的基本思想

EMD的核心思想特别朴素：任何复杂的非平稳信号，都可以看成若干个“固有模态函数”（IMF）叠加一个残余趋势项。所谓固有模态，你可以把它理解成一种局部对称的振荡模式，它的瞬时频率在整个时间轴上有物理意义。EMD的分解过程不需要任何先验知识，完全依赖信号自身的极值点分布，因此它被称为“自适应信号分解方法”。

打个比方：如果一个家庭里大人、孩子、宠物同时在活动，你要分析谁在什么时候做了什么。傅里叶变换就像是听了一整天之后统计“全家平均每小时走动次数”，小波变换像是根据你预想的“大人步伐频率”去筛选，而EMD更像是把录音逐帧拆开，根据每个人动作的高频、低频特征把各自的节奏自然分离出来。EMD输出的每一个IMF，携带的信号能量在时间尺度上逐渐降低：第一个IMF通常频率最高、包含细节最多，后面的IMF频率逐级降低，最后的残余项就是整体趋势。

1.3 典型应用场景和前提条件

我实际用过EMD的场景包括：轴承故障振动信号的特征提取、心电信号中的去噪和基线漂移校正、风速序列的分解预测、以及桥梁结构健康监测。生理信号领域用得尤其多，因为EEG、ECG、呼吸信号往往同时包含多种生理节律、噪声和伪迹，EMD能把它逐层剥开。

需要提醒的是，EMD并不是万能的。它对信号长度有一定要求，理论上至少得包含两个完整振荡周期，工程上建议采样点数至少要几百个点。另外它要求信号是实信号，复数信号需要特殊处理。如果信号靠一两个极值点硬撑，分解出的IMF可能是假模态，这一点新手很容易忽视。

1.4 学习EMD需要准备什么

使用MATLAB跑EMD，建议使用2018a或更高版本，因为从那个版本开始信号处理工具箱内置了emd函数，还有配套的hht函数可以直接画希尔伯特谱。如果你用的是老版本，或者没有安装信号处理工具箱，也没关系，后面我会给出一个手写教学版EMD，在纯MATLAB环境下就能运行。学习EMD不需要特别高深的数学底子，懂微积分基本概念、知道什么是插值和局部极值，就足够理解整条链路。

2. EMD算法的关键细节，不读懂这些后面全是坑

2.1 IMF的两个硬性判据

一个函数能不能被认定为固有模态函数，需要同时满足两个条件。

第一，在整个数据长度内，极值点（局部极大值加局部极小值）的个数必须与过零点的个数相等或最多相差一个。这条保证了IMF的每个振荡都有一个相对完整的“波峰-波谷”结构，频率不会忽快忽慢到失去意义。

第二，在任意一点上，由局部极大值确定的上包络线和由局部极小值确定的下包络线的均值必须为零。这里的“零”是工程意义上的近似，实际算法用阈值来判断。这条保证了IMF的波形关于时间轴局部对称，不会出现明显的直流漂移或不对称震荡。

之所以强调“局部”而不是整个信号，是因为EMD支持时变频率和时变振幅，每个IMF可以是一个调幅调频信号。它不像傅里叶分量那样必须是一个标准的正弦波，所以两个判据都是局部的。我在刚开始学的时候总想套用全局均值为零的理解，结果怎么也想不通，后来才明白这就是EMD自适应的精髓所在。

2.2 筛分流程的逐步拆解

EMD的分解过程有一个专门的名字：筛分（sifting）。整个流程可以拆成五个步骤。

第一步，找到当前信号的所有局部极大值和局部极小值点。第二步，用三次样条插值分别对极大值点和极小值点做包络，得到上包络线和下包络线。第三步，计算上下包络线的均值作为局部均值线。第四步，用原信号减去局部均值线，得到一个候选模态。第五步，检查候选模态是否满足IMF的两个判据，如果不满足，就把候选模态当作新的信号，重复前面四个步骤继续筛分；如果满足，就把它作为当前IMF提取出来，然后用原信号减去这个IMF，对剩余信号继续重复整个流程，直到剩余信号单调或小于预设阈值为止。

这里最容易让新手困惑的是“包络线”。我建议直接在MATLAB里把上下包络线和均值线画出来看一次，比读十遍文字描述都管用。上下包络线必须把信号夹在中间，均值线反映了信号局部的起伏趋势，每次减去均值线，就是把一个振荡周期里不对称的偏移剥掉，让波形逐步变得局部对称。这个过程就像把一根弯曲的铁丝一次次压直，但每次只矫正宏观弯曲，保留更细节的振荡。

2.3 停止准则和几个重要参数

筛分不能无限做下去，否则会把每个IMF都削成一个纯正弦波，反而失去物理意义。经典EMD中黄锷等人提出了标准差准则，用相邻两次筛分结果的相对差异来判断：[ SD = \sum_{t} \frac{|h_{k-1}(t) - h_k(t)|^2}{h_{k-1}^2(t)} ]当SD小于某个阈值（比如0.1到0.3之间），就认为筛分已经收敛，可以停止。实际操作中我习惯把阈值设在0.2，既能保证波形对称，又不至于过度筛分。

MATLAB内置emd函数里还有一些衍生参数值得关注。SiftMaxNumIterations控制找到每个IMF前最多进行多少次筛分，默认值好像是100，如果你发现分解出的IMF波形过于平缓或者丢失细节，可以试着把这个值调小。MaxNumIMF限制总共提取IMF的个数，实际信号分解出的IMF数量不会超过它。还有一个能量比阈值EnergyRatio，它控制剩余信号与原始信号能量之比，剩余信号能量占比足够小时算法也会提前停止。这些参数每一个都对应一种停止策略，理解它们的意义比死记默认值更重要。

2.4 端点效应和模态混叠，绕不开的两大顽疾

EMD用三次样条做包络时，信号首尾两端没有完整的极值点信息，样条插值容易在端点处大幅摆动，导致IMF的头尾出现严重失真，这就是端点效应。你可以把三次样条想象成一根有弹性的细木条，两端没有钉子钉住时它会随意翘起，让包络线和真实信号脱离。端点效应会在分解结果的头尾产生“假振荡”，直接影响瞬时频率计算。

模态混叠则是另一个更让人头疼的问题。当信号中存在间歇性高频成分（比如一段噪声脉冲夹杂在低频正弦波里）时，EMD可能把不同时间尺度的分量硬塞进同一个IMF，或者把同一个分量拆散到多个IMF里。反映在IMF上就是某个IMF里一会儿高频一会儿低频，物理含义变得含混不清。处理模态混叠的经典方法是集合经验模态分解（EEMD），我在后面实战部分会专门讲。

3. MATLAB环境下的实现：从手写函数到内置命令

3.1 工具箱与版本确认

如果你打算用官方内置函数，第一件事是确认自己的MATLAB版本在R2018a之后，并且安装了Signal Processing Toolbox。在命令窗口输入：

which emd

如果返回路径，说明可以直接使用；如果报错“未找到emd”，则需要检查工具箱或版本。我见过不少人卡在第一步，误以为emd是Deep Learning Toolbox里的函数，其实它属于信号处理工具箱。没有工具箱的也别急，接下来这个手写版本能帮你理解算法底层，同时也能在旧版本MATLAB上运行。

3.2 手写一个教学版EMD（完整代码）

为了讲清楚原理，我自己写过一版简化EMD，去掉了很多保护性判断，但核心逻辑完整。它使用三次样条插值构造上下包络线，用标准差准则控制筛分次数，适合用在教学和简单信号测试上。代码如下：

function [imfs, residue] = emd_manual(x, max_imf, sd_thresh) % 简化版经验模态分解 % 输入： % x - 实数列向量信号 % max_imf - 最多提取IMF个数，默认8 % sd_thresh - 筛分停止阈值，默认0.2 % 输出： % imfs - 行向量组成的矩阵，每一行是一个IMF % residue - 残余趋势项 if nargin < 2 || isempty(max_imf), max_imf = 8; end if nargin < 3 || isempty(sd_thresh), sd_thresh = 0.2; end x = x(:); n = length(x); imfs = zeros(max_imf, n); residue = x; for k = 1:max_imf h = residue; sd = 1; iter = 0; while sd > sd_thresh && iter < 200 iter = iter + 1; % 通过差分找局部极值点 d_h = diff(h); max_locs = find(diff(sign(d_h)) < 0) + 1; min_locs = find(diff(sign(d_h)) > 0) + 1; if length(max_locs) < 2 || length(min_locs) < 2 break; end % 护端：把首尾点加入极值序列，抑制端点摆动 ext_max = [1; max_locs; n]; ext_min = [1; min_locs; n]; val_max = [h(1); h(max_locs); h(end)]; val_min = [h(1); h(min_locs); h(end)]; % 三次样条包络 upper = interp1(ext_max, val_max, (1:n)', 'spline'); lower = interp1(ext_min, val_min, (1:n)', 'spline'); mean_env = (upper + lower) / 2; old_h = h; h = h - mean_env; % 标准差停止准则 sd = sum((old_h - h).^2) / sum(old_h.^2 + eps); end if length(max_locs) < 2 || length(min_locs) < 2 % 剩余信号没有足够极值点，不再继续 break; end imfs(k, :) = h; residue = residue - h; end % 去掉未用到的空行 imfs = imfs(any(imfs, 2), :); if size(imfs, 1) < max_imf imfs(max_imf, :) = 0; % 保持输出维度一致，可自行修改 end end

这个版本的代码刻意做了简化，有几个地方你如果直接拿去分析真实信号，可能会发现效果不如官方函数。比如极值点检测用差分法在信号极值点附近斜率变化很小时容易漏检，护端方式也只是简单地把首尾当作极值点，没有做镜像延拓。我建议把这段代码当作理解算法流程的“白盒”，搞懂之后再切换到官方函数做正式分析。

运行这段代码时，可以在脚本里这样调用：

fs = 1000; t = (0:999) / fs; x = sin(2*pi*50*t) + 0.5*sin(2*pi*5*t + 10*t.^2) + 0.3*randn(size(t)); [imfs, residue] = emd_manual(x, 6, 0.2);

用plot逐行画出IMF，你会发现第一个IMF基本是噪声和最高频成分，后面是调频信号，残余项接近直线的趋势。这个直观效果能帮你建立对EMD的信任感。

3.3 调用MATLAB自带的emd函数

手写版只适合教学，正式项目我推荐直接用内置函数。语法非常简洁：

[imf, residual] = emd(x);

内置函数返回多个IMF矩阵和残余向量。如果要限制IMF个数和筛分迭代次数，可以这样写：

[imf, residual] = emd(x, 'MaxNumIMF', 6, 'SiftMaxNumIterations', 150);

我更常用的是配合'Display'参数，它会在命令窗口输出每一步的筛选过程。例如：

[imf, residual] = emd(x, 'Display', 1);

这样可以直观看到算法在哪些位置停止筛分，也能反过来帮你理解参数对结果的影响。官方emd内置了端点效应处理，默认做镜像延拓，所以同样信号下它分解出的IMF通常比我的简化版更干净。

3.4 用hht画希尔伯特频谱

EMD的经典搭档是希尔伯特变换。对每个IMF做Hilbert变换，可以得到它的瞬时幅度和瞬时频率，把所有IMF的瞬时频率随时间变化画在同一张图上，就是Hilbert谱。MATLAB中一句话就能完成：

hht(imf, fs);

运行后会自动画出三个窗口：希尔伯特谱、瞬时频率曲线和边际谱。希尔伯特谱有色彩映射和时间轴，能清楚看到信号频率随时间的变化。边际谱则把所有时间点的瞬时幅度按频率累加，可以看成EMD版的“频谱”。我在实际项目中经常用边际谱替代普通FFT幅值谱，因为它由信号能量在不同瞬时频率上的分布构成，频率分辨率不再受到窗长限制，极端情况下甚至在理论上能够分辨两个间隔很窄的频率成分。

4. 实测案例：给一段非平稳混合信号做EMD

4.1 构造一个典型混合信号

为了直观感受EMD的价值，我构造了一段仿真信号，包含三个部分：一个稳定的50Hz正弦波、一个频率随时间漂移的调频分量，以及一个缓慢上升的线性趋势，最后加了少量白噪声。采样率设成1000Hz，时长1秒。这个信号在FFT下可能分不清调频分量和趋势，但EMD应该能把它们逐层拆开。

fs = 1000; t = (0:999)' / fs; f_trend = 1 + 3 * t; x = sin(2*pi*50*t) + sin(2*pi*(5 + 10*t).*t) + f_trend + 0.15*randn(1000,1);

这个信号里50Hz正弦是固定频率，调频分量的瞬时频率从5Hz线性增加到15Hz，趋势项一直在上升。噪声则是明显的毛刺。把x直接做EMD分解：

[imf, residual] = emd(x); disp(size(imf));

我这里跑出来的结果通常有3~4个IMF加一个残差，不同的MaxNumIMF设置会导致数量变化。

4.2 逐条解读IMF的含义

用subplot把IMF按顺序画出来，你会看到分层规律非常明显。

第一个IMF主要对应白噪声和部分高频细节，波形看起来毛糙，但幅度很小；第二个IMF对应50Hz正弦分量，波形非常干净，频率基本稳定；第三个IMF对应调频分量，你能在图上看到它的振荡周期从密变疏，因为瞬时频率从5Hz增加到了15Hz；最后剩下的残差则是一条近似线性上升的斜线，正好对应趋势项。

这就是EMD让人兴奋的地方：它并没有被提前告知信号里有哪几个分量，全凭信号自身的极值点结构，就自动把成分按频率从高到低剥离了出来。不过要注意，实际分解得到的IMF与原始分量不是一一精确对应，有时候频率相近的成分会发生轻微混合，比如噪声可能和50Hz谐波纠缠，这属于模态混叠的范畴，后文会讲对策。

4.3 用IMF做去噪和重构

既然噪声主要集中在第一个IMF，一个最简单的去噪策略就是丢弃第一个IMF，把剩下的IMF和残差相加重建信号。用代码表达就是：

x_recon = sum(imf(2:end, :), 1) + residual';

把重建信号与原始信号叠加画出来，能看到高频毛刺被大幅抑制，同时50Hz正弦和调频结构被保留。这一点在工程上非常实用：比如分析轴承振动数据时，第一个IMF往往对应高频噪声和干扰，去掉它之后故障冲击特征更容易暴露。

但我要提醒一句：不要盲目丢弃第一个IMF，因为有些信号的故障特征本身就是高频的，第一个IMF里可能藏着更重要的信息。正确做法是结合时域波形、瞬时幅度和频谱，判断某个IMF与目标特征是否相关，再决定是保留还是丢弃。我见过不少初学者把第一个IMF一律当作噪声拿掉，结果把真正的冲击信号也削掉了，这个坑需要留意。

4.4 EMD与FFT在同一信号上的直观对比

把同一段信号的FFT幅值谱和边际谱放在一起，对比效果非常明显。FFT的结果会显示在5~15Hz附近有一个带宽很宽的包，对应那个调频分量，因为普通FFT无法区分瞬时频率变化；还能看到50Hz的峰，以及一部分趋势项经Hanning窗泄漏后的低频隆起。边际谱则不同，它能在5~15Hz区域形成一个清晰平滑的带，准确反映瞬时频率的变化范围，50Hz处的峰也更尖锐。

另一个优势在于时频定位。FFT无法告诉你调频分量的频率是线性增加的，而Hilbert谱中你可以看到一条从左下角向右上角延展的亮线，瞬时频率随时间线性升高，一眼就能判断信号是调频信号。如果是机械故障产生的瞬态冲击，Hilbert谱上能看到在冲击时刻频率突然扩展的竖线，这在FFT里是绝对看不到的。

5. 实战经验总结：参数调优、边界处理和排错记录

5.1 端点效应的几种处理手段

就算用官方emd，端点效应也无法完全消除，只能减弱。最常见的处理是先对信号做镜像延拓：把信号左右两端各自向外翻转一段，当作极值点的“幻觉边界”，使得三次样条包络在原始信号两端有约束，而不是自由摆动。官方emd内部默认采用类似策略，但如果你手写版本，可以在包络前人为添加端点极值点。

我常用的另一种做法是“数据两端补一段预测值”，也就是用AR模型预测信号左右两侧的几个点，拼接到原始数据前后，做完分解后再截掉对应长度。这个方法在多步预测任务中效果很好，但引入的预测误差也会随分解层级传递。如果只是观察IMF形态，直接用截断法就行：分解后舍弃每个IMF首尾的5%~10%长度，不纳入后续计算。虽然浪费了一点数据，但能保证瞬时频率计算不发散。

5.2 模态混叠的补救方案：EEMD与CEEMDAN

当模态混叠严重时，标准EMD让人头疼。一种实用的做法是EEMD（集合经验模态分解），本质上是“加噪声后反复EMD再平均”：在原始信号上加入多组均值为零的白噪声，分别做EMD，再把所有分解结果按IMF序号取平均。因为白噪声在每次独立实现中都不同，随机噪声会相互抵消，而真实分量会被保留，从而抑制模态混叠。

MATLAB没有内置EEMD函数，但有第三方代码包，也可以在文件交换中心找到。我在没有外部工具箱时自己写了一段，核心就三句话：循环生成x + randn(size(x))*alpha，调用emd，把所有IMF累计后除以循环次数。加入白噪声的标准差通常是原始信号标准差的0.1~0.4倍，太小起不到作用，太大又会破坏真实信号。

CEEMDAN是全经验模态分解，它在EEMD基础上进一步优化，逐级对残差加入自适应噪声，IMF的完备性更好，计算量也更大。如果你只是做科研分析，建议优先试试CEEMDAN，分解结果通常比EEMD更干净。不过要注意，不论哪种方法，计算耗时都比单个EMD高出一个量级，处理长序列时要做好心理准备。

5.3 常见MATLAB报错与解决方案

我把在这类问题里遇到的典型报错整理成了表格，方便你遇到时快速对症处理。

有一个隐藏问题很难检查：当你把IMF矩阵单行传给plot时，如果矩阵行数是1，画出来的波形可能变成了按行显示的多色曲线，让我一度以为自己分解出错了。后来才发现这是因为MATLAB对向量和矩阵绘图行为不同，用plot(imf(1,:))就行。如果你发现某个IMF里出现了周期为1的“锯齿”状模式，先检查是不是绘图问题，再检查信号是否有零点偏移。

5.4 我踩过的坑和调参心得

最后说几个只有真正用过EMD才会遇到的细节。

第一个坑是采样率。EMD对极值点采样天然敏感，采样率太低会漏掉高频振荡的真实极值，导致第一个IMF和第二个IMF混在一起。我自己有个经验：目标最高频率成分至少需要20个采样点左右，如果信号带宽很高，先提高采样率再做分解，不要在欠采样数据上强行跑EMD。

第二个坑是筛分阈值。默认阈值0.2是在大量经验里磨合出来的，但它不代表对所有信号都最优。我处理很平滑的正弦信号时，阈值即使设成0.1也很快收敛；但处理带毛刺的振动信号时，如果阈值太小，筛分可能在局部振荡里来回震荡，收敛慢还产生不真实的IMF。遇到这种情况，我会把SiftMaxNumIterations调小到50~80，宁可牺牲一点对称性，也要保住物理意义。

第三个坑是不要把EMD结果当成标准数学分解。它没有唯一解，同一个信号用不同参数分解，得到的IMF可能不完全一样，这是算法本身的自适应特性。做研究时一定要固定参数并记录完整，否则复现性会很差。

第四个坑是关于趋势项。EMD最后的残差往往是趋势项，但它不是趋势的唯一估计。如果信号本身含有直流分量，残差会直接把它体现出来；如果信号没有趋势，残差会在零附近小幅波动。千万别把残差当作“噪声”，它在预测类任务里往往是重要特征。

我之前处理风电功率预测时，把一段很长的风速序列用EMD拆成8个IMF，然后分别送入LSTM，再叠加预测结果，效果比直接预测功率序列提升了大概15%。但也遇到过一个问题：分解后的IMF在数据边界存在严重端点效应，导致在线预测时最新的几个值完全失真。后来我用镜像延拓并对每个IMF同步外推，才解决了边界漂移的问题。这种细节只有亲自调试才能体会，所以我特别建议你拿到代码后，先拿仿真信号把每个环节都画图看一看，比用真实数据盲跑要有效率得多。

EMD不是一个银弹工具，但你一旦掌握它的脾性，就能在非平稳信号分析里多一把非常顺手的利器。希望这份从原理到实战的整理，能让你在MATLAB里把EMD用得明白、用出价值。