简介：这是一份面向MATLAB用户的EMD（经验模态分解）算法实现资源，核心程序emd.m完整覆盖包络线构造、IMF迭代提取及希尔伯特变换等关键步骤，适用于振动分析、故障诊断、生物医学信号处理等非线性非平稳信号研究场景。压缩包共44个文件，以38个m源文件为主，另有5个asv备份文件与1个fig图例文件，整体约90KB，便于对照学习与二次开发。除主程序外，资源还提供emd_Rilling_paracon.m等配套函数，并附有简单谐波、复杂谐波、间谐波、简单中断、复杂暂降等多个典型信号的分解示例，可直观验证算法在不同输入下的效果。该资源已有1314人浏览学习，适合需要快速上手EMD分解、理解IMF提取原理或开展信号处理实验的初学者与研究人员；通过研读源码与测试用例，读者可掌握完整分解流程，并直接移植到自己的项目中。

1. EMD matlab程序：从看不懂频谱到手写分解代码

某次拿设备振动信号找故障特征，频谱上只知道某个频段被噪声抬起来了，用滤波器来回切都觉得会把有效分量切掉。后来把这段信号交给EMD经验模态分解处理，跑一遍MATLAB程序，信号被自动拆成一组成分，高频噪声、冲击分量、趋势项分得干干净净，再画包络谱故障特征就非常明显了。这套EMD matlab程序不依赖复杂工具箱，核心代码就是极值点查找、三次样条包络和筛分循环三个部分，适合信号处理方向的学生，也适合做机械故障诊断的工程师拿来当预处理工具。下面从算法原理讲到代码怎么拆，到怎么避免端点效应那些坑，最后给一条改成EEMD的路线。

2. EMD分解原理与程序架构：先弄清楚筛分迭代再动手

2.1 EMD把信号拆成什么

EMD全称是经验模态分解，核心思路是把一段复杂的非线性非平稳信号，拆成一组分量的叠加。每个分量都是一个本征模函数，最后一个剩余分量是趋势项。本征模函数要满足两个条件：一是极值点数和过零点数相等或至多差1；二是上下包络的均值在局部意义上为0。这两个条件保证了每个IMF是单一尺度的振荡，能反映信号内在的物理成分。

拿振动信号举例，第一个IMF通常是最高频的噪声或冲击，后面的IMF对应轴承故障特征频率附近的分量，最后的残余项是整体的趋势。第二个条件里的“局部对称”在精确数学意义上很难量化，工程实现时用筛分阈值来近似，而不是严格计算包络均值零点。这个近似也是后面各种参数坑的来源。

2.2 筛分迭代过程

一次完整的筛分分四步：找出当前信号的全部局部极大值点和极小值点；用极大值点做一次三次样条插值得到上包络，用极小值点做三次样条得到下包络，把上下包络平均得到均值序列；让当前信号减去这个均值，把信号里的低频包络中心去掉，剩下的就是更高频的振荡候选；判断这个结果是否满足IMF条件，满足就输出，不满足就回到第一步继续筛。

内部反复执行上面四步是单个IMF的内循环。得到一个IMF之后，原信号减去这个IMF，剩余部分继续做同样的流程，直到剩余部分极值点不足两个、变成单调序列为止，这是外循环。整体思路不复杂，写成代码也就一百多行。真正的难点集中在极值点漏检和端点效应两个细节上，后面避坑章节会展开说。

2.3 程序结构和一个清晰的数据流

我不太喜欢把EMD做成一个黑匣子。拆这个项目时，我习惯把程序按职责拆成独立函数，这样后来换端点延拓策略或者换停止判据时，不用动主循环。拿到手的版本函数命名可能不一样，但内部划分大同小异。

程序的数据流是：emd_demo生成或读取信号后传给emd_decompose；主函数把当前待分解序列交给extrema拿极值点；包络函数用极值点做三次样条并计算上下包络均值；筛分循环把均值从原信号中剥离得到IMF；剥离后的残余项继续进入外循环，直到残余项极值不足。这个过程很直观，调试时可以在emd_decompose里加断点，看每一轮包络和均值长什么样。

2.4 停止条件：SD阈值与三阈值判据

筛分循环为什么不用固定迭代次数，而是用判据？因为筛分次数少了IMF不干净，多了会把信号筛成纯正弦，丢失调制信息。经典做法是计算SD，也就是前后两次筛分结果的能量相对变化：

