EEMD集合经验模态分解：MATLAB工具箱与模态混叠解决实战

简介：EEMD（集合经验模态分解）是针对EMD方法模态混叠缺陷提出的一种噪声辅助数据分析方法。这份MATLAB程序代码将白噪声均匀分布时频空间的自适应分解原理完整实现，提供可直接运行和调用的函数与脚本，面向信号处理、故障诊断、振动分析、生物医学数据分析等领域的研究者和工程师，适合需要应用HHT类时频分析方法的各类技术场景。资源包共46个文件，其中42个m文件构成算法主体，包括主分解程序、极值点查找、包络拟合、Hilbert变换、谱分析及显著性检验等完整模块，另有2个csv文件存放示例数据，1个txt文件为算法说明，1个p文件用于辅助功能；整个压缩包仅189KB，结构精简，便于快速下载和使用。目前已有5188人学习下载。通过这套代码，使用者能够直接执行EEMD分解得到各阶本征模态函数，结合附带的频谱分析、边际谱、置信区间等扩展程序，可进一步完成时频特征提取、信号重构与二次开发，是学习集合经验模态分解算法和推进科研工作的实用工具。

1. EEMD集合经验模态分解：一套解决模态混叠的MATLAB工具箱

EEMD（集合经验模态分解）这套MATLAB程序代码，是给被EMD模态混叠困扰的人准备的：信号里间歇出现的冲击成分，分解完会同时出现在两三个IMF里，特征频率根本数不出来。EEMD的处理思路很直接——给原始信号反复叠加白噪声、做多次EMD再取集合平均，用噪声的统计特性把混叠成分逼出来。代码包里有完整的EEMD主函数、底层EMD筛选循环、极值点查找工具和仿真验证脚本，拿到手改一行数据路径就能跑，适合做振动信号分析、故障诊断、去噪和时频分析，也适合接在机器学习特征提取流程的前端。

2. EEMD算法原理：集合平均为什么能把混叠压下去

2.1 从EMD的筛选过程看混叠从哪来

EMD的分解对象是信号的局部特征尺度，它不需要预设基函数，理论上能把任意信号拆成有限个固有模态函数（IMF）加一个残差。每个IMF要求极值点个数与过零点个数相等或至多相差一个，且上下包络的局部均值处处接近零。代码里的实现路径很固定：找到当前残差信号的所有极大值和极小值，用三次样条分别插值得到上包络与下包络，取上下包络均值作为局部均值线，把局部均值从信号中减去，重复这个过程直到筛分指标满足停止条件。每筛出一个IMF就从残差里扣掉，再对残差做同样的操作，直到残差只剩单调趋势。

这套逻辑在频率成分连续存在的信号上表现很好，比如两个正弦波叠在一起，分解结果能干净地按频率分层。问题出在间歇信号上：设想一个50Hz分量和一个只在很短时间窗内出现的80Hz分量叠加，80Hz分量在大多数时间点上没有极值点，三次样条包络在这些区间只能靠50Hz分量的极值形态来推测，结果包络被低频趋势带偏。筛分到最后，80Hz的高频能量被分摊到好几个IMF里，或者和50Hz挤在同一个IMF里。这是模态混叠最典型的症状——分解出的IMF数量和内容不稳定，物理意义模糊，这就是做故障诊断的人最头疼的问题。

2.2 白噪声为什么能帮忙：集合平均的统计机制

EEMD的改进思路很巧妙：既然高频分量的间歇性让极值点分布不均匀，那就主动往信号里添加宽带均匀的白噪声。白噪声的极值点在整个时间轴上密集且均匀，相当于给样条包络一个高频参考骨架，让每个局部的极值统计变得连续。对加了噪声的信号做一次EMD，结果是真实模态和噪声响应混在一起；重复NE次，每次使用独立的噪声实现，然后把所有分解结果按IMF序号取平均。白噪声和信号不相关，在平均中互相抵消，而真实模态成分因为每次都出现，会在平均里保留下来。

两个参数直接决定效果。NE是集合次数，实际使用中100到500在文献里最多见，日常调试先用50，最终跑结果用200起步。NE太小，平均后会有明显的噪声残留，IMF波形毛毛糙糙；NE太大只是线性增加计算量，精度提升会碰到天花板。Nstd是添加噪声的标准差相对于原始信号标准差的比例，最常用的值是0.2，范围在0.1到0.4之间。Nstd偏小，高频参考骨架强度不够，混叠压不住；偏大则信号被噪声改写，分解出的IMF保真度下降。参数没有普适最优值，碰到一套新数据，画几张频谱对比着调才是常态，这也是EEMD被一些人叫“玄学”的原因——但掌握了规律之后，参数其实很好定。

