简介：面向数字信号处理研究者与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.6 | 0.2 | 过小分解不充分，过大会产生伪模态 |
| MaxIter | 每个IMF最大筛选迭代次数 | 50 ~ 500 | 200 | 过小未分解完整，过大会过分解 |
| corrThr | 相关系数重构阈值 | 0.05 ~ 0.6 | 0.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.m：
function [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，配合批量去噪脚本即可形成完整工具链。
本文还有配套的精品资源，点击获取
