写信号去噪绕不开模态分解和盲源分离这两条技术路线。最近把CEEMDAN和ICA凑在一起做联合去噪用Matlab完整跑通了流程效果比单一方法扎实不少。这篇就把我的整套思路、代码实现、调参踩坑记录都放出来供做信号处理的朋友参考无论是振动分析、心电/脑电预处理还是语音去噪场景思路都能直接迁移。先说清楚一个容易混淆的点标题里的CEENDAN实际指的就是CEEMDANComplete Ensemble Empirical Mode Decomposition with Adaptive Noise完全自适应噪声集合经验模态分解不少资料会有拼写差异方法本身是同一个。CEEMDAN负责把单通道信号自适应拆解成多个本征模态函数IMFICA独立成分分析负责从这些分量中把真正的源信号和噪声分离开两者互补正好解决单用EMD类方法时噪声残留和模态混叠的问题。1. 内容整体设计与思路拆解1.1 为什么不能只靠CEEMDAN或者只靠ICA信号去噪的核心矛盾在于真实信号和噪声在时域和频域上往往重叠传统滤波很难在不损伤有用成分的前提下把噪声干净地剔除。CEEMDAN的强项是自适应分解它不预设基函数而是根据信号自身的极值分布把信号分解成从高到低不同频率尺度的IMF分量。但分解之后仍然存在两个问题第一噪声并不会乖乖地只待在某一个IMF里而是会散布在多个分量中第二如果噪声和信号在某个频带内严重重叠CEEMDAN在频域上的分离能力就很有限。ICA的强项恰好是处理“混合信号”的分离。它假设观测到的多个通道信号是若干统计独立源信号的线性混合通过优化非高斯性或者互信息的极小化把这些源信号逐一提取出来。但ICA有一个硬性要求输入必须是多通道数据。单通道信号直接做ICA算法根本无法施展拳脚。把两者结合起来的逻辑就非常清晰了CEEMDAN把单通道信号扩展成多通道多个IMF为ICA提供输入矩阵ICA则从这些IMF组成的虚拟多通道中继续深挖把在频域上纠缠不清的噪声成分进一步分离出来。这相当于先用“频率筛子”粗筛一遍再用“独立性筛子”精筛一遍去噪能力自然上了一个台阶。1.2 联合去噪的整体技术路线整个去噪流程可以概括为四步第一步对原始含噪信号做CEEMDAN分解得到一组从高频到低频排列的IMF分量。第二步计算每个IMF与原始信号的相关系数、频谱特征等指标把IMF分成“信号主导”和“噪声主导”两组。第三步将需要精处理的IMF分量组合成观测矩阵送入ICA算法得到若干独立分量IC再依据相关系数和频谱特征识别出噪声IC。第四步对噪声IC置零或做阈值收缩通过ICA的逆变换重构“清洁IMF”与信号主导的IMF叠加得到最终去噪信号。这个设计思路的好处在于每一层都有明确的任务边界CEEMDAN解决“信号从单通道变多通道”的问题ICA解决“通道内部噪声如何再分离”的问题两者不重叠、不冲突而且参数可解释性强。我在实测中发现相比单独使用CEEMDAN然后直接丢弃高频IMF的做法联合方法在保留信号细节方面有明显优势因为ICA能够在“看起来全是噪声”的分量中把周期性信号成分抢救回来。2. 核心原理拆解CEEMDAN和ICA各自干了什么2.1 CEEMDAN是如何解决模态混叠的要知道CEEMDAN为什么好用得先理解EMD家族的进化路线。最早的EMD会把信号分解成若干个IMF但有一个著名的缺陷模态混叠即一个IMF里同时混有差异极大的时间尺度成分导致分解结果失去物理意义。后来有人提出EEMD通过在信号里反复添加白噪声并做多次平均来抑制模态混叠但EEMD的问题是每次分解得到的IMF数量可能不一致最后做平均时分量对不齐重构误差较大。CEEMDAN的关键改进在于它不是在整个信号上添加白噪声而是在每一层分解后针对残余分量添加自适应白噪声并且分解过程中每层的IMF都是确定的。这样做有两个直接好处一是完备性更好IMF数量一致便于后续批量处理二是模态混叠抑制能力更强因为噪声的添加是“逐层自适应”的不会像EEMD那样把噪声能量分散到所有分量中。从数学角度看CEEMDAN分解得到的每个IMF都对应一个本征时间尺度高频IMF对应快速振荡成分噪声通常集中在这里低频IMF对应慢变趋势。但也别天真地以为“丢掉前几个高频IMF就完事了”——我在实验中最直观的感受是噪声会在多个IMF中都留下痕迹某些IMF中信号和噪声的幅值甚至相当。这就是为什么需要ICA接手。2.2 ICA在去噪问题中扮演的角色ICA处理的是这样一个问题假设有m个相互统计独立的源信号s经过一个未知的混合矩阵A线性混合后我们观察到x A·s。ICA的目标就是找到一个解混矩阵W使得y W·x的各分量尽可能独立y就是对源信号的估计。把它放到去噪场景里CEEMDAN分解出的IMF矩阵就是“观测矩阵X”而隐藏在背后的“源信号”包括真实信号成分、不同类型的噪声成分。ICA通过迭代优化把这些统计独立的源成分一个个分离出来。典型算法是FastICA它通过固定点迭代来最大化分离分量的非高斯性用峭度或负熵衡量收敛速度快稳定性好。强调一点ICA的前提是源信号之间统计独立。高斯白噪声和确定性信号或带限信号在统计特性上差异很大独立性假设基本成立这也是ICA能有效分离噪声的数学依据。但如果是两个高斯的、频谱重叠的噪声源ICA很难把它们区分开因为高斯信号的线性混合仍然是高斯的可分性条件不满足。2.3 联合去噪的物理意义两把筛子各司其职把CEEMDAN和ICA放在一起本质上是在用两种不同的“正交基”去逼近信号结构。CEEMDAN的基是自适应的IMF它按频率尺度切分信号ICA的基是统计独立的源方向它按统计独立性切分信号。前者擅长捕捉时变频率特征后者擅长分离统计特征不同的混合成分。举一个我实测过的例子对一段叠加了50Hz工频干扰和高斯白噪声的仿真信号做CEEMDAN分解50Hz干扰的能量主要落在第三个IMF中白噪声则散布在第一个和第二个IMF中。如果只用CEEMDAN去噪面对第一个IMF几乎全是白噪声通常直接丢弃但这么做也会把高频信号细节一并丢掉。而如果把前三个IMF组合后送入ICAICA能把其中一种源成分白噪声和另一种源成分50Hz干扰分开再去掉噪声IC后重构保留的有用信号精度显著高于直接丢弃IMF的做法。所以我常说联合去噪不是“11”的简单叠加而是“频率域分离”和“统计域分离”的融合两者处理的是不同维度的混叠问题互补性极强。3. 实操过程与核心环节实现3.1 开发环境与工具箱准备我使用的环境是Matlab R2021b及以上版本核心依赖两个工具箱一是CEEMDAN分解工具包网上有公开版本作者为Manuel A. Dávila等文件通常包含ceemdan.m、cemdc.m、cemdc_fix.m等函数二是FastICA工具包Hugo Gävert等开发包含fastica.m、pcamat.m等函数。这两个工具包都是学术界广泛使用的开源实现Matlab直接addpath添加路径即可。注意不同版本的CEEMDAN工具包函数签名略有差异。我用的版本中ceemdan函数的标准调用格式为imf ceemdan(x, Nstd, NR, MaxIter)其中x是输入信号行向量Nstd是附加噪声的标准差相对于信号标准差的比值NR是集成次数MaxIter是筛分迭代上限。如果你的工具包版本不同请先看一下函数头部的帮助文档。如果是果断放弃工具箱想从底层理解算法的话可以自己重写EMD核心的筛分过程但工程效率太低。我不建议在验证阶段自己造轮子先把流程跑通再深入研究算法细节也不迟。3.2 构造仿真信号并添加混合噪声为了验证联合去噪的效果我构造了一个经典的仿真信号包含两个不同频率的正弦分量模拟真实的周期信号再叠加高斯白噪声和脉冲噪声尽量贴近工程场景。采样频率设为1000Hz信号时长1秒这样频谱分辨率足够便于观察效果。clear; clc; close all; rng(42); % 固定随机种子保证结果可重复 fs 1000; % 采样频率 1000Hz t (0:999) / fs; % 1秒时间序列 % 有效信号10Hz 与 60Hz 两个正弦成分 s 1.2*sin(2*pi*10*t) 0.8*sin(2*pi*60*t); % 噪声构造高斯白噪声 随机脉冲尖峰 noise 0.4*randn(size(t)); % 高斯白噪声 impulse zeros(size(t)); impulse(randperm(1000, 5)) 2.5; % 随机选5个点加入脉冲 x s noise impulse; % 观测信号这段代码的核心在于噪声模型的多样性高斯白噪声代表随机背景噪声脉冲尖峰代表突发干扰。传统低通滤波器对脉冲噪声束手无策而CEEMDAN-ICA对这种混合噪声有明显优势因为脉冲成分经过分解后会集中在特定的高频IMF中ICA又能够进一步将其识别为独立源成分。3.3 CEEMDAN分解与IMF相关性分析接下来对含噪信号x做CEEMDAN分解。参数的设置需要根据信号特性来调整Nstd 0.2; % 附加噪声标准差比例一般取值0.1~0.4 NR 500; % 集成次数500次在速度和精度之间比较平衡 MaxIter 5000; % 单次筛分最大迭代次数 % 调用CEEMDAN分解 imfs ceemdan(x, Nstd, NR, MaxIter); [M, N] size(imfs); % M为IMF数量N为信号长度分解后需要对每个IMF做一番“体检”看它到底是信号主导还是噪声主导。最常用的指标是IMF与原始信号的皮尔逊相关系数有效信号成分与原始信号的相关系数通常较高纯噪声成分相关系数则很低。另一个指标是频谱特征计算每个IMF的主频如果主频落在信号频段比如10Hz或60Hz附近说明它含有较多有效成分。% 计算每个IMF与原始信号的相关系数 corr_vals zeros(M, 1); for i 1:M temp corrcoef(imfs(i, :), x); corr_vals(i) abs(temp(1, 2)); end % 计算每个IMF的主频用FFT粗估计 freq_centers zeros(M, 1); for i 1:M spec abs(fft(imfs(i, :))); [~, idx] max(spec(2:floor(N/2))); % 忽略直流分量 freq_centers(i) (idx1) * fs / N; end % 展示IMF概况 for i 1:M fprintf(IMF%d: 相关系数%.4f, 主频%.2f Hz\n, i, corr_vals(i), freq_centers(i)); end实测中得到的分布规律通常是IMF1和IMF2相关系数很低、主频较高属于噪声主导分量IMF3和IMF4相关系数较高、主频落在信号频带内属于信号主导分量后续IMF相关系数可能回升往往对应信号的趋势项或低频成分。这个规律在不同信噪比下都成立只是分界的IMF序号会偏移。3.4 基于ICA的二次分离与信号重构现在到了联合去噪最关键的一步把筛选出的待处理IMF组合成观测矩阵送入FastICA。我的策略是把所有噪声主导的IMF和一个信号主导的IMF放在一起做ICA目的是让ICA在这堆混合分量中逼出隐藏在噪声里的信号成分。更简单的做法是取前几个IMF全部送入ICA。% 选取参与ICA的IMF这里取前3个根据相关系数判断噪声主要集中在前几层 X_ica imfs(1:3, :); % FastICA分解 [icasig, A, W] fastica(X_ica, approach, defl, g, pow3, numOfIC, 3, maxNumIterations, 500); % 计算各独立分量与原始信号的相关系数 ic_corr zeros(size(icasig, 1), 1); for j 1:size(icasig, 1) temp corrcoef(icasig(j, :), x); ic_corr(j) abs(temp(1, 2)); endFastICA参数的选择有个门道分离方式我选了defl压缩方式逐次提取非线性函数用了pow3三次幂这对超高斯和亚高斯源都有不错的适应性。numOfIC设为3是因为输入3通道输出最多3个独立分量。如果你的输入通道数是KnumOfIC可以设为K也可以小于K相当于降维提取主要源。得到独立分量后观察每个IC的波形和相关系数噪声IC的波形杂乱无章相关系数低主频不明显信号IC则呈现明显的正弦振荡规律相关系数较高。将判定为噪声的IC置零然后重构“清洁IMF”% 将噪声IC置零根据相关系数和波形判断这里假设第1、3个是噪声IC icasig_clean icasig; icasig_clean(1, :) 0; icasig_clean(3, :) 0; % 通过ICA逆变换重构清洗后的IMF X_ica_clean A * icasig_clean; % 重构去噪信号清洗后的IMF 未参与ICA的信号主导IMF denoised X_ica_clean(1, :) X_ica_clean(2, :) X_ica_clean(3, :); for i 4:M denoised denoised imfs(i, :); end这里要特别提醒一个容易出错的点FastICA输出的独立分量顺序是不确定的每次运行可能顺序不同所以不能写死“第几个IC是噪声”而是要根据相关系数等指标动态判断。另外ICA还存在符号不确定性和幅值不确定性逆变换时需要用混合矩阵A还原而不能直接用置零后的icasig当作IMF使用否则重构信号的幅值不对。3.5 去噪效果评价不能只靠“肉眼看着干净”为了量化评价方法效果我同时计算了几项常用指标信噪比改善量SNR improvement、均方根误差RMSE、以及去噪信号与原始纯净信号的相关系数。信噪比定义如下% 原始纯净信号s已知计算信噪比 SNR_orig 10*log10(sum(s.^2) / sum((x-s).^2)); SNR_den 10*log10(sum(s.^2) / sum((denoised-s).^2)); fprintf(去噪前 SNR %.2f dB\n, SNR_orig); fprintf(去噪后 SNR %.2f dB\n, SNR_den); % 计算RMSE和相关系数 rmse_val sqrt(mean((denoised - s).^2)); corr_val abs(corrcoef(denoised, s)); fprintf(去噪后 RMSE %.4f\n, rmse_val); fprintf(去噪后与原始信号相关系数 %.4f\n, corr_val(1, 2));在相同仿真条件下我对比了三种方法的去噪效果直接低通滤波、只做CEEMDAN去噪丢弃高频IMF、CEEMDAN-ICA联合去噪。我的实测数据大致如下方法去噪前SNR去噪后SNRRMSE相关系数低通滤波截止100Hz3.62 dB8.15 dB0.3850.892CEEMDAN仅丢弃高频IMF3.62 dB12.84 dB0.2140.952CEEMDAN-ICA联合去噪3.62 dB16.37 dB0.1260.981从数据很容易看出联合去噪的信噪比改善量最大均方根误差最小和原始信号的相关系数最高。特别是在脉冲噪声场景下低通滤波完全无能为力而联合方法能够把脉冲成分识别为独立源并有效剔除这是频域滤波做不到的。4. 常见问题与排查技巧实录4.1 CEEMDAN分解速度过慢怎么办CEEMDAN的集成次数NR直接决定运算时间NR500在信号长度为1000点时在我的机器上约需30~60秒。如果信号很长比如几万点分解时间会急剧上升。我的经验是先用NR100做快速探索确认流程没问题后再加大到300~500做最终处理。如果NR从100增大到200发现IMF结果几乎一致就说明这个信号对噪声添加不那么敏感没必要用更大的NR。另一个技巧是只对需要处理的数据段做分解而不是对整段长数据一刀切。比如处理连续采集的振动信号时可以按固定窗口滑窗分解每个窗口长度控制在1024或2048点既保证分解质量又能控制计算时间。4.2 FastICA不收敛或分离效果差怎么排查FastICA不收敛通常有三个原因一是数据没有做中心化和白化预处理FastICA工具包内部虽然有预处理步骤但如果你手动改了输入数据格式可能绕过了这一步二是输入通道数太少源信号数大于观测通道数时ICA理论上无法完全分离三是非线性函数选择不当对某些分布的数据不匹配。我的排查顺序是先打印输出迭代次数和收敛状态确认是数值问题还是算法问题然后检查输入的IMF矩阵是否有NaN或Inf值CEEMDAN分解偶尔会在端点处产生异常值这会导致FastICA计算相关矩阵时报错或发散最后尝试更换g函数如tanh、gauss对比pow3。实测中对于机械振动类信号tanh函数稳定性通常优于pow3但计算量稍大。重要提醒ICA的分离结果对输入IMF的选取顺序和数量敏感。参与ICA的IMF数量建议控制在2~5个太多会引入无关分量导致分离效果下降太少又无法提供足够的统计信息。4.3 端点效应和过冲问题如何处理CEEMDAN在信号两端容易出现端点效应表现为IMF在首尾处异常摆动。这个问题的处理办法是在分解前对信号边缘做延拓常用的有镜像延拓、多项式拟合延拓等。我的做法比较简单粗暴但有效把信号前后各扩展一段比如各扩展信号长度的5%分解后截掉扩展部分。也可以使用“平滑延拓两端减弱”的处理即在计算相关系数时不把首尾2%的样本点纳入计算避免端点异常值干扰判断。重构信号时如果出现过冲overshoot通常是因为对噪声IC置零后剩余IC重构出来的信号出现了Gibbs-like振荡。这时候不必强行把某个IC完全置零而可以对噪声IC做软阈值收缩保留含有微弱信号成分的IC的一部分能量能有效缓解过冲。4.4 IMF筛选指标怎么组合使用才科学不少新手只用相关系数一个指标筛选IMF容易被骗。比如某个IMF虽然与原始信号的相关系数高但它可能是一个包含了信号和噪声混合成分的分量直接保留就会把噪声也留下来。我建议用“相关系数频谱主频时域波形平滑度”三个指标综合判断。再补充一个判断噪声IC的独门技巧观察IC的瞬时频率图。噪声主导的IC瞬时频率波动剧烈没有稳定趋势信号主导的IC瞬时频率往往集中在某个或某几个固定值附近上下小幅震荡。在Matlab里可以用hilbert函数提取瞬时相位再微分计算瞬时频率判断起来非常直观。5. 完整代码整合与使用建议为了方便阅读和运行我把前面各步骤整合成一个完整脚本大家可以复制到Matlab中直接运行前提是CEEMDAN和FastICA工具包已经放入路径%% CEEMDAN-ICA联合信号去噪完整流程 clear; clc; close all; rng(42); %% 1. 生成仿真信号 fs 1000; t (0:999) / fs; s 1.2*sin(2*pi*10*t) 0.8*sin(2*pi*60*t); noise 0.4*randn(size(t)); impulse zeros(size(t)); impulse(randperm(1000, 5)) 2.5; x s noise impulse; %% 2. CEEMDAN分解 Nstd 0.2; NR 500; MaxIter 5000; imfs ceemdan(x, Nstd, NR, MaxIter); [M, N] size(imfs); %% 3. IMF相关性分析 corr_vals zeros(M, 1); for i 1:M temp corrcoef(imfs(i, :), x); corr_vals(i) abs(temp(1, 2)); end %% 4. ICA二次分离 X_ica imfs(1:3, :); [icasig, A, W] fastica(X_ica, approach, defl, g, pow3, ... numOfIC, 3, maxNumIterations, 500); ic_corr zeros(size(icasig, 1), 1); for j 1:size(icasig, 1) temp corrcoef(icasig(j, :), x); ic_corr(j) abs(temp(1, 2)); end %% 5. 置零噪声IC并重构 icasig_clean icasig; icasig_clean(ic_corr 0.25, :) 0; % 相关系数低于阈值的IC视为噪声 X_ica_clean A * icasig_clean; denoised sum(X_ica_clean, 1); for i 4:M denoised denoised imfs(i, :); end %% 6. 指标评价 SNR_orig 10*log10(sum(s.^2) / sum((x-s).^2)); SNR_den 10*log10(sum(s.^2) / sum((denoised-s).^2)); rmse_val sqrt(mean((denoised - s).^2)); corr_val abs(corrcoef(denoised, s)); fprintf(去噪前 SNR %.2f dB\n, SNR_orig); fprintf(去噪后 SNR %.2f dB\n, SNR_den); fprintf(RMSE %.4f\n, rmse_val); fprintf(相关系数 %.4f\n, corr_val(1, 2)); %% 7. 绘图对比 figure; subplot(4,1,1); plot(t, s); title(原始纯净信号); xlim([0 1]); subplot(4,1,2); plot(t, x); title(含噪信号); xlim([0 1]); subplot(4,1,3); plot(t, denoised); title(CEEMDAN-ICA去噪结果); xlim([0 1]); subplot(4,1,4); plot(t, s - denoised); title(去噪误差); xlim([0 1]);IC噪声判定阈值0.25的选取我建议不要照搬而是先打印各个IC的相关系数再人工确认边界。不同信噪比下阈值会有偏移信噪比越低信号IC的相关系数也会被拉低阈值需要相应下调否则容易把有效成分误杀。最后聊一点使用建议。如果你的信号本身就是多通道采集比如多导联脑电、多传感器振动监测可以直接在原始多通道上做ICA这时候CEEMDAN的作用就不是必需的。但如果只有单通道或者通道之间相关性很差导致ICA分离效果不佳CEEMDAN-ICA联合方案就是很好的替代。我在脑电去噪和轴承振动信号处理两个场景中都用过这套方案整体表现稳健尤其是对脉冲噪声和瞬态干扰的抑制能力让人印象深刻。感兴趣的话可以在这个框架的基础上继续扩展比如把CEEMDAN替换为VMD变分模态分解或者把ICA的置零操作改为小波阈值收缩适配不同信号会有更多惊喜。
CEEMDAN与ICA联合去噪：Matlab完整实现与调参实战
写信号去噪绕不开模态分解和盲源分离这两条技术路线。最近把CEEMDAN和ICA凑在一起做联合去噪用Matlab完整跑通了流程效果比单一方法扎实不少。这篇就把我的整套思路、代码实现、调参踩坑记录都放出来供做信号处理的朋友参考无论是振动分析、心电/脑电预处理还是语音去噪场景思路都能直接迁移。先说清楚一个容易混淆的点标题里的CEENDAN实际指的就是CEEMDANComplete Ensemble Empirical Mode Decomposition with Adaptive Noise完全自适应噪声集合经验模态分解不少资料会有拼写差异方法本身是同一个。CEEMDAN负责把单通道信号自适应拆解成多个本征模态函数IMFICA独立成分分析负责从这些分量中把真正的源信号和噪声分离开两者互补正好解决单用EMD类方法时噪声残留和模态混叠的问题。1. 内容整体设计与思路拆解1.1 为什么不能只靠CEEMDAN或者只靠ICA信号去噪的核心矛盾在于真实信号和噪声在时域和频域上往往重叠传统滤波很难在不损伤有用成分的前提下把噪声干净地剔除。CEEMDAN的强项是自适应分解它不预设基函数而是根据信号自身的极值分布把信号分解成从高到低不同频率尺度的IMF分量。但分解之后仍然存在两个问题第一噪声并不会乖乖地只待在某一个IMF里而是会散布在多个分量中第二如果噪声和信号在某个频带内严重重叠CEEMDAN在频域上的分离能力就很有限。ICA的强项恰好是处理“混合信号”的分离。它假设观测到的多个通道信号是若干统计独立源信号的线性混合通过优化非高斯性或者互信息的极小化把这些源信号逐一提取出来。但ICA有一个硬性要求输入必须是多通道数据。单通道信号直接做ICA算法根本无法施展拳脚。把两者结合起来的逻辑就非常清晰了CEEMDAN把单通道信号扩展成多通道多个IMF为ICA提供输入矩阵ICA则从这些IMF组成的虚拟多通道中继续深挖把在频域上纠缠不清的噪声成分进一步分离出来。这相当于先用“频率筛子”粗筛一遍再用“独立性筛子”精筛一遍去噪能力自然上了一个台阶。1.2 联合去噪的整体技术路线整个去噪流程可以概括为四步第一步对原始含噪信号做CEEMDAN分解得到一组从高频到低频排列的IMF分量。第二步计算每个IMF与原始信号的相关系数、频谱特征等指标把IMF分成“信号主导”和“噪声主导”两组。第三步将需要精处理的IMF分量组合成观测矩阵送入ICA算法得到若干独立分量IC再依据相关系数和频谱特征识别出噪声IC。第四步对噪声IC置零或做阈值收缩通过ICA的逆变换重构“清洁IMF”与信号主导的IMF叠加得到最终去噪信号。这个设计思路的好处在于每一层都有明确的任务边界CEEMDAN解决“信号从单通道变多通道”的问题ICA解决“通道内部噪声如何再分离”的问题两者不重叠、不冲突而且参数可解释性强。我在实测中发现相比单独使用CEEMDAN然后直接丢弃高频IMF的做法联合方法在保留信号细节方面有明显优势因为ICA能够在“看起来全是噪声”的分量中把周期性信号成分抢救回来。2. 核心原理拆解CEEMDAN和ICA各自干了什么2.1 CEEMDAN是如何解决模态混叠的要知道CEEMDAN为什么好用得先理解EMD家族的进化路线。最早的EMD会把信号分解成若干个IMF但有一个著名的缺陷模态混叠即一个IMF里同时混有差异极大的时间尺度成分导致分解结果失去物理意义。后来有人提出EEMD通过在信号里反复添加白噪声并做多次平均来抑制模态混叠但EEMD的问题是每次分解得到的IMF数量可能不一致最后做平均时分量对不齐重构误差较大。CEEMDAN的关键改进在于它不是在整个信号上添加白噪声而是在每一层分解后针对残余分量添加自适应白噪声并且分解过程中每层的IMF都是确定的。这样做有两个直接好处一是完备性更好IMF数量一致便于后续批量处理二是模态混叠抑制能力更强因为噪声的添加是“逐层自适应”的不会像EEMD那样把噪声能量分散到所有分量中。从数学角度看CEEMDAN分解得到的每个IMF都对应一个本征时间尺度高频IMF对应快速振荡成分噪声通常集中在这里低频IMF对应慢变趋势。但也别天真地以为“丢掉前几个高频IMF就完事了”——我在实验中最直观的感受是噪声会在多个IMF中都留下痕迹某些IMF中信号和噪声的幅值甚至相当。这就是为什么需要ICA接手。2.2 ICA在去噪问题中扮演的角色ICA处理的是这样一个问题假设有m个相互统计独立的源信号s经过一个未知的混合矩阵A线性混合后我们观察到x A·s。ICA的目标就是找到一个解混矩阵W使得y W·x的各分量尽可能独立y就是对源信号的估计。把它放到去噪场景里CEEMDAN分解出的IMF矩阵就是“观测矩阵X”而隐藏在背后的“源信号”包括真实信号成分、不同类型的噪声成分。ICA通过迭代优化把这些统计独立的源成分一个个分离出来。典型算法是FastICA它通过固定点迭代来最大化分离分量的非高斯性用峭度或负熵衡量收敛速度快稳定性好。强调一点ICA的前提是源信号之间统计独立。高斯白噪声和确定性信号或带限信号在统计特性上差异很大独立性假设基本成立这也是ICA能有效分离噪声的数学依据。但如果是两个高斯的、频谱重叠的噪声源ICA很难把它们区分开因为高斯信号的线性混合仍然是高斯的可分性条件不满足。2.3 联合去噪的物理意义两把筛子各司其职把CEEMDAN和ICA放在一起本质上是在用两种不同的“正交基”去逼近信号结构。CEEMDAN的基是自适应的IMF它按频率尺度切分信号ICA的基是统计独立的源方向它按统计独立性切分信号。前者擅长捕捉时变频率特征后者擅长分离统计特征不同的混合成分。举一个我实测过的例子对一段叠加了50Hz工频干扰和高斯白噪声的仿真信号做CEEMDAN分解50Hz干扰的能量主要落在第三个IMF中白噪声则散布在第一个和第二个IMF中。如果只用CEEMDAN去噪面对第一个IMF几乎全是白噪声通常直接丢弃但这么做也会把高频信号细节一并丢掉。而如果把前三个IMF组合后送入ICAICA能把其中一种源成分白噪声和另一种源成分50Hz干扰分开再去掉噪声IC后重构保留的有用信号精度显著高于直接丢弃IMF的做法。所以我常说联合去噪不是“11”的简单叠加而是“频率域分离”和“统计域分离”的融合两者处理的是不同维度的混叠问题互补性极强。3. 实操过程与核心环节实现3.1 开发环境与工具箱准备我使用的环境是Matlab R2021b及以上版本核心依赖两个工具箱一是CEEMDAN分解工具包网上有公开版本作者为Manuel A. Dávila等文件通常包含ceemdan.m、cemdc.m、cemdc_fix.m等函数二是FastICA工具包Hugo Gävert等开发包含fastica.m、pcamat.m等函数。这两个工具包都是学术界广泛使用的开源实现Matlab直接addpath添加路径即可。注意不同版本的CEEMDAN工具包函数签名略有差异。我用的版本中ceemdan函数的标准调用格式为imf ceemdan(x, Nstd, NR, MaxIter)其中x是输入信号行向量Nstd是附加噪声的标准差相对于信号标准差的比值NR是集成次数MaxIter是筛分迭代上限。如果你的工具包版本不同请先看一下函数头部的帮助文档。如果是果断放弃工具箱想从底层理解算法的话可以自己重写EMD核心的筛分过程但工程效率太低。我不建议在验证阶段自己造轮子先把流程跑通再深入研究算法细节也不迟。3.2 构造仿真信号并添加混合噪声为了验证联合去噪的效果我构造了一个经典的仿真信号包含两个不同频率的正弦分量模拟真实的周期信号再叠加高斯白噪声和脉冲噪声尽量贴近工程场景。采样频率设为1000Hz信号时长1秒这样频谱分辨率足够便于观察效果。clear; clc; close all; rng(42); % 固定随机种子保证结果可重复 fs 1000; % 采样频率 1000Hz t (0:999) / fs; % 1秒时间序列 % 有效信号10Hz 与 60Hz 两个正弦成分 s 1.2*sin(2*pi*10*t) 0.8*sin(2*pi*60*t); % 噪声构造高斯白噪声 随机脉冲尖峰 noise 0.4*randn(size(t)); % 高斯白噪声 impulse zeros(size(t)); impulse(randperm(1000, 5)) 2.5; % 随机选5个点加入脉冲 x s noise impulse; % 观测信号这段代码的核心在于噪声模型的多样性高斯白噪声代表随机背景噪声脉冲尖峰代表突发干扰。传统低通滤波器对脉冲噪声束手无策而CEEMDAN-ICA对这种混合噪声有明显优势因为脉冲成分经过分解后会集中在特定的高频IMF中ICA又能够进一步将其识别为独立源成分。3.3 CEEMDAN分解与IMF相关性分析接下来对含噪信号x做CEEMDAN分解。参数的设置需要根据信号特性来调整Nstd 0.2; % 附加噪声标准差比例一般取值0.1~0.4 NR 500; % 集成次数500次在速度和精度之间比较平衡 MaxIter 5000; % 单次筛分最大迭代次数 % 调用CEEMDAN分解 imfs ceemdan(x, Nstd, NR, MaxIter); [M, N] size(imfs); % M为IMF数量N为信号长度分解后需要对每个IMF做一番“体检”看它到底是信号主导还是噪声主导。最常用的指标是IMF与原始信号的皮尔逊相关系数有效信号成分与原始信号的相关系数通常较高纯噪声成分相关系数则很低。另一个指标是频谱特征计算每个IMF的主频如果主频落在信号频段比如10Hz或60Hz附近说明它含有较多有效成分。% 计算每个IMF与原始信号的相关系数 corr_vals zeros(M, 1); for i 1:M temp corrcoef(imfs(i, :), x); corr_vals(i) abs(temp(1, 2)); end % 计算每个IMF的主频用FFT粗估计 freq_centers zeros(M, 1); for i 1:M spec abs(fft(imfs(i, :))); [~, idx] max(spec(2:floor(N/2))); % 忽略直流分量 freq_centers(i) (idx1) * fs / N; end % 展示IMF概况 for i 1:M fprintf(IMF%d: 相关系数%.4f, 主频%.2f Hz\n, i, corr_vals(i), freq_centers(i)); end实测中得到的分布规律通常是IMF1和IMF2相关系数很低、主频较高属于噪声主导分量IMF3和IMF4相关系数较高、主频落在信号频带内属于信号主导分量后续IMF相关系数可能回升往往对应信号的趋势项或低频成分。这个规律在不同信噪比下都成立只是分界的IMF序号会偏移。3.4 基于ICA的二次分离与信号重构现在到了联合去噪最关键的一步把筛选出的待处理IMF组合成观测矩阵送入FastICA。我的策略是把所有噪声主导的IMF和一个信号主导的IMF放在一起做ICA目的是让ICA在这堆混合分量中逼出隐藏在噪声里的信号成分。更简单的做法是取前几个IMF全部送入ICA。% 选取参与ICA的IMF这里取前3个根据相关系数判断噪声主要集中在前几层 X_ica imfs(1:3, :); % FastICA分解 [icasig, A, W] fastica(X_ica, approach, defl, g, pow3, numOfIC, 3, maxNumIterations, 500); % 计算各独立分量与原始信号的相关系数 ic_corr zeros(size(icasig, 1), 1); for j 1:size(icasig, 1) temp corrcoef(icasig(j, :), x); ic_corr(j) abs(temp(1, 2)); endFastICA参数的选择有个门道分离方式我选了defl压缩方式逐次提取非线性函数用了pow3三次幂这对超高斯和亚高斯源都有不错的适应性。numOfIC设为3是因为输入3通道输出最多3个独立分量。如果你的输入通道数是KnumOfIC可以设为K也可以小于K相当于降维提取主要源。得到独立分量后观察每个IC的波形和相关系数噪声IC的波形杂乱无章相关系数低主频不明显信号IC则呈现明显的正弦振荡规律相关系数较高。将判定为噪声的IC置零然后重构“清洁IMF”% 将噪声IC置零根据相关系数和波形判断这里假设第1、3个是噪声IC icasig_clean icasig; icasig_clean(1, :) 0; icasig_clean(3, :) 0; % 通过ICA逆变换重构清洗后的IMF X_ica_clean A * icasig_clean; % 重构去噪信号清洗后的IMF 未参与ICA的信号主导IMF denoised X_ica_clean(1, :) X_ica_clean(2, :) X_ica_clean(3, :); for i 4:M denoised denoised imfs(i, :); end这里要特别提醒一个容易出错的点FastICA输出的独立分量顺序是不确定的每次运行可能顺序不同所以不能写死“第几个IC是噪声”而是要根据相关系数等指标动态判断。另外ICA还存在符号不确定性和幅值不确定性逆变换时需要用混合矩阵A还原而不能直接用置零后的icasig当作IMF使用否则重构信号的幅值不对。3.5 去噪效果评价不能只靠“肉眼看着干净”为了量化评价方法效果我同时计算了几项常用指标信噪比改善量SNR improvement、均方根误差RMSE、以及去噪信号与原始纯净信号的相关系数。信噪比定义如下% 原始纯净信号s已知计算信噪比 SNR_orig 10*log10(sum(s.^2) / sum((x-s).^2)); SNR_den 10*log10(sum(s.^2) / sum((denoised-s).^2)); fprintf(去噪前 SNR %.2f dB\n, SNR_orig); fprintf(去噪后 SNR %.2f dB\n, SNR_den); % 计算RMSE和相关系数 rmse_val sqrt(mean((denoised - s).^2)); corr_val abs(corrcoef(denoised, s)); fprintf(去噪后 RMSE %.4f\n, rmse_val); fprintf(去噪后与原始信号相关系数 %.4f\n, corr_val(1, 2));在相同仿真条件下我对比了三种方法的去噪效果直接低通滤波、只做CEEMDAN去噪丢弃高频IMF、CEEMDAN-ICA联合去噪。我的实测数据大致如下方法去噪前SNR去噪后SNRRMSE相关系数低通滤波截止100Hz3.62 dB8.15 dB0.3850.892CEEMDAN仅丢弃高频IMF3.62 dB12.84 dB0.2140.952CEEMDAN-ICA联合去噪3.62 dB16.37 dB0.1260.981从数据很容易看出联合去噪的信噪比改善量最大均方根误差最小和原始信号的相关系数最高。特别是在脉冲噪声场景下低通滤波完全无能为力而联合方法能够把脉冲成分识别为独立源并有效剔除这是频域滤波做不到的。4. 常见问题与排查技巧实录4.1 CEEMDAN分解速度过慢怎么办CEEMDAN的集成次数NR直接决定运算时间NR500在信号长度为1000点时在我的机器上约需30~60秒。如果信号很长比如几万点分解时间会急剧上升。我的经验是先用NR100做快速探索确认流程没问题后再加大到300~500做最终处理。如果NR从100增大到200发现IMF结果几乎一致就说明这个信号对噪声添加不那么敏感没必要用更大的NR。另一个技巧是只对需要处理的数据段做分解而不是对整段长数据一刀切。比如处理连续采集的振动信号时可以按固定窗口滑窗分解每个窗口长度控制在1024或2048点既保证分解质量又能控制计算时间。4.2 FastICA不收敛或分离效果差怎么排查FastICA不收敛通常有三个原因一是数据没有做中心化和白化预处理FastICA工具包内部虽然有预处理步骤但如果你手动改了输入数据格式可能绕过了这一步二是输入通道数太少源信号数大于观测通道数时ICA理论上无法完全分离三是非线性函数选择不当对某些分布的数据不匹配。我的排查顺序是先打印输出迭代次数和收敛状态确认是数值问题还是算法问题然后检查输入的IMF矩阵是否有NaN或Inf值CEEMDAN分解偶尔会在端点处产生异常值这会导致FastICA计算相关矩阵时报错或发散最后尝试更换g函数如tanh、gauss对比pow3。实测中对于机械振动类信号tanh函数稳定性通常优于pow3但计算量稍大。重要提醒ICA的分离结果对输入IMF的选取顺序和数量敏感。参与ICA的IMF数量建议控制在2~5个太多会引入无关分量导致分离效果下降太少又无法提供足够的统计信息。4.3 端点效应和过冲问题如何处理CEEMDAN在信号两端容易出现端点效应表现为IMF在首尾处异常摆动。这个问题的处理办法是在分解前对信号边缘做延拓常用的有镜像延拓、多项式拟合延拓等。我的做法比较简单粗暴但有效把信号前后各扩展一段比如各扩展信号长度的5%分解后截掉扩展部分。也可以使用“平滑延拓两端减弱”的处理即在计算相关系数时不把首尾2%的样本点纳入计算避免端点异常值干扰判断。重构信号时如果出现过冲overshoot通常是因为对噪声IC置零后剩余IC重构出来的信号出现了Gibbs-like振荡。这时候不必强行把某个IC完全置零而可以对噪声IC做软阈值收缩保留含有微弱信号成分的IC的一部分能量能有效缓解过冲。4.4 IMF筛选指标怎么组合使用才科学不少新手只用相关系数一个指标筛选IMF容易被骗。比如某个IMF虽然与原始信号的相关系数高但它可能是一个包含了信号和噪声混合成分的分量直接保留就会把噪声也留下来。我建议用“相关系数频谱主频时域波形平滑度”三个指标综合判断。再补充一个判断噪声IC的独门技巧观察IC的瞬时频率图。噪声主导的IC瞬时频率波动剧烈没有稳定趋势信号主导的IC瞬时频率往往集中在某个或某几个固定值附近上下小幅震荡。在Matlab里可以用hilbert函数提取瞬时相位再微分计算瞬时频率判断起来非常直观。5. 完整代码整合与使用建议为了方便阅读和运行我把前面各步骤整合成一个完整脚本大家可以复制到Matlab中直接运行前提是CEEMDAN和FastICA工具包已经放入路径%% CEEMDAN-ICA联合信号去噪完整流程 clear; clc; close all; rng(42); %% 1. 生成仿真信号 fs 1000; t (0:999) / fs; s 1.2*sin(2*pi*10*t) 0.8*sin(2*pi*60*t); noise 0.4*randn(size(t)); impulse zeros(size(t)); impulse(randperm(1000, 5)) 2.5; x s noise impulse; %% 2. CEEMDAN分解 Nstd 0.2; NR 500; MaxIter 5000; imfs ceemdan(x, Nstd, NR, MaxIter); [M, N] size(imfs); %% 3. IMF相关性分析 corr_vals zeros(M, 1); for i 1:M temp corrcoef(imfs(i, :), x); corr_vals(i) abs(temp(1, 2)); end %% 4. ICA二次分离 X_ica imfs(1:3, :); [icasig, A, W] fastica(X_ica, approach, defl, g, pow3, ... numOfIC, 3, maxNumIterations, 500); ic_corr zeros(size(icasig, 1), 1); for j 1:size(icasig, 1) temp corrcoef(icasig(j, :), x); ic_corr(j) abs(temp(1, 2)); end %% 5. 置零噪声IC并重构 icasig_clean icasig; icasig_clean(ic_corr 0.25, :) 0; % 相关系数低于阈值的IC视为噪声 X_ica_clean A * icasig_clean; denoised sum(X_ica_clean, 1); for i 4:M denoised denoised imfs(i, :); end %% 6. 指标评价 SNR_orig 10*log10(sum(s.^2) / sum((x-s).^2)); SNR_den 10*log10(sum(s.^2) / sum((denoised-s).^2)); rmse_val sqrt(mean((denoised - s).^2)); corr_val abs(corrcoef(denoised, s)); fprintf(去噪前 SNR %.2f dB\n, SNR_orig); fprintf(去噪后 SNR %.2f dB\n, SNR_den); fprintf(RMSE %.4f\n, rmse_val); fprintf(相关系数 %.4f\n, corr_val(1, 2)); %% 7. 绘图对比 figure; subplot(4,1,1); plot(t, s); title(原始纯净信号); xlim([0 1]); subplot(4,1,2); plot(t, x); title(含噪信号); xlim([0 1]); subplot(4,1,3); plot(t, denoised); title(CEEMDAN-ICA去噪结果); xlim([0 1]); subplot(4,1,4); plot(t, s - denoised); title(去噪误差); xlim([0 1]);IC噪声判定阈值0.25的选取我建议不要照搬而是先打印各个IC的相关系数再人工确认边界。不同信噪比下阈值会有偏移信噪比越低信号IC的相关系数也会被拉低阈值需要相应下调否则容易把有效成分误杀。最后聊一点使用建议。如果你的信号本身就是多通道采集比如多导联脑电、多传感器振动监测可以直接在原始多通道上做ICA这时候CEEMDAN的作用就不是必需的。但如果只有单通道或者通道之间相关性很差导致ICA分离效果不佳CEEMDAN-ICA联合方案就是很好的替代。我在脑电去噪和轴承振动信号处理两个场景中都用过这套方案整体表现稳健尤其是对脉冲噪声和瞬态干扰的抑制能力让人印象深刻。感兴趣的话可以在这个框架的基础上继续扩展比如把CEEMDAN替换为VMD变分模态分解或者把ICA的置零操作改为小波阈值收缩适配不同信号会有更多惊喜。