2.3 EEMD与EMD的适用边界

把两者的差异整理成一张表，方便拿到新数据时快速决定用哪套：

对成分简单且各频段连续存在的信号，EMD又快又稳，没必要上EEMD。但对含冲击、调制、突变或者信噪比较低的工程信号，EEMD几乎成了必选项。还有一个容易被忽略的好处：EEMD引入白噪声后，不同样本之间分解出的IMF尺度范围更一致，这对后续批量提取特征、训练分类器很重要——EMD出来的特征在不同样本间对齐性差，EEMD要好很多。代价就是慢，一个10万点样本跑NE=200大概要一两分钟，属正常水平，优化方向是并行化或者减小NE。

3. 核心代码拆解：eemd.m、emd_core.m与极值点查找的完整链路

3.1 文件组成与调用关系

这套代码的结构很直接，没有复杂依赖，整个分解链路由三个函数组成：

调用关系是单向的：测试脚本 → eemd.m → emd_core.m → findExtrema.m。没有全局变量，没有跨函数的隐式状态，所有数据都是参数传递，这样的设计对二次开发最友好——想换噪声类型、改包络插值方式，只需要动对应函数的内部实现，其他层不用管。

3.2 eemd主函数：噪声注入、循环分解与集合平均

先看入口函数eemd.m的核心实现：

function [IMFs, residue] = eemd(signal, NE, Nstd) % EEMD 集合经验模态分解 % signal : 输入信号，行向量或列向量均可，内部统一转成行向量 % NE : 集成次数，建议 50~500，默认 100 % Nstd : 噪声相对标准差，建议 0.1~0.4，默认 0.2 % IMFs : 每行一个IMF，按频率从高到低排列 % residue: 最终残差，行向量 if nargin < 2 || isempty(NE), NE = 100; end if nargin < 3 || isempty(Nstd), Nstd = 0.2; end signal = signal(:)'; N = length(signal); noiseAmp = Nstd * std(signal); % 每次EMD产生的IMF个数可能不同，用cell暂存 allIMF = cell(NE, 1); maxRows = 0; for k = 1:NE noise = noiseAmp * randn(1, N); [imfK, resK] = emd_core(signal + noise); cur = [imfK; resK]; allIMF{k} = cur; maxRows = max(maxRows, size(cur, 1)); end % 按行对齐后做集合平均 IMFs = zeros(maxRows, N); for k = 1:NE rows = size(allIMF{k}, 1); IMFs(1:rows, :) = IMFs(1:rows, :) + allIMF{k}; end IMFs = IMFs / NE; residue = IMFs(end, :); % 平均后的最后一行是残差 IMFs = IMFs(1:end-1, :); % 其余行才是真正的IMF end

几个实现细节值得说。signal(:)'把任意形状的输入统一成行向量，后面所有处理都按行向量约定走，避免转置问题。noiseAmp = Nstd * std(signal)用信号标准差标定噪声幅度，这样不管输入信号量纲是毫米还是微米，噪声的相对强度都一致，正好对应参数定义。randn(1, N)每次循环生成一组独立的零均值单位标准差高斯白噪声，乘上noiseAmp就是本次要叠加的噪声序列。

集合平均的对齐策略是这套代码里最容易出错的地方。严格的EEMD概念要求按IMF的物理尺度对齐后平均，但在工程实现里，每次EMD分解出的IMF个数几乎必然不同。这里采用的是补零对齐：所有结果存进cell之后统计最大行数，把行数不足的后面行当零处理，最后统一除以NE。这种策略对前几个高频IMF效果很好，因为高频IMF的分解次序在每次试验中都很稳定；越往后的低频IMF，对齐误差会越大，这一点到避坑章节再细说。

3.3 底层EMD筛选：包络拟合、筛分停止条件与残差

emd_core.m是真正干重活的地方，EEMD的NE次循环里每一次都在跑这个函数。核心逻辑是两层循环：外层循环负责不断从残差里剥出IMF，内层循环负责对一个候选信号反复做“减去包络均值”的筛分操作，直到满足IMF条件。

