简介：经验模态分解（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 函数跑通第一组 IMF
3.1 解压后的路径处理：先确认你拿到的是哪套 API
zip 解压后第一件事不是双击运行，而是用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 会把噪声逐个筛成近似等幅振荡的伪分量；过小则包络均值没压到足够小。判断方式是观察残差：如果残差呈随机毛刺而非光滑趋势，说明迭代过度，调小迭代次数或换用Interpolation='pchip'。
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 三个频带的能量占比分别算出，作为分解质量与物理意义的双重复核。占异常集中的频带与已知故障特征频率吻合，结果才算闭环。
本文还有配套的精品资源，点击获取
