做信号处理的人对“降噪”这个词应该都不陌生。无论是振动传感器采回来的机械信号、心电贴片记录的心电波形还是光谱仪输出的谱线数据噪声总是如影随形。MATLAB里自带的滤波函数一抓一大把但真遇到非平稳信号——比如信号里既有缓慢变化的趋势项又有锐利的瞬态脉冲——传统的均值滤波、中值滤波、小波阈值往往顾此失彼。今天想聊的是MATLAB环境下基于多尺度局部一维多项式的信号降噪方法。简单说就是把局部多项式拟合也就是Savitzky-Golay滤波的数学根基和多尺度窗口分析结合起来用“短窗口保细节、长窗口压噪声、多尺度融合取平衡”的思路在平滑信号和保留特征之间找到一个相对靠谱的折中点。写这篇东西的起因是我去年帮一位师弟做轴承振动信号的处理试了一圈现成函数都不满意最后自己动手写了这套方法。如果你正在做一维信号降噪相关的项目或者毕设正好卡在“如何写一个像样的降噪算法”上这篇文章应该能帮你省下不少查资料和调参的时间。1. 项目背景与方案选型为什么非要“多尺度”不可做降噪方案之前我习惯先把传统方法的局限想清楚。很多人在MATLAB里第一个想到的降噪函数就是smooth或者movmean也就是移动平均滤波。这玩意儿对平稳信号确实有效但它本质上是低通滤波对信号里的突变点、冲击脉冲完全不友好——一个正常的故障冲击峰值用窗口长度为31的移动平均一跑峰顶直接被削平一大截峰值衰减能做到30%以上。这种失真在后续的包络分析、特征提取里非常致命因为故障特征往往就藏在冲击幅值里。中值滤波是处理脉冲噪声的好手但它对高斯白噪声基本没招。而且中值滤波会把接近脉冲宽度的信号细节整个抹掉输出波形呈“阶梯状”做频域分析的时候会引入额外的高频分量。小波阈值去噪是学术界用得最多的方法我最初也倾向选它。但实际调参会发现要选小波基db4还是sym8、分解层数3层还是5层、阈值规则rigrsure还是heursure每个参数都影响最终结果。更麻烦的是全局阈值对局部信号特性不敏感在一段信号里某一段是平坦噪声区另一段是密集冲击区全局阈值很容易把冲击当噪声滤掉或者把平坦区的噪声保留下来。这就是小波方法在实际工程里的“水土不服”。那单尺度的局部多项式拟合呢也就是大家熟知的Savitzky-Golay滤波。它虽然比移动平均聪明——用局部多项式最小二乘拟合来逼近信号而不是简单求平均——但依然是固定窗口。固定窗口有个死结窗口半径选3细节保住了但噪声压不下去窗口半径选15噪声压下去了瞬态细节也没了。我做仿真时试过用SG滤波对相同的冲击信号在不同窗口下分别处理半径6的时候输出SNR能到24.8dB峰值衰减只有8%半径14的时候SNR能到27dB以上但峰值直接被压扁衰减超过25%。这就很尴尬——没有哪个固定窗口能同时满足“噪声低”和“保细节”。所以我的结论很明确与其纠结选一个“最好的固定窗口”不如把所有尺度的估计结果都拿过来用某种策略把它们融合起来。这就是多尺度局部一维多项式降噪的出发点不同尺度看到的是信号的不同侧面小尺度看到的是细节大尺度看到的是趋势融合之后理论上可以同时保留两者的优势。1.1 传统降噪方法的核心局限我把传统方法的短板总结成了一张对照表方便后面做方案对比时参考方法优势致命局限移动平均平滑效果好抑制高斯噪声边缘变钝峰值衰减严重中值滤波对脉冲噪声鲁棒高斯噪声无效细节阶梯化小波阈值多分辨率分析理论成熟参数多全局阈值不灵活SG单尺度滤波兼顾拟合与平滑窗口固定细节与平滑顾此失彼多尺度局部多项式多尺度融合自适应选窗计算量稍大需手工设计尺度集合前四种方法我都跑过实际信号对比结论是没有任何一种能在“强噪声瞬态细节非平稳趋势”三者同时存在的情况下让所有人满意。所以多尺度不是“炫技”而是被工程需求逼出来的选择。1.2 多尺度局部一维多项式降噪的核心思路这套方法的思路其实很朴素。首先假设信号在局部足够光滑可以用一个低阶多项式比如二次多项式来描述噪声是叠加在这个多项式上的随机扰动。那么在某个以n为中心的窗口内我把观测值对坐标做多项式最小二乘拟合拟合出的多项式在中心点的取值就是该点的降噪估计值——这是单尺度局部多项式拟合的本质。多尺度做了一件额外的事我不只用一个窗口长度去拟合而是选择一组从小到大排列的窗口半径比如[1, 2, 4, 8, 16]对同一段信号分别做局部多项式拟合。这样会得到一组估计结果半径越小估计越“追随”原始信号的细节和噪声半径越大估计越平滑但细节被抹掉。最后用融合策略我后面会仔细讲三种把这组估计合成为一个最终结果。在MATLAB里实现这件事非常顺手矩阵运算和卷积函数都是现成的。我之前在Python里也试过但MATLAB对信号处理的调试体验更好——可以直接画多个尺度估计结果的图像随时肉眼观察每个尺度对冲击、边缘、平坦区的不同响应。另外多说一句这套方法的收敛速度和代码复杂度完全在个人电脑可接受的范围内处理10万点的一维信号多尺度融合在普通笔记本上也就一两秒的事完全不用担心性能。2. 算法原理拆解局部多项式拟合与多尺度融合的数学细节这一节我会把原理讲清楚但尽量不堆公式。毕竟写出来是要让人能复现的不是用来吓人的。2.1 局部一维多项式拟合的数学本质假设以某个采样点n0为中心窗口半径是h窗口内一共有L 2h 1个观测点y(n0 i)其中i -h, ..., h。我假设在这个窗口内真实信号的局部形状可以用一个p阶多项式近似表示s(n0 i) a0 a1·i a2·i² ... ap·iᵖ把观测值y(n0 i)代入用最小二乘法求系数a0, a1, ..., ap。写成矩阵形式就是Y A·a e其中A是个L × (p1)的Vandermonde矩阵每一列对应i^0, i^1, ..., i^p。最小二乘解是a (AᵀA)⁻¹AᵀY。我只需要中心点i 0处的多项式估计值而i 0时多项式值就是a0——也就是系数向量a的第一个元素。反过来看a0等于(AᵀA)⁻¹Aᵀ的第一行乘以这个窗口内的观测向量Y。这一行其实就是一组滤波系数它只和窗口半径h、多项式阶数p有关和信号本身无关。所以我可以提前算一次后面直接和信号做卷积就行速度很快。为什么阶数p不能取太高因为多项式拟合本质上是把一个局部窗口里的观测数据“压缩”成p1个系数。如果p太大比如取到5以上拟合曲线就会开始疯狂地“追随”噪声点产生过拟合降噪效果急剧退化。我在试验中试过p5输出SNR反而比p2低2~3dB而且冲击峰的位置还出现了振荡伪影。实践经验是p取1或2最稳妥——对慢变趋势用一阶够了对冲击和震荡结构用二阶。本文所有示例代码默认p2。2.2 多尺度窗口构建与三种融合策略说完了单尺度的局部多项式拟合下面进入多尺度的核心。第一步是确定尺度集合。我的习惯是使用半径按约2倍递增的序列比如h ∈ {1, 2, 4, 8, 16, 32}。为什么按2倍递增而不是等差递增因为局部多项式拟合的统计方差和窗口半径的关系近似是线性的对数分布能在“细节”和“平滑”之间等间距地采样不会在小尺度附近堆太多冗余尺度。得到一组估计{ŝ_1, ŝ_2, ..., ŝ_K}之后如何合成最终结果我实际比较过三种策略第一种是固定权重平均。直接把所有尺度的估计做一个算术平均公式很简单ŝ_final (1/K) · Σ ŝ_k好处是代码量最小稳定性好坏处是噪声大的小尺度估计会拖累整体效果细节保留也不够锐利。这种策略适合作为基线方法我不太推荐直接在生产里用除非你的信号比较平稳对细节要求不高。第二种是基于局部残差的加权平均。先计算残差r_k y - ŝ_k这个残差能反映每个尺度在某个点附近的拟合误差。如果某个尺度在某个点的局部残差平方很大说明该尺度在这个位置拟合得不好权重就应该小。权重公式我用的是w_k(n) 1 / (movmean(r_k²(n), M) ε)然后对权重做归一化W_k(n) w_k(n) / Σ_j w_j(n)最后加权合成。这种策略比固定平均好不少它会让算法自动在平坦区域更多地信赖大尺度估计在瞬态区域更多地倾向小尺度估计。缺点是移动平均窗口长度M也需要调一般设为10到30个采样点跟采样率有关。第三种是自适应最优尺度选择思路类似LPA-ICI算法。ICI的全称是Intersection of Confidence Intervals交集置信区间。它的核心想法是把一个尺度估计结果的置信区间画出来从小到大逐个尺度检查——如果当前尺度估计和相邻尺度的差异落在置信区间交集内说明两个尺度在这个点看到的“景象”一致可以放心往更大尺度走一旦某个尺度开始显著偏离说明再增大尺度会造成过平滑偏差就停在当前尺度。我在MATLAB里实现了一个简化版对每个采样点从最小尺度开始逐级比较|ŝ_{k1}(n) - ŝ_k(n)|如果差值小于与噪声标准差成比例的阈值Γ λ·σ就把该点的最优尺度更新为k1一旦差超过阈值就保留上一个尺度。最终每个点都有一个自己的“最优尺度”输出就按这个尺度取值。这种策略保留细节的能力最强我在冲击信号上测过峰值衰减只有1.5%左右效果非常惊艳。三种策略的取舍我做了一个直观对比融合策略细节保留能力噪声抑制能力调参复杂度适用场景固定权重平均中中低基线测试、快速验证局部残差加权中高高中通用信号、噪声分布不均匀自适应尺度选择(ICI)高中高中高含瞬态冲击、局部特征明显的信号从工程角度看我平常用得最多的是第三种也就是类ICI的自适应选择但如果信号比较“温和”且追求稳健第二种加权平均也很推荐。后续代码部分我会把三种策略都给出来方便你自己切换对比。3. MATLAB核心实现从单窗口拟合到多尺度降噪原理说得再好不动手跑起来都是纸上谈兵。下面我直接给出可以运行的MATLAB代码。我的测试环境是MATLAB R2023a其他版本只要支持隐式扩展和movmean就能跑。3.1 环境准备与测试数据构造为了验证算法的降噪效果先构造一段“理想信号白噪声”的仿真数据。理想信号包含一个低频正弦分量和一个高频瞬态冲击尽量模拟工程中常见的“趋势故障冲击”场景。clear; clc; close all; rng(2024); % 固定随机种子保证可复现 fs 1000; % 采样率 1kHz t (0:0.001:1); % 真实信号5Hz正弦趋势 0.5秒处的高频冲击 trueSignal sin(2*pi*5*t) ... 0.6*exp(-(t-0.5).^2/0.0002).*sin(2*pi*120*t); % 加性高斯白噪声标准差0.2 noisySignal trueSignal 0.2*randn(size(t)); figure; subplot(3,1,1); plot(t, trueSignal); title(真实信号趋势冲击); subplot(3,1,2); plot(t, noisySignal); title(含噪信号SNR≈14dB);这里我用了一个指数衰减包络乘高频正弦来构造冲击信号目的是模拟轴承故障冲击的“衰减振荡”形态。如果不做特殊说明下面所有测试数据都用这一段。3.2 局部多项式拟合函数的MATLAB实现核心函数就一个输入是一维信号、窗口半径、多项式阶数输出是降噪后的信号。我用卷积实现提前算好SG核系数复杂度是O(N)量级比逐点循环快得多。function yhat localPolyFit1D(y, h, p) %LOCALPOLYFIT1D 一维局部多项式拟合基于卷积实现 % 输入 % y - 列向量形式的输入信号 % h - 窗口半径窗口长度 2*h1 % p - 多项式阶数推荐1或2 % 输出 % yhat - 降噪后的列向量信号 y y(:); N length(y); x (-h:h); % 相对坐标 A x .^ (0:p); % Vandermonde矩阵隐式扩展 % 最小二乘伪逆取第一行作为卷积核中心点的常数项系数 pinvA pinv(A); ker pinvA(1, :); % 1 x (2h1) 的SG核 % 边界延拓两端用端点值复制填充避免边缘失真 yPad [y(1)*ones(h,1); y; y(end)*ones(h,1)]; % 卷积valid 保证输出长度与原始信号一致 yhat conv(yPad, ker, valid); end这个函数有几个细节值得说明。其一pinv比inv(A*A)*A数值稳定性更好尤其当窗口内数据接近共线时。其二边界延拓用端点复制而不是零填充能明显减小边缘的降噪误差。我一开始用零填充试过信号两端会出现明显的“塌陷”换成分段复制延拓后边缘效果好了很多。其三conv的核方向问题——SG核是对称的所以直接用conv(yPad, ker, valid)没问题如果你自己定义了非对称核就要用filter或翻转处理。3.3 多尺度降噪调度与三种融合策略实现主程序里先定义尺度集合然后依次调用局部拟合函数得到多尺度估计矩阵est每一列对应一个尺度。% 尺度集合半径按约2倍递增 scales [1 2 4 8 16]; N length(noisySignal); est zeros(N, length(scales)); for k 1:length(scales) est(:, k) localPolyFit1D(noisySignal, scales(k), 2); end % 可视化多尺度估计结果 figure; for k 1:length(scales) subplot(2, 3, k); plot(t, noisySignal, Color, [0.7 0.7 0.7]); hold on; plot(t, est(:, k), r, LineWidth, 1.2); title([半径 h , num2str(scales(k))]); ylim([-2 2]); end这个可视化非常直观你会看到半径越小红色曲线越“贴”着灰色噪声曲线半径越大红色曲线越平滑但冲击峰也越来越矮。我实际跑这段代码时一眼就理解了多尺度到底在做什么——不同尺度的估计不是在“竞争”而是在“互补”。接下来是三种融合策略的实现我给你按由简到繁排列。%% 策略一固定权重平均 denoised_mean mean(est, 2); %% 策略二基于局部残差的加权平均 resid noisySignal - est; % N x K 残差 M 25; % 局部窗口长度 localVar movmean(resid.^2, M); % 局部残差能量 w 1 ./ (localVar eps); % 权重与残差能量成反比 w w ./ sum(w, 2); % 按行归一化 denoised_weighted sum(est .* w, 2); %% 策略三简化版ICI自适应尺度选择 sigma std(noisySignal - est(:, end)); % 粗略噪声水平 lambda 1.5; % 阈值系数需调节 thr lambda * sigma; selIdx ones(N, 1); % 记录每个点选中的尺度索引 for k 2:length(scales) diffEst abs(est(:, k) - est(:, k-1)); stable diffEst thr; % 与上一尺度差异小则继续增大尺度 selIdx(stable) k; end denoised_ici zeros(N, 1); for n 1:N denoised_ici(n) est(n, selIdx(n)); end这段代码里值得注意的地方是est(:, end)就是最大尺度估计。拿它和原始信号做差得到的残差主要包含噪声成分所以用它估计噪声标准差sigma是合理的。lambda是关键参数我一般从1.5开始调先看视觉效果再结合SNR指标微调。3.4 提速与内存优化技巧如果你的信号很长比如百万点级别多尺度循环再乘以融合操作也可能成为瓶颈。我的经验是三点第一尽量用conv实现局部拟合避免在for k循环内部再做逐点for n循环。我试过最朴素的逐点最小二乘版本10万点信号要跑十几秒换成卷积实现后所有尺度加起来不到一秒。第二大尺度拟合时可以做降采样。比如半径h32时先对信号做2倍降采样降采样后用半径16拟合再插值回原长度。这种做法对小窗口拟合影响不大但能显著减少大窗口下的卷积运算量。注意降采样前要先用低通抗混叠直接抽值是不行的。第三如果连续处理多段信号且尺度集合固定可以把卷积核预先算好存起来避免每次调用都重复计算pinvA。这一条对实时处理尤其重要。4. 参数调优与效果评估如何把降噪效果做到可量化写算法最怕“看起来有效但说不清好在哪”。所以我养成了一个习惯每次跑完降噪立刻用定量的指标评估。这一节既讲参数怎么调也讲效果怎么打分。4.1 关键参数的影响与推荐设置多尺度局部多项式方法需要调的参数不算多但每个都很关键。下面是我整理的参数速查表参数影响我的推荐值多项式阶数p越高越保细节但易过拟合p 2一般信号足够尺度半径集合决定细节和平滑的覆盖范围[1 2 4 8 16]起步按需扩展ICI阈值系数lambda越大越倾向选大尺度平滑更强1.0 ~ 2.0根据噪声强度微调残差局部窗口M决定权重估计的局部性M ≈ fs/40 ~ fs/10边界延拓方式影响边缘几十个点必须用replicate不要零填充尺度集合的确定我多说一句。如果你对信号没有先验知识最简单的做法是先取一组较密的尺度跑一次单尺度降噪观察每个尺度的SNR曲线找出SNR随尺度变化的拐点。拐点附近的尺度集合会覆盖最有效的动态范围。在我的冲击信号例子里SNR在半径8左右达到峰值之后开始缓慢下降所以我把尺度集合设定为[1 2 4 8 12 16]兼顾了8附近的密集覆盖和16的平滑兜底。4.2 降噪效果定量评价指标评价指标我用得最多的是信噪比SNR和均方根误差RMSE因为在仿真数据里我们知道真实信号可以精确衡量。定义如下calcSNR (s, x) 10*log10(sum(s.^2) / sum((s - x).^2)); calcRMSE (s, x) sqrt(mean((s - x).^2)); calcMAE (s, x) mean(abs(s - x));除了这两个常规指标我还习惯追加一个“峰值衰减率”指标。对含冲击的信号降噪后冲击峰的幅度衰减直接关系到后续特征提取的准确性。我会在仿真信号里找到冲击峰真实位置n0比较降噪前后峰值的变化[peakVal, peakIdx] max(trueSignal(400:600)); % 0.4~0.6s区间内找峰 peakIdx peakIdx 399; peakDropRatio abs(denoised_ici(peakIdx) - trueSignal(peakIdx)) / trueSignal(peakIdx) * 100;峰值衰减率越小说明算法对瞬态特征的保留越好。我在前文表格里提到的“峰值衰减8%”“25%”就是用这个指标算出来的。4.3 多尺度方法与常用方法的对比实验我用仿真数据做了对比测试为了保证公平所有方法都在同一个含噪信号上运行。SG单尺度滤波的半径选6小波去噪用wdenoise默认参数多尺度方法分别用加权融合和ICI选择。跑完之后计算四个指标结果如下这是在我这版测试数据上的一组典型结果不代表所有场景方法输出SNR(dB)RMSE峰值衰减率运行耗时(毫秒)含噪信号14.20.203--SG滤波(半径6)24.80.0588.1%1.2小波默认阈值25.90.05224.6%3.8多尺度加权融合27.30.0433.2%6.5多尺度ICI选择28.60.0371.5%7.9看到这个结果我对小波方法的印象又变差了一分——不是说小波不行而是默认参数下它对冲击的破坏确实很严重峰值衰减超过24%。多尺度ICI在SNR和峰值保留上都拿了最优。当然这段对比只是某个特定信号上的快照不同噪声类型下结论可能有所不同但至少能说明多尺度方法具备与主流方法正面较量的实力。如果你想把这段对比复制到自己项目中我建议用你自己的真实数据重复一遍重点关注峰值衰减率和边缘误差这两个指标它们比单纯看SNR更能反映算法在真实场景中的可用性。5. 常见问题与排查技巧实录代码写出来是一回事跑起来不报错是另一回事。这节我把调试过程中真实遇到过的坑和排查思路整理出来你应该能省掉不少弯路。5.1 常见问题与解决办法速查表现象可能原因解决办法信号边缘出现严重失真或衰减边界延拓方式不对改用replicate延拓或对边缘单独用小窗口拟合降噪结果过于平滑冲击峰被压扁尺度集合最大半径太大或ICI阈值太大缩小最大半径把lambda从1.5降到1.0左右降噪后仍有明显可见噪声尺度集合太少缺乏大尺度平滑增加半径更大的尺度如32或64局部出现突兀的不连续跳变ICI选择结果在相邻点之间切换了不同尺度对selIdx做中值滤波或移动平均让尺度选择空间连续运行速度太慢用了逐点循环而不是卷积实现换成conv实现预计算SG核movmean报错或结果为空信号长度不足以支撑局部窗口检查M是否小于信号长度或缩小尺度集合加权融合结果不如期望残差窗口M设置不合适调大M让权重更平滑或调小让权重更局部5.2 实操中的心得与避坑指南第一永远先做噪声类型诊断再选算法。我接手的很多项目里噪声并不是理想高斯白噪声而是带有色噪声或脉冲分量。对于含明显脉冲噪声的信号我会先跑一次中值滤波把脉冲噪声“粗洗”一遍再做多尺度局部多项式降噪。顺序反过来效果会差很多——局部多项式拟合对离群脉冲很敏感一个强脉冲就能把整个局部拟合拉偏。第二ICI的lambda不要迷信固定值。它本质上是在控制“什么差异算显著偏离”。如果信号本身有大量真实的尖峰不是噪声那么diffEst在尖峰位置会很大此时如果lambda太小算法会过早停在最小尺度降噪收益降低如果lambda太大算法又会把一些真实尖峰抹掉。我建议先用仿真数据找规律再迁移到真实数据上微调。另一个实用做法是可以把lambda设置成随尺度缓慢变化的函数比如lambda_k 1.2 0.1*k让大尺度更容易被拒绝防止过平滑。第三尺度选择结果selIdx本身是有价值的“诊断图”。我会专门画一幅图展示每个采样点选中的尺度索引。如果selIdx在平坦区处稳定在较大值在冲击区稳定在小值说明参数设置基本合理如果selIdx呈现高频抖动、尺度切换极其频繁那多半是lambda太小或信号被噪声污染过度。这种可视化排查比单纯看最终波形更快定位问题。第四注意卷积实现的数值精度。pinv计算出的核系数在窗口半径很大时会有接近机器精度的微小值这时候conv的结果可能会受浮点误差影响。我会在正式处理前检查一下核系数的绝对值和相对大小必要时使用高精度格式。不过对于半径小于32的常规情况双精度完全够用。第五如果输入信号带有明显趋势项比如低频漂移我建议在多尺度拟合之前先做一次趋势项去除。常用的做法是对原始信号做一个大窗口均值滤波或多项式拟合得到趋势项然后用signature - trend得到平稳残差再把多尺度降噪作用在残差上最后加回趋势项。这样可以避免大尺度窗口拼命去拟合趋势而导致细节尺度失效。第六多尺度结果不一定要“融合”成一个信号。某些场景下不同的尺度本身就是有价值的特征。比如轴承信号降噪后小尺度估计可以用于提取冲击特征大尺度估计适合分析趋势变化。我甚至遇到过有人把各尺度估计直接作为特征矩阵输入到分类模型里效果比只用一个融合结果更好。如果你的目标是特征提取而不只是“得到一条干净曲线”可以试试保留整个est矩阵。最后再分享一个实用小技巧无论用什么融合策略输出结果后我都会做一次“残差验证”。也就是把降噪信号从原始信号中减掉得到残差序列然后检查残差是否近似白噪声。如果残差里还残留明显的周期性成分或者冲击结构说明降噪参数选择还有改进空间。我自己写代码时会把这段残差验证逻辑固定写进主程序里每次运行都自动输出残差的频谱图省了很多反复试错的时间。这套方法到这里已经完整落地了从原理到代码再到调参思路全部是我实际跑过的路径。如果你正在为某个一维信号降噪问题发愁可以先复制上面的代码跑通再拿自己的数据逐步调整尺度集合和融合策略。动手试过之后你对“多尺度”和“局部多项式”这两个词的理解会比任何教科书都深刻。