function [IMF, residue] = emd_core(x) % 经典EMD筛选：输入一行信号，输出IMFs与残差 IMF = []; residue = x(:)'; while true % 残差极值点不足时，认为剩下的就是趋势项 [maxP, ~, minP, ~] = findExtrema(residue); if numel(maxP) < 2 || numel(minP) < 2 break; end h = residue; SD = 1; iter = 0; while SD > 0.2 && iter < 50 % 筛分停止：SD阈值0.2，迭代上限50 [maxP, maxV, minP, minV] = findExtrema(h); if numel(maxP) < 2 || numel(minP) < 2 break; % 极值点不足，包络无法构造 end t = 1:numel(h); upper = spline(maxP, maxV, t); % 三次样条上包络 lower = spline(minP, minV, t); % 三次样条下包络 meanEnv = (upper + lower) / 2; % 局部均值线 hNew = h - meanEnv; SD = sum((h - hNew).^2) / (sum(h.^2) + eps); h = hNew; iter = iter + 1; end IMF(end+1, :) = h; %#ok<AGROW> residue = residue - h; end end

包络均值用的是MATLAB内置的spline，以极值点位置为节点、极值数值为节点值，在一整条时间轴上进行三次样条插值。SD是经典筛分停止指标，定义是相邻两次筛分结果的能量相对变化量，阈值取0.2，这是经典EMD方法里推荐的区间下限，实际代码里取0.2到0.3都很常见，阈值越小IMF的局部对称性越好，但计算也更慢。iter < 50是一个保护壳，防止某些病态信号筛分不收敛时无限循环。

外层循环的退出条件是残差极值点少于两个，意味着剩下的信号要么单调要么是常数，已经是趋势项。每次residue = residue - h之后，没有单独对h做残差判定，而是把判定放到下一轮外层循环的开头，这样写的好处是逻辑清晰：筛出一个IMF就扣掉一个，扣到残差没有可分解的极值为止。

3.4 极值点查找与端点约束

findExtrema.m是整个分解的底层依赖，它的质量直接决定包络是否靠谱。实现上用相邻三点比较，不用先差分再找过零，避免差分放大量化噪声：

function [maxPos, maxVal, minPos, minVal] = findExtrema(x) % 相邻三点比较法找局部极值，端点在循环里自然被排除 x = x(:)'; N = numel(x); maxPos = []; maxVal = []; minPos = []; minVal = []; for i = 2:N-1 if x(i) >= x(i-1) && x(i) > x(i+1) maxPos(end+1) = i; %#ok<AGROW> maxVal(end+1) = x(i); %#ok<AGROW> elseif x(i) < x(i-1) && x(i) <= x(i+1) minPos(end+1) = i; %#ok<AGROW> minVal(end+1) = x(i); %#ok<AGROW> end end end

这个实现故意把端点排除在极值候选之外，是因为三次样条在端点处行为不稳定，把端点当极值会引入包络过冲。>=和>的不对称比较是为了应付平坦的极值平台：对一段连续相同的数值，平台最后一个点被记为极大值，平台前不会重复计数，避免同一平台产生多个极值。这套查找逻辑对采样率足够高的信号很稳，如果信号本身有较多高频毛刺，极值会偏密，包络会跟着毛刺走，建议先做低通滤波再进EEMD。

4. 跑通与验证：用已知频率的仿真信号检验分解质量

4.1 准备工作与调用方式

这套代码不需要额外工具箱，核心函数spline、randn、corrcoef、fft都是MATLAB基础函数，放在一个目录下即可运行。建议把eemd.m、emd_core.m、findExtrema.m三个文件放进项目根目录，然后用一个测试脚本调用。以下脚本是标准的验证流程：构造一个成分已知的仿真信号，调用EEMD分解，再检查分解结果是否能还原出原始成分。

clear; clc; close all; % 仿真信号：50Hz正弦 + 120Hz正弦 + 白噪声 fs = 1000; % 采样率 1000 Hz t = 0:1/fs:1-1/fs; % 1 秒时长 x = 1.0*sin(2*pi*50*t) + 0.6*sin(2*pi*120*t) + 0.3*randn(1, fs); rng(2024); % 固定随机种子，保证EEMD可复现 [IMFs, residue] = eemd(x, 100, 0.2); % 可视化：原信号 + 各IMF + 残差 figure; subplot(size(IMFs,1)+2, 1, 1); plot(t, x); title('原始信号'); grid on; for k = 1:size(IMFs,1) subplot(size(IMFs,1)+2, 1, k+1); plot(t, IMFs(k,:)); title(['IMF' num2str(k)]); grid on; end subplot(size(IMFs,1)+2, 1, size(IMFs,1)+2); plot(t, residue); title('残差'); grid on;

