简介:围绕CEEMD完全集合经验模态分解、IMF分量相关系数筛选与信息熵特征提取的一套MATLAB实验材料,面向从事非线性非平稳信号处理、故障诊断与特征工程的学习者和研究人员,帮助打通从信号分解、分量筛选、熵特征量化到特征向量矩阵构建的完整分析链路。压缩包共4个文件,约60KB,以3个.m脚本为主体,承载分解、极值点计算与信息熵特征提取等核心流程,另附1个.xlsx样本数据用于直接跑通实验;体积轻量,便于快速查阅与二次修改。目前已有5248人学习下载,具备一定的社区验证基础。读者可据此复现CEEMD分解得到各阶IMF、用相关系数剔除冗余分量、逐分量计算信息熵并拼装成可供分类器使用的特征矩阵,同时从脚本组织方式中体会MATLAB在振动信号特征提取上的实现习惯,适合作为课程实验、课题预研或算法对比的起步模板。
1. 一条非平稳信号的两难:为什么 IMF 越多不等于信息越多
把一段非平稳信号直接送进分类器,准确率常常卡在六七成再难提高。很多人的第一反应是换模型、调超参,但问题往往出在特征侧:这段信号里同时存在趋势、周期性冲击和高频噪声,它们的尺度相差好几个数量级,模型看到的只是被搅在一起的数字序列。用 MATLAB 搭建的这条CEEMD 分解 → IMF 分量相关系数筛选 → 信息熵特征流水线,解决的正是这件事:先把混叠信号拆成若干本征模态分量,再判断哪些分量真正承载状态信息,最后把每个分量的波形差异压缩成几个标量。
反直觉的地方在于,IMF 数量越多并不等于信息越多。CEEMD 会把一段信号分解成八到十二个分量,但真正进入特征向量的通常只有三到五个——决定是哪几个的,是相关系数,并非能量大小。搞清楚这条链路每一环在做什么、为什么这么做,比背下一段调用代码有用得多。
2. CEEMD 的噪声对消:成对加噪为什么比平均 N 次更干净
2.1 筛分过程:IMF 到底是怎么被剥出来的
EMD 的全部机制就藏在一个循环里。对当前待处理序列 h(t),找出所有局部极大值点,用三次样条插值出上包络;极小值点插出下包络;两者取平均得到 m(t);然后 h ← h − m。这一步叫一次筛分,反复执行到上下包络的均值近似为零、过零点数与极值点数相差不超过一个,就得到一个 IMF。原信号减掉这个 IMF,残差继续筛,直到残差只剩下两三个极值点或已经单调为止。
function [h, nIter] = sifting(x, sdThr, maxIter) % 单次筛分:反复剥掉上下包络均值,直到满足停止准则 h = x(:); for nIter = 1:maxIter idxMax = find(h(2:end-1) > h(1:end-2) & h(2:end-1) >= h(3:end)) + 1; idxMin = find(h(2:end-1) < h(1:end-2) & h(2:end-1) <= h(3:end)) + 1; if numel(idxMax) < 2 || numel(idxMin) < 2, break; end % 极值不足,停止 t = (1:numel(h))'; upper = interp1(t(idxMax), h(idxMax), t, 'spline'); % 上包络 lower = interp1(t(idxMin), h(idxMin), t, 'spline'); % 下包络 m = (upper + lower) / 2; % 包络均值 hPrev = h; h = h - m; % 剥掉均值 sd = sum((hPrev - h).^2) / max(sum(hPrev.^2), eps); if sd < sdThr, break; end % 标准差停机准则 end这里有两个参数会直接改变分解结果的松紧程度。sdThr是标准差停机准则的阈值,常用 0.2~0.3:取值太小,筛分次数过多,IMF 的幅值会被反复平均压扁,物理意义变模糊;取值太大,分量之间又会互相串扰,本应分开的尺度混在同一个 IMF 里。maxIter是兜底上限,多数实现取 10 次左右,防止在噪声段陷入死循环。代码里的极值判定用了严格大于和大于等于的组合,这是为了避免平台区(相邻采样点数值相同)被重复计成极值,导致样条插值出现奇异点。
2.2 模态混叠与 EEMD 的残余噪声
筛分机制有一个固有缺陷:包络靠极值点插值,极值点的分布决定了包络的形状。当信号里出现间歇性冲击——某段时间极值密集,某段时间几乎平坦——样条在平坦区只能拉出一条长弧去连接相隔很远的两个极值点,包络均值被这条长弧严重污染,高频成分就被“搬”到了低频 IMF 里。这就是模态混叠,是 EMD 的血泪经验里绕不开的第一个坑。
EEMD 的思路很直接:往信号里叠加零均值白噪声再做分解,重复 N 次,把对应阶的 IMF 取平均。白噪声在整个频带上均匀分布,等于给每个尺度都补上了极值点,样条不再需要靠长弧撑场面,尺度分离自然稳定。因为噪声独立同分布、均值为零,理论上平均之后会抵消。但抵消是统计意义上的,残余噪声量级大约是 ε/√N(ε 为噪声幅值)。想让残噪降低一个数量级,重复次数得加一百倍——这就是很多人跑了一整夜仍然觉得“不够干净”的根源,也是 EEMD 在工程上最尴尬的地方。
2.3 互补对消:+n 与 −n 两次分解的机制
CEEMD 把随机性从“靠次数堆”改成了“靠结构对消”。它不再独立生成 N 段噪声,而是生成一段噪声 w,同时构造 x + w 和 x − w 两个输入,各自分解,再把两组 IMF 相加取半。加性噪声的一阶项在配对求和时直接抵消,残余噪声的衰减速度远快于 1/√N,通常几十对就能达到 EEMD 几百次的效果。
function [IMFs, res] = ceemd(x, Npair, ratio, seed) % x 原始信号,列向量 % Npair 互补噪声对数,常用 50~100 % ratio 噪声幅值系数,相对 std(x),常用 0.1~0.3 % seed 随机种子,固定后结果可复现 rng(seed); x = x(:); sd = std(x); N = numel(x); pool = {}; % 收集每一次分解的 IMF 矩阵 for k = 1:Npair w = ratio * sd * randn(N,1); % 一段噪声,正负各用一次 [imfP, ~] = emd_decompose(x + w); % 加正噪声分解 [imfN, ~] = emd_decompose(x - w); % 加负噪声分解 pool{end+1} = imfP; pool{end+1} = imfN; end K = max(cellfun(@(A) size(A,2), pool)); % 各次阶数不齐,取最大阶 IMFs = zeros(N, K); for i = 1:numel(pool) A = pool{i}; IMFs(:, 1:size(A,2)) = IMFs(:, 1:size(A,2)) + A; % 缺的高阶补零 end IMFs = IMFs / numel(pool); % 2*Npair 次分解的集成平均 res = x - sum(IMFs, 2); % 平均后不再精确重构,残差留用ratio控制噪声幅值,取 0.1~0.3 是经验区间:太小起不到补极值点的作用,太大会让低阶 IMF 被噪声淹没,反而降低信噪比。Npair取 50~100 已经能拿到稳定结果,再往上收益递减,而计算量是线性增长的。代码里各次分解的 IMF 阶数常常不一致,直接截断到最小阶会丢掉高频细节,所以按最大阶对齐;缺的分量补零只是权宜之计,更严谨的做法是用该次分解的残差补齐缺失的高阶项,保证每一组都满足 ΣIMF + res = x。严格来说,EMD 不是线性算子,E(x+w) + E(x−w) ≠ 2E(x),所以这里的对消只是近似。CEEMDAN 改成了逐级加噪、每级只取第一阶 IMF 的集成均值,并把残差显式保留下来,重构是精确成立的——这是它比 CEEMD 更值得优先考虑的地方。
3. 相关系数筛选:IMF 和原信号的那个数到底在度量什么
3.1 正交性假设下,相关系数其实是贡献占比
IMF 之间近似正交,这一点是理解筛选逻辑的前提。如果各分量彼此正交,那么原信号的能量约等于各分量能量之和,而分量与原信号的相关系数平方之和也近似为 1。换句话说,ρ_i 反映的是第 i 个尺度成分在原信号里占多大比重、与原信号的整体走势有多强的线性关联。ρ 接近 1,说明这个分量几乎是原信号的缩影,通常出现在趋势项或者主导频率所在的那一阶;ρ 接近 0,说明这个分量几乎不携带原信号的结构,多半是分解过程自己造出来的伪分量。
伪分量的判别还有一个容易忽略的维度:除了看 IMF 与原信号的相关性,还要看 IMF 与相邻 IMF 的相关性。如果某个分量与原信号几乎不相关,却和前一个分量高度相关,那它大概率不是独立的物理尺度,而是筛分停止判据不严格留下的残留。工程上常见的处理是直接剔除相关系数低于阈值的分量,但如果这些分量还要参加重构,就得先做这一步判别,否则会把伪分量当成有效信息一起加回去。
3.2 阈值怎么定:固定值、自适应还是找拐点
相关系数阈值是这条流水线上最容易被随意设定的参数。常见的三种策略差别很大:固定阈值(0.1 或 0.2)实现简单,但不同信号的相关性分布差异很大,一个在振动信号上工作良好的阈值换到脑电信号上可能一个分量都筛不出来;相对阈值(ρ_i > α·max(ρ))随信号自适应,α 取 0.1~0.2,能保证至少保留主导分量;拐点法把 ρ 序列按阶次画出来,找下降最陡的那个位置作为截断点,物理含义最清楚,代价是需要人工看一次图。
我更倾向双下限的组合:max(0.1, 0.1*max(rho))。绝对下限挡住那些几乎不相关的伪分量,相对下限保证在信号整体相关性偏低时不至于把所有分量都扔掉。还有一个必须单独处理的对象是最后的残差项——它是趋势,与原信号的相关性往往很高,但它不包含任何振荡结构,计算样本熵或者排列熵没有意义,必须显式排除,不要让它混进特征向量。
function [rho, keep] = imf_select(x, IMFs, thr) % IMFs 为 N×K 矩阵,每列一个 IMF,不含残差项 xc = x(:) - mean(x); K = size(IMFs, 2); rho = zeros(1, K); for k = 1:K ck = IMFs(:,k) - mean(IMFs(:,k)); rho(k) = (xc' * ck) / (norm(xc) * norm(ck) + eps); % 皮尔逊相关系数 end if nargin < 3 || isempty(thr) thr = max(0.1, 0.1 * max(rho)); % 绝对下限 + 相对下限 end keep = find(rho >= thr); % 保留的分量序号 end分母加eps是为了防止某个分量在极短窗内被分解成全零序列时出现除零。keep返回的是序号而不是矩阵本身,后面组装特征时按序号取列,可以保持阶次信息——这一点很重要,因为熵值和阶次强相关,打乱顺序会让特征失去可解释性。
3.3 一个反直觉的取舍:IMF1 不要急着扔
很多教程会说“IMF1 是高频噪声,直接剔除”。这个结论在平稳、高信噪比的场景下成立,但在故障诊断这类任务里经常翻车:轴承局部剥落、齿轮断齿产生的周期性冲击,其特征频带恰恰落在最高阶的 IMF 里,把 IMF1 一律丢掉等于把最有判别力的信息扔了。合理的做法是让相关系数说话——如果 IMF1 与原信号的相关性确实很低,说明它主要是噪声;如果相关性可观,那它承载的就是冲击结构,应该保留。这条经验没法写进固定阈值,只能靠对具体信号的观察积累。
4. 信息熵特征:怎么把一段波形压成能分类的向量
4.1 熵在度量什么,为什么对 IMF 有效
熵度量的是一段序列的不确定性或者不规则程度。同一个系统在正常状态和故障状态下,振动或者电生理信号的规则性会发生系统性变化:出现周期性冲击时,波形重复性增强,熵值下降;出现宽带噪声或者混沌行为时,熵值上升。把这种变化直接量化成标量,比让分类器自己去卷积核里找要省事得多。
关键在于,熵必须施加在单一尺度的分量上才有判别力。原始信号里混着趋势、基频、谐波和噪声,整体算出来的熵是所有尺度复杂度的加权混合,不同故障模式之间的差异容易被相互抵消。CEEMD 把尺度拆开之后,每个 IMF 的熵各自反映那一层的规则性,形成一组天然的多维特征。这也是这套流水线的核心价值:利用分解提供尺度分辨能力,熵只负责在单个尺度上量化。
4.2 样本熵与排列熵:参数比公式更影响结果
样本熵的计算逻辑是:把序列切成若干长度为 m 的片段,统计两两之间距离小于容限 r 的匹配对数 B;再统计长度 m+1 的匹配对数 A;样本熵取 −ln(A/B)。它的优点是抗短时干扰、对数据长度要求不高,缺点是对 r 极其敏感。
function se = sampen(x, m, r) x = x(:)'; N = numel(x); B = 0; A = 0; for i = 1 : N-m for j = i+1 : N-m if max(abs(x(i:i+m-1) - x(j:j+m-1))) <= r, B = B + 1; end % m 维匹配 if max(abs(x(i:i+m) - x(j:j+m) )) <= r, A = A + 1; end % m+1 维匹配 end end se = -log(A / max(B, eps)); % 两个条件概率之比取对数 end距离用的是切比雪夫范数(各维差值的最大值),这是样本熵的标准定义。m 一般取 2;r 取 0.1~0.25 倍标准差,但这个标准差必须是当前这个 IMF 自己的标准差,不能用全局统一值。这一点是踩坑率最高的参数问题,下文还会展开。
排列熵走的是另一条路:不比较数值距离,只比较大小顺序。长度为 m 的片段里,元素的相对大小构成一种排列模式,一共 m! 种;统计所有模式出现的频率,求香农熵再除以 ln(m!) 归一化到 [0,1]。
function pe = permen(x, m, tau) x = x(:)'; N = numel(x); n = N - (m-1)*tau; P = zeros(n, m); for i = 1:n P(i,:) = x(i : tau : i+(m-1)*tau); % 按延迟 tau 取 m 个点 end [~, ord] = sort(P, 2); % 升序名次即排列模式 code = (ord-1) * (m.^(0:m-1))'; % 排列编码成整数便于计数 cnt = accumarray(code+1, 1, [m^m 1]); p = cnt(cnt>0) / n; pe = -sum(p .* log(p)) / log(factorial(m));% 归一化到 [0,1] endm 取 3~7,τ 通常取 1。m 太小时最多只有 6 种模式,分辨率不够;m 太大则需要的数据长度急剧上升,而且容易产生大量未曾出现的模式,导致熵值被系统性高估。排列熵的计算量比样本熵小一个量级,对冲击类信号尤其敏感,很适合放在滑窗实时计算的场合。这两种熵并不冲突,同时使用能覆盖不同的复杂度侧面。
4.3 能量熵与特征向量组装
第三个维度是能量熵。先算每个 IMF 的能量 E_i = Σ c_i²,归一化成概率 p_i = E_i / ΣE,再取 −Σ p_i ln p_i 并除以 ln K 归一化。它描述的是能量在各尺度上分布得有多均匀:单一频率主导时能量集中在某一阶,熵值低;宽带冲击时能量分散,熵值高。计算代价几乎为零,但信息量和前两种熵不重叠,值得一起放进去。
E = sum(IMFs.^2, 1); % 各保留分量的能量 p = E / sum(E); H = -sum(p .* log(p)) / log(numel(p)); % 归一化能量熵,范围 [0,1]完整的特征向量按 [样本熵 × K, 排列熵 × K, 能量熵 × 1] 拼装,维度为 2K+1。K 是筛选后保留的分量数,通常在 3~5 之间,对应特征维数 7~11 维,这个量级对 SVM、随机森林这类分类器非常友好,不容易过拟合。
单条信号只能产出一个向量,要构成训练集必须分帧:
win = 1024; hop = 256; % 窗长与步长 nSeg = floor((numel(x) - win) / hop) + 1; F = zeros(nSeg, 2*K + 1); for s = 1:nSeg seg = x((s-1)*hop + (1:win)); [IMFs, ~] = ceemd(seg, 50, 0.2, 42); % 每帧独立分解,种子固定 [~, keep] = imf_select(seg, IMFs); % 逐分量算样本熵、排列熵,再算能量熵,填入 F(s,:) end窗长 1024、步长 256 对应 75% 重叠,是为了让相邻样本之间有足够的相关性又不至于完全重复。这里有个绕不开的代价:逐帧做 CEEMD,计算量是整段分解的 nSeg 倍。如果帧数上千,跑一整晚毫不夸张。折中方案是整段分解一次得到全局 IMF,再把这些 IMF 按同样的窗切分算熵——代价是每帧的 IMF 不再是该帧自身的自适应分解,边界处会有轻微失配,但对缓变状态识别通常可以接受。
5. 避坑与排查:这条流水线最容易翻车的五个环节
5.1 端点发散,相关系数整体偏高
现象:分解出来的每个 IMF 在首尾各几十个点幅度异常放大,相关系数普遍虚高,样本熵要么算出来接近 0,要么直接报错。
原因:三次样条在数据端点没有约束,只能靠最外侧的极值点外推,误差按指数向外放大,这是 EMD 系方法的固有缺陷,不是实现 bug。
解决:分解前做镜像延拓或者极值延拓,把信号两端各延长一段再插值;或者更简单粗暴——分解后统一裁掉两端各 10% 的数据,只在中段计算相关系数和熵。后者实现成本最低,代价是有效样本变短,窗长需要相应加大。
5.2 各帧 IMF 阶数不一致,特征矩阵拼不起来
现象:循环跑到中途报“维度不匹配”,或者强行拼出来的矩阵大量NaN。
原因:每帧信号长度、噪声实现都不同,筛分停止的时机不同,分解出的 IMF 个数自然不一样,有时是 8 阶,有时是 10 阶。
解决:先固定保留前 K 阶(K 取所有帧中最小阶数,或者直接取 4),其余阶数与残差合并;也可以在 CEEMD 内部按最大阶对齐、缺项用该次残差补齐。绝不能用cell2mat直接拼,那是必翻车。
5.3 样本熵的 r 用了全局固定值
现象:低阶 IMF 的样本熵正常,高阶 IMF 的样本熵全部饱和到接近 0 或者接近上界,特征向量区分度极差。
原因:IMF 的幅值随阶次下降,高频分量幅值可能是低频分量的几十倍。用同一个 r,对幅值小的分量来说容限过宽,几乎所有片段都匹配,熵值被压到接近 0。
解决:每个 IMF 单独计算r = 0.2 * std(imf_k),逐分量自适应。这是这条流水线上最值得优先检查的一条。
5.4 随机种子没固定,同一段信号跑两次结果不一样
现象:重复运行同一份数据,特征值有百分之几的浮动,交叉验证的准确率忽高忽低。
原因:CEEMD 的噪声实现是随机的,randn每次调用都不一样,分解结果因此带有随机性,下游的相关系数、熵全部跟着变化。
解决:在每次分解入口处显式rng(seed),并且在论文或者报告里写清楚种子值。同时把Npair提到 50 以上,把随机性压到可忽略的量级。可复现性不是洁癖,是排查问题的前提。
5.5 归一化时机错误,把熵信息抹掉了
现象:模型在训练集上表现很好,换一批数据准确率断崖式下跌。
原因:对整个数据集(含测试集)算均值和方差做归一化,测试集信息泄漏进训练过程;或者对每一个窗口单独做零均值单位方差归一化,把窗口之间的幅值差异抹平,而能量熵恰恰依赖这个差异。
解决:先划分训练集和测试集,归一化参数只用训练集统计量;对依赖幅值信息的特征(能量熵、能量占比)不做逐窗标准化,只在必要时做全局缩放。
6. 设计取舍:特征维数、计算代价和可解释性怎么平衡
这套流水线的每一个参数都在三者之间做交换,没有全局最优解,只有场景匹配。
| 配置项 | 取小值的后果 | 取大值的后果 | 工程建议 |
|---|---|---|---|
| 噪声对数 Npair | 残噪大、结果不可复现 | 耗时线性上升 | 50~100,种子固定 |
| 噪声幅值 ratio | 补极值不足、仍有混叠 | 低阶 IMF 被淹没 | 0.1~0.3,先看 IMF1 是否被淹 |
| 保留分量数 K | 特征维数低但可能丢判别信息 | 维数膨胀、混入伪分量 | 3~5,用相关系数定 |
| 熵的种类 | 信息面窄 | 维数翻倍、冗余 | 样本熵 + 排列熵 + 能量熵 |
在验证方法上,我建议不要一上来就训分类器看准确率——那个数字受模型和划分方式影响太大,看不清特征本身的贡献。先做一个增量消融:只用原始信号的统计量(均值、方差、峭度)做一组基线;只做 CEEMD 不筛选、全分量算熵做第二组;筛选后再算熵做第三组。三组的准确率差值,才是分解和筛选各自的真实增益。如果第二组和第三组差不多,说明筛选没起作用,问题多半在阈值;如果第一组和第三组差别很小,那这套流水线对这个信号本身就没多少价值,换思路比赛调参更划算。
另一个快速筛选办法是算类间可分性,比如 Fisher 判别比(类间方差除以类内方差之和),它不需要训练任何模型,几秒钟就能告诉你某个特征维度值不值得留。维数上去之后再做一次相关性检查,两列特征相关系数超过 0.9 的,留一个就够,另一列纯属给分类器添乱。
计算代价的大头永远在 CEEMD:2×Npair 次筛分,每次筛分内部还要迭代十几次样条插值。熵的计算相比之下可以忽略。所以真要提速,方向是减少 Npair、缓存分解结果、或者把逐帧分解改成整段分解后切窗,而不是去优化熵函数里的那个双重循环。
从可解释性角度看,这套方法的优势恰恰在于每一维都有明确含义:第 k 阶 IMF 的样本熵高,意味着那一层尺度的波形更不规则;能量熵高,意味着能量在尺度上分布更分散。这个性质在需要给出判断依据的场合非常值钱,也是它比端到端深度模型更适合小样本场景的原因。
从那以后我每次搭这条流水线,都会先把分解结果画出来看一眼:IMF 两端有没有发散、阶数是不是逐帧一致、相关系数序列的拐点在哪一阶。这三张图看完,后面的参数基本就不用猜了。希望这些踩过的坑能帮你少走一些弯路。
本文还有配套的精品资源,点击获取