CEEMD完整MATLAB代码注释：从原理到参数调优

简介：一套带完整注释的CEEMD（完全集合经验模态分解）MATLAB实现代码，适合信号处理领域的研究人员、工程技术人员以及需要学习EMD/CEEMD算法原理的初学者。包内共4个文件，包含可运行的主程序（.m）、源码压缩包、说明文档（.txt）及示例压缩包，整体大小约1.03MB，结构紧凑便于快速部署，在MATLAB环境中可直接运行。已有1753人关注学习，应用场景覆盖地震数据分析、生物信号处理和机械故障诊断等非线性非平稳信号分析任务。代码中逐段注释了初始化、噪声注入、EMD迭代提取、残差计算与噪声分离等关键流程，读者可对照注释理解算法每一步的数学逻辑，并通过可视化IMF分量直接观察分解效果，从而快速上手将该方法迁移到自己的研究或工程项目中；压缩包内还附带可运行的示例包与使用说明，方便直接验证和二次开发。

1. 为什么你需要的是一份带完整注释的CEEMD的MATLAB代码

CEEMD（互补集合经验模态分解）的MATLAB代码，网上流传的版本不少，带完整注释的却很少。很多“完整程序”下载下来，跑通一次，换一个信号就崩，或者分解出来的IMF和论文里的图对不上，最后只能一句一句猜代码在干什么。这篇笔记把一份能直接复制、逐段注释的CEEMD的MATLAB完整程序拆开讲：主函数怎么组织、EMD核心怎么写、正负噪声结果怎么合并、Nstd和NE这些参数怎么调，以及最容易翻车的几个地方。适合刚接触自适应分解、或者手上已有EMD/EEMD代码但想换成CEEMD的人，新手能照着跑通，熟手也能拿来当参数排查清单。

2. CEEMD在信号分解里的位置：原理、选型与边界

CEEMD属于经验模态分解（EMD）家族。EMD把一个非平稳信号拆成若干本征模态函数（IMF）和一个残差趋势项，前提是每个IMF的上下包络均值趋近于零。EMD本身不需要先验基函数，看起来很自由，实际做起来有两个老问题：模态混叠和端点效应。模态混叠最典型的画面，是一个IMF里同时出现高频振荡和低频起伏，明明是两个时间尺度，却被塞进同一个分量里。EEMD用加白噪声再平均的方式缓解模态混叠，CEEMD又在EEMD的基础上，把单侧随机噪声改成正负成对噪声，进一步压低集成平均后留下的噪声尾巴。下面把这条演进线讲清楚，顺便把CEEMD和CEEMDAN这两个经常被混为一谈的名字掰开。

2.1 从EMD到EEMD再到CEEMD：单侧噪声为什么残留，成对噪声为什么能抵消

EMD的筛分过程可以简化成一句话：反复求上下包络的均值并从原信号里减去，直到剩余波形局部对称。这个操作对间歇性的高频分量非常敏感——一段高速振动中间夹着几拍低速摆动，包络会被低速摆动拽偏，高频分量就混进了低频IMF。EEMD的思路是，给信号加多次随机白噪声再做EMD，因为噪声在不同试验里是随机的，集成平均后噪声项互相抵消，模态混叠的毛刺被“抹平”了。代价也很直接：白噪声是单侧加进去的，有限次试验平均后不可能完全消干净，总有一部分随机残留在IMFs里。残留的方差大致和集成次数成反比，想压得越低，NE就要取得越大，计算时间跟着线性上涨。

CEEMD把每个噪声试验拆成正负一对：同一份噪声，一次加到信号上，一次从信号上减掉。两个试验分别做EMD，得到两组IMFs，再按层配对取平均。因为正负噪声完全对称，配对平均时噪声项的对消效率比EEMD高一个数量级，同样的集成次数下，IMF的噪声毛刺明显更小。计算成本则很透明：每做一次集成，CEEMD要跑两次EMD，NE次集成就是2×NE次EMD，比EEMD贵一倍。这个成本换来的好处，是更干净的IMF和更接近原信号的残差。

2.2 CEEMD与CEEMDAN的区别：拿到的代码到底是哪一个

网上搜“CEEMD matlab”，很大概率会搜到CEEMDAN实现，因为这类文章标题经常混写。CEEMDAN的全称是“完全自适应噪声集合经验模态分解”，它每一层分解之前，先计算当前残差的标准差，再把自适应系数乘上白噪声加到残差里，一层一层提取IMF。CEEMD相对简单：它始终对原始信号加减成对噪声，不逐层调整噪声强度。判断一份代码到底是哪一个，不需要看标题，看循环里有没有“计算残差标准差再乘系数”这一步就行。如果有，那是CEEMDAN；如果只是每轮生成一份噪声、加正的和加负的分别做EMD，才是这里要讲的CEEMD。

