简介这套CEEMD完全集合经验模态分解的MATLAB实现代码带有逐段中文注释主要面向需要分析非线性、非平稳信号的科研人员、算法工程师及高年级本科生也适合希望弄清EMD与CEEMD区别、并将分解算法用于实际数据的读者。资源整体包含4个文件包含1份.m源代码、1份可运行的RAR演示包、1份ZIP备份及1份TXT说明文档总大小约1.03MB。目前已有1753人下载学习。代码注释完整贯穿算法流程从初始化、白噪声注入、上下包络线计算到IMF的迭代提取、残差更新和噪声分离均配有通俗说明同时解释了IMF定义中局部极值点数量与零点数量相等或相差1等要点方便对照程序理解。压缩包内附可运行示例运行后可直接查看原始信号与各IMF分量的分解结果可进一步迁移到地震数据、生物医学信号、机械设备故障诊断等场景。借助这份带注释的程序零基础学习者可以降低入门门槛有经验的开发者也能在此基础上快速修改参数或嵌入自己的项目。1. CEEMD 到底是什么写给卡在模态混叠和噪声残留的人做振动信号分解时直接用 EMD 往往会被模态混叠气得够呛——两个频率相近的成分挤进同一个 IMF 里换 EEMD 加白噪声混叠解决了可重构信号总差一点因为噪声残留进去了。CEEMD 在这条路上多走了一步给信号同时加一对正负白噪声让噪声在集成平均中互相抵消既保留 EEMD 抑制模态混叠的能力又让分解结果可以接近完美重构。本文拆的这份带注释的 MATLAB 完整程序就是干这个的。适合那些做轴承故障诊断、振动分析、非线性时间序列分解的 MATLAB 使用者尤其是被 EMD/EEMD 折腾过、想换更稳定方案的工程师。2. 从 EMD 到 CEEMD这条迭代路线解决了两件事2.1 EMD 的两个步骤筛分 IMF再逐层剥出残差EMD 的思想一句话讲完任何信号都能被分解成若干本征模态函数IMF加上一个残差。IMF 有两个硬性要求过零点数和极值点数最多相差一个上下包络的均值局部为零。算法分内外两层循环内层是筛分外层是剥层。筛分循环的动作是这样的找到信号的局部极大值和局部极小值用三次样条分别拟合出上包络和下包络取两条包络的平均值从信号里减掉这个平均包络重复这个过程直到剩余信号满足 IMF 条件。剥层循环则在每次筛出一个 IMF 后把它从原信号中扣除对剩余部分继续筛下一个 IMF直到剩余信号极值点太少或者趋于单调算法自然结束。这个设计有个天然问题局部极值点的分布直接决定了包络线形态。如果信号里混着间歇性高频成分比如低频正弦上叠加了一个冲击脉冲包络线就会跟着剧烈摆动筛出来的第一个 IMF 往往把不同频带的成分搅在一起这就是模态混叠。模态混叠的典型表现是某个 IMF 的时间轴上先是一段 50 Hz 振荡过一会儿变成 120 Hz 振荡没有办法直接用它提取故障特征频率。2.2 EEMD 用白噪声解决了混叠但引入了噪声残留EEMD 的思路相当巧妙既然单次分解容易混叠那就往信号里反复加白噪声做多次 EMD。高频白噪声会在整个频带上密布极值点让不同尺度的信号成分在筛分时自动分离到不同 IMF 中。每次分解得到一组 IMF最后把所有组结果做集成平均随机白噪声因为统计上不相关而被大幅抵消。但这个方案有一个骗得过眼睛、骗不过数字的问题白噪声统计期望虽然为零有限次集成却不可能完全抵消干净。尤其当噪声幅值取到常见的 0.1 倍信号标准差时平均后残留的噪声会并进各个 IMF最终把重构信号与原始信号的相对误差推到百分之一量级。对需要把各分量精确叠加恢复原始信号的场景比如故障信号幅值本身就是诊断依据这个误差显然没法接受。2.3 CEEMD 的正负噪声对让噪声自己抵消自己CEEMD 的改进在数学上非常干净每次迭代同时生成一条正噪声和一条等幅负噪声分别加到原信号上做两次互相独立的 EMD然后把两个结果平均。表面上看只是多算了一次但这个对称平均让单次迭代的正负噪声在平均时直接抵消——因为正噪声和负噪声符号相反、幅值相同二者相加再除以二随机噪声项迅速收缩。这个方案带来的直接收益是重构误差从 EEMD 的百分之一量级下降到千分之一量级同时 IMF 数量也变得更加稳定。我实际使用时的经验区间是噪声幅值比 Nstd 取 0.05 到 0.2集成次数 NE 取 50 到 200。工程里多数信号在 NE 取 100 左右就能得到很稳定的结果再往上加次数收益很小只会白白拖慢运行时间。这里要特别区分一个概念CEEMD 和 CEEMDAN 不是同一个算法。CEEMD 只做正负噪声对和集成平均CEEMDAN 还会额外对噪声做 EMD在每一层残差上继续加自适应噪声收敛速度更快、模态混叠更少但代码复杂度完全不是一个量级。标题写的是 CEEMD拿到代码包里如果看到有个函数叫 ceemdan先看清楚再调参数这两种方法的 Nstd 和 NE 取值范围和含义不完全一样。3. 核心代码拆解CEEMD 的 MATLAB 实现要点拿到代码包先别急着跑先分清三个文件各自干什么主函数负责噪声对生成和整体集成平均内核函数负责单次 EMD 分解调用脚本负责构造信号、调参和画图。下面按这个顺序逐个拆开讲。3.1 CEEMD 主函数噪声对生成与集成平均主函数是整个程序的骨架核心逻辑集中在三处按信号标准差缩放生成噪声对、单次迭代内部先做一次正负平均、全部迭代结束后再做一次集成平均。参考实现如下function [IMFs, residue] ceemd(data, Nstd, Ne, MaxIter) % CEEMD互补集合经验模态分解 % 输入 % data 待分解信号1xN 行向量或 Nx1 列向量均可 % Nstd 白噪声幅值比推荐 0.05~0.2 % Ne 集合次数推荐 50~200 % MaxIter 筛分最大迭代次数防止死循环默认 500 % 输出 % IMFs 各阶本征模态函数矩阵每行一个 IMF % residue 最后残余分量 if nargin 2 || isempty(Nstd), Nstd 0.1; end if nargin 3 || isempty(Ne), Ne 100; end if nargin 4 || isempty(MaxIter), MaxIter 500; end data data(:).; % 强制转成行向量统一后续处理 N length(data); allImfs cell(Ne, 1); % 每轮结果独立存放 for iter 1:Ne noise randn(1, N) * Nstd * std(data); % 按信号标准差缩放噪声 [imfP, rP] emdOnce(data noise, MaxIter); % 正噪声分解 [imfM, rM] emdOnce(data - noise, MaxIter); % 负噪声分解 nRow min(size(imfP, 1), size(imfM, 1)); % 对齐IMF层数 imfAvg (imfP(1:nRow, :) imfM(1:nRow, :)) / 2; rAvg (rP rM) / 2; allImfs{iter} [imfAvg; rAvg]; % 残差跟在最后 end minRows min(cellfun((x) size(x, 1), allImfs)); % 取所有轮次最小层数 IMFs zeros(minRows, N); for iter 1:Ne IMFs IMFs allImfs{iter}(1:minRows, :); % 集成累加 end IMFs IMFs / Ne; % 集成平均 residue IMFs(end, :); % 最后一层是残差 IMFs IMFs(1:end-1, :); % 去掉残差行 end这段代码最容易被新手忽略的地方是min对齐操作。不同噪声扰动下同一信号分解出来的 IMF 层数可能不一致比如正噪声分解出 7 层负噪声分解出 8 层直接相加会报维度错误。取两者最小层数再对齐是工程上最稳妥的做法代价是偶尔会丢掉多出来的那层实际使用中影响很小。std(data)这行是关键噪声幅值不是固定数字而是跟随信号标准差缩放。这样无论你分析的是幅值 0.01 的采集信号还是幅值 100 的加速度信号Nstd 取 0.1 的含义都是一致的参数不会因为量纲不同而失效。3.2 内层 EMD 筛分循环极值点、样条包络和停止准则emdOnce 是每轮迭代里真正跑 EMD 的内核函数。它的结构是两层循环外层循环剥取一个又一个 IMF内层循环负责让信号满足 IMF 条件。function [IMFs, residue] emdOnce(x, MaxIter) % 单次EMD分解内核 % 输入 % x 待分解信号行向量 % MaxIter 筛分最大迭代次数 % 输出 % IMFs 分解得到的各阶IMF每行一个 % residue 剩余信号 IMFs []; residue x(:).; while true % 剩余信号极大值或极小值不足2个时无法构造包络停止 pksMax findpeaks(residue); pksMin findpeaks(-residue); if length(pksMax) 2 || length(pksMin) 2 break; end h residue; sd 1; cnt 0; % 内层筛分反复减去包络均值直到标准差小于0.2 while sd 0.2 cnt MaxIter cnt cnt 1; % 提取局部极大值和极小值的位置与数值 [vm, lm] findpeaks(h); [vn, ln] findpeaks(-h); vn -vn; if length(vm) 2 || length(vn) 2 break; end % 用三次样条拟合上下包络 t 1:numel(h); upper spline(lm, vm, t); lower spline(ln, vn, t); envMean (upper lower) / 2; hNew h - envMean; % 标准差判据相邻两次筛分结果的差异百分比 sd sum((h - hNew).^2) / sum(h.^2); h hNew; end IMFs [IMFs; h]; residue residue - h; if cnt MaxIter warning(筛分达到最大迭代次数结果可能未收敛); end end end三处细节值得多看一眼。第一findpeaks返回的位置 lm 和 ln 是索引值对应信号长度spline用这些离散点做插值但如果信号的第一个点或最后一个点不是极值点样条包络在边界处会失去约束出现大幅甩尾这就是后面要讲的端点效应。第二停止准则用的是相邻两次筛分的标准差相对变化量阈值 0.2 是行业通用经验值。取太大会过早停止得到的 IMF 不完整频率混叠没除干净取太小会大幅增加迭代次数信号到几千个点时跑起来相当慢。第三如果信号含有直流分量筛分出来的第一个 IMF 会整体偏高这是因为包络均值没有完全消除趋势项。严格的 EMD 实现会在内层循环中额外扣除局部均值漂移但这个精简版靠后续残差把趋势项带走对多数诊断场景够用。3.3 参数怎么选Nstd、NE 的取值区间与调参顺序我每次调参都遵循一个固定顺序先用 Nstd 0.1、Ne 100 跑一版然后看两个指标——重构相对误差和 IMF 波形形态。重构误差小于千分之一说明噪声对抵消得好IMF 波形没有断档、没有突变说明参数合适可以直接用。如果重构误差偏大优先把 Ne 加到 200 再试而不是动 Nstd。原因是 Ne 只增加计算量不改分解的特性Nstd 一动分解结果的形态可能整体变化还得重新验证。Nstd 低于 0.05 时白噪声密度不够模态混叠会回来高于 0.2 时分解结果里能明显看到噪声留下的细碎波纹IMF 变得毛糙。这两个边界我都在真实信号上踩过宁可在 0.05 到 0.15 之间选点。4. 完整实战跑通一次 CEEMD 并看懂结果4.1 构造仿真信号与主脚本为了验证程序正确性和观察分解效果我习惯先造一个多分量叠加信号。这个信号包含一个 50 Hz 主频和一个 120 Hz 的高频成分再加一点白噪声模拟轴承早期故障的振动特征。%% 主脚本构造仿真信号并调用CEEMD clear; clc; close all; rng(42); % 固定随机种子保证结果可复现 Fs 1000; % 采样率 1000 Hz t 0:1/Fs:1; % 时长 1 秒 N length(t); % 仿真信号50Hz 主成分 120Hz 高频成分 白噪声 sig 4*sin(2*pi*50*t) 2*sin(2*pi*120*t) 0.5*randn(1, N); % 调用CEEMD分解Nstd0.1Ne100 [IMFs, residue] ceemd(sig, 0.1, 100, 500); fprintf(分解得到 %d 个IMF\n, size(IMFs, 1)); % 重构检验 recon sum(IMFs, 1) residue; relErr norm(sig - recon) / norm(sig); fprintf(重构相对误差: %.6f\n, relErr);rng(42)这一行很关键。由于 CEEMD 涉及随机噪声如果不固定随机种子每次运行得到的 IMF 都会有些许不同初学者容易被误导以为程序有问题。固定种子后任何人都能复现同样的结果。重构误差的计算放在这里是为了第一时间验证程序基本正确如果 relErr 在千分之一量级说明正负噪声对抵消效果到位如果误差到了百分之一量级先查是不是少加了一行 IMF 或残差再查参数。4.2 IMF 输出格式与结果可视化CEEMD 的返回值 IMFs 是一个二维矩阵每行对应一阶 IMF每列对应一个时间点。第一行是频率最高、能量最强的分量之后逐行频率降低最后一行是残差代表信号的整体趋势。%% 可视化各阶IMF figure(Color, w, Position, [100 100 560 700]); nIMF size(IMFs, 1); for k 1:nIMF subplot(nIMF 1, 1, k); plot(t, IMFs(k, :), LineWidth, 1); ylabel(sprintf(IMF%d, k)); xlim([0 1]); if k 1 title(CEEMD分解结果); end end subplot(nIMF 1, 1, nIMF 1); plot(t, sig, r, LineWidth, 1); xlabel(时间/s); ylabel(原始信号);正常结果应该有这样的特征IMF1 对应 120 Hz 的高频分量IMF2 对应 50 Hz 的主分量IMF3 之后幅值很小、频率极低残差接近直线。如果 IMF1 同时出现 50 Hz 和 120 Hz 的频率切换说明 Nstd 太小模态混叠没有压住要把 Nstd 从 0.1 调到 0.15 再试。我一般不会只画这一张图还会用pspectrum对每个 IMF 单独做频谱图确认各阶 IMF 的中心频率清楚分开。这一步能快速暴露参数问题比只看时域波形直观得多。4.3 重构检验与分解效果评价重构检验除了看相对误差还有一个严格判定把分解得到的全部 IMF 加上残差逐点相加与原始信号做逐点差值查看最大绝对误差出现在哪个位置。端点位置误差偏大属正常但如果在信号中间出现明显尖峰说明某个 IMF 在筛分时卡到了噪声上。另一个常用评价指标是正交性指数。EMD 理论认为各 IMF 之间应相互正交即任意两阶 IMF 的内积接近零。%% 计算正交性指标 io 0; for i 1:size(IMFs, 1) for j i1:size(IMFs, 1) io io abs(sum(IMFs(i, :) .* IMFs(j, :))); end end io io / sum(sum(IMFs.^2)); fprintf(正交性指标: %.4f\n, io);这个指标的值越小越好一般在 0.01 到 0.1 之间。如果超过 0.2说明各 IMF 之间有明显的能量泄漏可能是停止准则太松也可能是噪声幅值取得不合适。我拿到一个分解结果惯例是先跑这个脚本再根据 relErr 和 io 两个数字决定要不要调参。5. 避坑手册CEEMD 最常见的五个翻车现场5.1 端点效应让 IMF 两头炸开现象分解结果在信号开头和结尾几十个点明显发散波形向外甩开一个大尾巴甚至幅度比信号本身还大。原因三次样条包络需要在端点处给出边界约束而信号端点通常不是极值点样条插值缺少约束后会在端点外推时失控这种误差再经过筛分循环放大就形成端点效应。解决我常用的方案是镜像延拓——把信号左右两端各延伸一个窗口长度做镜像对称后参与极值点和包络计算分解完成后再裁剪掉延拓部分。代码里可以在调用 emdOnce 之前对 x 做一次延拓分解后裁剪回来。如果不想改代码也可以在分解前对信号两端各加一个汉宁窗的 5% 衰减减少端点突变。5.2 IMF 数量不一致导致平均错位现象Ne 次集成中有的轮次分解出 7 个 IMF有的分解出 9 个最终平均后某几阶 IMF 波形明显失真幅值比预期小很多。原因CEEMD 主函数用了 min 对齐策略多出的层被丢弃。不同噪声扰动下筛分终止条件分裂层数波动正常但如果波动太大说明 Nstd 取值偏高。解决先查看各轮次的 IMF 层数分布如果大多数在 8 层但偶尔出现 6 层把 Nstd 降低到 0.08 再试。另一个做法是保留最大层数不足层数的轮次用零补齐再平均但这样会在最后几层引入虚假的低频分量我一般不推荐。5.3 噪声幅值一个值打天下就是翻车现象同样的 Nstd 0.1分解一组幅值很大的信号效果很好换一组微弱信号就完全失效IMF 里全是噪声。原因Nstd 是按信号标准差缩放的相对值但它没有考虑信噪比本身。如果原始信号的信噪比极低加入的噪声反而会盖过有效信号筛分出来的第一层可能全是噪声。解决先用snr(sig)或snr(sig, noise)估算信噪比低于 10 dB 时把 Nstd 调到 0.05 左右同时把 Ne 提高到 200靠更多次集成来压噪声。我处理弱故障信号时踩过这个坑调完参数后分解结果明显干净很多。5.4 运行时间太长等不起怎么办现象信号长度 10 万点Ne 取 100跑一次要十几分钟二次调参根本没法做。原因EMD 内核是双重循环里层筛分最坏情况要迭代几百次每次还要做一次样条插值本身计算量就大Ne 取 100 等于是把这一整套重复 200 次。解决三个手段按顺序用——先降采样到 1 万点左右一次性调参确定参数后再全量分解再把 MaxIter 从 500 剪到 200多数信号 50 次内收敛剪短后能省掉一大半计算最后用parfor把 for iter 1:Ne 改成并行循环并行池开满能加速 4 到 6 倍。这三个一起上10 万点的信号能压到两三分钟。5.5 重构误差偏大先查这两个地方现象重构相对误差不是千分之一量级而是到了百分之几明显不对劲。原因多数情况是两层出了问题——要么输出的 IMFs 矩阵被截断了残差行被误删要么调用时返回的其实是去趋势后的信号原始直流分量没参与重构。还有一类原因是代码里 allImfs 每轮存的行数不一致最终平均时 minRows 裁掉了一些有效层。解决先手动执行recon sum(IMFs,1) residue; plot(t, sig - recon)画出逐点误差波形。误差在两端大是端点效应在中间大是层数截断整体是一个常数值则是直流分量丢失。确定类型后再针对性处理别盲改参数。6. 进阶用法用三个指标给分解结果打分在把 CEEMD 结果用于故障诊断之前我总会强制自己花两分钟做一次量化验证。除了前面讲的重构相对误差和正交性指标我还会加一个瞬时频率包络检查对每个 IMF 做希尔伯特变换求瞬时频率曲线看它是否在一个合理频带内平稳波动。%% 检查IMF1的瞬时频率 h1 hilbert(IMFs(1, :)); instFreq Fs * diff(unwrap(angle(h1))) / (2*pi); fprintf(IMF1瞬时频率均值: %.2f Hz\n, mean(instFreq));把三个数字放在一起看重构误差在千分之一以内、正交性指标在 0.1 以下、瞬时频率均值与理论值偏差不超过 5%这个分解结果才值得拿去做特征提取。三个指标里有一个不达标就回炉调参这是我踩了不少坑之后才养成的检查习惯。整套流程走下来CEEMD 的分解质量基本不会再有玄学成分希望帮到你。本文还有配套的精品资源点击获取