简介:这份资源是面向信号处理与数据分析方向的科研人员和工程技术人员整理的CEEMDAN MATLAB仿真程序源码,用于解决非线性、非平稳时间序列分解中EMD易出现模态混合与虚假频率的问题。CEEMDAN通过向信号中注入自适应噪声并迭代平均,能获得更纯净的内在模态函数,适用于工程振动、地球科学、生物医学等领域的复杂信号解析。资源包共6个文件,以4个p文件(核心函数与FFT相关实现)、1个m主程序脚本和1个docx代码说明文档为主,压缩包约210KB,结构紧凑,便于直接调用与二次开发。目前已有2152人学习下载,说明其在同类算法实现中具有较高的参考价值。读者可借助主程序快速跑通分解流程,结合说明文档理解各函数职责与参数含义,并在此基础上替换自己的数据开展实验,掌握从数据预处理、噪声注入、EMD迭代到IMF结果分析的完整思路,提升信号分解与特征提取的实践能力。
1. CEEMDAN MATLAB仿真程序源码:从模态混叠到可复现的分解链路
做非平稳信号处理的人,迟早会撞上模态混叠这堵墙。经验模态分解(EMD)在轴承故障、心电、地震波、语音这些场景里被反复使用,但一旦信号里出现间歇性成分,本征模态函数(IMF)就会串味——同一个时间尺度被拆到两条IMF里,或者一条IMF里混着两个尺度。集合经验模态分解(EEMD)靠加白噪声多次平均缓解了这个问题,代价是计算量和残留噪声。CEEMDAN(Complete Ensemble EMD with Adaptive Noise)在每一阶分解时自适应地加入有限幅值的噪声,并把这一阶的残差作为下一阶的输入,理论上能做到「重构误差几乎为零」且模态对齐。这篇笔记围绕 CEEMDAN MATLAB仿真程序源码展开,讲清楚这套源码能解决什么、核心函数怎么搭、参数怎么调、跑不通时先看哪里。适合已经用过 EMD/EEMD、准备把 CEEMDAN 落到自己数据上的工程师,也适合刚拿到一份 CEEMDAN MATLAB仿真程序源码却不知道从哪下手的新手。
2. CEEMDAN 的算法骨架与源码模块划分
拿到一份 CEEMDAN MATLAB仿真程序源码,第一件事不是急着 run,而是先看清它由哪几块组成。绝大多数实现都遵循同一套骨架:噪声生成与幅值控制、单次 EMD 分解、IMF 的集合平均、残差更新、停止判据。把这五块认全,后面调参和排错才有抓手。
2.1 为什么 CEEMDAN 要在每一阶单独加噪声
EEMD 的做法是「一次性加噪声 → 分解 → 多次平均」,噪声在整个过程中只加一次,平均之后残留噪声会留在低阶 IMF 里。CEEMDAN 换了个思路:第 k 阶分解时,往当前残差里加入的是「经过 EMD 分解后的噪声分量」,而不是原始白噪声。这样每一阶的噪声都是可控的、可被后续阶次抵消的。
用公式描述,设原始信号为 x(t),第 k 阶残差为 r_k(t),第 k 阶要加入的噪声分量记为 E_1(w^i),则第 k 阶的输入是:
r_k(t) + beta_k * E_1(w^i(t))其中 beta_k 是第 k 阶的噪声幅值系数,i 是集合次数索引。对每个 i 做一次 EMD,取第一条 IMF,再对所有 i 求平均,就得到第 k 阶的 IMF。这个「取第一条 IMF 再平均」是 CEEMDAN 和 EEMD 最本质的区别,也是它重构误差小的原因。
源码里通常用一个循环实现这件事,核心变量是Nstd(噪声标准差比例)和NR(集合次数)。这两个参数直接决定分解质量和耗时,后面单独讲。
2.2 源码目录与函数职责
一份结构清晰的 CEEMDAN MATLAB仿真程序源码,一般会拆成下面几个文件。不同作者的命名习惯不同,但职责划分大同小异:
| 文件/函数 | 职责 | 输入 | 输出 |
|---|---|---|---|
ceemdan.m | 主分解函数,控制阶次循环 | 信号、Nstd、NR、MaxIter | IMF 矩阵、残差 |
emd.m | 单次经验模态分解 | 信号、停止判据 | IMF 矩阵 |
getImf.m | 取指定阶 IMF 的辅助函数 | IMF 矩阵、阶次 | 单条 IMF |
plotImf.m | 分解结果可视化 | IMF 矩阵、时间轴 | 图形句柄 |
demo_ceemdan.m | 示例脚本,构造仿真信号并调用主函数 | 无 | 图形与指标 |
主函数ceemdan.m的典型签名是:
function [imf, residue, its] = ceemdan(x, Nstd, NR, MaxIter, SNRFlag) % x : 输入信号,列向量 % Nstd : 噪声标准差相对于信号标准差的比值,常用 0.2 % NR : 集合平均次数,常用 100~500 % MaxIter : 单次 EMD 的最大迭代次数,防止死循环 % SNRFlag : 是否启用 SNR 判据,1 启用 0 关闭 % imf : 分解得到的 IMF 矩阵,每列一条 IMF % residue : 最终残差 % its : 实际分解出的 IMF 条数这里要强调一点:Nstd不是绝对幅值,而是相对信号标准差的比值。很多人第一次改这个参数时直接填 0.2 以为是幅值,结果信号量纲一变就崩。源码里一般会有一行x_std = std(x); noise = Nstd * x_std * randn(...),看到这行就明白它是比例了。
2.3 用仿真信号验证分解链路
在动真实数据之前,先用一个已知成分的仿真信号跑通链路,这是我最推荐的做法。构造一个由三个不同频率分量叠加的信号,看 CEEMDAN 能不能把它们分开:
% demo_ceemdan.m 片段:构造仿真信号并分解 fs = 1000; % 采样率 1000 Hz t = (0:1/fs:1-1/fs)'; % 1 秒时间轴 f1 = 5; f2 = 50; f3 = 120; x = sin(2*pi*f1*t) + 0.6*sin(2*pi*f2*t) + 0.3*sin(2*pi*f3*t); x = x + 0.05*randn(size(t)); % 加一点观测噪声 Nstd = 0.2; % 噪声比例 NR = 200; % 集合次数 MaxIter = 500; [imf, res, its] = ceemdan(x, Nstd, NR, MaxIter, 0); % 重构误差检查 x_rec = sum(imf, 2) + res; err = norm(x - x_rec) / norm(x); fprintf('相对重构误差 = %.3e\n', err);这段代码的逻辑是:先造一个成分已知的信号,调用主函数分解,再把所有 IMF 加残差重构回原信号,用相对误差衡量分解的完备性。正常情况下err应该在 1e-12 到 1e-15 量级,如果跑到 1e-3 以上,说明源码里的残差更新或平均环节有问题。参数上,Nstd=0.2是文献里最常用的起点,NR=200是精度和耗时的折中,MaxIter=500是防止某次 EMD 不收敛。
跑完这一步,再用plotImf或直接plot(imf)看每条 IMF 的波形,前几条应该对应高频分量,后面几条对应低频。如果发现某条 IMF 里明显混着两个频率,那就是模态混叠没压住,需要回头调Nstd和NR。
3. 关键参数怎么设:Nstd、NR 与停止判据
CEEMDAN 的分解质量几乎全压在三个参数上:噪声比例Nstd、集合次数NR、以及单次 EMD 的停止判据。这三个参数没有万能值,但有明确的取舍逻辑。下面把每个参数的物理含义、推荐区间和调整方向讲透。
3.1 Nstd 的取值区间与信号量纲的关系
Nstd控制加入噪声的幅值。太小,噪声不足以触发极值点,模态混叠压不住;太大,噪声会污染分解结果,低阶 IMF 里出现明显噪声成分。经验区间是 0.1 到 0.3,默认 0.2。
| Nstd | 效果 | 适用场景 |
|---|---|---|
| 0.05~0.1 | 噪声残留少,但对间歇成分抑制弱 | 信噪比高、成分平稳的信号 |
| 0.2 | 通用起点,平衡抑制与残留 | 大多数振动、生理信号 |
| 0.3~0.5 | 抑制强,但低阶 IMF 噪声明显 | 强间歇、强非平稳信号 |
调整方向很直接:如果分解后前两条 IMF 里能看到明显的高频毛刺,说明Nstd偏大;如果某条 IMF 里两个尺度还是缠在一起,说明Nstd偏小。我一般会先固定NR=200,在 0.1、0.2、0.3 三个值上各跑一次,对比 IMF 的频谱再定。
有一点容易被忽略:Nstd是相对信号标准差的比值,所以如果输入信号做过归一化,Nstd的绝对效果会变。源码里如果先做了x = x / std(x),那Nstd就相当于绝对幅值了。看源码时务必确认这一步有没有做,否则换一份数据结果就对不上。
3.2 NR 的收敛判断:什么时候 100 次就够了
NR是集合平均次数。理论上越大越好,因为噪声的平均效果随次数增加而提升,但耗时是线性增长的。实际经验是:NR从 100 增加到 200,重构误差和模态对齐改善明显;从 200 到 500,改善有限;超过 500 基本是浪费算力。
判断NR够不够,有个简单办法:固定其他参数,跑NR=50, 100, 200, 400四组,看第 2 条 IMF 的标准差变化。如果从 200 到 400 变化小于 1%,就说明 200 已经收敛。
% NR 收敛性检查 NR_list = [50, 100, 200, 400]; for k = 1:length(NR_list) [imf_tmp, ~, ~] = ceemdan(x, 0.2, NR_list(k), 500, 0); std_imf2(k) = std(imf_tmp(:,2)); end disp(std_imf2);这段代码跑完,如果std_imf2在 200 和 400 之间几乎不变,就可以把NR定在 200。注意每次调用ceemdan内部都有随机噪声,严格来说应该固定随机种子再比较,否则会有小幅波动。源码里如果用了rng控制种子,记得在循环外设一次。
3.3 停止判据:SD 阈值与 MaxIter 的配合
单次 EMD 的停止判据决定一条 IMF 什么时候算「筛干净了」。最常见的是标准差判据(SD):当相邻两次筛分结果的相对标准差小于阈值时停止。阈值一般取 0.2 到 0.3,太小会导致过度筛分、IMF 条数暴增,太大则筛不干净。
源码里通常写成:
% EMD 内部的 SD 停止判据 SD = sum((h_prev - h_curr).^2) / sum(h_prev.^2); if SD < SD_threshold || iter > MaxIter break; endSD_threshold默认 0.2,MaxIter默认 500。这两个要配合看:如果MaxIter经常被触发(可以在源码里加个计数器打印),说明SD_threshold太严或者信号本身极值点太多,这时候要么放宽阈值到 0.3,要么接受更多迭代次数。反过来,如果分解出的 IMF 条数明显多于信号实际成分数,就是阈值太松,往 0.15 收一收。
提示:改停止判据之前,先把
MaxIter调大观察实际迭代次数分布,再决定是改阈值还是改上限,不要两个一起动。
4. 从源码到结果:分解、重构与频谱验证
参数定下来之后,接下来是把分解结果验证一遍。CEEMDAN 的输出是 IMF 矩阵加残差,验证分三步:重构完备性、模态对齐、频谱归属。这三步走完,才能说这次分解是可信的。
4.1 重构误差与 IMF 条数的合理性检查
重构完备性是 CEEMDAN 的招牌,源码正确的话误差应该在机器精度量级。检查方法前面已经给过,这里补充 IMF 条数的合理性判断。
IMF 条数没有固定标准,但有个经验:对于采样率 fs、信号长度 N 的信号,IMF 条数大致在 log2(N) 附近。比如 N=1000,log2(1000)≈10,分解出 8 到 12 条都算正常。如果只有 3 条,说明停止判据太松,很多细节被合并了;如果有 30 条,说明过度筛分,后面对每条 IMF 做分析时会被噪声条数拖累。
% IMF 条数与能量分布检查 [imf, res, its] = ceemdan(x, 0.2, 200, 500, 0); energy = sum(imf.^2, 1); % 每条 IMF 的能量 energy_ratio = energy / sum(energy); bar(energy_ratio); xlabel('IMF 序号'); ylabel('能量占比'); title('各 IMF 能量分布');正常情况下,能量应该集中在前几条 IMF,后面的 IMF 能量递减。如果能量分布很平均,或者某条中间 IMF 能量异常高,往往意味着模态混叠还在。
4.2 用希尔伯特谱确认模态对齐
模态对齐是 CEEMDAN 相对 EEMD 的核心优势,验证方法是做希尔伯特谱。对每条 IMF 做希尔伯特变换,得到瞬时频率,再画时频图。如果模态对齐好,每条 IMF 的瞬时频率应该是一条相对平滑的曲线,不会出现同一时刻多个频率跳变。
% 希尔伯特谱验证模态对齐 fs = 1000; N = size(imf, 1); t = (0:N-1)' / fs; for k = 1:size(imf, 2) h = hilbert(imf(:,k)); inst_freq = fs/(2*pi) * diff(unwrap(angle(h))); subplot(size(imf,2),1,k); plot(t(2:end), inst_freq); ylabel(['IMF', num2str(k)]); end xlabel('时间 (s)');这段代码逐条画瞬时频率。看的时候重点看两点:一是频率曲线是否平滑,二是不同 IMF 的频率区间是否分离。如果两条 IMF 的瞬时频率在某段时间重叠,说明模态对齐没做好,回到第 3 章调Nstd。
4.3 频谱对比:CEEMDAN 与 EEMD 的差异怎么看
如果手头同时有 EEMD 的实现,做个对比会更直观。同样的信号,分别用 EEMD 和 CEEMDAN 分解,对同序号 IMF 做 FFT,看频谱纯度。
| 指标 | EEMD | CEEMDAN |
|---|---|---|
| 重构误差 | 1e-3 ~ 1e-5 | 1e-12 ~ 1e-15 |
| 模态混叠 | 低阶 IMF 偶有残留 | 明显改善 |
| 计算耗时 | 较低 | 约 1.5~2 倍 |
| 参数敏感度 | Nstd、NR 都敏感 | Nstd 敏感,NR 相对宽容 |
这张表不是让你背,而是给你一个预期:CEEMDAN 的代价就是耗时,换来的是重构精度和模态对齐。如果你的应用对实时性要求高、信号又比较平稳,EEMD 可能更划算;如果做的是故障诊断、生物医学信号这类对模态纯度要求高的场景,CEEMDAN 值得多花那点时间。
5. 避坑与排查:CEEMDAN 源码跑不通时先看这几处
这一章是我自己踩过的坑,也是帮别人看 CEEMDAN MATLAB仿真程序源码时最常遇到的五类问题。每条按「现象 → 原因 → 解决」写,照着排查基本能覆盖八成故障。
5.1 报错「矩阵维度不一致」
现象:调用ceemdan时报维度不匹配,或者sum(imf,2)报错。
原因:输入信号是行向量,而源码内部按列向量处理。MATLAB 里randn(size(x))会跟着行向量生成行噪声,后续拼接时维度对不上。
解决:在调用前强制转列向量x = x(:);,或者在主函数入口加一行x = x(:);。这是最高频的坑,没有之一。
5.2 分解结果只有两三条 IMF
现象:跑完发现its只有 2 或 3,明显少于预期。
原因:停止判据SD_threshold设得太大,或者MaxIter太小导致每次 EMD 提前退出。
解决:先把SD_threshold从 0.3 降到 0.2,再把MaxIter从 100 提到 500,重跑观察。如果还是少,检查信号本身是不是太短——N 小于 100 时,极值点不够,分解不出多条 IMF 是正常的。
5.3 重构误差大得离谱
现象:norm(x - x_rec)/norm(x)跑到 1e-2 以上。
原因:残差更新环节写错,常见的是把「当前残差减去当前 IMF」写成了「减去所有 IMF 之和」,或者平均时忘了除以NR。
解决:打开主函数,找到残差更新那一行,确认是residue = residue - imf_k;而不是residue = x - sum(imf,2);。再检查集合平均是不是imf_k = sum(imf_k_all, 2) / NR;。
5.4 每次运行结果都不一样
现象:同样的输入和参数,两次运行 IMF 波形有肉眼可见差异。
原因:源码内部用了randn生成噪声,但没有固定随机种子。
解决:在调用ceemdan之前加rng(1);(或任意固定整数)。如果源码内部自己调了rng('shuffle'),把它删掉或改成固定值。做对比实验时这一步必须做,否则结论不可复现。
5.5 大数据量下内存爆掉
现象:信号长度到十万点以上时,MATLAB 报内存不足。
原因:源码把每次集合的 IMF 都存下来再平均,内存占用是N × NR × 阶数。
解决:改成增量平均,每算完一次集合就累加,不保存全部中间结果。把imf_all = zeros(N, NR);这类预分配改成累加器imf_sum = imf_sum + imf_tmp;,最后除以NR。这样内存占用降到N × 阶数,代价是代码稍复杂一点。
注意:改增量平均时,残差更新也要跟着改,确保每一阶用的是当前阶的平均 IMF,而不是某一次集合的 IMF。
6. 把 CEEMDAN 接到自己的数据上:一个可复用的封装习惯
源码跑通、参数调好之后,真正的工作量在于把它接到自己的数据流里。我自己的习惯是写一层薄封装,把「读数据 → 预处理 → 分解 → 特征提取 → 保存」串成一条链,主函数只负责分解,其余逻辑放在外面。这样换数据时只改封装层,不动核心源码。
封装层的关键是统一输入输出格式。我一般约定:输入是N×1的 double 列向量,输出是一个结构体,包含imf、residue、params、metrics四个字段。params记录本次用的Nstd、NR、MaxIter,metrics记录重构误差和 IMF 条数。这样后面做批量实验时,每个样本的分解配置都可追溯,不会出现「这组结果当时用的什么参数」这种后悔药都买不到的情况。
function result = run_ceemdan_pipeline(x, cfg) % 统一封装:预处理 + 分解 + 指标 x = x(:); x = detrend(x); % 去趋势,避免残差里混入线性项 x = x / std(x); % 归一化,让 Nstd 效果稳定 [imf, res, its] = ceemdan(x, cfg.Nstd, cfg.NR, cfg.MaxIter, 0); x_rec = sum(imf, 2) + res; result.imf = imf; result.residue = res; result.params = cfg; result.metrics.recon_err = norm(x - x_rec) / norm(x); result.metrics.num_imf = its; end这段封装里有两个动作值得说。一是detrend,很多振动信号带线性趋势,不去掉的话趋势会跑到最后一条 IMF 或残差里,干扰后续分析。二是归一化,归一化之后Nstd的相对含义变成绝对含义,换数据时参数不用重调,这是省事的关键。代价是分解结果作用在归一化信号上,如果后续要还原物理量纲,记得乘回std(x)。
调用时这样写:
cfg.Nstd = 0.2; cfg.NR = 200; cfg.MaxIter = 500; result = run_ceemdan_pipeline(x, cfg); fprintf('IMF 条数 = %d, 重构误差 = %.3e\n', ... result.metrics.num_imf, result.metrics.recon_err);批量跑的时候,把cfg存成表格的一行,结果结构体存成.mat,后面做特征工程直接读。这套习惯让我在换数据集时基本不用改核心代码,省下的时间都花在特征和模型上。
最后一个技巧:如果你要对比不同Nstd下的分解效果,别每次手动改参数重跑,写个网格循环,把Nstd和NR的组合遍历一遍,用重构误差和 IMF 条数做筛选条件,自动挑出最优组合。这样一次跑完,比手动试快得多,也避免了我早期那种「试到第五组忘了第一组参数」的翻车。希望帮到你。
本文还有配套的精品资源,点击获取