这两个方法没有绝对的优劣。CEEMDAN收敛通常更快，噪声残留也更小，但多了一套逐层噪声系数，调参时解释起来很麻烦，不同论文对系数序列的取法也不一致。CEEMD概念更直白：正负噪声对消、集成平均、残差由重构关系兜底，每一步都能在代码里指出来。如果你只是想把手头的EEMD升级成结果更干净、注释更完整的方案，CEEMD更合适。

2.3 什么信号适合CEEMD：选型判断与不适合的三种情况

CEEMD适合非平稳、非线性、且怀疑有模态混叠的信号。常见落地场景包括机械振动与轴承故障特征提取、心电和脑电等生理信号分析、风速和地震动数据处理，以及把IMF的能量占比作为特征输入机器学习模型。这类信号有一个共同点：没有固定的周期结构，用傅里叶拆会拆出一大片连续的频带，用CEEMD拆出来的IMF才有物理意义。

不适合的情况我也遇到过不少。第一种是信号太短，只有几十个点，包络插值根本撑不起来，分解出一堆伪IMF。第二种是强实时场景，CEEMD的2×NE次EMD串行计算太慢，NE取100时千点信号就得跑几十秒，做在线监测要考虑降级方案。第三种是分量频率极度重叠，两个分量的瞬时频率几乎贴合在一起，单靠EMD家族很难切开，这时候该先做解调或者直接换同步挤压变换。选型判断放在最前面，能省掉后面一堆调参时间。

3. 加注释完整程序：主函数ceemd.m与三层循环结构

这一章直接给代码。整个程序按“一个函数文件ceemd.m”来组织，主函数在最上面，局部函数跟在后面。调用方式只有一行：[imf, residual] = ceemd(x, Nstd, NE, MaxIter, seed)。主函数要干的事，是把CEEMD的三层循环理清楚：外层遍历集成次数，中层做正负两个噪声试验，内层是每轮试验内部的EMD筛分。层数对齐的问题也在主函数里处理，不然不同试验出来的IMF数量不一致，后面根本没法求平均。

3.1 主函数ceemd.m的结构：三层循环与输入输出

function [imf, residual] = ceemd(x, Nstd, NE, MaxIter, seed) % CEEMD 互补集合经验模态分解（Complementary Ensemble Empirical Mode Decomposition） % 输入： % x - 待分解信号，行向量或列向量均可 % Nstd - 辅助白噪声幅值，相对原始信号标准差的比例，默认 0.2 % NE - 集成次数，正负噪声各跑 NE 次，默认 100 % MaxIter - 单次 EMD 的最大筛分次数，默认 100 % seed - 随机种子，传入后结果可复现，默认随机 % 输出： % imf - K x N 矩阵，K 为 IMF 层数，每一行是一个 IMF % residual - 1 x N 残差，满足 x = sum(imf,1) + residual x = x(:).'; % 统一成行向量，避免行列向量的歧义 N = length(x); if nargin < 2 || isempty(Nstd), Nstd = 0.2; end if nargin < 3 || isempty(NE), NE = 100; end if nargin < 4 || isempty(MaxIter), MaxIter = 100; end if nargin < 5 || isempty(seed) rng('shuffle'); % 不传种子，每次结果不同 else rng(seed); % 固定种子，保证可复现 end % 先用无噪信号跑一次 EMD，得到基准层数 K0 [imf0, ~] = emd_core(x, MaxIter); K0 = size(imf0, 1); imf_avg = zeros(K0, N); % 积分量，最后除以 NE for k = 1:NE % 生成一份标准正态白噪声，再乘幅值 noise = Nstd * std(x) * randn(1, N); % 正噪声试验：x + noise [imf_p, res_p] = emd_core(x + noise, MaxIter); % 负噪声试验：x - noise [imf_n, res_n] = emd_core(x - noise, MaxIter); % 把正负试验的 IMF 都对齐到 K0 层 imf_p = align_to_K(imf_p, res_p, K0); imf_n = align_to_K(imf_n, res_n, K0); % 正负对消平均，再累加到总积分里 imf_avg = imf_avg + (imf_p + imf_n) / 2; end imf = imf_avg / NE; % 跨试验集成平均 % 残差直接由重构关系得到，避免多轮平均把误差带进残差 residual = x - sum(imf, 1); end