sd = sum((h_{k-1} - h_k).^2) / sum(h_{k-1}.^2)

当sd小于0.25左右时认为本轮筛分已经稳定，停止内循环。阈值太大IMF粗糙，太小则过筛分。后来研究里又提出了更严格的三阈值判据，用包络幅值做归一化，判断每一时刻的均值相对幅值是否足够小，同时允许少量时刻越界。三阈值判据计算复杂，但对间歇性信号更稳。如果你的信号是机械振动，我建议先用SD等于0.25跑通，再去研究更严格的判据。

换句话说，EMD的代码难点不在主循环本身，而在极值点漏检和端点效应，这两块恰恰是后面所有坑的来源。

3. 核心代码拆解：极值点、包络均值与筛分主循环

3.1 极值点查找与平坦段处理

很多EMD实现的第一个坑在extrema函数。直接用find加diff找极值点，遇到信号中间有平台段就会漏检，包络插值跟着变形。我一般这样写：

function [indmin, indmax] = extrema(x) % 输入：x 为列向量信号 % 输出：indmin 极小值索引；indmax 极大值索引 x = x(:)'; % 一阶差分符号：正表示上升，负表示下降，0 表示平坦 dx = sign(diff(x)); % 平坦段处理：连续 0 都继承前一个方向，让极值点落在平台的末端 for k = 2:length(dx) if dx(k) == 0 dx(k) = dx(k - 1); end end % 符号由 +1 变成 -1 是局部极大值，由 -1 变成 +1 是局部极小值 % diff 会让序列短一位，索引加 1 找回原位置 change = diff(dx); indmax = find(change == -2) + 1; indmin = find(change == +2) + 1; % 极值点数量不足时把两端补进来，避免样条插值没法做 if numel(indmax) < 2 || numel(indmin) < 2 indmax = unique([indmax, 1, numel(x)]); indmin = unique([indmin, 1, numel(x)]); end end

这段代码里平坦段处理是核心。信号从上升到平台再到下降，差分符号会是+1、0、0、-1，替换后变成+1、+1、+1、-1，极值点就被定位在平台末端，符合EMD对极值点的定义。change等于-2对应符号从+1切到-1，等于+2对应从-1切到+1。最后那次边界补充很关键，样条插值至少需要两个点，极值点不足时程序会直接报错。

3.2 三次样条包络与均值序列

拿到极值点后，上包络由极大值点插值构造，下包络由极小值点插值构造。常见做法是用csape，它对非均匀采样的时间轴更友好，比内置spline的边界条件更可控：

function [m, a] = envelope_mean(t, x, indmin, indmax) % 输入：t 时间轴；x 信号；indmin、indmax 极值点索引 % 输出：m 上下包络均值；a 包络半幅值 up = csape(t(indmax), x(indmax), 'natural'); u = fnval(up, t); lo = csape(t(indmin), x(indmin), 'natural'); l = fnval(lo, t); m = (u + l) / 2; a = (u - l) / 2; end

自然边界条件会让包络在两端更平滑，但端点飞翼问题主要出在极值点稀疏的区域，样条外插值不受控。把a一并返回是我自己的习惯，SD判据用不到它，但后面改三阈值判据时直接拿a做归一化，不用再改函数签名。如果出现上下包络交叉这种异常，优先排查extrema返回的位置，多半是极值点漏检。

3.3 筛分主循环

主循环是EMD程序的骨架，内外两层迭代都在这里：

function [imf, r] = emd_decompose(x, max_iter, sd_tol) % EMD 分解主入口 % 输入：x 列向量；max_iter 内循环最大筛分次数；sd_tol 停止阈值 % 输出：imf，每行一个IMF，最后一行是残余项；r 单独返回残余项 x = x(:)'; t = 1:numel(x); imf = []; r = x; stop = false; while ~stop h = r; iter = 0; while iter < max_iter iter = iter + 1; [indmin, indmax] = extrema(h); if numel(indmin) < 2 || numel(indmax) < 2 % 极值点不足，继续筛只会得到无意义振荡 break; end m = envelope_mean(t, h, indmin, indmax); h_prev = h; h = h - m; % SD 判据，前后两次迭代的能量相对变化量 sd = sum((h_prev - h).^2) / sum(h_prev.^2); if sd < sd_tol break; end end imf = [imf; h]; r = r - h; % 外循环停止：残余项极值点不超过2个，就认为是趋势 [indmin, indmax] = extrema(r); if numel(indmin) + numel(indmax) < 3 stop = true; end end imf = [imf; r]; end

