CEEMD分解、IMF相关系数筛选与信息熵特征提取实战

简介这份资料包围绕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 到 CEEMDIMF 的筛分条件与互补加噪的重构恒等式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 到 10sdTh就是下面要说的 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 榨成纯调频信号sdThSD 停止阈值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每行一个 IMFx: 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) ρ(k1) 的位置之前的保留之后的剔除。这个策略在机械振动和脑电信号上表现都不错因为它自动适应了信号的衰减速度。缺陷是当 ρ 序列单调递减时不存在局部极小得退回到固定阈值兜底。累积贡献率从能量角度切入把 ρ 从大到小排序后累加达到 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(k1) 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再统计 m1 维的匹配概率 A^{m1}SampEn -ln(A/B)。值越小说明序列越规律、越可预测越大说明越接近随机。参数上m 通常取 2r 取 0.1 到 0.25 倍标准差。r 的敏感度很高同一个 IMFr 从 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: 全部 IMFkeepIdx: 筛选后保留的索引 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.01siftMax限制在 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 阶数和特征向量的差异应该在数值误差量级。第三特征有效性只用训练集做特征选择测试集完全不参与避免信息泄漏把精度虚抬上去。从那以后我每次调这条链路都强制先跑一遍重构误差和可复现性检查确认分解本身没问题再去动熵和分类器。跳过这步直接调模型十有八九是在给一个错误的特征矩阵调参。希望帮到你。本文还有配套的精品资源点击获取