简介：这份资料包围绕CEEMD完全集合经验模态分解展开，面向从事非线性、非平稳振动信号分析与特征提取的MATLAB学习者，帮助完成从IMF分量筛选到信息熵特征构建的完整实验流程。压缩包共4个文件，含3个m脚本与1个xlsx数据表，体积约60KB；脚本对应分解、极值处理与特征计算，表格提供示例信号数据，便于直接运行和对照。已有5248人学习，说明该流程在信号处理练习中具有一定参考价值。读者可获得CEEMD分解与IMF相关系数筛选的实现思路，理解用相关系数剔除冗余分量、再计算信息熵以刻画复杂度与随机性的做法，并看到如何将筛选后的分量与熵特征组织成特征向量矩阵，供后续分类模型使用。整体结构紧凑，适合具备MATLAB基础、希望把分解、特征工程与机器学习预处理串起来的学习者参考。
1. 加噪声反而更稳：CEEMD 分解后的三个数字才值钱
往原始信号里主动掺一段白噪声再分解，结果反而更稳定、更可复现——CEEMD（互补集合经验模态分解）的第一行 MATLAB 代码就在做这件反直觉的事。它要治的是 EMD 的老毛病：模态混叠。间歇性的高频冲击会让极值点分布突变，包络随之变形，同一频率成分被劈到两个 IMF 里，或者一个 IMF 里塞进两个尺度。
但分解只是前戏。真正决定这套特征能不能用的，是后面两步：用相关系数把噪声主导的 IMF 分量剔掉，再用信息熵把剩下的分量压成分类器能吃的向量。整条链路「CEEMD 分解 → IMF 相关系数筛选 → 信息熵特征」在 MATLAB 里不到百行，可每一个阈值都在悄悄决定特征里还剩多少信息量。
2. EMD 到 CEEMD：IMF 的筛分条件与互补加噪的重构恒等式
2.1 本征模态函数的两个硬条件
IMF 不是随便切一段信号就叫的，它必须同时满足两个条件。第一，整段数据内局部极值点数目（极大值加极小值）与过零点数目之差不超过 1。第二，任意时刻，由局部极大值点拟合出的上包络与局部极小值点拟合出的下包络，其均值为零。
第一个条件是在近似保证窄带性。过零点刻画振荡本身，极值点刻画包络的起伏，两者数量匹配意味着信号围绕零轴做单一尺度的振荡，不会出现「大周期里套着小周期」的复合结构。第二个条件是在保证局部零均值，把骑在慢变趋势上的偏置剥掉——筛分过程反复减均值，做的就是这件事。
这两个条件在实际实现里都带容差，不可能严格成立。判断函数通常写成「极值点数与过零点数之差 ≤ 1 且包络均值绝对值小于某个阈值」，阈值定得松，IMF 里残留趋势；定得紧，筛分次数暴涨。
2.2 三次样条包络与停止准则
一次完整的筛分只有五步，但每一步都藏着实现细节。先找出当前序列的全部局部极大值和局部极小值，再用三次样条分别插值成上、下包络，取两者均值为 m(t)，用 h(t) = x(t) - m(t) 得到候选分量；判断 h 是否满足 IMF 条件，不满足就把 h 当成新的输入再筛一遍。
function [imf, res] = emd_core(x, maxIMF, siftMax, sdTh) % 极简筛分核心，只表达流程，端点延拓与极值点保护从略 res = x(:)'; imf = zeros(maxIMF, numel(x)); for k = 1:maxIMF h = res; for s = 1:siftMax [idxMax, idxMin] = find_extrema(h); % 局部极大/极小值位置 upper = spline(idxMax, h(idxMax), 1:numel(h)); % 三次样条上包络 lower = spline(idxMin, h(idxMin), 1:numel(h)); % 三次样条下包络 m = (upper + lower) / 2; hNew = h - m; sd = sum((h - hNew).^2) / (sum(h.^2) + eps); % Huang 的 SD 准则 h = hNew; if sd < sdTh, break; end % 满足阈值提前退出 end imf(k, :) = h; res = res - h; % 残余分量继续分解 endmaxIMF一般取floor(log2(N)) - 1量级，再多就会分解出没有物理意义的极低频分量；siftMax是单次筛分的次数上限，典型值 8 到 10；sdTh就是下面要说的 SD 阈值。三个参数里最容易踩的是sdTh。
停止准则主流有三种。Huang 提出的 SD 准则用相邻两次筛分结果的归一化能量差做判据，SD 取 0.2 到 0.3 比较常见；S 次准则更直观，连续 S 次（典型 3 到 5 次）筛分结果都满足「极值点数等于过零点数」就停；还有一种是干脆固定筛分次数，实现最简单但一刀切。
筛分过度是个隐蔽的翻车点。继续筛下去，IMF 会被榨成纯粹的调频信号，幅值几乎恒定，包络均值恒为零，数学上完美满足 IMF 条件，物理意义却没了——它不再对应任何真实的振动模式。这种分量算出来的能量熵会异常偏小，看起来「很干净」，实际上是在给分类器喂假特征。
端点效应是另一个必须处理的问题。三次样条需要极值点作为插值节点，端点外侧没有节点，样条会向外发散，形成首尾剧烈摆动的「飞边」。这个误差不会只停留在第一层，它会进入残余分量，逐层向内传播。常见的补救是镜像延拓、极值点延拓，或者用 AR 模型预测端点外的走势，目的都是给样条补出合理的边界节点。
2.3 模态混叠到底是怎么发生的
设想一段振动信号，绝大部分时段只有低频成分，但在 t = 0.3s 附近出现一次高频冲击。做筛分时，高频冲击贡献的极值点只出现在那一小段窗口里，上包络在其余时段完全由低频极值点决定。减均值的时候，高频冲击要么被包络「拉平」，要么被劈成两半分到相邻两个 IMF 里。
根源在于 EMD 的尺度分离完全依赖局部极值点分布，没有任何全局频率基准。极值点分布一变，包络就变，筛分路径就变。同一个信号加一点微小扰动，分解出来的 IMF 阶数都可能不同——很多人觉得 EMD 有点玄学，玄学感就是从这来的。EEMD 的思路是：既然极值点分布不稳，那就人为把它「抹匀」。
2.4 EEMD 的残余噪声与 CEEMD 的互补对
EEMD 每次加入不同 realization 的白噪声，各自分解，再把同阶 IMF 做集成平均。噪声的作用是给每个时间点都补充极值点，让包络在统计意义上更均匀。代价是平均之后噪声不会完全消失，残余幅值大约是 ε/√M（ε 为噪声幅值，M 为集成次数）。M = 100 时残余还有 10%，想把误差压到 1% 需要一万次分解，计算上不现实。
CEEMD 换了个思路：不靠增加次数，靠对称性。每次构造一对信号 x + εn 和 x - εn，分别分解，再把同阶 IMF 取平均。噪声在正负两个版本里相位相反，相加时奇次项直接抵消，偶次项被平均削弱。同等集成次数下，重构误差比 EEMD 低一到两个数量级。
function [imfs, res] = ceemd(x, NE, noiseAmp, maxIMF, siftMax, sdTh) % NE: 互补对数，每对包含 +n 和 -n 两次独立分解 N = numel(x); imfSum = zeros(maxIMF, N); sigma = std(x); for e = 1:NE n = randn(1, N) * noiseAmp * sigma; % 噪声幅值按信号标准差缩放 [imfP, ~] = emd_core(x + n, maxIMF, siftMax, sdTh); [imfM, ~] = emd_core(x - n, maxIMF, siftMax, sdTh); imfSum = imfSum + imfP + imfM; % 正负对同阶相加，噪声奇次项抵消 end imfs = imfSum / (2 * NE); res = x - sum(imfs, 1); % 余项理论上只剩趋势noiseAmp是相对噪声幅值系数，典型 0.1 到 0.2 倍标准差。太小（低于 0.005 倍）改变不了极值点分布，CEEMD 直接退化成 EMD；太大会引入虚假分量，低阶 IMF 里全是噪声残余。NE是互补对数，50 到 200 之间，计算时间随它线性增长，因为每对要做两次完整 EMD。
| 参数 | 含义 | 典型取值 | 调错的后果 |
|---|---|---|---|
| noiseAmp | 噪声幅值 / 信号标准差 | 0.1 ~ 0.2 | 太小退化为 EMD，太大产生虚假分量 |
| NE | 互补对数 | 50 ~ 200 | 太小残余噪声不消，太大计算爆炸 |
| maxIMF | 最大 IMF 层数 | floor(log2(N))-1 | 过多产生无意义低频分量 |
| siftMax | 单次筛分上限 | 8 ~ 10 | 过大把 IMF 榨成纯调频信号 |
| sdTh | SD 停止阈值 | 0.2 ~ 0.3 | 过小导致过筛分 |
提示：
randn不固定种子的话，每次运行结果都不一样，调试阶段很难判断是参数问题还是随机性。正式跑之前先rng(0)固定种子，保证可复现。
3. 相关系数筛 IMF：噪声主导分量的判据与阈值
3.1 相关系数凭什么能区分信号与噪声
Pearson 相关系数的定义是 ρ(a,b) = Σ(a-ā)(b-b̄) / √(Σ(a-ā)² · Σ(b-b̄)²)。拿它来筛 IMF，背后有三条性质撑腰。
第一，EMD 筛分保证各阶 IMF 在频带上近似分离，分量之间近似正交，所以单个 IMF 与原信号的相关性大致反映它在原信号中的能量贡献。第二，白噪声与任何确定性信号的相关期望为零，样本相关系数随长度 N 增大以 1/√N 量级收敛到零，因此噪声主导的高阶 IMF 天然相关偏低。第三，相关系数做了归一化，对不同幅值的分量可以放在同一把尺子上比较，这点比直接用能量更公平。
但要留个心眼：低频 IMF 的相关系数也会偏低。一个周期可能比整个样本长度还长，样本内根本看不到完整的相关结构，算出来的 ρ 自然小。这就是后面要讲的坑——一刀切会把有效的趋势分量误删。
3.2 MATLAB 里的手写实现
corr函数能用，但它内部会重复去均值、做类型检查，批量处理成千上万个样本时开销不小。手写内积版本更快，也更清楚每一步在干什么。
function rho = imf_corr(imfs, x) % imfs: K×N，每行一个 IMF；x: 1×N 原始信号；返回 K×1 相关系数 K = size(imfs, 1); rho = zeros(K, 1); xc = x(:) - mean(x); nx = norm(xc) + eps; % 加 eps 防止全零信号除零 for k = 1:K ic = imfs(k, :)' - mean(imfs(k, :)); rho(k) = abs(ic' * xc) / (norm(ic) * nx); end end分母是去均值后的二范数乘积，等价于标准差乘积，所以结果就是标准的 Pearson 系数。取绝对值是因为我们只关心相关强度，不关心正相关还是负相关——一个分量与原始信号反相，同样说明它承载了结构信息。eps保护针对的是端点飞边把某个 IMF 压成近似常数的极端情况，不加的话会直接返回 NaN，后面整条特征链全废。
3.3 阈值怎么定：固定值、局部极小还是累积贡献
固定阈值最省事。ρ 小于 0.1 或 0.2 判定为噪声主导，直接剔除。优点是批量处理同源信号时行为一致，缺点是面对不同信噪比的信号毫无适应性——高噪信号可能把所有分量都筛掉，低噪信号可能一个都筛不掉。血泪经验是：固定阈值只适合在数据源单一、已经验证过的流水线里用。
第一个局部极小点更稳一些。相关系数序列通常呈现「前几个大、中间下降、尾部略回升」的形状，从 k = 2 开始找第一个满足 ρ(k) < ρ(k-1) 且 ρ(k) < ρ(k+1) 的位置，之前的保留，之后的剔除。这个策略在机械振动和脑电信号上表现都不错，因为它自动适应了信号的衰减速度。缺陷是当 ρ 序列单调递减时不存在局部极小，得退回到固定阈值兜底。
累积贡献率从能量角度切入：把 ρ 从大到小排序后累加，达到 0.9 到 0.95 就截止。物理意义清晰，但阈值仍然要人工定，而且排序会打乱 IMF 的阶数顺序，后续拼特征矩阵时得记住原始索引。
function keepIdx = select_imf(rho, mode, thr) K = numel(rho); switch mode case 'fixed' keepIdx = find(rho >= thr)'; % thr 典型 0.1~0.2 case 'firstmin' keepIdx = 1:K; for k = 2:K-1 if rho(k) < rho(k-1) && rho(k) < rho(k+1) keepIdx = 1:k-1; % 局部极小点之前全保留 break; end end if numel(keepIdx) == K && all(rho < thr) % 单调递减时兜底 keepIdx = find(rho >= thr)'; end case 'cumenergy' [~, ord] = sort(rho, 'descend'); c = cumsum(rho(ord)) / sum(rho); n = find(c >= thr, 1); % thr 典型 0.9~0.95 keepIdx = sort(ord(1:n))'; % 恢复原始阶数顺序 end end筛选之后一定要做重构验证：把保留的 IMF 加上残余分量求和，与原始信号算 RMSE。如果 RMSE 相比筛选前突然跳升一个量级，说明筛过头了，把有效分量也扔了。这一步几乎不花时间，却能挡住大部分阈值设错的情况。
4. 信息熵特征：能量熵、样本熵与特征矩阵拼装
4.1 能量熵：分量能量分布的香农熵
能量熵刻画的是「能量在各个 IMF 上摊得有多开」。先算每个 IMF 的能量 E_k = Σ_t imf_k(t)²，归一化成占比 p_k = E_k / ΣE，再套香农熵 H = -Σ p_k ln p_k。
function H = energy_entropy(imfs) % imfs: K×N，每行一个 IMF；返回标量熵值（单位 nat） E = sum(imfs.^2, 2); % 各分量能量，K×1 p = E / sum(E); % 能量占比 p(p <= 0) = eps; % 避免 log(0) H = -sum(p .* log(p)); end熵值范围是 [0, ln K]，K 个分量能量完全均分时取最大值。熵大说明能量分散、信号成分复杂；熵小说明能量集中在少数几个分量上，结构简单。用log2则单位变成 bit，量级上更适合和别的特征拼在一起。
这里有个口径问题容易被忽略：到底对全部 IMF 算，还是只对筛选后的 IMF 算。两种做法都有人用，但含义完全不同——对全部算，噪声分量的能量占比会稀释熵值，等于把噪声水平也编码进去了；对筛后集合算，熵值只反映有效分量的能量结构。选哪种都行，但整条流水线必须统一，否则训练集和测试集的特征含义对不上，模型精度会莫名其妙地掉。
除了熵，单个分量的能量占比 p_k 本身也是好特征，粒度比熵更细。我一般把「熵 + 前几个分量的占比」一起拼进去，让分类器自己决定用哪个。
4.2 样本熵与排列熵：衡量分量的复杂度
样本熵衡量时间序列的规律性。给定嵌入维 m 和容限 r，统计长度为 m 的向量对在 Chebyshev 距离下小于 r 的匹配概率 B^m，再统计 m+1 维的匹配概率 A^{m+1}，SampEn = -ln(A/B)。值越小说明序列越规律、越可预测，越大说明越接近随机。
参数上，m 通常取 2，r 取 0.1 到 0.25 倍标准差。r 的敏感度很高：同一个 IMF，r 从 0.1σ 调到 0.2σ，样本熵能差出一倍。工程上的做法是先在一小批样本上扫一遍 r，挑类间区分度最大的那个值固定下来。
排列熵走的是另一条路。对 IMF 做相空间重构（嵌入维 m 取 3 到 7，延迟 τ 取 1），每个窗口内的 m 个值按大小排序得到一个序数模式，总共 m! 种。统计每种模式的出现概率 p_j，熵值 H = -Σ p_j ln p_j，再除以 ln(m!) 做归一化。
function pe = perm_entropy(x, m, tau) % m: 嵌入维；tau: 延迟；返回归一化排列熵，范围 [0,1] N = numel(x); n = N - (m - 1) * tau; if n <= 1, pe = 0; return; end patterns = zeros(1, n); for i = 1:n seg = x(i : tau : i + (m-1)*tau); [~, ord] = sort(seg); % 排序得到序数模式 patterns(i) = sum((ord - 1) .* factorial(m-1:-1:0)) + 1; % 编码成整数 end cnt = histcounts(patterns, 1:(factorial(m)+1)); p = cnt / sum(cnt); p(p <= 0) = eps; pe = -sum(p .* log(p)) / log(factorial(m)); % 归一化到 [0,1] end排列熵的优势是计算快、抗噪、对长序列友好，代价是丢掉了幅值信息，只保留序关系。如果两类信号的区别主要体现在幅值而非波形顺序上，排列熵就区分不开。
| 熵类型 | 计算复杂度 | 参数敏感度 | 适合场景 |
|---|---|---|---|
| 能量熵 | O(N) | 无参数 | 分量能量分布的整体描述 |
| 样本熵 | O(N²) | 高，r 极敏感 | 短序列的规律性度量 |
| 排列熵 | O(N·m log m) | 中，m 影响有限 | 长序列、实时场景 |
4.3 特征矩阵拼装与归一化
单个样本走完前面三步，得到的是 L 个筛后 IMF，每个 IMF 算一到两个熵，拼成 1×L 或 1×2L 的向量。所有样本按行堆起来，就是分类器的输入矩阵。
function feat = build_feature(imfs, keepIdx) % imfs: 全部 IMF；keepIdx: 筛选后保留的索引 L = numel(keepIdx); feat = zeros(1, L * 2); for j = 1:L seg = imfs(keepIdx(j), :); feat(2*j - 1) = energy_entropy(seg); feat(2*j) = perm_entropy(seg, 4, 1); end end归一化必须整批做，不能每个样本各算各的。正确顺序是：先把所有样本的特征拼成 M×D 矩阵，用训练集的均值和标准差做 z-score，再把同一组统计量套到测试集上。如果每个样本单独归一化，等于把样本间的相对差异抹平了，测试集精度会莫名其妙地崩。
注意：能量熵量级在 1 附近，排列熵归一化后在 0 到 1 之间，样本熵可能到 2 以上。不做归一化直接送进 SVM，核函数会被大量级特征主导，小量级特征等于没加。
5. 常见问题排查：这条链路上最容易翻车的五个点
5.1 端点飞边让 IMF 首尾剧烈发散
现象很直观：画出来的 IMF 在首尾两端振幅远大于中段，像被人拎起来抖了两下。相关系数会因此虚高——飞边那一段幅值大，内积贡献大，算出来的 ρ 不能反映真实的分量强度。
原因在于三次样条的外推行为。极值点只存在于数据内部，端点外侧没有插值节点，样条多项式会按自己的曲率向外延伸，形成发散。这个误差还会进入残余分量，下一层分解继续放大。
解决办法优先级从高到低：镜像延拓（把信号首尾各取一段翻转拼接，人为造出边界极值点）、极值点延拓（在端点外补几个虚拟极值点）、直接截断（丢掉首尾各 5% 的样本再算特征）。前两种改动小、信息损失少，第三种最粗暴但最省事，适合样本长度充裕的场景。
5.2 筛分次数过多，IMF 变成纯调频信号
现象是 IMF 波形看起来「太干净」：幅值几乎恒定，频率单一，包络平得像一条直线。这种分量算出来的能量熵异常小，排列熵接近 0，在特征空间里显得特别突出，但它对应的其实不是任何真实的物理模式。
成因是 SD 阈值设得太小，或者 S 次准则要求太严。筛分本质上是反复减均值，减到极致就只剩下纯调频成分，数学上完美满足 IMF 的两个条件，物理意义却丢光了。
解决上，SD 取 0.2 到 0.3 而不是 0.01，siftMax限制在 8 到 10 次，两个约束同时生效。判断是否过筛分有个土办法：看 IMF 的瞬时幅值方差，方差接近零就说明被榨干了。
5.3 噪声幅值和集成次数不匹配，CEEMD 退化成 EMD
现象是每次运行结果都不一样，IMF 阶数飘忽，重构误差 RMSE 在 5% 以上下不来。把noiseAmp从 0.1 调到 0.15 再跑，阶数又变了——看起来像是算法本身不稳定。
根因有两个。noiseAmp太小（低于 0.005 倍标准差）时，加进去的噪声改变不了极值点分布，EEMD 那套「抹匀」机制完全失效，等价于直接跑 EMD。NE太小则残余噪声消不掉，平均之后仍然有明显扰动。
解决：noiseAmp用 0.1 到 0.2 倍标准差，NE不低于 50，并且rng固定种子。验证方法很直接——把重构信号与原始信号做差，残差的标准差应该小于原始信号标准差的 1%。
5.4 相关系数阈值拍脑袋，把有效低频分量误删
现象是筛选后重构信号丢了整体趋势，只剩高频振荡，分类精度不升反降。查相关系数序列会发现，最低频那一两个 IMF 的 ρ 值确实很小，被阈值一刀切掉了。
原因是低频分量的周期可能比样本长度还长，样本内看不到完整的相关结构，Pearson 系数天然偏低。这跟「噪声主导」是两回事，但数值上表现一样。
解决：对最低频的一到两个分量单独处理，改用能量占比判断——如果它的能量占总量超过 5%，就保留，不管相关系数多小。或者干脆对低频段不做筛选，只筛中高频。这个改动很小，但对趋势敏感的信号（比如温度、位移）效果差别很大。
5.5 熵特征维度不统一，训练测试分布错位
现象是训练集上精度 98%，测试集只有 60%，交叉验证的折间方差也大得离谱。查特征矩阵会发现，不同样本的列数对不上，或者同一列的量级差异悬殊。
维度对不上通常是因为筛选阶段保留的 IMF 个数随样本变化。有的样本筛出 5 个，有的筛出 7 个，拼出来的特征向量长度不一致，只能截断或者补零，信息就乱了。
解决：筛选后固定保留 L 个分量，L 取所有样本中出现次数最多的那个值，不够的用零填充并在末尾加一个掩码位。归一化统一在拼完矩阵之后做，统计量只从训练集算，测试集直接套用。这两条听起来是常识，但赶进度的时候最容易跳过，翻车也最惨。
6. 设计取舍：熵特征选型与验证方法
这条链路上没有哪一步存在唯一正确答案，全是取舍。
分解层数上，多分一层能多拿一个尺度的信息，但也多一份端点误差和过筛分风险。我的习惯是maxIMF取floor(log2(N)) - 1，宁可少分一层，也不让最底层那个几乎全是数值噪声的分量混进特征。
熵的选择上，能量熵几乎免费（O(N)），排列熵中等，样本熵最贵（O(N²)）且对 r 极其敏感。如果样本量上万，样本熵基本可以排除；如果信号短、规律性差异是主要判别依据，样本熵反而最有效。我一般先用能量熵加排列熵打底，跑一轮交叉验证看混淆矩阵，哪两类分不开再针对性补样本熵。
相关系数筛选和「全量喂给分类器」也是个取舍。理论上让分类器自己学权重更优雅，实践中 IMF 维度一高，小样本训练必然过拟合，先做物理意义上的筛选反而更稳。但筛选阈值不要定太死，留一两个边缘分量进去，让模型有一点冗余空间。
验证方法上，我固定跑三件事。第一，重构误差：保留的 IMF 加残差与原始信号算 RMSE，超过原始信号标准差的 5% 就回头查分解参数。第二，可复现性：固定rng连跑三次，IMF 阶数和特征向量的差异应该在数值误差量级。第三，特征有效性：只用训练集做特征选择，测试集完全不参与，避免信息泄漏把精度虚抬上去。
从那以后我每次调这条链路，都强制先跑一遍重构误差和可复现性检查，确认分解本身没问题，再去动熵和分类器。跳过这步直接调模型，十有八九是在给一个错误的特征矩阵调参。希望帮到你。
本文还有配套的精品资源，点击获取