rng(2024)这行在调用eemd之前执行，它把随机数生成器的种子固定下来，这样eemd内部每次randn生成的噪声序列都是确定的，分解结果完全可复现。如果去掉这行，第二次运行同一段代码，IMF的波形会略有不同，这不代表程序有bug，而是EEMD的随机性使然。

这段仿真信号包含三个成分：幅度1.0的50Hz正弦、幅度0.6的120Hz正弦、幅度0.3的高斯白噪声。理想情况下分解结果应该是：IMF1是120Hz成分加少量噪声残留，IMF2是50Hz成分，余下的IMF幅度很小，残差接近一条直线。跑通后先用眼睛看图，能看出IMF1和IMF2的周期分别对应120Hz和50Hz，就说明分解方向对了。

4.2 用重构和频谱做量化验证

肉眼看图能发现大问题，但细节要靠量化指标。第一个必做的检查是重构验证：把所有IMF加残差，跟原始信号算相关系数，这个值应该无限接近1。

recon = sum(IMFs, 1) + residue; cc = corrcoef(x(:), recon(:)); fprintf('原始信号与重构信号相关系数: %.6f\n', cc(1,2));

如果相关系数明显小于0.999，说明某一步的IMF筛选有问题，常见原因是SD阈值设得太大导致IMF欠筛分，或者信号末端有异常延拓。第二个检查是频谱对照，把IMF1和IMF2做FFT，看峰值频率是否落在预设的120Hz和50Hz上：

Nfft = 2^nextpow2(fs); fAxis = (0:Nfft-1) * fs / Nfft; Y1 = abs(fft(IMFs(1,:), Nfft)); Y2 = abs(fft(IMFs(2,:), Nfft)); [~, i1] = max(Y1(1:floor(Nfft/2))); [~, i2] = max(Y2(1:floor(Nfft/2))); fprintf('IMF1峰值频率: %.1f Hz\n', fAxis(i1)); fprintf('IMF2峰值频率: %.1f Hz\n', fAxis(i2));

正常情况下输出应该是120.0 Hz和50.0 Hz，误差在频率分辨率以内。这两个检查做完，基本可以确认这套EEMD代码在机器上跑通了，接下来换真实数据才有底气。

4.3 换真实数据时要注意的输入差异

仿真信号验证通过后，把数据源换成实测信号有两种常见方式。如果数据是Excel或CSV文件，用readmatrix读进来，取其中一列作为信号向量；如果是MATLAB工作区现成的变量，直接传给eemd即可。很多实测信号带有直流偏置，输入前最好先减均值，因为直流分量会占用一个IMF位置，让分解结果多出一行本可以避免的趋势项。

实测数据的采样率和时长影响参数选择。采样率决定能分解到的最高频率，采样定理在这里同样适用，EEMD处理不了高于采样率一半的频率成分。数据长度方面，EEMD对短样本（几百个点）的分解效果不太稳定，端点效应占比太大，建议至少准备2000点以上。如果实测信号特别长，几十万点，直接跑NE=200会非常慢，先降采样或者分段处理是更现实的选择。

5. EEMD实战避坑：五个高频翻车点与排查记录

5.1 同一份数据两次运行结果对不上

现象：脚本什么都没改，只是重新跑了一遍，IMF1的波形明显变了，幅度和相位都不一致。 原因：EEMD内部用randn生成随机噪声，两次运行的噪声序列不同，集合平均在NE不够大时无法完全抹平随机差异。这不是bug，是算法自带特性。 解决：需要可复现结果的场景（论文图表、对比实验），在调用eemd之前执行rng(固定值)，比如rng(2024)，这样randn的序列每次都一样。如果更看重统计稳定性而不是可复现性，就把NE提到200以上，让平均更充分。

5.2 Nstd设得小，加了噪声还是混叠

现象：Nstd取0.01，分解结果和直接跑EMD差不多，IMF1频谱里同时出现两个峰值。 原因：噪声幅度过小，白噪声提供的极值参考骨架强度不够，间歇性高频分量对包络的干扰没有被压制住，集合平均救不回来。 解决：把Nstd提到0.1以上，多数振动信号从0.2起步调试。调参时盯着IMF频谱看，目标频率能单独落在某个IMF里且没有明显串扰，就算到位。注意不要盲目加大到0.5以上，噪声过大会把IMF波形撑变形，反而引入虚假成分。

