简介面向非线性、非平稳信号处理学习者的CEEMDAN算法实现与示例包聚焦经验模式分解中模态混叠问题的改进思路适合已了解EMD、EEMD基础并希望进一步掌握自适应噪声集成方法的研究生、工程技术人员参考。压缩包共5个文件以4个m脚本和1个mat数据文件为主整体约32KB脚本承担算法实现与调用示例mat文件提供可直接载入的实验信号便于对照运行和观察分解结果。该资源已有5935人学习下载说明其在相关课程与课题中具有一定参考价值。读者可从中获得CEEMDAN的迭代分解流程、自适应噪声添加与集成平均的实现脉络并借助心电图等示例数据理解IMF分离、残差计算及参数调整对分解稳定性的影响为地震、生物医学、金融时序或机械故障诊断等场景的算法复现与实验分析提供可操作的起点。1. 模态混叠这道坎为什么 EMD 越「自适应」越难复现同一段轴承振动信号换一台机器、换一份代码、甚至只是把采样起点挪了几百个点EMD 分解出来的 IMF 个数就从 7 个变成 9 个第 3 阶分量时域波形看着还挺像频谱一画却跑到了另一个频带。这种「结果随实现漂移」的现象是 EMD 系算法最劝退的地方也是 EEMD、CEEMDAN 这一整条改进链存在的理由。搞懂这条链等于搞懂「自适应分解」这类方法的能力上限和代价它不依赖预设基函数因此必须用某种方式自己定义「什么是一个模态」而这个定义一旦落到极值点、包络、停止准则这些离散实现上就天然带着敏感性。下面的推演从筛分循环的底层动作开始一路走到 CEEMDAN 的自适应噪声注入把每一步「为什么这么设计」拆开。2. 筛分循环极值点、三次样条与 IMF 判据如何决定分解结果EMD 的全部魔法都压在一个叫筛分sifting的循环里。给定一段实信号 $x(t)$算法不做任何频域变换只反复执行三步找出所有局部极大值和极小值用插值分别连成上包络 $u(t)$ 和下包络 $l(t)$减去两者的均值 $m(t)(u(t)l(t))/2$得到新的 $h(t)x(t)-m(t)$。如此迭代若干次直到 $h(t)$ 满足「过零点数与极值点数相差不超过 1且上下包络均值处处为零」就把这个 $h(t)$ 认定为第一阶 IMF从原信号里减掉对残差重复整个过程。IMF 的定义只有两条极值点和过零点交替出现、上下包络关于零轴对称。它刻意避开了「频率」这个概念因为自适应分解的前提就是不预设频率。2.1 极值点与包络三次样条为什么是「必须」而不是「可选」包络插值的选择直接决定了分解质量。线性插值会让包络出现折角均值曲线跟着抖动筛分次数被迫增加三次样条保证二阶导数连续包络光滑均值曲线也更贴近真实趋势。代价是三次样条在数据两端会做外推端点附近没有极值点约束多项式很容易发散这就是端点效应第一阶 IMF 的两端常常出现大幅摆动越往后分解摆动越往内污染。常见的工程处理是延拓镜像延拓把端点附近的波形对称翻折出去极值延拓在两端各补一个极大值和一个极小值再在延拓后的序列上做样条最后裁掉多余部分。延拓长度一般取一个到两个极值周期太短压不住发散太长会把虚假的周期结构带进包络。这段逻辑在实现里通常长这样import numpy as np from scipy.interpolate import CubicSpline def find_extrema(x): d np.diff(x) # 差分由正转非正记极大由负转非负记极小等号归到前一点避免平台段漏检 idx_max np.where((d[:-1] 0) (d[1:] 0))[0] 1 idx_min np.where((d[:-1] 0) (d[1:] 0))[0] 1 return idx_max, idx_min def envelope_mean(x, n_ext2): idx_max, idx_min find_extrema(x) if len(idx_max) 2 or len(idx_min) 2: return None # 极值不足筛分无法继续 t np.arange(len(x)) # 端点镜像延拓把首尾各 n_ext 个采样点翻折出去抑制样条外推发散 x_ext np.r_[x[n_ext:0:-1], x, x[-2:-n_ext-2:-1]] t_ext np.arange(-n_ext, len(x) n_ext) imax_ext, imin_ext find_extrema(x_ext) upper CubicSpline(t_ext[imax_ext], x_ext[imax_ext])(t) lower CubicSpline(t_ext[imin_ext], x_ext[imin_ext])(t) return (upper lower) / 2.0find_extrema用差分符号判断极值把等于零的情况归到前一个点是为了避免平坦段被重复计数。n_ext控制镜像延拓的深度取值太小端点仍会翘取值太大则会把边界处的虚假极值当成真实结构。envelope_mean返回的是均值曲线筛分循环拿它去减而不是直接拿包络。这里有个容易忽略的细节上下包络必须在同一组时间点上求值插值后的包络才可逐点相减。2.2 筛分停止准则SD 阈值和 S 数到底在约束什么筛分不能无限做下去。迭代次数过多IMF 会退化成调幅的常值序列失去物理意义迭代太少包络均值没压到零模态混叠就留在分量里。早期实现用标准差判据$$SD_k\sum_{t}\frac{[h_{k-1}(t)-h_k(t)]^2}{h_{k-1}^2(t)}$$当 $SD_k$ 落在 0.2 到 0.3 之间就停止。这个阈值的物理含义是「相邻两次筛分结果的相对变化已经足够小」但它对幅值接近零的采样点极其敏感分母一小$SD$ 就爆掉所以后来又有了 S 数判据连续 $S$ 次筛分中过零点数和极值点数都满足「相差不超过 1」就认为分量已经稳定。工程实现里常见做法是两者取与S 数卡住结构SD 卡住幅度任一不满足就继续迭代同时设一个硬上限比如 20 次防止死循环。停止准则的差异是「同一段信号跑出不同 IMF 数」的头号来源。SD 阈值从 0.2 调到 0.3某些分量就会提前停止筛分残余的均值曲线被推到下一阶后面的分量跟着整体偏移。这也是为什么比较两套 EMD 结果之前必须先对齐停止准则否则比较的其实是两种算法。2.3 残差终止什么时候该停手外层循环的终止条件同样影响结果。常见判据是残差单调化当残差里最多只剩一个极值点或者极值点数少于 2 个无法再构造包络分解结束。另一种做法是限制 IMF 总数或者当残差能量低于原始信号的一个比例时停止。限制层数看起来省事但它掩盖了「残差里还有没有可分解结构」这个问题如果残差本身就是一段趋势项强行多分一层只会得到一条没有物理意义的低频摆动。残差处理还有个隐性约定很多实现把最后一次筛分的结果直接当作趋势项 $r(t)$而不是继续当作 IMF。判断标准是残差是否单调或者其极值点数量是否已经不足以支撑一次完整筛分。这个约定不写进论文但直接决定了「重构信号 各 IMF 之和 残差」这条等式能不能成立。完备性是后面所有改进算法都要反复校验的指标而它的第一道关口就在残差终止这里。3. EEMD 的噪声赌局加多少白噪声、集多少次才不是碰运气EEMD 的思路听起来有点反直觉往信号里加白噪声反而能把模态混叠压下去。做法是构造 $N$ 组独立的噪声样本 $n_i(t)$每组按 $x_i(t)x(t)\varepsilon n_i(t)$ 各自做一次完整 EMD得到各自的 IMF 集合再按阶次取平均。噪声在每一组里打乱了信号的极值点分布让原本会挤在同一阶里的不同尺度成分被随机分配到不同的筛分路径上集合平均又把噪声本身的贡献抵消掉留下的就是更稳定的分量。噪声在这里不是干扰而是一种「随机探针」。3.1 白噪声为什么能把混叠「摊平」等效滤波器组的解释白噪声的功率谱在整個频带内平坦意味着它的能量均匀铺在所有尺度上。当它叠加到信号上信号的极值点间隔被随机扰动原本因为极值点间距跳变而混在一起的成分被迫在不同的筛分轮次里分开。更严格的一种解释来自等效滤波器组白噪声在 EMD 下会被分解成一组近似二进的带通分量加了噪声的信号相当于被投影到这组滤波器上每个 IMF 对应一个频带频带之间重叠越小混叠越少。这个解释有个重要推论噪声幅值必须足够大到能改变极值点分布但又不能大到淹没信号本身。幅值太小极值点几乎不变集合平均等于白做幅值太大信号的低频结构被噪声主导IMF 会退化成噪声的分解结果。经验取值是噪声标准差取信号标准差的 0.2 倍左右对于以高频冲击为主的信号可以降到 0.1对于低频缓变信号可以升到 0.3。这个系数没有理论最优解只能按信号的主频分布试。3.2 两个参数定生死噪声幅值与集合次数EEMD 的结果由两个参数主导噪声幅值 $\varepsilon$ 和集合次数 $N$。二者的关系可以用一个近似式描述$$e\frac{\varepsilon}{\sqrt{N}}$$$e$ 是残余噪声的相对误差。想让残余噪声减半集合次数要翻四倍。这解释了为什么 EEMD 慢$N$ 取 100 是常态取 400 也不稀奇每次都要跑一遍完整 EMD计算量线性增长。实际实现里通常固定随机种子让结果可复现否则同一组参数跑两遍IMF 的末尾几位小数都会不一样。import numpy as np def eemd(x, noise_std0.2, trials100, max_imf10, seed42): rng np.random.default_rng(seed) # 固定种子保证集合可复现 sigma noise_std * np.std(x) # 噪声幅值按信号标准差缩放 all_imfs [] for _ in range(trials): xn x sigma * rng.standard_normal(len(x)) imfs emd_once(xn, max_imf) # 单次 EMD返回 IMF 列表 all_imfs.append(imfs) # 按阶对齐后取平均阶数不足的用零补齐避免形状不一致报错 k min(len(imfs) for imfs in all_imfs) return [np.mean([imfs[j] for imfs in all_imfs], axis0) for j in range(k)]noise_std是相对幅值乘的是信号标准差而不是方差避免幅值量纲混淆。trials直接决定残余噪声水平也直接决定耗时。max_imf限制单次分解的层数防止个别噪声样本分解出过多分量导致对齐困难。seed不是可有可无的装饰没有它「同一段信号两次结果不同」会被误判成算法不稳定。3.3 残噪与模态数不一致EEMD 的固有代价EEMD 解决了混叠但引入了两个新问题。第一是残余噪声集合平均只能把噪声压到 $\varepsilon/\sqrt{N}$压不到零低频 IMF 上会残留一层细密抖动做包络谱分析时这层抖动会抬高底噪。第二是模态数不一致不同噪声样本分解出的 IMF 个数可能不同按阶平均时只能取最小阶数多出来的分量被丢弃等于损失了一部分信号能量完备性因此被破坏。这两个问题指向同一个根源EEMD 是在「原始信号」这一层加噪噪声和信号从头到尾绑在一起平均只能消掉噪声的均值消不掉噪声对分解路径的影响。要让噪声只在需要的地方出现、并且能被精确扣除就得把加噪的位置从输入端挪到每一阶残差上这正是 CEEMDAN 做的事。4. CEEMDAN 的逐级扣除把噪声分批注入残差为什么能收敛CEEMDAN 的核心改动只有一句话不再对原始信号加噪而是对每一阶残差加噪并且加进去的噪声以「IMF 分量」的形式出现可以被逐级减掉。它定义了一个算子 $E_k(\cdot)$含义是对给定信号做 EMD 并取第 $k$ 阶 IMF。第一阶 IMF 由 $x\beta_0 E_1(n_i)$ 分解得到减去之后得到残差 $r_1$第二阶 IMF 由 $r_1\beta_1 E_2(n_i)$ 分解得到再减依此类推。噪声不再是一次性铺满全频带而是按阶注入、按阶扣除残余噪声因此被压到远低于 EEMD 的水平。4.1 逐级残差注入CEEMDAN 的递推结构把递推写清楚设原始信号为 $x$第 $k$ 阶残差为 $r_k$$n_i$ 是第 $i$ 组标准高斯白噪声$N$ 是集合次数。第一阶$$IMF_1\frac{1}{N}\sum_{i1}^{N}E_1\left(x\beta_0 E_1(n_i)\right)$$$$r_1x-IMF_1$$第 $k$ 阶$k\ge 2$$$IMF_k\frac{1}{N}\sum_{i1}^{N}E_1\left(r_{k-1}\beta_{k-1}E_k(n_i)\right)$$$$r_kr_{k-1}-IMF_k$$关键在于算子用的是 $E_1$对「残差加上第 $k$ 阶噪声分量」这个新信号只取它的第一阶 IMF而不是继续分解到第 $k$ 阶。这样每一轮只需要一次「分解到第一层就停」的操作计算量比 EEMD 的完整分解低同时噪声分量 $E_k(n_i)$ 的频带与当前残差的主频带对齐注入的扰动更有针对性。def ceemdan(x, noise_std0.2, trials100, max_imf10, seed42): rng np.random.default_rng(seed) sigma noise_std * np.std(x) imfs, r [], x.copy() for k in range(1, max_imf 1): # 每阶重新抽噪声并取该噪声的第 k 阶 IMF 作为注入分量 noise_imfs [emd_once(sigma * rng.standard_normal(len(x)), k)[k-1] for _ in range(trials)] # 对残差加噪后只分解第一层取平均得到当前阶 IMF candidates [emd_once(r nk, 1)[0] for nk in noise_imfs] imf_k np.mean(candidates, axis0) imfs.append(imf_k) r r - imf_k if np.count_nonzero(np.diff(np.sign(np.diff(r)))) 2: break # 残差已单调无法继续筛分 return imfs, rnoise_imfs里对每组噪声单独做 EMD 取第 $k$ 阶这一步是 CEEMDAN 与 EEMD 的分水岭噪声不是直接加而是先被分解成与目标阶同频带的分量。emd_once(r nk, 1)只取第一层是控制计算量的关键。break条件用残差的极值点数量判断和前面 EMD 的残差终止逻辑保持一致。实际实现里β_k常取 $\varepsilon_0\cdot\mathrm{std}(r_{k-1})$让注入强度随残差能量自适应衰减noise_std在这里扮演的是 $\varepsilon_0$ 的角色。4.2 自适应系数 β_k 与算子 E_k$\beta_k$ 的取法直接影响两件事注入的噪声能不能改变残差的极值点分布以及扣除之后会不会留下可见的残余。固定 $\beta$ 的问题在于残差能量随阶数下降固定幅值的噪声到后面会相对过大把低频分量搅乱自适应取 $\beta_k\varepsilon_0\cdot\mathrm{std}(r_{k-1})$ 让噪声强度跟着残差走每一阶的信噪比大致保持一致。算子 $E_k$ 的实现有个坑它要求对噪声样本做至少 $k$ 层分解而不同噪声样本分解出的层数可能不足 $k$。常见处理是取不到就补零或者对该样本重新抽噪声直到层数够。这个细节不写进论文但会让两套实现的第 5、6 阶 IMF 出现肉眼可见的差异。另外$E_k(n_i)$ 本身也是随机的集合次数 $N$ 不够时这个分量的估计误差会直接传进 $IMF_k$表现为高阶 IMF 上残留一层低频起伏。4.3 完备性从哪来残余噪声又去了哪CEEMDAN 的卖点之一是完备性把各阶 IMF 和最终残差加起来能几乎无损地重构原信号。原因在于递推结构本身——每一阶都是从残差里减去一个明确计算出来的分量减法是精确的重构误差只来自浮点运算量级在 $10^{-12}$ 以下。EEMD 做不到这一点因为它的平均过程和残差丢弃会损失能量。残余噪声被压低的机制也不难理解EEMD 的噪声在第一次分解时就混进了所有阶平均只能消均值CEEMDAN 的噪声在每一阶单独注入且注入的是与当前残差同频带的 IMF 分量该分量在下一轮会被当作信号的一部分再分解、再扣除。噪声的脚印被限制在当阶不会跨阶累积。代价是计算结构更复杂每阶都要对 $N$ 组噪声做一次 EMD 取第 $k$ 阶再对 $N$ 个候选做一次 EMD 取第一阶实际耗时不一定比 EEMD 低尤其在 $N$ 取大值时。变体 ICEEMDAN 换了个思路不算 IMF 的平均改算局部均值的差用「残差减去局部均值」代替「残差减 IMF」进一步减少残余噪声和虚假模态。它的实现复杂度和 CEEMDAN 接近参数含义也类似选哪个更多取决于信号本身的频带分布而不是理论优劣。5. 复现排查同一段信号跑出两套 IMF 的五个原因分解结果对不上先别怀疑算法。EMD 系算法的输出是「信号 实现细节」的联合函数任何一处实现差异都会被放大成 IMF 数目或波形的偏移。下面按排查顺序列五条每条都是实际项目里翻过车的。5.1 先固定随机性把「玄学」变成可复现EEMD 和 CEEMDAN 都带随机采样随机种子不固定两次运行的结果在小数点后就不一样集合次数越少差异越明显。排查第一步永远是固定种子并把种子写进配置而不是留在代码默认值里。固定之后再比较参数否则比较的是噪声样本不是算法。第二步是对齐停止准则。SD 阈值、S 数、最大筛分次数、残差终止条件这四项里任何一项不同IMF 的阶数和波形都会变。跨实现对比时先把这四项打印出来核对再谈结果差异。某实验室做过一次对比两套代码的 IMF 数目差了 2 个最后定位到一方把 SD 阈值设成了 0.3、另一方设成了 0.2跟算法本身没关系。第三步是检查采样率和信号长度。端点效应的污染范围与极值周期数相关采样率一变同样的物理时长对应的采样点数变了延拓深度按点数算就会失配。5.2 五条典型踩坑记录现象第一阶 IMF 两端出现大幅摆动越靠近边界幅度越大。原因三次样条在端点外推发散没有做延拓或延拓太浅。解决把镜像延拓深度从 1 提到 2 到 3 个极值周期或者在两端各补一个极值点后再插值然后裁掉延拓段。现象分解出的 IMF 数目比预期少最后一阶看起来像原始信号的趋势项。原因残差终止条件太宽松残差里还有可分解的振荡结构就被判定为单调。解决把终止判据从「极值点少于 2 个」改成「极值点数与过零点数之差连续两次不满足 IMF 条件且残差能量低于阈值」同时检查max_imf是否设得太小。现象EEMD 低频 IMF 上有一层细密抖动包络谱底噪偏高。原因集合次数不足残余噪声停在 $\varepsilon/\sqrt{N}$ 的水平。解决把trials从 50 提到 200 以上或者改用 CEEMDAN 把噪声按阶扣除如果信号本身低频能量弱同时把noise_std降到 0.1。现象CEEMDAN 第 3 阶之后的 IMF 幅值突然变小重构误差却仍然很小。原因β_k用了固定值残差能量下降后噪声相对过大部分信号能量被当成噪声扣掉但因为减法完备重构误差看不出来。解决改成自适应 $\beta_k\varepsilon_0\cdot\mathrm{std}(r_{k-1})$并检查噪声样本分解层数不足时是否补了零。现象同一份数据在两台机器上跑IMF 波形在小数点后第 6 位开始分叉累积到高阶后肉眼可见。原因浮点运算顺序不同并行求和、BLAS 实现差异加上筛分循环对初值敏感。解决固定求和顺序集合平均用手写的顺序累加而不是并行 reduce比较结果时用容差而不是逐位相等。这五条的共同点是它们都不在算法的数学定义里而在实现的边角上。EMD 系算法像半个黑匣子输入和输出之间隔着一串实现选择排查的本质是把这些选择一个个显式化。6. 选型与验证从 EMD 到 CEEMDAN 的取舍边界在哪三档方法不是「越新越好」而是各自对应一类信号条件。EMD 适合极值点分布均匀、端点效应可控的短信号胜在快和可解释EEMD 适合模态混叠严重、但对残余噪声不敏感的场景比如只看前几阶高频 IMF 的故障冲击检测CEEMDAN 适合需要完备重构、或者要做全阶频谱分析的场景代价是参数更多、实现更复杂。维度EMDEEMDCEEMDAN模态混叠明显显著缓解进一步缓解残余噪声无约 $\varepsilon/\sqrt{N}$远低于 EEMD完备性依赖残差处理平均与丢阶会损失能量重构误差接近浮点精度主要参数停止准则、延拓方式噪声幅值、集合次数噪声幅值、集合次数、$\beta_k$ 策略计算量单次分解$N$ 次完整分解每阶 $N$ 次浅分解可复现性高确定性依赖随机种子依赖随机种子与算子实现验证一套分解结果是否可信不能只看波形「像不像」。至少跑三项重构误差把各阶 IMF 加残差与原信号相减范数比应在 $10^{-10}$ 以下正交性指标计算 $\sum_{ij}\langle IMF_i,IMF_j\rangle/\langle x,x\rangle$健康分解通常在 0.05 以下偏大说明阶间泄漏已知分量校验构造一个由两个不同频率正弦加趋势项组成的合成信号看分解出的分量中心频率是否落在预设值附近频带是否重叠。三项都过再拿去分析真实信号否则后面的物理解释都是空中楼阁。选型上还有一个常被忽略的取舍集合次数和噪声幅值是一对耦合参数单独调一个没有意义。噪声幅值决定「能不能打散混叠」集合次数决定「打散之后能不能平均回来」前者不够后者再多也救不了前者过大后者要成倍增加才能压住残余。我一般的做法是先固定 $\varepsilon0.2$、$N100$ 跑一遍看 IMF 数目是否稳定再以 0.05 为步长调 $\varepsilon$最后按残余噪声目标反推 $N$。有一次做一组非平稳冲击信号EMD 和 EEMD 的前两阶几乎一致第 3 阶开始分叉当时以为是算法问题查了两天才发现是延拓深度在两次实验里不一样端点效应把低频分量污染了。从那以后我每次换实现、换参数都强制先跑一遍合成信号的重构误差和正交性指标两项都过再碰真实数据这个习惯省下的时间比调参本身多得多。这份资料把 EMD、EEMD、CEEMDAN 的推导、伪代码和参数说明串在一条线上适合对着源码逐段核对上面这些实现细节需要的话可以取来对照自己的实现逐项排查。希望帮到你。本文还有配套的精品资源点击获取