非平稳信号特征提取：CEEMD、IMF相关系数与信息熵实战

简介面向信号处理与机器学习特征提取学习者的一套MATLAB实践材料围绕CEEMD分解、IMF分量相关系数筛选与信息熵特征构建展开可用于非线性非平稳振动信号的分解、冗余分量剔除及分类特征矩阵生成等场景。压缩包共4个文件含3个m脚本与1个xlsx数据表整体约60KB脚本对应分解与特征计算流程表格可作为实验样本数据便于直接运行和对照理解。已有5248人学习下载说明该流程在相关课程与课题中具有一定参考价值。读者可据此掌握从信号分解、分量相关性分析到信息熵量化与特征向量矩阵组织的完整思路并借助MATLAB实现和可视化结果加深对IMF时间尺度特性、分布均匀性及随机性的认识适合作为课程实验、论文复现或特征工程入门的轻量素材。1. 一段非平稳信号摆在面前CEEMD、IMF相关系数、信息熵这三样东西怎么串成一条特征流水线手上有一批振动信号、脑电片段或者轴承加速度数据想送进分类器第一步就卡在特征怎么构造上。直接上时域统计量均值、峭度、RMS信息量太薄上 FFT 又默认信号是平稳的碰上冲击、突变、工况切换就散架。这时候常见的解法是走「自适应分解 分量筛选 熵特征」这条路先用 CEEMD互补集合经验模态分解把信号拆成若干个本征模态函数IMF再用 IMF 与原始信号的相关系数把伪分量剔掉最后对保留下来的分量逐个算信息熵拼成一张维度固定、物理意义清楚的特征表。这套组合在 MATLAB 里能完整落地不需要深度学习工具箱也不依赖大样本。下面按「分解怎么调参 → 分量怎么筛 → 熵怎么算 → 哪里会翻车 → 怎么验证」的顺序把这条流水线拆开讲。2. CEEMD分解的MATLAB落地噪声幅值Nstd与集成次数NE怎么定才不白算2.1 EMD的模态混叠与CEEMD用正负噪声对补了什么EMD 的核心是筛分把信号上下包络的均值一层层剥掉剩下的震荡部分就是 IMF。它的问题出在「筛分」这个动作完全由数据自己决定没有基函数约束一旦信号里出现间断、脉冲或者尺度相近的两个成分同一个 IMF 里就会混进两个物理上不相干的振荡这就是模态混叠。混叠一旦发生后面算的相关系数和熵全都建立在错误的分量上特征表看着有数实际是噪声。EEMD 的思路是往信号里加白噪声再分解利用白噪声在时频平面上的均匀分布把极值点「铺平」让不同尺度落到各自的 IMF 上最后对多次集成结果取平均。代价是白噪声没法完全抵消重构残差里留了一层噪声。CEEMD 的改动很小但很关键每次集成加一对正负相反的白噪声分解完再一起平均。因为正负噪声的期望为零残余噪声在成对平均时相互抵消重构误差比 EEMD 小一个量级。理解这一点很重要——CEEMD 不是新算法它是 EEMD 的配对采样版本参数含义几乎完全继承。在 MATLAB 里落地要注意版本问题。R2018b 之后emd函数进了 Signal Processing Toolbox可以直接调更早的版本或者没有该工具箱的环境得用第三方实现社区里流传较广的是基于筛分循环自己写的版本函数签名会不一样。CEEMD 本身没有官方函数得自己在emd外面套一层集成循环。下面这份封装是我常用的写法逻辑就是把配对加噪、分解、补零对齐、平均这四步写死。function [IMF, RES] ceemd2(x, Nstd, NE, K, seed) % CEEMD正负成对高斯白噪声辅助的EMD集成 % x 输入一维信号 % Nstd 噪声标准差相对系数常用 0.1~0.2 % NE 集成次数对数总分解次数为 2*NE % K 统一输出的IMF层数用来对齐不同次分解的维度 % seed 随机种子固定后结果可复现 if nargin 5, seed 20240101; end rng(seed, twister); % 关键不固定种子两次跑出来不一样 x x(:); N numel(x); acc zeros(N, K); cnt 0; for k 1:NE w Nstd * std(x) * randn(N, 1); % 噪声幅值按信号标准差缩放 imp emd(x w, MaxNumIMF, K, MaxIteration, 100); imn emd(x - w, MaxNumIMF, K, MaxIteration, 100); acc acc padcols(imp, K) padcols(imn, K); cnt cnt 2; end IMF acc / cnt; % 成对平均残余噪声互相抵消 RES x - sum(IMF, 2); % 残余项单独留下 end function Y padcols(X, K) % 不同次分解出的IMF个数可能不一致统一补零到K列 [n, m] size(X); Y zeros(n, K); Y(:, 1:min(m, K)) X(:, 1:min(m, K)); end这段代码里padcols是必须的。emd每次返回的 IMF 列数受信号局部极值点数影响同一段信号加不同的噪声层数可能差一两层。如果不做对齐矩阵拼接会直接报维度错误或者更隐蔽地——把某次分解的第二层错位接到另一次的第一层上平均出来的 IMF 物理意义全乱。补零对齐虽然会引入一点点幅值偏差零填充本身不参与后续重构只影响循环内的累加但它保证了每一列对应的是同一尺度的分量这是后面算相关系数的前提。2.2 参数取值Nstd、NE、K分别卡在什么区间Nstd控制噪声幅值相对信号标准差的比例。取太小比如 0.01噪声不足以铺平极值点模态混叠照旧取太大0.5 以上噪声本身会引入大量虚假分量而且残余噪声抵消不干净。工程上 0.10.2 是覆盖大多数振动和生理信号的区间信号本身信噪比高、冲击成分明显时偏 0.1信号平稳、尺度接近时偏 0.2。判别办法不是拍脑袋同一段信号取 0.1、0.15、0.2 各跑一次看第一层 IMF 的主频是否稳定如果三次结果的主频漂移超过 10%说明噪声幅值选得不合适。NE是集成对数总分解次数是它的两倍。这个参数直接决定耗时NE50意味着 100 次 EMD一段 2000 点的信号大概几秒到十几秒看机器。从误差收敛的角度NE从 50 往上加重构误差下降很快超过 200 之后收益就很平了。我的习惯是先按NE50把整条流水线跑通、确认特征表维度对得上再上NE100或200出最终结果。中途调参时用 50别一上来就 200否则改一次参数等两分钟调试节奏会被拖死。K是统一输出的层数取值要跟信号长度挂钩。经验公式是K fix(log2(N)) - 1N1024时约 9 层N2048时约 10 层。取多了后面几层基本是数值噪声还会让样本熵算出来全是 0取少了会切掉低频趋势重构残差偏大。一个实用的检查是看RES的幅值如果残余项的标准差超过原信号标准差的 5%说明层数给少了或者信号里有明显的趋势项没剥干净。2.3 分解结果的自检三件事必须在算相关系数之前做完分解跑完不要急着往下走先做三项自检。第一sum(IMF,2) RES与原始信号的误差应该在 1e-10 量级超了说明平均环节出了错。第二第一层 IMF 的主频应该落在信号最高频段如果它比第二层还低多半是筛分停止条件设得太松MaxIteration太小。第三画一张 IMF 的时域叠图肉眼确认没有某一层突然出现跟其他层完全不搭的冲击——那基本就是模态混叠留下的痕迹得回头调Nstd。这三步花不了两分钟但能挡掉后面 80% 的返工。我自己就吃过亏有一次直接跳过自检算完一整套熵特征去做分类准确率死活上不去回头一画图才发现第五层 IMF 里混进了工频干扰相关系数还偏偏不低硬筛根本没筛掉。3. IMF分量筛选相关系数阈值定在0.1还是0.3扔掉的分量到底是不是噪声3.1 相关系数在IMF上度量的是什么IMF 与原始信号的相关系数本质上是度量这个分量里含有多少原始信号的线性成分。CEEMD 分解出来的前几层通常是高频噪声和冲击成分中间几层承载主要能量最后几层是低频趋势。相关系数的典型形态是一条先升后降的曲线峰值出现在能量最集中的那一两层两头低。这里有个容易踩的认知误区相关系数低不等于该分量没用。有些故障特征恰恰藏在能量很小的高频 IMF 里它的相关系数可能只有 0.05但对分类的判别力很强。反过来相关系数高的分量也可能只是把工频或者趋势项原样搬了一遍。所以相关系数筛选的作用是「剔除明显的伪分量」不是「挑出最有用的分量」这两件事别混。计算上要用绝对值。IMF 与原始信号可能出现负相关相位相反带符号的相关系数会把这种情况算成负值如果阈值是正数这些分量会被误删。手写实现比调corr稳不依赖 Statistics Toolboxfunction r imfCorr(IMF, x) % 逐列计算IMF与原始信号的皮尔逊相关系数绝对值 x x(:); xc x - mean(x); K size(IMF, 2); r zeros(K, 1); for i 1:K yi IMF(:, i) - mean(IMF(:, i)); r(i) abs(sum(xc .* yi) / (norm(xc) * norm(yi) eps)); end endeps是防零除的norm直接算二范数比std更直观。这段代码在N2000、K10时耗时可以忽略不用向量化优化。3.2 三种阈值策略的实算对比固定阈值、相对阈值、累积贡献率是三种最常见的筛法效果差别不小。固定阈值就是r 0.1或0.2写起来简单但它对信号幅值和长度的变化不敏感。同一台设备换个工况相关系数整体抬高或压低固定阈值就失灵了。相对阈值取r alpha * max(r)alpha常见 0.1它跟着信号自身走鲁棒性好得多。累积贡献率是先按相关系数降序排取累积和达到总和的某个比例比如 0.85为止适合分量数少、需要保留更多信息的情况。function idx selectIMF(r, mode, param) % mode: fixed 固定阈值 | relative 相对最大值的比例 | cumrate 累积贡献率 switch mode case fixed idx find(r param); case relative idx find(r param * max(r)); case cumrate [rs, ord] sort(r, descend); c cumsum(rs) / sum(rs); idx sort(ord(1:find(c param, 1))); end idx idx(:); end用同一段轴承冲击信号实测fixed(0.2)留下 4 层relative(0.1)留下 6 层cumrate(0.85)留下 5 层。三种结果里相对阈值那组在后续分类上的表现最稳因为它自动适应了这段信号相关系数整体偏低的特性。我一般默认用relative(0.1)只有在样本量特别小、需要尽量压缩特征维度时才切到cumrate。筛完之后有两种用法硬筛选是只保留选中的分量其余直接丢掉软加权是给每个分量乘上它的相关系数再求和构造一个加权复合分量。硬筛选维度低、解释性好软加权保留了更多信息但维度不降反升。分类任务用硬筛选趋势监测或者异常检测用软加权更合适。3.3 相关系数之外的两个补充判据相关系数有它的盲区尤其是当伪分量和有效分量能量接近的时候。补两个判据能让筛选更稳。第一个是过零率。噪声型 IMF 的过零率接近 0.5每点都变号而有效分量因为有一定周期性过零率明显偏低。算出来之后过零率超过 0.45 的分量可以直接标记为高频噪声不管它的相关系数是多少。第二个是能量占比。E_i / sum(E)小于 1% 的分量即使相关系数过线它对最终特征的贡献也可以忽略保留下来只是徒增维度。这两个判据跟相关系数取交集能挡掉绝大多数「相关系数够高但其实没用」的分量。function [keep, zcr, ep] refineIMF(IMF, r, alpha) % 相关系数 过零率 能量占比三重判据 [~, K] size(IMF); zcr zeros(1, K); ep zeros(1, K); for i 1:K s IMF(:, i); zcr(i) sum(abs(diff(sign(s))) 0) / numel(s); end E sum(IMF.^2, 1); ep E / sum(E); keep find(r alpha * max(r) zcr 0.45 ep 0.01); keep keep(:); end注意diff(sign(s))遇到零值会漏计一次过零对结果影响很小不用专门处理。真正要注意的是ep 0.01这个门槛信号长度越长、分量越多单层能量占比越低0.01在小样本上可能筛得过狠样本少于 500 点时降到0.005。4. 信息熵特征提取样本熵、排列熵、能量熵在IMF上的参数与实算4.1 三种熵各自刻画什么别混着用熵在特征工程里的作用是描述复杂度但不同熵的侧重点差得很远混着用会让特征表失去物理意义。样本熵Sample Entropy度量的是时间序列中相似模式重复出现的概率值越小说明信号越规律、越可预测值越大说明越随机。它对参数r很敏感r取固定常数在不同幅值的 IMF 上完全不可比。排列熵Permutation Entropy只看序关系不看具体数值对幅值变化免疫计算快适合捕捉突变但它对噪声也免疫——噪声和真实突变在排列熵上长得一样。能量熵Energy Entropy最简单就是各 IMF 能量占比的香农熵反映能量在不同尺度上的分散程度值大说明能量分布均匀值小说明能量集中在少数几个分量上。我的组合方式是对每个保留的 IMF 分别算样本熵和排列熵再对整个分量集合算一个能量熵。这样特征维度是2 * K_keep 1既有逐尺度的细节也有全局的分布描述。4.2 样本熵手写实现绕开O(N²)内存的那点事样本熵的定义不复杂把序列构造成长度为m的模板向量统计距离小于r的模板对数B再对m1长度重复一次得到A最后S -ln(A/B)。问题在于标准实现要构造N×N的距离矩阵N5000时就是 2000 万个元素的矩阵内存直接爆掉。所以必须写成双重循环按需计算虽然慢一点但内存恒定。function S sampen(x, m, r) % 样本熵双循环实现避免构造 N×N 距离矩阵 x x(:); N length(x); if N m 2, S NaN; return; end B 0; A 0; for i 1:N - m xi x(i : i m - 1); for j i 1 : N - m % 切比雪夫距离模板各点差值的最大绝对值 if max(abs(xi - x(j : j m - 1))) r B B 1; end end end for i 1:N - m xi x(i : i m); for j i 1 : N - m if max(abs(xi - x(j : j m))) r A A 1; end end end S -log((A eps) / (B eps)); end参数m一般取 2这是文献里最常用的值取 3 会让模板匹配率骤降、统计不稳定。r取0.1到0.25倍该分量自身的标准差注意是「自身」每个 IMF 单独算不能全用原始信号的 std。这一条是血泪经验早期我图省事用统一的r结果高频分量的熵值系统性偏高低频分量偏低特征表里这两类分量的熵根本不在一个量纲上标准化都救不回来。A和B加eps是为了防止某段信号特别规律、A0时出现-Inf。真遇到A0说明这段信号在该参数下完全可预测熵值应该记 0 而不是无穷大可以在调用前先判断if A 0, S 0; return; end。4.3 排列熵与能量熵嵌入维数、延迟与归一化排列熵的两个参数是嵌入维数m和延迟tau。m太小2区分度不够太大8 以上时m!种排列在短序列上根本采不全反而失真。工程上m取 4 到 6 最稳tau一般取 1。归一化要除以log(factorial(m))这样熵值落在[0,1]区间不同m之间的结果才可比。function H permen(x, m, tau) % 排列熵归一化到 [0,1] x x(:); N length(x); M N - (m - 1) * tau; if M 1, H NaN; return; end P zeros(M, m); for i 1:M P(i, :) x(i : tau : i (m - 1) * tau); end [~, ord] sort(P, 2); % 每行得到一个秩排列 code zeros(M, 1); for i 1:M c 0; for j 1:m c c * m (ord(i, j) - 1); end code(i) c; end p accumarray(code 1, 1, [m^m, 1]) / M; % 给定上界避免维度不一致 p p(p 0); H -sum(p .* log(p)) / log(factorial(m)); endaccumarray的第三个参数给定了输出维度[m^m, 1]这是必须的。不给的话某些排列没出现时输出长度会变跟M对不上。m6时m^m是 46656accumarray生成这么长的向量完全没压力m再加就不建议了。能量熵就三行但要记得归一化function H energyen(IMF) % 能量熵基于各IMF能量占比归一化到 [0,1] E sum(IMF.^2, 1); p E / sum(E); p p(p 0); H -sum(p .* log(p)) / log(numel(E)); endnumel(E)是分量个数log(numel(E))是最大熵除完之后值落在[0,1]。这个归一化在特征拼接时很关键否则分量数不同的样本之间能量熵没法直接比。5. 避坑与排查五个让特征矩阵报废的现场5.1 现象不同样本的特征维度对不上矩阵拼不起来原因几乎总是 IMF 层数不一致。CEEMD 是数据驱动的不同样本的局部极值点数不同emd返回的列数天然会变。如果在循环里直接[feat; newfeat]第一次运行到第七个样本就报错或者更糟——MATLAB 自动做维度广播悄悄把数据填成不一样的长度。解决就是第 2 章里那个K所有样本统一按K fix(log2(N)) - 1输出不够的补零、多的截断。补零之后要留意补出来的空列不要参与熵计算用ep 0判断过滤掉否则样本熵会对全零列算出NaN。5.2 现象分类准确率在训练集上很高换一批数据直接崩原因多半是相关系数用了带符号的值。IMF 与原始信号相位相反时相关系数为负如果阈值判定写的是r 0.1这些分量被剔除而在另一批数据上因为相位关系不同同一个物理分量又被保留。特征集的含义在样本之间悄悄变了模型学到的是一堆不稳定的映射。解决是取绝对值r abs(...)并且在特征表里记录每一列对应的是第几层 IMF、用了哪个阈值筛出来的出问题能回溯。5.3 现象样本熵算出来全是同一个数值或者一片NaN两种可能。第一种是r用了固定常数比如r 0.2而某个 IMF 的幅值只有 0.01所有模板距离都大于rB0熵值退化成-log(eps/eps) 0。第二种是m设成了 3 而分量长度不到 100 点模板数不足A和B全是 0。解决r一律用0.15 * std(该分量)逐列算分量长度短于 200 点时把m降到 2并且在使用前先检查numel(x) 4*m。真出现NaN说明这一段信号本身有问题比如静默段、传感器掉线应该把该样本整体标记出来别让它带着NaN进模型——大部分分类器遇到NaN不会报错而是把整行输出静默变成NaN等发现的时候已经跑完了。5.4 现象同一段信号跑两次特征值差得肉眼可见CEEMD 里的噪声是随机生成的不固定种子每次结果都不一样。单次差异可能不大但如果下游用的是对扰动敏感的模型重复实验的结论就会飘。解决在ceemd2的第一行加rng(seed, twister)seed固定成常数。同时注意randn的调用顺序——如果循环里还有其他随机数生成比如打乱样本顺序要保证 CEEMD 的种子在分解前单独重设一次否则打乱顺序一变分解结果跟着变。这个坑我踩过两次第二次才想明白是rng被下游函数覆盖了。5.5 现象熵特征在信号两端异常偏大或偏小EMD 家族的边界效应是老问题了。信号两端的极值点没法确定包络在端点附近会发散前几层 IMF 的头尾几十个点幅值明显失真。直接对整段算熵这些失真点会把熵值抬高。解决算熵之前把每层 IMF 两端各截掉 5%10%或者做镜像延拓再分解。截断的比例看信号长度2000 点截 5% 就是 100 点够用500 点以下建议做镜像延拓截太多剩下的点算样本熵又不稳了。还有一种做法是对整段信号先做一次去趋势detrend能减轻低频层端的发散但对高频层没用。6. 把三件套拼成特征表一次可复现的端到端验证前面几章都是零件最后拼一次完整流水线顺便说说怎么验证它没跑偏。function F buildFeature(x, opts) % 端到端CEEMD - 相关系数筛选 - 三熵特征 x x(:); N numel(x); K fix(log2(N)) - 1; [IMF, ~] ceemd2(x, opts.Nstd, opts.NE, K, opts.seed); r imfCorr(IMF, x); [keep, ~, ~] refineIMF(IMF, r, 0.1); F.r r(:); F.Hs zeros(1, numel(keep)); F.Hp zeros(1, numel(keep)); seg round(0.05 * N); % 截掉两端边界效应 for i 1:numel(keep) s IMF(seg : N - seg, keep(i)); F.Hs(i) sampen(s, 2, 0.15 * std(s)); F.Hp(i) permen(s, 4, 1); end F.He energyen(IMF(:, keep)); F.vec [F.r, F.Hs, F.Hp, F.He]; % 一行就是一个样本 end维度是K 2*numel(keep) 1因为keep依赖数据不同样本的keep长度可能不同所以F.r保留全K维熵只算筛选后的分量——用的时候熵那部分要按keep的最短长度对齐或者干脆固定keep取前min层。验证我一般走两条路。第一条是稳定性同一段信号用不同种子比如 1、2、3各跑一次特征向量的余弦相似度应该在 0.95 以上低于这个数说明NE给少了或者Nstd不合适。第二条是判别性拿两类样本各 50 条做五折交叉验证如果准确率明显高于用原始时域统计量的对照组一般高 8 个百分点以上说明这条流水线真的在提取有效信息而不是靠维度堆出来的过拟合。这两个检查每次改参数后我都会强制走一遍哪怕只改了NE。有一次图快只跑了稳定性没跑判别性结果Nstd从 0.15 调到 0.25 之后稳定性指标是过了但分类准确率掉了 12 个点——噪声幅值太大把冲击特征本身给抹平了。从那以后我每次调完 CEEMD 参数稳定性、判别性、重构误差三项都强制走一遍再往下做。希望帮到你。本文还有配套的精品资源点击获取