这段代码的核心在第20行到第30行的循环里。每次集成先生成一份噪声，正负各做一次EMD，然后把两层结果平均。配对平均这一步就是CEEMD消除噪声残留的关键：正试验和负试验用的是同一份噪声，符号相反，叠加后噪声项在理论上被对消，剩下的才是信号本身的分量。最后一行残差不是直接取平均残差，而是用x - sum(imf,1)倒推。这样无论中间对齐怎么折腾，重构误差永远受限，不会越积越大。

需要提醒的是，nargin用来处理省略参数的情况，前四个参数都有默认值。实际使用时，我一般建议至少把seed显式传进去，方便横向对比不同信号的分解结果。

3.2 Nstd、NE、MaxIter如何设：参数表与调试顺序

调试顺序比具体数值更重要。第一次跑通，用Nstd=0.2, NE=50, MaxIter=100快速看结果，确认分解层数和大致形态没问题，再把NE调到100到200之间用于正式实验。NE超过200之后，噪声对消带来的改善会越来越不明显，我做过实验对比，NE从100加到400，IMF波形差异肉眼几乎看不出来，耗时却翻了四倍。MaxIter要结合信号长度看：信号越长，单次筛分需要更多迭代才能把包络均值压到阈值以下，给100基本够用；如果发现分解出的IMF有明显“振荡锯齿”，可以把MaxIter提到150再观察。

3.3 align_to_K：层数不一致时，多出来的IMF去哪了

function imf_out = align_to_K(imf_in, res_in, K0) % 将一次 EMD 的 IMF 矩阵对齐到基准层数 K0 % 超过 K0 的分量叠加回残差，不足 K0 的用残差补齐 imf_out = imf_in; N = size(imf_out, 2); K = size(imf_out, 1); if K > K0 % 多出来的分量并回残差，保证能量不凭空消失 res_in = res_in + sum(imf_out(K0+1:end, :), 1); imf_out = imf_out(1:K0, :); elseif K < K0 % 缺的层用这次试验的残差填充 imf_out(end+1:K0, :) = repmat(res_in, K0-K, 1); end end

这个局部函数解决的是CEEMD实现里绕不开的痛点：每次EMD做出来的IMF数量不一样。加了噪声的试验，可能比无噪基准多拆出一层，也可能少拆一层。如果不处理，正负试验的IMF矩阵行数对不上，配对平均直接报错。对齐策略是“多退少补”：多出来的IMFs叠加到残差里，缺的层用残差填充。为什么要把多出来的叠回残差而不直接扔掉？因为CEEMD最终要求sum(imf,1) + residual = x，扔了任何分量都会破坏这个关系，叠回残差至少保证能量还在。补齐的逻辑虽然简单，但它是整个程序能稳定运行的前提，没有这层处理，NE稍微大一点就会遇到矩阵维度不匹配的报错。

4. EMD核心子函数与正负噪声结果的层数对齐：从包络到残差

主函数再完整，EMD核心没写对，CEEMD也是白搭。CEEMD里的每一次试验，本质上都要调用一次完整的EMD分解，所以emd_core和envelope这两个局部函数决定了分解质量。这一章把它们单独拎出来讲，顺便交代怎么把上一章的主函数和这一章的子函数拼成一个能直接运行的文件。

4.1 emd_core.m：一次EMD的筛分循环与停止准则

function [imf, residual] = emd_core(x, MaxIter) % 单次 EMD 分解核心，供 ceemd 内部调用 x = x(:).'; imf = []; residual = x; while true h = residual; % 当前待筛分的信号 for iter = 1:MaxIter h_prev = h; [upper, lower] = envelope(h); % 上下包络 m = (upper + lower) / 2; % 包络均值 h = h - m; % 减去均值，向局部对称逼近 % 标准差停止准则：两次筛分之间的变化足够小就停 sd = sum((h_prev - h).^2) / (sum(h_prev.^2) + eps); if sd < 0.2 break; end end imf(end+1, :) = h; % 当前筛分结果作为一个 IMF residual = residual - h; % 扣掉已经提取的分量 % 残差极值点太少，无法继续构造包络，停止分解 if sum(islocalmax(residual)) + sum(islocalmin(residual)) < 2 break; end end end

筛分循环里发生的事，可以理解成把波形一层层“削平”：每次算出上下包络的均值，从当前信号里减掉，让剩余波形越来越接近零均值、对称的形状。sd是衡量两次筛分结果差异的指标，阈值0.2是EMD文献里最常见的经验值。这里用eps防止除零，信号静默时不会报错。外层while循环每次提取一个IMF，然后看残差还剩几个极值点。极值点少于2，说明剩下来的基本是单调趋势项，没有办法再做包络了，分解就到此为止。

