简介:面向非线性、非平稳信号处理学习者的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_ext=2): 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_std=0.2, trials=100, max_imf=10, seed=42): 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], axis=0) 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_{i=1}^{N}E_1\left(x+\beta_0 E_1(n_i)\right)$$
$$r_1=x-IMF_1$$
第 $k$ 阶($k\ge 2$):
$$IMF_k=\frac{1}{N}\sum_{i=1}^{N}E_1\left(r_{k-1}+\beta_{k-1}E_k(n_i)\right)$$
$$r_k=r_{k-1}-IMF_k$$
关键在于算子用的是 $E_1$:对「残差加上第 $k$ 阶噪声分量」这个新信号,只取它的第一阶 IMF,而不是继续分解到第 $k$ 阶。这样每一轮只需要一次「分解到第一层就停」的操作,计算量比 EEMD 的完整分解低,同时噪声分量 $E_k(n_i)$ 的频带与当前残差的主频带对齐,注入的扰动更有针对性。
def ceemdan(x, noise_std=0.2, trials=100, max_imf=10, seed=42): 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, axis=0) 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 适合需要完备重构、或者要做全阶频谱分析的场景,代价是参数更多、实现更复杂。
| 维度 | EMD | EEMD | CEEMDAN |
|---|---|---|---|
| 模态混叠 | 明显 | 显著缓解 | 进一步缓解 |
| 残余噪声 | 无 | 约 $\varepsilon/\sqrt{N}$ | 远低于 EEMD |
| 完备性 | 依赖残差处理 | 平均与丢阶会损失能量 | 重构误差接近浮点精度 |
| 主要参数 | 停止准则、延拓方式 | 噪声幅值、集合次数 | 噪声幅值、集合次数、$\beta_k$ 策略 |
| 计算量 | 单次分解 | $N$ 次完整分解 | 每阶 $N$ 次浅分解 |
| 可复现性 | 高(确定性) | 依赖随机种子 | 依赖随机种子与算子实现 |
验证一套分解结果是否可信,不能只看波形「像不像」。至少跑三项:重构误差,把各阶 IMF 加残差与原信号相减,范数比应在 $10^{-10}$ 以下;正交性指标,计算 $\sum_{i<j}\langle IMF_i,IMF_j\rangle/\langle x,x\rangle$,健康分解通常在 0.05 以下,偏大说明阶间泄漏;已知分量校验,构造一个由两个不同频率正弦加趋势项组成的合成信号,看分解出的分量中心频率是否落在预设值附近,频带是否重叠。三项都过,再拿去分析真实信号,否则后面的物理解释都是空中楼阁。
选型上还有一个常被忽略的取舍:集合次数和噪声幅值是一对耦合参数,单独调一个没有意义。噪声幅值决定「能不能打散混叠」,集合次数决定「打散之后能不能平均回来」,前者不够,后者再多也救不了;前者过大,后者要成倍增加才能压住残余。我一般的做法是先固定 $\varepsilon=0.2$、$N=100$ 跑一遍看 IMF 数目是否稳定,再以 0.05 为步长调 $\varepsilon$,最后按残余噪声目标反推 $N$。
有一次做一组非平稳冲击信号,EMD 和 EEMD 的前两阶几乎一致,第 3 阶开始分叉,当时以为是算法问题,查了两天才发现是延拓深度在两次实验里不一样,端点效应把低频分量污染了。从那以后我每次换实现、换参数,都强制先跑一遍合成信号的重构误差和正交性指标,两项都过再碰真实数据,这个习惯省下的时间比调参本身多得多。这份资料把 EMD、EEMD、CEEMDAN 的推导、伪代码和参数说明串在一条线上,适合对着源码逐段核对上面这些实现细节,需要的话可以取来对照自己的实现逐项排查。希望帮到你。
本文还有配套的精品资源,点击获取