简介: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 |
|---|---|---|
| 模态混叠 | 间歇信号下容易发生 | 加噪平均后显著抑制 |
| 分解稳定性 | 同段数据多次分解结果一致 | 依赖噪声随机性,固定种子才可复现 |
| 计算开销 | 单次分解,速度快 | 需要NE次EMD,速度慢NE倍 |
| 参数依赖 | 基本无参数 | 需调NE和Nstd |
| 适用场景 | 成分简单、连续信号 | 振动、故障、生物医学等复杂信号 |
对成分简单且各频段连续存在的信号,EMD又快又稳,没必要上EEMD。但对含冲击、调制、突变或者信噪比较低的工程信号,EEMD几乎成了必选项。还有一个容易被忽略的好处:EEMD引入白噪声后,不同样本之间分解出的IMF尺度范围更一致,这对后续批量提取特征、训练分类器很重要——EMD出来的特征在不同样本间对齐性差,EEMD要好很多。代价就是慢,一个10万点样本跑NE=200大概要一两分钟,属正常水平,优化方向是并行化或者减小NE。
3. 核心代码拆解:eemd.m、emd_core.m与极值点查找的完整链路
3.1 文件组成与调用关系
这套代码的结构很直接,没有复杂依赖,整个分解链路由三个函数组成:
| 文件 | 职责 | 被谁调用 |
|---|---|---|
| eemd.m | 噪声注入、循环调用EMD、集合平均 | 用户入口 |
| emd_core.m | 单次经典EMD筛选,输出IMFs和残差 | eemd.m |
| findExtrema.m | 查找局部极大值和极小值的位置与数值 | emd_core.m |
调用关系是单向的:测试脚本 → 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参数跑偏的翻车概率降掉大半。希望帮到你。
本文还有配套的精品资源,点击获取