islocalmax要求MATLAB R2017b以上版本。如果用的是老版本，可以自己写一个极值查找：idx = find(diff(sign(diff(h))) == -2) + 1;效果接近，但斜率和重复点较多的信号要额外处理。我一般优先用islocalmax，新版本环境基本都支持。

4.2 envelope.m：三次样条包络与端点补偿

function [upper, lower] = envelope(h) % 用三次样条分别拟合上下包络，并做端点补偿 N = length(h); pMax = find(islocalmax(h)); pMin = find(islocalmin(h)); if numel(pMax) < 2 || numel(pMin) < 2 % 极值点不够，无法插值，直接返回原信号，筛分会自动停止 upper = h; lower = h; return; end % 端点补偿：把两端点作为伪极值加入，抑制三次样条外推产生的飞翼 pMax = [1, pMax, N]; pMin = [1, pMin, N]; % 三次样条插值，得到连续包络 UP = spline(pMax, h(pMax), 1:N); LOW = spline(pMin, h(pMin), 1:N); % 防止数值误差导致上下包络交叉 upper = max(UP, LOW); lower = min(UP, LOW); end

包络是EMD的“地基”。每个IMF筛分步骤都要先求上下包络，包络扭曲，整条IMF就废了。MATLAB的spline函数在数据点之间做三次样条插值，点数足够时效果好，但数据点之外的外推行为不可控，这就是端点飞翼的根源。代码里的端点补偿，是把序列首尾两个点直接当成伪极值加进插值点集合，让包络在两端被“按”住，而不是由最外面的极值点斜率任意外推。代价是端点附近的包络会被拉向信号本身的值，可能轻微影响首尾几个数据点的IMF幅值，但整体比飞翼可控得多。

upper = max(UP, LOW)这行不是多余的。样条插值在某些畸变点附近可能出现上下包络交叉，交叉后包络均值会偏离真实中心，后面筛分就会越算越偏。做一次取大取小，至少保证上包络不低于下包络。

4.3 合并成一个ceemd.m：文件组织与调用示例

把第3章的主函数、align_to_K，以及这一章的emd_core、envelope按顺序放进同一个ceemd.m文件里，主函数在最上面，其他作为局部函数跟在后面，保存后就能直接调用。局部函数不会暴露到命令行，也不会和别的同名函数冲突。调用方式如下：

fs = 1000; t = (0:999) / fs; x = sin(2*pi*40*t) + 0.6*sin(2*pi*5*t) + 0.2*randn(1, 1000); [imf, res] = ceemd(x, 0.2, 20, 100, 2024); % 查看分解层数 disp(size(imf, 1)); % 验证重构误差 err = norm(x - (sum(imf, 1) + res)) / norm(x); fprintf('重构误差: %.3e\n', err);

先用NE=20跑通，确认文件、调用、重构误差都正常，再加大NE做正式实验。重构误差如果明显大于机器精度，说明某个环节的层数对齐或残差计算写错了，这时候别急着调参数，先回头查代码。

5. CEEMD参数调优与常见问题排查：五个踩坑记录

这一章写给已经跑通程序、开始认真调参数的人。下面五条都是实际使用中反复出现的坑，每条按“现象、原因、解决”的顺序写，方便你在出问题时快速定位。

5.1 噪声幅值和集成次数不匹配：结果全是噪声毛刺

现象：IMF第一层全是细碎的毛刺，看不到清晰的振荡包络，把NE从50加到200，毛刺只是变短一点，并没有消失。原因是Nstd太大，或者NE太小。Nstd取到0.4以上时，白噪声在配对平均后仍会留下肉眼可见的痕迹，尤其高频第一层IMF对噪声最敏感。解决方法是把Nstd压回0.1到0.2之间，同时保证NE至少50。调试顺序我一般固定为：先设NE=50，Nstd从0.1开始往上加，每换一个Nstd都看第一层IMF的毛刺程度；找到一个临界点后，再调NE把残余噪声进一步抹平。不要一上来就NE=500，那不是调参，是烧机器。

5.2 端点飞翼：三次样条包络的边界失控

现象：分解结果两端剧烈张开，像一个甩开的翅膀，中段看起来正常，两端完全不能用。原因是三次样条在最外侧数据点以外的外推不受控制，尤其在信号端点恰好不是极值点时，包络会按最后一段的斜率冲出去。解决方法是上一章envelope里写的端点补偿，把首尾两点作为伪极值加入插值。如果信号两端摆动依然明显，可以把端点补偿换成镜像延拓：在信号两端各镜像一段局部极值，再做样条插值，这样首尾极值更自然。镜像延拓会增加一点代码量，但处理强非平稳信号时值得一试。

