
简介经验模态分解EMD及其改进算法是处理非线性、非平稳信号的常用工具。这份MATLAB代码包面向信号处理学习者与工程研究人员集中提供了EMD、EEMD、CEEMD、CEEMDAN四种方法的完整实现并附带示例音频便于从实际信号中观察分解效果。压缩包共8个文件以6个m脚本为主涵盖核心算法、极值点查找、主程序及可视化等模块整体大小228KB。目前已有964人学习下载。代码结构清晰配有主绘图程序用户可借助自带wav数据快速运行直观对比不同算法在模态混叠抑制和分解稳定性上的差异。对于想深入理解希尔伯特-黄变换原理、或需要在项目中应用自适应噪声完备集合经验模态分解的读者这套代码能够提供可直接修改的参考起点帮助理解迭代分解流程、边界处理与瞬时频率提取等关键细节。1. 从压缩包到 IMF各种 EMD 代码 matlab.zip 能在项目里干什么把「各种EMD代码matlab.zip」解压出来的工程师多半不是来逛代码博物馆的而是手里正压着一组非线性非平稳信号轴承振动、脑电、风速、设备温度曲线或者某只股票的分钟线。EMD经验模态分解不预设正弦基或小波基而是按信号自身的极值分布把数据一层一层筛成本征模态函数IMF再配合 Hilbert 变换得到瞬时频率谱。这个 zip 之所以叫「各种」因为里面通常不止一个 emd.m还会带上 eemd、ceemdan 变体、测试数据和画 HHT 谱的脚本——三套思路停止条件不同适用信号不同。这份压缩包不是给你读的是给你跑的。它帮你省掉了从论文到可执行代码的翻译成本剩下的问题全在调用、参数和结果判定上。下面按「sifting 到底在筛什么 → 解压后怎么跑通第一组 IMF → 模态混叠和报错怎么调 → 批量和有效性检验怎么做」四步展开命令直接抄。2. EMD 的 sifting 递推先分清 emd、eemd 与 ceemdan 再动手2.1 IMF 的两个硬条件为什么「包络均值归零」是收敛方向任何一个 IMF 都必须同时满足两个条件全段时间内过零点数目与极值点数目相等或至多相差 1任意位置由局部极大值拟合的上包络与局部极小值拟合的下包络的均值趋近于 0。第一个条件保证分量是窄带振荡第二个条件保证它关于时间轴局部对称。这两个条件不是空泛定义它们直接决定了 sifting 的收敛方向。原始信号往往叠着趋势项、噪声和多个振荡模态直接做 Hilbert 变换得不到有物理意义的瞬时频率。sifting 要做的事情就是反复减去上下包络的均值把藏在信号里的局部非对称性一层层抽走直到剩余波形满足上述两个条件。抽出来的每一层就是一个 IMF最后剩下的残余项通常是一条单调趋势或极低频缓变曲线。理解这个递推过程比背出「自适应分解」「后处理基」这类词有用得多。因为 zip 包里所有变体本质都是对「如何减均值、减到什么时候停」这两个环节做修改。先看清楚基础版本后面调 eemd 参数时才不会盲调。2.2 sifting 循环的 matlab 教学实现一次迭代在算什么读一个陌生的 matlab 代码包最忌讳上来就看几十个文件的互相调用。先把 sifting 的最小闭环写出来再回头读源码五分钟就能对上号。下面是一个教学级单次 sifting 实现省略了端点延拓等工程细节但完整保留了主循环骨架。function imf sifting_once(x, maxIter, sdTol) % 单次 sifting从 x 中筛出一个 IMF教学简化版 x x(:); % 强制列向量 h x; N numel(x); t (1:N); for iter 1:maxIter % 找局部极大值与极小值的位置 [~, locMax] findpeaks(h); [~, locMin] findpeaks(-h); % 对 -h 找峰等价于找 h 的谷 if numel(locMax) 2 || numel(locMin) 2 break; % 极值点不足包络无法继续拟合 end % 三次样条拟合上下包络端点直接带首尾点简化 up spline([1; locMax; N], [h(1); h(locMax); h(N)], t); lo spline([1; locMin; N], [h(1); h(locMin); h(N)], t); m (up lo) / 2; % 上下包络均值 hNew h - m; % 减去均值得到更接近 IMF 的波形 sd sum((h - hNew).^2) / sum(h.^2); % 标准 SD 停止准则 h hNew; if sd sdTol break; end end imf h; end这段代码的逻辑是每次迭代先找峰值和谷值用三次样条各拟合一条包络取平均后从当前信号里减掉。findpeaks需要信号处理工具箱对-h做峰值检测是为了复用同一个函数找谷值。sd是相邻两次迭代结果之间归一化能量差低于阈值就认为包络均值已经足够接近零停止本级 sifting。两个参数直接影响结果maxIter过小会提前结束IMF 里还残留可见的包络不对称sdTol设到 1e-5 以下时迭代会明显变慢但对多数工程信号收益很小常见做法是设在 0.05 到 0.001 之间。注意包络端点直接用首尾点参与样条这是刻意简化真实代码包里会在这里做镜像延拓或端点极值外插正是后面要讲的端点效应来源。2.3 代码包里的三种变体emd、eemd 与 ceemdan 怎么选zip 包里常见的emd.m、eemd.m、ceemdan.m最大的区别不在外层调用而在 sifting 里如何处理模态混叠。模态混叠指单个 IMF 里混进了不同时间尺度的成分典型表现是高频间歇信号把低频连续振荡「撕」成几段。应对思路有三种对应三个文件。变体核心思想主要代价适用场景EMD直接对原信号逐级 sifting结果不稳定、易模态混叠信号干净、单次定性分析、教学验证EEMD多次给原信号加白噪声后分别 EMD再对 IMF 集合取平均计算量大、重构误差非零含冲击、间断成分的振动与生物电信号CEEMDAN每级分解时自适应加入噪声分量边分解边消噪最慢但重构残差极小需要精确复现各 IMF 能量占比的定量分析选型依据不复杂信号平滑且只关心大概分层直接用 EMD信号里有明显间歇冲击先试 EEMD要拿 IMF 能量做健康指标或者训练模型优先 CEEMDAN。EEMD 加噪声的本质是用噪声填满信号缺口处的极值分布因此噪声幅值和集合次数是两个必调参数。但这里先不展开不同代码包的参数名和默认值不同第三节用内置函数跑通后第四节再用具体命令讲怎么配对调。3. 解压到出图用 matlab 内置 emd 函数跑通第一组 IMF3.1 解压后的路径处理先确认你拿到的是哪套 APIzip 解压后第一件事不是双击运行而是用which emd确认当前命令行会命中哪个文件。较新版本的 matlab 在信号处理工具箱里已内置emd语法是imf emd(x)右侧的 Name-Value 参数与 Flandrin 系离线包完全不同。离线包通常要求把整个目录加入 path而且不同作者写的签名差异很大有的返回imf有的返回[imf, residual, info]有的把停止准则写在parament结构体里。which emd addpath(D:\work\emd_toolbox); which emd第一次which看到的是内置函数路径addpath之后再次检查如果路径仍是内置目录说明离线包的函数名与内置函数重名且优先级更低。常见做法是把离线包目录放到 path 最前面或者直接给离线包里的emd.m改名例如emd_rilling.m再把所有内部对emd的调用同步替换。这一步不做干净后面所有基于该包的案例都会跑出「看似正常、实际是另一套代码在干活」的结果。3.2 合成信号上的最小可运行调用从构造数据到画图不引入任何真实传感器数据先用一段合成信号把整条链路打通。下面的信号叠加了 10 Hz 正弦、40 Hz 正弦、线性趋势和弱高斯白噪声结构简单到可以肉眼核对分解结果。fs 1000; dt 1/fs; t (0:dt:1-dt); x sin(2*pi*10*t) 0.5*sin(2*pi*40*t) 0.05*t 0.02*randn(size(t)); imf emd(x, MaxNumIMF, 4, SiftMaxIterations, 200); residual x - sum(imf, 2); % 残余趋势 原信号减所有 IMF tiledlayout(size(imf,2) 1, 1); for k 1:size(imf,2) nexttile; plot(t, imf(:,k)); ylabel(sprintf(IMF %d, k)); end nexttile; plot(t, residual, k); ylabel(residual);返回的imf是N x m矩阵每列是一个 IMF按频率从高到低排列m由MaxNumIMF限制。x - sum(imf,2)得到残差就是那个0.05*t线性趋势。绘图用纵向平铺每个子图里的振荡频率应当随阶数升高而降低如果 IMF1 里同时出现 10 Hz 和 40 Hz 的形态说明出现了轻微混叠通常是噪声或间歇分量的影响第四节专门处理。3.3 emd() 的四个常调 Name-Value 参数内置emd是生产级实现端点处理、停止准则和包络插值都已经工程化大多数情况不需要改默认值。真正值得手动干预的只有下表中四个参数。参数名作用使用建议MaxNumIMF限制最大分解层数只关心高频故障分量时设为 3~5能显著提速SiftMaxIterations单次 sifting 最大迭代次数默认值对多数信号够用分解结果出现等幅振荡时可调小Interpolation包络插值方式可选spline或pchip信号毛刺多时选pchip可减少过冲Display是否打印 sifting 进度批量处理时设为 0避免刷屏拖慢速度调用方式统一为imf emd(x, MaxNumIMF, 4)这样的 Name-Value 写法。Display参数在封装成函数批量跑时很有用默认的逐行输出写进日志文件后会让“看运行结果”变成了“翻日志文本”。3.4 从 IMF 到 HHT 谱一张图判断分量是否可解释对分解出的 IMF 做 Hilbert 变换得到瞬时幅值和瞬时频率把所有 IMF 的结果拼成时频图就是 Hilbert 谱。matlab 内置hht直接吃imf矩阵和采样率。fs 1000; [hspec, f, ti] hht(imf, fs, FrequencyLimits, [0 200]); imagesc(ti, f, hspec); set(gca, YDir, normal); xlabel(Time (s)); ylabel(Frequency (Hz)); colorbar;hspec是频率乘时间的幅值矩阵颜色越亮表示该时刻该频率的能量越强。健康信号在 HHT 谱上表现为几条稳定的亮线如果亮线断裂或在相邻频率间来回跳动通常对应两种可能信号本身非平稳或者分解中残留了混叠分量。到这里链路是通的接下来进入所有 EMD 类代码包最核心的调参环节。4. 模态混叠与端点飞翼EEMD 参数和三个高频报错排查4.1 模态混叠为什么会出现问题根源在极值分布被破坏模态混叠的本质是一个 IMF 的 sifting 过程中上下包络被来自其他时间尺度的极值「拉扯」。最典型的是幅值小的连续高频信号叠加在幅值大的低频信号上且高频信号中间存在幅值接近零的间隙。间隙处没有高频极值样条包络在该区间直接跟随低频形态走几个周期后高频分量在间隙段被合并进相邻 IMF 或从残余趋势中消失。主流解法是 EEMD对原信号加多次白噪声后分别做完整 EMD再对同一阶 IMF 取集合平均。白噪声在每一段都提供均匀分布的额外极值填补了间隙使每次分解的包络都遵循同一组统计规律。代价是速度成倍下降且集合平均只能让随机噪声项互相抵消不能保证各阶平均 IMF 严格满足 IMF 定义。4.2 eemd 的两个必调参数Nstd 与 NE 怎么配离线 eemd 包常见调用签名是imf eemd(x, Nstd, NE)。Nstd是附加白噪声标准差与原始信号标准差之比NE是集合平均次数也就是重复 EMD 的次数。Nstd 0.1; % 噪声幅值典型范围 0.02 ~ 0.4 NE 200; % 集合次数常用 100 ~ 500 imf_e eemd(x, Nstd, NE);Nstd太小压不住模态混叠太大则会把噪声本身分解成若干内在模态污染低频 IMF。经验上含冲击振动的机械信号用 0.1 到 0.2生物电信号用 0.2 到 0.3比较平滑的振荡信号用 0.02 到 0.05。NE决定集合平均的统计精度误差随1/sqrt(NE)下降从 100 次加到 200 次误差约降到 70%从 200 加到 800 次只再降一半而耗时线性增长。常见做法是先按NN 100跑一遍看混叠是否消失不够再翻倍而不是一上来配 1000 次。提示eemd 内置了随机数生成同一段信号每次运行结果会有微小差异。做对比实验时先执行rng(2024)固定种子再调用 eemd否则两次结果之间的差异无法归因于参数变化。4.3 端点效应与 sifting 停止准则飞翼从哪来端点飞翼指 IMF 两端出现幅值突然放大的振荡尾巴。原因是样条包络在端点附近缺少外部极值点约束包络走向被首尾几个点强行主导。内置emd已做延拓但信号两端如果起止值差异很大飞翼仍然存在。处理手段是截掉两端振荡段分解完成后按瞬时频率不发散的范围裁剪数据通常直接丢弃每个 IMF 首尾各 1% 到 2% 的采样点再用findpeaks统计 IMF 周期时才不会把飞翼的假极值算进去。停止准则同样影响端点质量。SiftMaxIterations过大时sifting 会把噪声逐个筛成近似等幅振荡的伪分量过小则包络均值没压到足够小。判断方式是观察残差如果残差呈随机毛刺而非光滑趋势说明迭代过度调小迭代次数或换用Interpolationpchip。4.4 三个高频报错NaN、维度与函数找不到运行 zip 包里的脚本报错点通常集中在数据入口和工具箱依赖。先跑一遍下面的自检能把八成问题定位到具体一行。if ~isreal(x) || any(~isfinite(x)) error(信号需要是实数且不含 NaN/Inf); end if size(x, 1) 1 x x(:); % 行向量转列向量 end which emd自检之后对照常见报错现象原因处理Input must be real and finite数据含 NaN、Inf 或复数用rmmissing清洗检查传感器标定是否输出复数格式findpeaks/spline报错或未定义缺少信号处理工具箱或第三方包没进 path先ver(signal)查工具箱再检查which命中文件分解耗时肉眼可见地长信号太长或SiftMaxIterations设置过大先降采样到感兴趣频带的 5~10 倍限制MaxNumIMF同一信号每次结果不一样在线包内部用了随机噪声rng固定种子后重跑最容易被忽略的是第三行第三方 eemd 代码里常嵌套子函数子文件名与函数定义名不一致时addpath能通过主函数入口但内部调用仍会失败。解决方式是用depfun列主脚本全部依赖逐一确认依赖文件确实在解压目录里。5. 让结果可信批量分解与基于白噪声的 IMF 有效性检验5.1 批量处理一个目录下的 csv 信号真实项目里不会有单条信号跑半天的场景更多是把一整批 csv 灌进脚本得到每段数据的分量文件。批量脚本的骨架如下。files dir(samples/*.csv); rng(2024); for i 1:numel(files) data readmatrix(fullfile(files(i).folder, files(i).name)); vec rmmissing(data(:, 2)); % 取第二列并丢弃缺失值 imf emd(vec, MaxNumIMF, 4); save(sprintf(imf_%02d.mat, i), imf, vec); end每次循环先去掉缺失值再统一分解参数。save同时存原始序列和分量后续计算 HHT 谱或统计频带能量时不需要重读 csv。批量跑完先检查各文件的size(imf,2)是否一致如果某些样本只分解出 2 层说明这段信号本身平稳性较好属正常现象。5.2 白噪声统计检验怎么判断 IMF 不是随机噪声分解结果是否可信最便宜的验证是用白噪声生成一段测试序列重复分解并统计各阶 IMF 的能量密度与平均周期的关系。对白噪声做 EMD 后各 IMF 近似为带通滤波输出其能量密度与平均周期的乘积应接近常数。代码实现如下。n 4096; rng(1); wn randn(n, 1); wimf emd(wn); E sum(wimf.^2, 1) / n; % 每个 IMF 的能量密度 Tk zeros(1, size(wimf, 2)); for k 1:size(wimf, 2) [~, loc] findpeaks(wimf(:, k)); Tk(k) 2 * n / max(numel(loc), 1);% 平均周期 2N / 极值点数 end disp([E .* Tk]); % 应近似为常数E .* Tk数值沿各阶 IMF 几乎不变时说明该代码包对随机噪声没有系统性倾向后续对真实信号分解出的低频 IMF 才有统计意义。若乘积出现明显的逐阶漂移优先检查是否漏设了MaxNumIMF导致残差里藏了噪声成分。5.3 边际谱与频带能量统计一键把分解结果变成结论Hilbert 谱在时间维上积分得到边际谱横轴是频率纵轴是该频率在所有时刻的累积能量。这是从「分解出几个分量」到「哪个频带贡献了多半能量」最快的通道。[hspec, f] hht(imf, fs); marginal sum(hspec, 2); % 按频率行求和 figure; plot(f, marginal); xlabel(Frequency (Hz)); ylabel(Marginal amplitude);需要定量指标时在边际谱上分频带求和例如机械故障诊断里把 0~50 Hz、50~200 Hz、200~500 Hz 三个频带的能量占比分别算出作为分解质量与物理意义的双重复核。占异常集中的频带与已知故障特征频率吻合结果才算闭环。本文还有配套的精品资源点击获取