5.3 端点处IMF剧烈甩尾

现象：每个IMF的中间段很正常，但首尾几十个点出现大幅摆动，幅度甚至超过信号本身。 原因：三次样条在数据域两端没有足够的外部极值点约束，包络外推失控，这就是经典的端点效应。EEMD叠加的白噪声在端点处同样无约束，会让甩尾更明显。 解决：最省事的办法是在分解前对信号做镜像延拓，也就是把两端的一段数据镜像翻转接上去，分解完再按原始长度截断。这个延拓数据不参与结果输出，只用来给样条包络提供端点附近的极值点。另一个办法是强制让上下包络在端点处等于信号端点值，把包络“钉”在端点上，能缓解但不能根治。我一般用镜像延拓，效果最直接。

5.4 NE拉满但运行时间不可接受

现象：NE设到500，数据长度10万点，跑很长时间都没出结果。 原因：EEMD的计算量等于NE乘以单次EMD的计算量，而单次EMD内部的筛选迭代次数取决于数据的复杂程度。长信号加高NE，计算量很容易膨胀到无法接受。 解决：先把NE压到30到50，确认分解效果正常再逐步加大。大多数工程场景NE=100到200已经够用，不要一上来就追求500。若数据真的很长，考虑降采样后再分解，因为提高采样率并不会增加可分解的频率范围上限，只会成倍增加计算时间。

5.5 集合平均后高阶IMF形状发虚

现象：IMF1、IMF2很清晰，但IMF3往后幅度变小、波形模糊，和单次EMD结果对不上。 原因：每次EMD分解出的IMF个数不完全相同，集合平均按行数对齐时，高频IMF序号稳定所以对齐正确，低频IMF在不同次试验里可能多一个少一个，强行按序号平均就把不对应的模态混在一起了。 解决：如果只关心前几个高频IMF，补零平均完全够用；要完整使用所有IMF，就需要改为按相关性对齐——计算每次分解的各IMF与原始信号的相关系数，把相关性相近的项归为一组再平均。工程上做故障诊断和特征提取，几乎都是取前3到5个IMF，这个坑其实可控。

6. 二次开发：把EEMD接进特征提取流水线

6.1 从IMF矩阵快速提取特征

EEMD分解本身不是终点，它更常用的角色是特征提取的前置步骤。以故障诊断为例，分解出IMF后我会算三个量：每个IMF的能量占比、与原始信号的相关系数、峭度。三个量拼成一个特征向量，可以直接喂给分类器。

function feat = imfFeatures(x, IMFs) % 逐IMF计算能量占比、相关系数、峭度 nImf = size(IMFs, 1); energy = sum(IMFs.^2, 2); feat = zeros(nImf, 3); for k = 1:nImf feat(k, 1) = energy(k) / sum(energy); % 能量占比 c = corrcoef(x(:), IMFs(k, :)); feat(k, 2) = c(1, 2); % 与原始信号相关性 centered = IMFs(k, :) - mean(IMFs(k, :)); feat(k, 3) = mean(centered.^4) / (mean(centered.^2)^2 + eps); % 峭度 end end

能量占比刻画每个模态在信号中的贡献权重，相关系数衡量IMF与原始信号的线性相似度，峭度则对冲击成分敏感，轴承点蚀这类故障会在特定IMF上表现出高峭度。需要时频图时，对选中的IMF调用MATLAB的hilbert函数取解析信号，再对瞬时相位求导得到瞬时频率，这一步是基础功能，不用额外工具箱。

6.2 参数自适应与批量处理的改造

批量处理时参数不能每段手工调，我一般按信噪比粗估Nstd：信噪比低的信号本身比较脏，白噪声要加足才能压住混叠，Nstd取大值；干净信号取小值。NE固定100到200，待处理数据逐段循环分解、计算特征，攒成特征矩阵交给分类器。

这一套参数自适应我踩过坑：曾经有一批数据信噪比参差不齐，固定Nstd=0.2时分解质量时好时坏，改成逐段估计噪声水平再映射Nstd之后，效果才稳定下来。从那以后，我每次新接一批数据，都会强制走一遍快速体检：固定随机种子、NE压到50、重构对比相关系数、看前两个IMF频谱，四项检查全过才把NE放大做正式分解。这套流程加起来不到一分钟，却能把EEMD参数跑偏的翻车概率降掉大半。希望帮到你。

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