传感器数据的降噪我入行那会儿第一个正式任务就是处理振动信号里的高频毛刺。当时图省事直接套滑动平均结果曲线确实光滑了可峰值被压扁、阶跃沿糊成一团被做测试的同事拿去对照原始数据一眼就看出了问题。后来接触到基于局部一维多项式的平滑方法才意识到降噪本质上不是把线画平滑而是要在保留结构的前提下去掉随机成分。再多走一步把单一窗宽换成多尺度组合就变成了既能照顾平坦段、又能守住尖峰和阶跃的方案。这篇文章就围绕MATLAB环境下的多尺度局部一维多项式信号降噪聊聊原理、完整实现、实验效果以及我在实际调试中踩过的一些坑。适合刚接触信号处理又不想只看理论公式的工程师也说给那些在MATLAB里试过sgolayfilt却不知道怎么调参数的人听。1. 为什么选择多尺度局部一维多项式三个理由1.1 局部多项式比滑动平均“聪明”在哪滑动平均的问题很直观它假设窗口里所有点的重要性一样所以遇到阶跃、尖峰这种局部结构时会把突变平均掉。局部多项式拟合的思路不同它在一个局部窗口内用低阶多项式去逼近真实信号再用这个多项式在中心点的值作为估计值。如果把多项式阶数设为0得到的其实是滑动平均阶数设为1相当于拟合一条直线能保住线性变化的趋势阶数设为2则能拟合曲率对慢慢弯曲的波形更友好。换句话说局部多项式是滑动平均的一般化版本多出来的阶数带来了“信号在局部是光滑函数”这样一个建模能力。这个区别在实际波形上非常明显。我做过一个标准测试方波加上正弦波再叠一点白噪声。滑动平均会把方波的上升沿变成一个斜坡而一阶局部多项式即使窗宽相同上升沿也要陡得多。原因是窗口中包含阶跃点时最小二乘拟合出的直线仍然以阶跃为中心线穿过给出的中心估计值更接近真实跳变不会只取周围点的简单平均。1.2 单一尺度解决不了的难题局部多项式要发挥作用必须选一个窗宽。窗宽小局部窗口里的信号变化可以被低阶多项式抓住但样本少噪声没被充分平均估计方差大。窗宽大平坦段的噪声能压得很低但遇到尖峰、阶跃、快速振荡这类结构多项式会在窗口内产生明显偏差。更麻烦的是真实信号往往不是只有一种结构一段长长的平缓波动以后突然来一个脉冲然后又进入比较稳定的平台。这种时候不管固定窗宽取多少都会顾此失彼。小窗宽能把脉冲守住但平缓段毛刺多大窗宽能换来平缓段的光滑脉冲却被削平。单一尺度的局部多项式本质上只有一个“分辨率”这是它的结构性局限不是调参数能完全弥补的。所以我想如果能同时跑好几个尺度的局部多项式让它们各管一块再把结果按某种规则拼起来是不是比单尺度更稳这就是多尺度思路的出发点。1.3 多尺度组合的本质偏差和方差的折衷任何平滑方法都要面对一对矛盾窗宽越大方差越小但偏差可能变大窗宽越小偏差小但方差大。多尺度方法想做的事情不是去猜一个“最优窗宽”而是承认不同位置的最优窗宽不同通过加权融合把多个尺度的优点合并起来。最常见的直观做法是让每个尺度中残差较小的区域拥有更高的权重。某个尺度如果拟合得好残差里基本是噪声能量低如果拟合得不好残差里会混入信号成分能量偏高。于是残差能量倒数可以当作一个粗略的置信度。更严格的做法是引入交叉验证对每个尺度打分但工程上先跑一个加权版往往已经能拿到不错的效果。而且多尺度输出的差异本身也携带信息如果大尺度和小尺度输出差异很小说明局部结构简单这时应该信任大尺度差异很大说明存在小尺度结构要信任小尺度。把这个逻辑做成权重等于把“结构复杂度”也纳入了降噪判断。2. 算法拆解局部拟合、尺度构造与融合策略2.1 局部一维多项式拟合的数学形式在时刻n取一个以其为中心的窗口窗口半径是h窗宽是w2h1。假设窗口内的真实信号可以用k阶多项式近似f(t) ≈ a0 a1(t - n) a2(t - n)² ... ak(t - n)^k用最小二乘求解系数然后把tn代入多项式得到中心点的估计值。因为(t-n)在中心为0所以估计值其实就是常数项 a0。这个特性让局部多项式对中心点的估计非常稳定不受坐标平移影响。如果只做一维信号MATLAB里最现成的实现是Savitzky-Golay滤波器也就是sgolayfilt(x, order, framelen)。它的本质就是给每个样本算一组固定的卷积系数而这组系数只由多项式阶数和窗宽决定与信号内容无关。所以它运行非常快适合批量处理。在MATLAB里使用时要注意框架长度framelen必须是奇数而且必须满足framelen order否则矩阵会变奇异拟合结果直接爆掉。这个约束在参数扫描时很容易被忽略我后面会专门讲。2.2 窗口宽、阶数、步长三个参数的桥梁作用多尺度局部多项式里核心是窗宽序列。窗宽决定局部范围也就是“时间分辨率”。阶数决定局部建模能力也就是“多边形复杂度”。步长在这里不常用因为我们是逐点估计不是从信号里抽样。但如果你处理的是超长信号也可以每隔几个点估计一次再插值牺牲一点分辨率换计算速度。窗宽序列怎么给我的经验是从小到大按倍数取比如[5, 11, 21, 41, 81]。最小窗宽对应信号里最小的有意义结构通常取你能容忍的最小时间尺度最大窗宽对应平坦段允许的最大平均范围一般不要超过信号长度的十分之一。阶数推荐从2开始。阶数太低只能建模分段线性的信号阶数太高局部窗口里数据量不够拟合会开始跟着噪声走。在窗宽较小的时候尤其危险所以我的默认组合是阶数2、窗宽序列[5, 11, 21, 41]。如果信号有明显分段线性特征比如方波、脉冲串我会把阶数降到1窗宽序列不变阶跃沿保存得更好。2.3 多尺度输出的加权融合方式多尺度输出融合可以有三种思路。第一种是简单平均所有尺度等权。这么做实现最方便但会把大尺度在小结构上的偏差也平均进去效果一般。第二种是全局权重对整个信号内的每个尺度计算一个残差RMS取RMS倒数的平方作为权重。这种方式能把明显失真的尺度压下去但对非平稳信号来说全局权重太粗糙前半段效果好的尺度在后半段可能完全不管用。第三种是局部权重也就是逐点或分块计算置信度。比如在每个尺度滤波后用滑动窗口统计残差平方的局部均值取其倒数作为这个采样点的权重然后再做一次时间平滑避免权重在相邻点之间跳变。这个方案相对复杂一些但对非平稳噪声和结构突变的适应性明显更好。我下面的完整实现会把第三种方式写出来同时保留全局权重的版本方便你自己对比。3. MATLAB完整实现从函数到完整脚本3.1 基础函数实现全局加权版先看第一个版本逻辑清晰适合理解整体流程function y msmp_denoise_global(x, order, wins) % MSMP_DENOISE_GLOBAL 多尺度局部多项式降噪全局权重融合 % 输入: % x 一维信号可以是行向量或列向量 % order 局部多项式阶数推荐 1 或 2 % wins 窗宽向量每个元素为奇数且大于 order % 输出: % y 长度和 x 相同的降噪信号 s size(x); x x(:); % 统一为列向量 N length(x); L length(wins); yAll zeros(N, L); rAll zeros(N, L); for i 1:L w wins(i); % 核心一行Savitzky-Golay平滑 ys sgolayfilt(x, order, w); yAll(:, i) ys; rAll(:, i) x - ys; % 残差 end % 每个尺度的残差RMS rmsAll sqrt(mean(rAll.^2, 1)); % 权重与残差RMS的平方成反比 wArr 1 ./ max(rmsAll(:).^2, 1e-12); wArr wArr / sum(wArr); % 加权融合 y yAll * wArr(:); y reshape(y, s); end这段代码非常短但值得解释几个点。残差RMS是“这个尺度从信号里拿掉了多少东西”。如果某个尺度把信号里的有效成分也当噪声滤掉了残差里就会出现明显的信号结构RMS会变大权重自然下降。这样等于自动给尺度打分。当然它没法区分“残差里的信号结构”和“残差里的噪声”所以这个方法保守如果两个尺度效果差不多它会偏向RMS更小的也就是偏向更平滑的结果。这是它的一个已知缺陷后面我们用局部权重版本弥补。3.2 逐点局部自适应多尺度版为了处理非平稳信号我把权重做成逐点变化的。核心变化是不再对整段信号只算一个权重而是对每个采样点算一组权重。权重的来源是局部残差能量在每个尺度滤波后用滑动窗统计残差平方的均值再取倒数。为了防止相邻点权重剧烈跳变融合后还要对权重序列做一次平滑。function y msmp_denoise_adapt(x, order, wins) % MSMP_DENOISE_ADAPT 多尺度局部多项式降噪逐点局部权重融合 % 输入输出格式同 msmp_denoise_global s size(x); x x(:); N length(x); L length(wins); yAll zeros(N, L); powerAll zeros(N, L); for i 1:L w wins(i); ys sgolayfilt(x, order, w); yAll(:, i) ys; r x - ys; % 局部残差能量窗口宽度取当前尺度的一半减少对远端数据的依赖 localWin max(3, floor(w/2)); powerAll(:, i) movmean(r.^2, localWin, Endpoints, fill); end % 权重 局部残差能量的倒数先防止除零 W 1 ./ max(powerAll, 1e-12); % 对权重做时间平滑抑制权重本身的高频抖动 smoothWin max(5, 2 * floor(length(x)/200) 1); W movmean(W, smoothWin, Endpoints, fill); % 归一化每一行的权重 Wsum sum(W, 2); W W ./ max(Wsum, 1e-12); % 融合 y sum(yAll .* W, 2); y(isnan(y)) x(isnan(y)); % 如果有NaN保留原始值 y reshape(y, s); end这个版本相对更实用。movmean的Endpoints, fill会让信号两端出现NaN导致最后的求和结果也有NaN所以我专门做了NaN兜底如果某个点因为端点效应得到NaN就直接保留原始值。这段代码在MATLAB R2020a之后都能跑不需要额外工具箱只要Signal Processing Toolbox提供sgolayfilt即可。3.3 演示脚本非平稳噪声环境下的对比光看函数不够我搭一个模拟场景来验证效果。信号包含正弦振荡、方波和尖峰噪声方差在后半段突然加大用来模拟传感器跑一段时间后漂移或干扰变强的情况。fs 1000; t (0:1/fs:4-1/fs); N length(t); % 干净信号5Hz正弦 1Hz方波 一个尖峰 xTrue 2 * sin(2*pi*5*t) ... 1.2 * square(2*pi*1*t, 35) ... 2 * exp(-((t - 2.5).^2) / (2 * 0.005^2)); % 非平稳噪声后2秒噪声幅度变大 sigma 0.15 0.45 * (t 2); xNoisy xTrue sigma .* randn(N, 1); % 一组较宽的窗 wins [5, 11, 21, 41]; y1 msmp_denoise_global(xNoisy, 2, wins); y2 msmp_denoise_adapt(xNoisy, 2, wins); y3 sgolayfilt(xNoisy, 2, 11); % 单尺度对比 % 指标 snr_true (y) 20 * log10(norm(xTrue) / norm(xTrue - y)); fprintf(全局多尺度 SNR: %.2f dB\n, snr_true(y1)); fprintf(局部多尺度 SNR: %.2f dB\n, snr_true(y2)); fprintf(单尺度 SNR: %.2f dB\n, snr_true(y3));在我的一次运行里局部多尺度版SNR比单尺度11点窗高1.8 dB左右全局多尺度版居中。重点不是这个具体数值而是趋势噪声非平稳时单尺度几乎没有调整空间局部多尺度至少能在噪声大的区段自动降低对大窗的依赖保住更多细节。如果你想深挖可以把wins换成更密的序列比如[5,7,9,11,15,21,29,41]结果会更稳定但计算时间涨得也快。我后面会给出比较实用的参数模板。4. 实验效果与参数选择建议4.1 定量指标SNR、RMSE、MAE怎么算评估降噪效果我习惯同时看三个指标。SNR关注整体能量比信号能量和误差能量相比有多大。RMSE是均方根误差把误差平方后取平均再开方对大误差很敏感。MAE是绝对误差均值对大误差没那么敏感更适合看算法是否在某些点上出现严重失真。rmse sqrt(mean((xTrue - y).^2)); mae mean(abs(xTrue - y)); snr 20 * log10(norm(xTrue) / sqrt(N * rmse^2));实际使用时如果手头没有干净的xTrueSNR算不了那就只能靠残差分析和目检。残差如果还有明显的波形结构说明信号被滤掉了残差如果是白噪声的样子说明降噪比较彻底。我会把残差频谱也画出来看看有没有固定频率的泄漏。4.2 阶数、窗宽、融合方式对结果的影响阶数的影响最直接。order1时局部模型是直线对方波和快变沿友好但慢变曲线的拟合残差会偏大。order2时能拟合弯曲是最通用的选择。order3以上我不太建议除非源信号本身采自一个已知的高阶动力学系统否则多项式阶数越高越容易把局部窗口内的细节错误放大。窗宽的影响更微妙。最小窗宽设得太小比如等于order1拟合几乎就是插值噪声没有机会被平均最小窗宽设得太大尖峰就会被削。最大窗宽设得太大如果信号里有缓慢漂移多尺度融合会把漂移当噪声处理导致基线被拉平。我见过最典型的问题就是有人把窗宽设到信号长度的一半以上结果整段信号几乎被压成一条直线这是误区。融合方式上全局权重适合“整段信号噪声水平基本一致”的情况。逐点局部权重适合“噪声方差变化、信号结构密度变化”的情况。代价是逐点权重对参数更敏感权重平滑窗如果太大局部权重就失去局部性退化成全局权重如果太小则可能出现权重频繁切换的噪声。4.3 我常用的几组参数模板我整理了三组实际项目中用得最多的参数直接抄作业基本可行信号类型orderwins融合方式温度、压力等缓慢物理量2[15, 31, 61, 121]全局权重振动、音频、瞬态信号1或2[5, 11, 21, 41]局部权重生理信号如ECG、肌电2[7, 13, 25, 51]局部权重外加阶跃保护这些模板不是铁律但能省掉前期一多半的试错。更重要的是拿到新数据的第一件事永远是画直方图、画频谱搞清楚噪声到底是什么颜色再决定要不要用局部多项式。如果噪声本身有强窄带干扰比如50Hz工频我会先在频域做陷波再做多尺度局部多项式降噪。5. 实踩的坑与排查技巧5.1 五个高频问题和对应解法局部多项式方法看着简单实际调试时问题并不少。我整理了一个对照表都是自己或同事真遇到过的现象可能原因解决办法输出仍有高频毛刺最小窗宽还是太大或order太低把最小窗宽降到5或把order升到2阶跃沿被磨圆大尺度权重过高局部权重未生效改用局部权重版检查权重平滑窗是否过大尖峰衰减明显窗宽序列中的最小值不够小加入w3或5的窗并提高该尺度权重上限端点翘起局部多项式在边界外拟合失真先镜像延拓40个点滤波后裁掉输出比输入更“抖”阶数过高局部过拟合降低order检查wins里小窗是否过少端点翘起这个问题特别容易被忽视。sgolayfilt在处理边界时会用单边多项式外推如果信号边界有明显趋势外推就会翘得厉害。我在处理长记录时通常先做镜像延拓padLen 50; xp [flipud(x(1:padLen)); x; flipud(x(end-padLen1:end))]; % 滤波 yp msmp_denoise_adapt(xp, order, wins); % 裁掉 y yp(padLen1:end-padLen);这个方法不用改函数内部通用性很强。镜像延拓必须基于信号本身的连续趋势如果是纯随机噪声延拓效果不明显。5.2 几个容易被忽略的实现细节第一个是sgolayfilt的输入维度问题。它默认对每一列做滤波如果你把行向量传进去输出也是行向量但这和很多其他函数不太一致容易在后续拼接时出错。我在函数开头统一x x(:)就是为了避免这种坑。第二个是窗宽必须严格大于阶数且为奇数。如果你用for循环扫描窗宽偶尔传入偶数MATLAB不会报错但滤波结果会悄悄改变很难排查。建议在函数开头加一个校验if any(mod(wins, 2) 0) || any(wins order) error(窗宽必须全为奇数且必须大于多项式阶数); end第三个是NaN。如果原始信号里有一段缺失值直接跑sgolayfilt会把NaN扩散到邻近点。稳妥做法是先插值或填充再用一个遮罩把降噪结果中对应位置替换回NaN避免污染统计指标。上面局部自适应版里的NaN兜底只做了一层防护更严格的数据应该单独做缺失值处理。第四个是计算量。逐点局部权重需要把每个尺度都保存一份完整结果如果信号很长、尺度很多内存占用会非常可观。1e7长度的数据配8个尺度大约需要640MB内存用于纯矩阵存储加上残留和权重很容易逼近1GB。我建议超过500万点就分块处理或者减少尺度数量到4个以内。5.3 我遇到的一次典型翻车一次做超声回波信号降噪信号本身有尖锐的起振我贪心把窗宽序列设成了[3, 7, 15, 31, 63]order用了3。结果局部权重版本把起振点附近的权重全部给了最小窗因为最小窗残差能量最低导致输出几乎没降噪起振后还是有大量高频毛刺。后来我把order降到1并且强制了一个权重比例上限任何单一尺度的权重不超过总权重的50%。这个限制很简单W min(W, 0.5);在归一化之前加这一行就避免了小窗独霸权重。这个改动不明显但对瞬态信号效果立竿见影也让结果对窗宽选择不那么敏感了。6. 扩展思路把多尺度局部多项式用到更复杂的地方6.1 从一维到二维和多通道一维局部多项式可以直接推广到二维图像也就是局部二元多项式曲面拟合。MATLAB里没有直接对应的sgolayfilt二维版但可以用分块处理每一行分别用一维多尺度局部多项式降噪结果只解决“行向”噪声对于各向同性噪声更好的做法是构造二维Savitzky-Golay核不过代码复杂度会高很多。如果只想快速验证二维扩展可以把图像按行展开成多个通道用同一个函数处理再拼回来至少在处理扫描图像纹理时是有效的。多通道信号比如三轴加速度计或EEG多导联通道间往往存在相关性。多尺度局部多项式目前没有利用这个信息但可以先对每个通道独立降噪再做通道间校准。若通道间会存在相位差延迟补偿得多做一步否则降噪结果里会出现通道间不自然的错位。6.2 与自适应、鲁棒化结合的升级方向局部多项式对异常值很敏感因为最小二乘没有抗野值能力。如果信号里有尖峰干扰而且你不想保留它直接滤波会把尖峰摊到邻近点上。更稳的方式是在迭代中引入鲁棒权重先算一次局部多项式残差对残差大的点降权再做一次局部多项式类似迭代加权最小二乘。这个思路可以方便地嵌入多尺度框架中每个尺度各自做鲁棒迭代最后融合。另一个升级方向是让窗宽本身也变成局部的每个点估计一个最合适的窗宽而不是简单地在多个固定窗宽间插值。这个思路更彻底但需要额外的窗宽选择准则计算量会明显上升。多尺度融合实际上是一个折衷它的优点是能保留多个尺度各自的信息不会发生“选错唯一窗宽”这种灾难性错误。最后的扩展方向是跟频域方法结合。局部多项式是时域方法对非平稳结构和突变捕捉好频域滤波对窄带噪声抑制强。多尺度局部多项式输出后再用频谱残留分析做一次后验检查观察哪个频段的噪声没有压下去就可以补充一个陷波器或小波阈值步骤。我通常把这套组合叫作“时频双通道降噪”在工程上比单独用任何一种方法都稳。在实际使用中我的体会是多尺度局部多项式最大的价值不是某一次调试跑出多高的SNR而是它给降噪问题提供了一个“由粗到细”的观察框架。你看到的输出不是单一滤镜的结果而是多个分辨率下的模型在互相博弈后达成的共识。这种思路迁移到数据拟合、异常检测、趋势提取上同样能带来启发。如果真的要对一个全新的降噪任务做方案选型我会先用这个多尺度框架跑一遍从多尺度的残差和权重分布里快速读出信号的局部结构特征再决定要不要上更复杂的深度网络或非局部算法。这个前置诊断价值比那一点SNR提升更值得你花时间去掌握。