代码逻辑顺序是：外循环负责剥离IMF，内循环负责把当前IMF筛干净。内循环里每次都要重新找极值点，因为h每次减掉均值后极值位置会变，必须重新识别。为什么用前后两次h的能量差而不用包络均值直接判断？因为包络均值理论上为0，实际计算总有误差，能量差能更稳定地表达这个近似程度。

两个参数需要说明一下。max_iter默认取100，如果信号含强噪声，筛分上百次仍不收敛，这个上限能强制出结果，避免死循环。sd_tol默认取0.25，这是经验值。实测阈值从0.2放到0.3，对前两个IMF影响不大，对后面的IMF会明显改变分解结果，调参时优先改这个数。

3.4 跑起来：仿真信号验证

拿到代码先不要拿真实数据试，建议先用已知组成的仿真信号验证分解正确性。下面的demo把50Hz和120Hz两个正弦叠加，再加一个线性趋势项：

Fs = 1000; t = (0:1/Fs:1)'; % 两个正弦 + 趋势项 + 小幅噪声 x = sin(2*pi*50*t) + 0.5*sin(2*pi*120*t) + 0.2*t + 0.05*randn(size(t)); [imf, r] = emd_decompose(x, 100, 0.25); % 重构误差检查，EMD 完备性验证 x_rec = sum(imf, 1); fprintf('重构误差：%.3e\n', max(abs(x_rec - x'))); for k = 1:size(imf, 1) subplot(size(imf, 1) + 1, 1, k); plot(t, imf(k, :)); ylabel(sprintf('IMF%d', k)); end subplot(size(imf, 1)+1, 1, size(imf,1)+1); plot(t, x, 'k'); % 原始信号对比

运行后你应该看到：第一个IMF大体上是120Hz正弦，第二个IMF大体上是50Hz正弦，最后一个IMF是接近直线的趋势项，重构误差在1e-14量级。如果第一个IMF同时混入50Hz和120Hz，说明两个正弦频率太接近或者噪声太大，把频率拉开到50Hz和300Hz再试。这个验证过程很重要，它确认代码本身没毛病，后面遇到真实数据时，你才有底气判断问题是出在程序还是数据。

提示：如果demo结果和预期差太多，先检查t和x的维度，EMD相关函数多数默认输入为列向量，传成行向量会让imf和后续绘图对不上。

4. 避坑指南：EMD分解中常见的问题与排查办法

4.1 端点飞翼：分解出来的IMF两端向外甩

现象是画第一个IMF时，两端幅度突然放大，像海鸥翅膀一样翘起来，严重的时候整个两端都被污染，只有中间还是正常的。原因是信号两端附近往往没有极值点，样条在端点外部没有锚点，外插值不受控制，这跟数据本身无关，是插值算法的固有缺陷。

解决的办法是跑EMD前先做端点延拓。我常用的方式是把信号在左右两端做镜像对称，镜像长度取信号长度的十分之一左右，补出虚拟极值点后再做包络。端点处有了极值锚定，飞翼会明显缩小。另外，后面做特征提取时只用中间稳定的段，两端直接截掉，也是一个实用的兜底办法。

4.2 模态混叠：一个IMF里混进了两个尺度

现象是明明该是单一频率的分量，画出来明显有两段不同频率的振荡。原因是信号中间歇性出现的高频分量干扰了筛选过程。这个问题在机械振动里特别常见，轴承故障产生的冲击是间歇性的，EMD会把它和连续振动分量搅在一起。

解决的办法是换EEMD，通过加入成对白噪声让不同尺度的信号自动分离，多次平均后噪声抵消，模态混叠能被明显压制。如果你只是快速看看结果，也可以先把冲击段从信号里切掉再分解，但这样做会破坏信号的连续性，只是临时方案。

4.3 过筛分：IMF被筛成一个纯正弦

现象是IMF看起来太完美了，和正弦波一样，但和原始信号的包络调制关系对不上。原因是内循环迭代次数过多或者SD阈值太严格，导致每个IMF都被强制筛成局部对称的纯调频形态。这是新手最容易忽略的坑，因为结果看起来比真实信号还干净，反而会误以为分解很成功。