5.3 集成次数盲目加大：算力燃烧，收益趋近于零

现象：NE从200加到800，跑了几个小时，画出来的IMF图和NE=200时几乎一模一样。原因是CEEMD的成对噪声已经把随机残留压得很低，这个阶段的主要误差来源是单次EMD的包络不稳定和模态混叠，不是白噪声残留，加NE解决不了。解决方法是区分“随机噪声残留”和“系统性分解误差”。前者表现为IMF粗糙、带毛刺，可以用NE压；后者表现为IMF波形稳定但明显混叠，比如一个IMF里同时有快慢两个节奏，这时候应该调Nstd或MaxIter，而不是无脑加NE。我现在的习惯是先NE=100看结果，不满意就回到Nstd和MaxIter找问题。

5.4 信号太短还想分解：极值点不够，结果不可信

现象：一段长度不到100的信号，CEEMD跑出了七八层IMF，每一层都弯弯曲曲，看着很丰富，实际上毫无物理意义。原因是极值点太少，三次样条在稀疏极值点之间强行插值，造出很多伪包络，筛分过程也受影响，于是层层分解出现大量虚假模态。解决方法是分解前先做个快速检查：sum(islocalmax(x)) + sum(islocalmin(x))如果小于6，直接放弃CEEMD，别硬跑。至少要有几百个数据点、十几个极值点，分解结果才谈得上可信。如果数据本质上就短，更合理的做法是接受限制，改用更适合短序列的方法，而不是把CEEMD的阈值调低强行分解。

5.5 结果不可复现：随机种子没管住

现象：同一段信号，同一个参数，跑两次的IMF不完全一致，画图对比时能看到细微差异，写报告时不知道该信哪一次。原因是randn产生的随机噪声每次都不同，CEEMD的集成平均过程本质上是随机算法，不固定种子就没有可复现性。解决方法是传seed参数，并且在每次实验前保存再恢复随机状态：

s = rng; % 保存当前随机状态 [imf, res] = ceemd(x, 0.2, 100, 100, 2024); rng(s); % 恢复，避免影响后续代码

这里有个容易忽略的细节：我的ceemd内部会调用rng改全局随机状态。如果不恢复，同一脚本里后续其他随机操作就会被这次分解强行改写，间接导致后面一步的蒙特卡洛结果复现不了。所以不只是传seed进去，还要在调用后恢复现场。固定seed后，同一参数同一输入得到同一结果，写论文和调bug都会省心很多。

6. 进阶验证：三招确认CEEMD分解结果真的可用

跑通程序只是第一步，真正要命的是“跑出来了，但不知道对不对”。下面三个验证手法，我每次做CEEMD都会走一遍，能挡住大部分隐性错误。

第一招，用合成信号检验重构误差和正交性。构造一个已知频率成分的信号，分解后先算重构误差，再算IMF之间的正交性指数。重构误差反映整个程序的自洽性，正交性指数反映IMF是否真的把不同时间尺度分开：

fs = 1000; t = (0:999) / fs; x = 1.2*sin(2*pi*50*t) + 0.8*sin(2*pi*5*t) + 0.2*randn(1, 1000); [imf, res] = ceemd(x, 0.2, 100, 100, 2024); % 重构误差：正常应该接近机器精度 err = norm(x - (sum(imf,1) + res)) / norm(x); % 正交性指数：IMF之间相关性越低越好 K = size(imf, 1); io = 0; for i = 1:K for j = i+1:K io = io + abs(sum(imf(i,:) .* imf(j,:))); end end io = io / (sum(sum(imf.^2, 2)) + eps); fprintf('重构误差: %.3e\n正交性指数: %.3f\n', err, io);

第二招，用EEMD对照看噪声残留。同一信号分别用EEMD和CEEMD分解，对比第一层IMF的功率谱：CEEMD第一层里应该少一些宽频抬升的噪底，频谱更干净。第三招，用Hilbert检查最高频IMF的瞬时频率是否平滑。瞬时频率剧烈跳变，说明模态分离不彻底，里面还混着别的成分。

我自己最早用CEEMD时，只盯着IMF图看形态，结果有一版代码残差计算错误，分解出来的IMFs看起来挺好，重构误差却到了0.1量级，画图根本看不出来。后来养成习惯：每次跑完先看重构误差，再看第一层IMF的频谱，最后把seed连同参数一起存进.mat文件。这三个动作用不了两分钟，却救了我好几次。希望帮到你。

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