★ FEATURED ARTICLE 简介这份资源面向信号处理与机器学习方向的学习者与研究人员聚焦非线性、非平稳振动信号的分解与特征提取。包内共4个文件以3个MATLAB脚本.m和1个Excel数据表.xlsx为主压缩包约60KB脚本可直接运行数据表提供实测样本便于复现与二次开发。资源围绕CEEMD完全集合经验模态分解展开先由分解脚本得到一系列IMF分量再通过相关系数筛选出最具代表性的分量去除冗余与无关成分随后计算各分量的信息熵量化其分布均匀性与随机性最终组合成特征向量矩阵为后续分类或模式识别任务提供输入。已有5245人学习下载说明该思路在工程实践中具有一定参考价值。对于希望掌握从信号分解、特征筛选到特征构建完整流程的读者这份代码与数据可作为可直接上手的实践起点。1. 从一段非平稳振动信号说起这套 CEEMD-相关系数-信息熵特征到底能干什么手里有一段风机轴承的振动信号采样率 12.8 kHz采样 10 秒画出来是一团毛刺直接扔进 SVM 或神经网络准确率大概率在 60% 上下晃。问题不在分类器在于原始时序信号维度太高、噪声太重特征没提干净。这套「CEEMD 分解 → IMF 分量相关系数筛选 → 信息熵特征提取」的组合就是解决这个问题的经典链路先用互补集合经验模态分解把非平稳信号拆成若干个本征模态函数再用相关系数把混进来的伪分量剔掉最后对保留下来的分量算信息熵拼成一个低维、有物理意义的特征向量。它适合做旋转机械故障诊断、生理信号分类、地震波分析这类「信号非平稳、样本量不大、又不想上深度学习」的场景。MATLAB 是这条链路最顺手的工具因为 EMD 工具箱和熵计算都能几行搞定不用自己从零写筛分过程。2. CEEMD 分解为什么不用 EMDMATLAB 里怎么落地2.1 EMD 的模态混叠与 CEEMD 的补集思路经验模态分解EMD的核心是「筛分」找信号所有局部极大值和极小值分别拟合上、下包络取均值得到一条低频趋势原信号减去这条趋势就得到一个 IMF重复直到剩余分量单调。听起来干净但实际信号里一旦出现间歇性高频成分同一个 IMF 里会混进尺度差异很大的振荡这就是模态混叠。混叠一旦发生后面算出来的熵值就不可信了。EEMD 的思路是往信号里加白噪声利用噪声的均匀分布特性把不同尺度的成分「推开」多次分解取平均。但 EEMD 加的是独立噪声重构时残留噪声会污染分量。CEEMD 改成正负成对加噪每次加一对幅值相同、符号相反的噪声分别分解再对应 IMF 求平均。这样噪声在平均过程中相互抵消重构误差接近零同时保留了 EEMD 抗混叠的能力。代价是计算量翻倍但对 10 秒量级的信号MATLAB 里跑一次也就几秒到几十秒完全可接受。2.2 MATLAB 实现从信号读入到 IMF 矩阵MATLAB 官方从 R2018b 起在 Signal Processing Toolbox 里提供了emd函数但 CEEMD 需要自己封装正负成对加噪逻辑。下面这段是我常用的实现输入是一维列向量输出是 IMF 矩阵每列一个分量和残差。function [imfs, residual] ceemd_custom(x, num_imfs, num_ensembles, noise_std) % x 输入信号列向量 % num_imfs 期望分解出的 IMF 个数 % num_ensembles 正负成对加噪的组数常用 50~100 % noise_std 噪声标准差一般取 0.1~0.2 倍信号标准差 x x(:); N length(x); imfs_sum zeros(N, num_imfs); for k 1:num_ensembles % 生成一对正负噪声 noise noise_std * std(x) * randn(N, 1); x_pos x noise; x_neg x - noise; % 分别做 EMD取前 num_imfs 个分量 imf_pos emd(x_pos, MaxNumIMF, num_imfs); imf_neg emd(x_neg, MaxNumIMF, num_imfs); % 对齐长度后累加 imfs_sum(:, 1:size(imf_pos,2)) imfs_sum(:, 1:size(imf_pos,2)) imf_pos; imfs_sum(:, 1:size(imf_neg,2)) imfs_sum(:, 1:size(imf_neg,2)) imf_neg; end imfs imfs_sum / (2 * num_ensembles); residual x - sum(imfs, 2); end逻辑说明循环里每次生成一对符号相反的噪声分别叠加到原信号上做 EMD把两组 IMF 累加。循环结束后除以总分解次数2 倍组数得到平均 IMF。residual是原信号减去所有 IMF 之和正常应该接近零如果残差幅值超过信号幅值的 5%说明num_imfs设小了或者噪声幅值不合适。参数说明num_imfs一般取fix(log2(N)) - 1左右N 是信号长度比如 N12800 时取 12 到 13num_ensembles取 50 起步要求高就上 100再高收益递减noise_std是关键参数取 0.1 到 0.2 倍信号标准差太小抗混叠效果弱太大引入伪分量。我一般先跑一次看 IMF 的频谱分布再回头调这三个参数。2.3 分解结果怎么判断靠不靠谱跑完 CEEMD 别急着往下走先看三件事。第一各 IMF 的过零点数是否随阶数递减正常应该从高频到低频单调变化如果中间某个分量过零点数突然反弹说明混叠没压住。第二残差幅值是否足够小前面说了 5% 是红线。第三把前三个 IMF 画出来看有没有明显的端点飞翼——EMD 类方法在信号两端容易发散如果飞翼严重要么加窗要么用镜像延拓预处理。这三步花不了两分钟但能省掉后面反复调特征的功夫。3. 相关系数筛选把伪分量剔掉再算熵3.1 为什么不能对所有 IMF 直接算熵CEEMD 分解出来的分量里前几个高频 IMF 往往是噪声主导最后几个低频 IMF 可能是趋势项或分解伪影。如果把这些分量一股脑拿去算信息熵噪声分量的熵值接近随机分布的理论上限趋势分量的熵值又异常低两类极端值混进特征向量分类器会被带偏。常见做法是算每个 IMF 与原信号的皮尔逊相关系数保留相关系数大于某个阈值的分量。这个阈值的选取有讲究固定阈值比如 0.1在部分信号上会漏掉有用分量我一般用「相关系数均值 标准差」作为自适应阈值或者画相关系数曲线找拐点。3.2 相关系数计算与分量筛选代码% 假设 imfs 是上一步得到的 IMF 矩阵x 是原信号 num_imfs size(imfs, 2); corr_coef zeros(num_imfs, 1); for i 1:num_imfs % 皮尔逊相关系数 r corrcoef(imfs(:, i), x); corr_coef(i) abs(r(1, 2)); end % 自适应阈值均值 标准差 threshold mean(corr_coef) std(corr_coef); selected_idx find(corr_coef threshold); % 如果筛选后分量太少退化为保留前 60% if length(selected_idx) 3 [~, sort_idx] sort(corr_coef, descend); selected_idx sort_idx(1:max(3, round(0.6 * num_imfs))); end imfs_selected imfs(:, selected_idx);逻辑说明先逐个算 IMF 与原信号的相关系数绝对值再用均值加标准差做阈值。这个自适应规则的好处是不依赖人工设定固定值对不同类型的信号都能给出一个相对合理的筛选结果。兜底逻辑是防止阈值过高导致只剩一两个分量那样特征维度太低分类器学不到东西。参数说明corrcoef返回 2×2 矩阵取(1,2)位置就是两个向量的相关系数。阈值规则可以换成「均值 0.5 倍标准差」让筛选更宽松或者「均值 1.5 倍标准差」更严格。如果信号本身信噪比很低建议用后者宁可少留几个分量也别把噪声放进来。3.3 筛选后的分量对齐与重构验证筛选完记得做一次重构验证把选中的 IMF 加起来和原信号对比看主要波形结构有没有丢失。如果重构信号和原信号的相关系数低于 0.9说明筛选太狠把承载主要能量的分量也扔了。这时候要么放宽阈值要么检查 CEEMD 的噪声参数是不是设得太大导致能量分散。这一步很多人跳过结果后面分类效果差回头查半天才发现是筛选阶段就把有用信息丢了。4. 信息熵特征提取排列熵、样本熵、近似熵怎么选怎么算4.1 三种熵的适用边界信息熵在故障诊断里常用的有排列熵、样本熵、近似熵。排列熵计算最快对突变敏感适合在线监测场景但它只看序数关系对幅值变化不敏感。样本熵抗噪能力好适合信噪比低的信号缺点是参数敏感嵌入维数和相似容限设不好结果波动大。近似熵和样本熵思路接近但自匹配问题导致偏差现在用得少了。我一般对筛选后的每个 IMF 同时算排列熵和样本熵拼成二维特征比单用一种熵更稳。4.2 排列熵的 MATLAB 实现function pe permutation_entropy(x, m, tau) % x 输入分量 % m 嵌入维数常用 3~7 % tau 延迟时间常用 1 x x(:); N length(x); % 重构相空间 num_vectors N - (m - 1) * tau; patterns zeros(num_vectors, m); for i 1:num_vectors for j 1:m patterns(i, j) x(i (j - 1) * tau); end end % 每个向量转成序数模式 [~, order] sort(patterns, 2); % 把序数模式编码成唯一整数 pattern_idx zeros(num_vectors, 1); for i 1:num_vectors code 0; for j 1:m code code * m (order(i, j) - 1); end pattern_idx(i) code; end % 统计各模式出现概率 [counts, ~] histcounts(pattern_idx, BinMethod, integers); p counts / sum(counts); p p(p 0); % 归一化排列熵 pe -sum(p .* log(p)) / log(factorial(m)); end逻辑说明先把一维序列重构成num_vectors × m的相空间矩阵每行是一个嵌入向量。对每行做排序得到序数排列模式再把排列编码成整数方便统计。最后按香农熵公式计算除以log(factorial(m))做归一化让熵值落在 0 到 1 之间。参数说明m取 3 到 7太小区分度不够太大计算量涨且对噪声敏感我一般从 4 开始试。tau取 1 对大多数振动信号够用如果信号采样率极高、相邻点高度相关可以取 2 或 3 做降相关。归一化是可选的但做多分量特征拼接时建议归一化避免不同分量的熵值量纲不一致。4.3 样本熵的 MATLAB 实现与参数陷阱function se sample_entropy(x, m, r) % x 输入分量 % m 嵌入维数常用 2 % r 相似容限常用 0.1~0.25 倍标准差 x x(:); N length(x); r r * std(x); % 计算 m 维和 m1 维的匹配数 B count_matches(x, m, r, N); A count_matches(x, m 1, r, N); if B 0 se NaN; return; end se -log(A / B); end function count count_matches(x, m, r, N) count 0; for i 1:N - m for j i 1:N - m if max(abs(x(i:im-1) - x(j:jm-1))) r count count 1; end end end end逻辑说明样本熵的核心是统计在容限r下m 维向量匹配的对数 B 和 m1 维向量匹配的对数 A取负对数比。双重循环是 O(N²) 复杂度N 上万时在 MATLAB 里会明显变慢实际用的时候建议对分量做降采样或者用pdist配合阈值统计加速。参数说明m通常取 2r取 0.1 到 0.25 倍标准差。这里有个血泪经验r设太小匹配数接近零熵值不稳定甚至返回 NaN设太大所有向量都匹配熵值趋近零特征失去区分度。我一般先算分量标准差取 0.15 倍作为起点跑完看熵值分布如果方差太大就微调。另外样本熵对数据长度敏感N 低于 500 时结果不可靠信号短的话要么换排列熵要么做数据增强。5. 避坑与排查这套链路最容易翻车的五个地方5.1 现象CEEMD 跑完 IMF 个数每次不一样原因emd函数默认的筛分停止准则和最大迭代次数没有固定加上每次加的噪声是随机的导致分解层数波动。解决在调用emd时显式指定MaxNumIMF并且在 CEEMD 主循环开始前用rng(42)固定随机种子。这样同一段信号每次跑出来的 IMF 矩阵完全一致方便复现和调参。5.2 现象相关系数筛选后只剩一个分量原因信号信噪比极低或者 CEEMD 的noise_std设得太大导致能量被分散到很多伪分量上每个分量和原信号的相关系数都不高。解决先把noise_std降到 0.05 倍标准差重跑如果还是不行改用固定阈值 0.05 而不是自适应阈值保证至少留下几个分量。另外检查信号有没有做去均值处理带直流偏置的信号会让相关系数计算失真。5.3 现象样本熵返回 NaN原因r设得太小m 维匹配数 B 为零代码里直接返回 NaN。解决在count_matches里加一个判断如果 B 小于 10自动把r放大 1.5 倍重新计算最多迭代三次。或者干脆换排列熵排列熵不存在匹配数为零的问题鲁棒性更好。5.4 现象特征向量拼好后分类准确率还不如原始信号原因熵特征没有做归一化不同分量的熵值量纲差异大分类器被大数值特征主导。解决对每个特征维度做 z-score 标准化即减去均值除以标准差。如果用的是 SVM标准化后还要检查核函数参数RBF 核的gamma对特征尺度很敏感标准化前后需要重新调参。5.5 现象整个流程跑一次要十几分钟原因CEEMD 的num_ensembles设得太大加上样本熵的双重循环计算量叠加。解决num_ensembles从 100 降到 30 到 50对大多数信号精度损失在 2% 以内。样本熵部分把分量降采样到 2000 点以内再算或者用parfor把不同分量的熵计算并行化。MATLAB 的并行工具箱在这个场景下提速很明显前提是循环之间没有数据依赖。6. 进阶技巧把熵特征和时域统计量拼成混合特征向量单靠熵特征分类器能拿到的信息有限。我一般会在熵特征基础上再拼几个时域统计量均方根、峭度、峰值因子。这三个量计算快对冲击性故障敏感和熵特征形成互补。具体做法是对筛选后的每个 IMF 算排列熵、样本熵、均方根、峭度、峰值因子五个值拼成一个行向量多个分量就拼成矩阵最后展平成一维特征向量。feature_vec []; for i 1:size(imfs_selected, 2) comp imfs_selected(:, i); pe permutation_entropy(comp, 4, 1); se sample_entropy(comp, 2, 0.15); rms_val rms(comp); kurt_val kurtosis(comp); crest_val max(abs(comp)) / rms_val; feature_vec [feature_vec, pe, se, rms_val, kurt_val, crest_val]; end % z-score 标准化 feature_vec (feature_vec - mean(feature_vec)) / std(feature_vec);这段代码里rms、kurtosis都是 MATLAB 内置函数直接调用即可。crest_val是峰值因子反映信号中的冲击成分。标准化放在最后一步对所有特征统一做避免量纲问题。验证这套特征有没有用我习惯用两个指标一是用fitcsvm跑五折交叉验证看准确率是否比原始信号直接分类高 15 个百分点以上二是画特征向量的 t-SNE 散点图看不同类别的样本有没有聚成簇。如果散点图混在一起说明特征区分度不够回头检查 IMF 筛选和熵参数。从那以后我每次搭这套链路都会先跑一遍 t-SNE 看一眼特征分布再往下走省得后面调分类器调半天发现是特征的问题。希望帮到你。本文还有配套的精品资源点击获取
阅读完成 · 觉得有帮助？