解决的办法是把max_iter从100降到10到20，sd_tol从0.25放到0.3。我把筛选迭代上限写进函数参数而不是写死，就是为了排查这种问题。如果某个IMF的物理意义不明，优先怀疑过筛分，不要先怀疑程序。

4.4 平坦段漏检极值点

现象是对含有阶跃或平台特征的信号做分解，第一个IMF严重变形，但重构误差仍然很小。原因是extrema函数遇到连续相等的数据点没有正确处理，把平台中间误判为极值，导致包络形状错误，筛分过程顺着错误的包络跑偏了。

解决的办法是用3.1里的平坦段处理逻辑，把连续0符号继承前一个方向，让极值点定位在平台末端。如果跑别的开源版本遇到结果不同，先看它的extrema是不是直接find加diff，很多时候分歧就出在这里。两个版本跑出不同IMF个数，不代表哪个程序坏了，而是极值定义不一样。

4.5 分解结果不可直接对比

现象是同一段信号，换一台机器或者换一个EMD实现的包，跑出来的IMF个数和形状都不一样。原因是端点延拓方法、极值定义、包络函数边界条件这些细节不同，都会导致结果偏差，EMD本身没有绝对唯一解。

解决的办法是固定程序版本，记录每次运行时的信号长度、延拓方式和关键参数。比较不同版本结果时，我主要对比IMF的个数和重构误差，而不是逐点对比波形。在工程上，只要分解结果能稳定支撑后续特征提取，就不必纠结是不是最优分解。

5. 二次开发：把基本EMD改成EEMD并接上故障特征提取

5.1 加噪声平均改造成EEMD

改进思路很简单：对原始信号加一组正负成对的白噪声，分别做EMD，最后把所有结果按IMF序号取平均。白噪声会让不同尺度的信号自动落到相应的IMF里，多次平均后噪声相互抵消：

function imf = eemd(x, ensemble, noise_amp) % x 输入信号；ensemble 总体次数；noise_amp 噪声幅值 imf_sum = []; for k = 1:ensemble x_noisy = x + noise_amp * randn(size(x)) * std(x); [imf_k, ~] = emd_decompose(x_noisy, 100, 0.25); if k == 1 imf_sum = imf_k; else % 对齐IMF个数：少的继续累加已有的 n = min(size(imf_sum, 1), size(imf_k, 1)); imf_sum(1:n, :) = imf_sum(1:n, :) + imf_k(1:n, :); end end imf = imf_sum / ensemble; end

噪声幅值一般取信号标准差的0.2倍到0.4倍，总体次数取100到500。取多了计算开销大，取少了噪声消不掉。这段代码为了简短做了对齐处理，实际工程里更稳妥的做法是记录每次IMF的个数，用插值对齐后再平均。

5.2 用包络谱找故障特征频率

EMD的典型工程用法不是直接看IMF波形，而是选出一两个包含故障特征的IMF做包络谱。常见做法是计算每个IMF和原始信号的相关系数，挑相关系数最高的前几个，做Hilbert变换得到解析信号，再对解析信号幅值做FFT：

% 挑选与原始信号相关性最高的 IMF c = zeros(size(imf, 1), 1); for k = 1:size(imf, 1) R = corrcoef(x(:), imf(k, :).'); c(k) = abs(R(1, 2)); end [~, idx] = sort(c, 'descend'); sel = idx(1:min(3, size(imf, 1))); % 对选出的 IMF 做包络谱 env = abs(hilbert(imf(sel(1), :))); spec = abs(fft(env));

包络谱里观察故障特征频率及其倍频，是旋转机械诊断的常规流程。注意横轴频率要按采样率换算，否则倍频对不上。

5.3 参数接口和工程落地

把程序留三个可调参数就好：延拓方式、SD阈值、最大筛分次数。延拓方式影响端点，SD阈值影响分解细腻度，筛分次数影响过筛分。我通常把这三个参数做成函数入参，而不是埋在代码里写死，这样批量处理数据时可以在循环里挨个调。

从那以后，我每次拿到一批新的振动信号，都会先用默认参数跑一遍，看IMF个数和重构误差，再检查前几个IMF波形是否有飞翼和混叠，最后才进入特征提取。这个流程虽然多花五分钟，但能挡掉后面大多数分析翻车的可能。希望帮到你。

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