简介面向数字信号处理研究者与Matlab学习者这款基于龙格库塔优化算法RUN改进ICEEMDAN的去噪源码包针对含噪信号提供从分解、重构到性能评估的完整Matlab实现适合处理非平稳、非线性信号。压缩包共含34个文件包括16个m源码、14个csv测试数据、3张结果效果图与1个txt说明文档整体仅124KB轻量便携。其中m文件覆盖主程序、ICEEMDAN与CEEMDAN分解函数、RUN优化模块以及MSE、PSNR、信噪比等评价指标csv数据包含不同信噪比样本png图片可直接查看去噪效果。目前已有61人学习零基础用户只需替换数据即可运行代码经过实测能直观对比去噪前后波形帮助深入理解RUN-ICEEMDAN的参数寻优与信号重构机制是数字信号处理课程设计或算法对比实验的有力工具。1. RUN-ICEEMDAN信号去噪把参数搜索交给龙格库塔优化器数字信号去噪不能只依赖低通滤波。对于振动、生物电、水声这类非平稳信号ICEEMDAN能把信号自适应地分解成若干固有模态函数再从噪声主导的分量中把有效成分挑出来。但ICEEMDAN本身有两个麻烦一是噪声幅值、最大筛选次数要人定二是不同信号之间参数完全不通用。RUN-ICEEMDAN的思路是用龙格库塔优化器RUN去搜索这些参数目标函数定义为重构信号与带噪信号的误差或噪声估计让算法自己找到当前信号最优的分解参数组合。这样得到的去噪结果在信噪比和波形保真度上通常优于固定参数也避免了大量手动试参。这篇内容适合处理实测信号、又希望去噪流程自动化的工程师和科研人员原理、代码、参数边界和坑都会讲透。2. RUN-ICEEMDAN去噪的原理分解、优化与目标函数2.1 ICEEMDAN的分解逻辑与去噪切入点先回顾一下为什么需要ICEEMDAN。传统EMD能把非平稳信号按局部极值包络分解成若干IMF但存在模态混叠即一个IMF里混入频率差异很大的成分。EEMD通过加入白噪声平均来改善代价是分解结果残留噪声CEEMDAN把噪声添加从原始信号挪到每次筛余后的残差上降低了早期IMF的噪声残留但部分阶段仍会生成伪模态。ICEEMDAN在CEEMDAN基础上改用特殊的自适应噪声并在筛余过程中估计局部均值而不是直接使用含噪信号残差因此分解余量更干净伪模态更少。对去噪任务来说ICEEMDAN能按频率尺度把噪声和信号拆开。越靠前的IMF频率越高越靠后的IMF频率越低高频噪声通常集中在IMF1到IMF3但具体落在哪个IMF取决于采样率和信号内容。直接丢掉前几个IMF的做法很危险因为有用信号的高频成分也可能出现在这些位置。工程常见做法是计算每个IMF与原始信号的相关系数、能量密度、排列熵等统计量再根据阈值决定保留哪些IMF。其中相关系数阈值连续可调正好适合作为优化变量。2.2 RUN优化器如何实现参数寻优RUN的全称是Runge Kutta optimizer灵感来自四阶龙格库塔方法。它把当前解、历史最优解和全局最优解的差异当作“斜率”用RK积分方式更新位置。和粒子群、差分进化相比RUN的更新方向带有数值积分的平滑效果在低维连续参数空间上收敛快不容易在后期震荡。位置更新的核心可以概括成x_new x_i (1/6)(k1 2k2 2k3 k4)其中k1到k4来自不同解之间的线性组合。在MATLAB实现里每个粒子保存当前解每次迭代评估适应度再按RK系数更新位置。比起标准PSO的“速度位置”RUN多了一步类似积分的中值处理因此它对边界附近粒子的拉回也更平缓。RUN适合优化连续参数。ICEEMDAN参数中噪声相对幅值Nstd、最大筛选迭代数MaxIter、集合次数NR、相关系数阈值corrThreshold其中NR通常为整数且对结果影响不大可以固定为100到200。我们需要让RUN处理的是Nstd、MaxIter、corrThreshold三个连续或近似连续变量MaxIter虽然是整数但可以在目标函数里round取整不影响优化器主循环。2.3 适应度函数直接最小化重构误差不理想很多初学RUN-ICEEMDAN的人会把“去噪后的信号与原始含噪信号之差的均方误差MSE”作为目标函数让RUN去最小化它。这等于鼓励优化器把所有IMF全部保留因为全保留时重构信号等于原始带噪信号误差为零。显然这不是我们要的去噪效果。更合理的做法是估计去噪后的信噪比SNR让RUN最大化SNR。没有干净信号时可以用IMF1的噪声主导特性做鲁棒估计。常见公式是noiseStd median(abs(IMF1)) / 0.6745; snrEst 10 * log10( var(x) / noiseStd^2 );其中0.6745来自标准正态分布的中位数绝对偏差与标准差的关系。IMF1通常包含最多的高频噪声用这个方式估计噪声标准差在语音增强和机械振动里都有应用。估计出的SNR不精确但足以比较不同参数组合的优劣。为了避免优化器为了噪声指标乱选IMF我一般会在目标函数里加入一个惩罚项比如fval -snrEst 5 * mean(selMask)selMask表示被选中的IMF比例。加入惩罚后优化器会更倾向于用尽量少的IMF获得较高的信噪比这符合去噪应当“去除噪声、不损失有效成分”的直觉。3. MATLAB环境下RUN-ICEEMDAN去噪的完整实现3.1 主流程脚本与调用关系在MATLAB中实现RUN-ICEEMDAN去噪通常需要三层结构最外层是主脚本负责读取信号、设置参数边界、调用优化器、重构和绘图中间层是RUN优化器最里层是ICEEMDAN分解函数。下面的主脚本在MATLAB R2019b及以上版本可运行前提是你已经有一个可用的iceemdan_mine函数这个函数可以从你拿到的源码包中提取也可以自行实现。% demo_run_iceemdan.m clc; clear; close all; rng(2025); fs 1000; t (0:1/fs:1); signal 2*sin(2*pi*5*t) 1.2*sin(2*pi*50*t) 0.8*sin(2*pi*120*t); noise randn(1, length(t)) * 0.4; x signal noise; lb [0.1, 50, 0.05]; % Nstd下限, MaxIter下限, corrTh下限 ub [0.6, 500, 0.6]; % Nstd上限, MaxIter上限, corrTh上限 options.NP 20; % 种群数量 options.MaxIter 30; % 优化迭代次数 [bestP, bestF] run_optimizer_mine(x, lb, ub, options); fprintf(最优参数: Nstd%.3f, MaxIter%d, corrThr%.3f\n, ... bestP(1), round(bestP(2)), bestP(3)); % 用最优参数做最终分解 [imf, R] iceemdan_mine(x, bestP(1), 100, round(bestP(2)), 2); % 计算每个IMF与原始信号的相关系数 corrs zeros(size(imf, 1), 1); for k 1:size(imf, 1) cc corrcoef(imf(k, :), x); corrs(k) cc(1, 2); end % 相关系数大于阈值则保留重构时始终加上趋势项R mask corrs bestP(3); y_rec R sum(imf(mask, :), 1); figure; subplot(3, 1, 1); plot(t, x); ylim([-4 4]); title(含噪信号); subplot(3, 1, 2); plot(t, y_rec); ylim([-4 4]); title(RUN-ICEEMDAN去噪结果); subplot(3, 1, 3); plot(t, signal - y_rec); title(去噪误差);脚本逻辑并不复杂。第一步生成由三个正弦叠加的模拟信号加高斯白噪声。第二步设置参数边界lb和ub分别对应噪声幅值、最大筛选迭代数和相关系数阈值。第三步调用run_optimizer_mine返回最优参数。第四步用最优参数调用iceemdan_mine完成分解。第五步基于相关系数阈值选择IMF并重构。这里的关键点是相关系数阈值越高保留的IMF越少去噪越激进阈值越低保留的IMF越多波形更完整。3.2 一个可用的RUN优化器实现下面是run_optimizer_mine的完整实现。为了便于理解我保留RUN算法的核心RK更新去掉了论文里的一些随机缩放因子并没有照搬完整版本而是提供一个可复现的精简形式。重点是展示“如何用RK增量更新候选解”。function [best_x, best_f] run_optimizer_mine(sig, lb, ub, opt) % RUN优化器精简实现 % 输入: % sig : 待去噪信号 % lb : 参数下边界1*dim % ub : 参数上边界1*dim % opt : 结构体包含NP和MaxIter % 输出: % best_x : 最优参数 % best_f : 最优适应度 NP opt.NP; T opt.MaxIter; dim length(lb); % 初始化种群 X repmat(lb, NP, 1) rand(NP, dim) .* repmat(ub - lb, NP, 1); f zeros(NP, 1); for i 1:NP f(i) obj_iceemdan(X(i, :), sig); end [best_f, idx] min(f); best_x X(idx, :); Xbest best_x; for t 1:T for i 1:NP % 四个RK斜率分量a向量保证每个维度有不同的缩放 a rand(1, dim); xm (Xbest X(i, :)) / 2; k1 a .* (Xbest - X(i, :)); k2 a .* (xm - X(i, :)) rand(1, dim) .* (Xbest - xm); k3 a .* (xm - X(i, :)); k4 a .* (Xbest - X(i, :)) rand(1, dim) .* (X(i, :) - Xbest); Xnew Xbest (1 / 6) * (k1 2 * k2 2 * k3 k4); % 限制在边界内 Xnew max(min(Xnew, ub), lb); fnew obj_iceemdan(Xnew, sig); if fnew f(i) X(i, :) Xnew; f(i) fnew; end if fnew best_f best_f fnew; best_x Xnew; Xbest Xnew; end end end end这段代码里的斜率组合参考了四阶RK的加权思路。实际RUN论文的更新公式还会引入局部候选解和随机缩放因子这里为了确保可读性做了简化。每轮迭代里所有粒子都用同一个全局最优Xbest参与计算粒子之间没有直接通信靠历史最优和全局最优驱动收敛。边界修正使用max(min(Xnew, ub), lb)防止Nstd出现负数。目标函数obj_iceemdan在3.3中单独说明。它负责调用ICEEMDAN分解计算相关系数和噪声估计。3.3 目标函数与ICEEMDAN封装下面这段函数是本方案的核心目标函数。它不直接计算重构误差而是估计信噪比并加入选择惩罚避免优化器全保留IMF。function fval obj_iceemdan(params, sig) % RUN优化的目标函数 % params: [Nstd, MaxIter, corrThr] Nstd params(1); maxIter round(params(2)); corrThr params(3); % 调用ICEEMDAN分解100为集合次数2为噪声类型 [imf, ~] iceemdan_mine(sig, Nstd, 100, maxIter, 2); % 为空或分解失败时返回极大值 if isempty(imf) fval 1e10; return; end nImf size(imf, 1); corrs zeros(1, nImf); for i 1:nImf c corrcoef(imf(i, :), sig); corrs(i) c(1, 2); end % 选择相关系数大于阈值的IMF selMask corrs corrThr; % 用IMF1估计噪声标准差计算信噪比估计值 noiseStd median(abs(imf(1, :))) / 0.6745; powerTotal var(sig); snrEst 10 * log10(powerTotal / noiseStd^2 eps); % 加入选择比例惩罚让优化器避免全保留 penalty 5 * mean(selMask); % 最大化SNR等价于最小化负SNR fval -snrEst penalty; end需要说明的是这里调用iceemdan_mine函数输入参数依次是原始信号、噪声幅值、集合次数、最大筛选迭代数和噪声类型。集合次数固定为100在离线处理中足够了。噪声类型2表示使用ICEEMDAN标准噪声模式一般不要改。用IMF1估计噪声标准差有一个例外如果原信号本身信噪比很高IMF1可能不全是噪声而包含高频信号这样估计出的噪声功率偏大导致SNR估计偏低。但在RUN优化过程中所有参数组合共享同一个估计标准排序仍然有效。想要更精确可以使用IMF1和IMF2的噪声功率中位数作为噪声估计。主脚本、优化器和目标函数组合起来就是一个可运行的RUN-ICEEMDAN去噪流程。下面一个章节讲参数怎么调、怎么评判、以及最容易出现的坑。4. RUN-ICEEMDAN参数设置、评价指标与常见坑4.1 三个关键参数的影响范围RUN-ICEEMDAN实际需要调的核心参数量不大。多数情况下优化器搜索的是三项Nstd、MaxIter和corrThr。下表给出我在不同信号测试中常用的范围和经验初始值供参考。参数物理含义推荐搜索范围初始值影响Nstd添加的噪声相对幅值0.1 ~ 0.60.2过小分解不充分过大会产生伪模态MaxIter每个IMF最大筛选迭代次数50 ~ 500200过小未分解完整过大会过分解corrThr相关系数重构阈值0.05 ~ 0.60.3阈值越高保留IMF越少去噪越强但容易损失信号Nstd直接影响分解质量。ICEEMDAN需要在信号中加入有限幅值噪声来激发极值点Nstd太小无法改变极值分布分解结果和普通EMD接近Nstd太大则注入过多噪声低频IMF也会被噪声污染。从实际处理看振动信号取0.2到0.3生物电信号取0.1到0.2强噪声情况下可以到0.4以上。MaxIter影响分解完善程度。ICEEMDAN每个IMF的提取都通过迭代筛选完成MaxIter太小会导致IMF未完全分离相邻模态混在一起太大则筛选过度IMF趋向纯正弦成分破坏信号包络。一般50到200够用如果信号本身很长可以适当提高到300以上。corrThr决定重构阈值。这个参数在优化变量里特别重要因为它直接控制哪些IMF进入重构。实验中发现阈值在0.2到0.4之间变化对去噪结果影响不大但超过0.5后波形会明显变瘦低于0.1时噪声基本没有被去除。因此我会把搜索范围设置成0.05到0.6让优化器在两端之外找不到更优解。4.2 评价指标怎么算去噪效果需要用多个指标衡量不能只看时域波形。MATLAB中常见的三个指标如下。% 假设clean为干净参考信号denoised为去噪后的信号 SNR 10 * log10( sum(clean.^2) / sum((clean - denoised).^2) ); RMSE sqrt(mean((clean - denoised).^2)); CC corrcoef(clean, denoised); CC CC(1, 2);信噪比SNR反映整体能量误差RMSE反映逐点误差相关系数CC反映波形形态相似度。在只有含噪信号没有干净参考时可以用4.1的噪声估计法算SNR估计值。源码包里一般有比较函数通常是cal_SNR、cal_RMSE、cal_PRD三个文件具体名称要看压缩包内的命名。4.3 容易踩的几个坑坑一优化目标选成全保留IMF。前面已经提到直接用重构信号与带噪信号的MSE做目标会导致优化器选择全部IMF去噪变成无操作。解决方法是使用噪声估计信噪比加稀疏惩罚或者在目标函数中强制加入“保留IMF比例不超过60%”的硬约束。坑二边界效应。ICEEMDAN分解时信号两端容易受包络拟合影响产生端点发散。去噪后在首尾能看到明显波动。常见做法是在优化和最终分解之前先对信号做对称延拓分解完再裁剪。MATLAB中可用padarray实现x_ext padarray(x, [0, length(x)/10], symmetric); % 分解得到imf_ext后截取中间原始长度 imf imf_ext(:, length(x)/10 1 : length(x) length(x)/10);坑三趋势项误丢弃。有些低频趋势项不是噪声而是信号本身的一部分。重构时如果只把相关系数高的IMF相加趋势项R就会被丢掉导致信号均值漂移。正确做法是在重构时始终加上残差R就像主脚本中的y_rec R sum(imf(mask,:), 1)那样。坑四优化迭代中分解函数报错。RUN每次迭代都要调用ICEEMDAN不同参数组合可能导致分解失败。在obj_iceemdan里务必要加isempty判断失败时返回一个很大的fval避免优化器把无效解当成最优解。5. 进阶把RUN-ICEEMDAN封装成可复用去噪函数实际项目里不会只处理一条信号。把优化流程封装成标准函数输入一段信号输出去噪结果和最优参数能省下不少重复劳动。下面这段函数头设计可以保存为run_iceemdan_denoise.mfunction [denoised, bestP, info] run_iceemdan_denoise(x, fs, varargin) % 通用RUN-ICEEMDAN去噪函数 % 输入 % x : 待去噪信号行向量 % fs : 采样率用于后续可视化暂未使用 % varargin : NP、MaxIter、lb、ub % 输出 % denoised : 去噪后的信号 % bestP : 优化得到的最优参数 % info : 结构体包含SNR估计、选取IMF序号、迭代曲线 p inputParser; addParameter(p, NP, 20); addParameter(p, MaxIter, 30); addParameter(p, lb, [0.1, 50, 0.05]); addParameter(p, ub, [0.6, 500, 0.6]); parse(p, varargin{:}); opt.NP p.Results.NP; opt.MaxIter p.Results.MaxIter; [bestP, ~] run_optimizer_mine(x, p.Results.lb, p.Results.ub, opt); [imf, R] iceemdan_mine(x, bestP(1), 100, round(bestP(2)), 2); corrs zeros(size(imf, 1), 1); for k 1:size(imf, 1) c corrcoef(imf(k, :), x); corrs(k) c(1, 2); end mask corrs bestP(3); denoised R sum(imf(mask, :), 1); info.bestP bestP; info.selIMF find(mask); info.snrEst 10 * log10(var(x) / (median(abs(imf(1, :))) / 0.6745)^2); end调用方式很直接y run_iceemdan_denoise(x, fs, NP, 30, MaxIter, 50);封装好之后可以用pwelch对比去噪前后的频谱确认高频噪声是否被抑制同时观察有效频段是否保留。另一个快速的验证技巧是计算残差信号的短时自相关去噪效果好的话残差应该接近白噪声。在线场景里RUN优化30次迭代仍然偏慢一个实用技巧是先在不同SNR条件下离线做几次优化把得到的参数保存为查找表。实际运行时直接按信号预估SNR查表省去在线优化过程。每隔一段时间再用RUN跑一次离线优化更新查找表。MATLAB中把优化记录写入表格文件可以用writetable读取用readtable配合批量去噪脚本即可形成完整工具链。本文还有配套的精品资源点击获取