有段时间我在处理振动传感器的采样数据目标是识别轴承磨损瞬间产生的脉冲尖峰。信号里确实有个很明显的峰值但考虑到需要先做平滑降噪我直接用了最常用的滑动平均moving average。结果跑完一看脉冲尖峰基本消失了剩下一个极不起眼的鼓包峰值幅度掉了快一半。问题不是平滑本身而是平滑的方式不对。后来换成线性加权滑移平均算法同样的窗口长度下局部小峰值被清晰地保留下来噪声也压下去了。这篇文章就聊聊这个算法在Matlab里的实现思路、参数选择以及我实测中踩到的一些坑给做信号处理、时序分析、传感器数据预处理的同学做个参考。1. 为什么滑动平均总会抹掉小峰值权重摊薄的数学解释1.1 简单滑移平均的本质是一个矩形窗滑动平均是最朴素的平滑方法把最近L个点相加再除以L。从滤波器角度看它就是一个长度为L的矩形窗窗口内每个点的权重都是1/L。也就是说不管数据是刚刚发生的最新值还是已经过去很久的旧值对当前输出的贡献完全相同。问题恰恰出在这里。如果你想捕捉的局部小峰值刚刚出现在窗口最末端它只占总权重的一小部分。举个具体例子窗口长度L10时即使信号里这个脉冲的真实幅度是基线的3倍经过10点滑动平均后这个脉冲对输出的贡献最多只有原始幅度的10%。如果窗口里还有其他小幅波动这个小峰值会被彻底淹没。再加上窗口持续滑动峰值会在连续几个采样点里被反复摊薄视觉上就是一个矮胖的鼓包。这就是平滑后峰值消失的直观原因。1.2 线性加权滑移平均的权重设计与公式线性加权滑移平均的思路很直接既然新数据比旧数据更值得关注就按数据的年龄分配权重。窗口内最旧的数据权重为1越新的数据权重越大最新数据权重为L。权重序列可以写成 [L, L-1, ..., 1]最新在前也可以写成 [1, 2, ..., L]最旧在前本质相同。归一化之后它仍然是一个有效的加权平均所有权重之和等于 L(L1)/2。对当前时刻t和窗口长度L计算公式是y(t) [L·x(t) (L-1)·x(t-1) ... 1·x(t-L1)] / [L(L1)/2]从直观理解它还是平均只不过对最近的数据更偏心。从滤波器角度看它把矩形窗换成了一条递减的斜线窗频域上低通滤波的特性仍然保留但过渡带和纹波特性发生了变化。更重要的是这种设计让新信息在输出中拥有更高话语权这正是捕捉局部小峰值所需要的特性。1.3 两个指标直观对比峰值保留率与噪声衰减为了说明小峰值保留效果我构造两个极端场景算一笔账。场景一信号里有一个孤立的单点脉冲幅度A其他位置为零。在脉冲刚出现的那一刻简单滑动平均的输出是 A/L线性加权滑移平均的输出是A·L / [L(L1)/2] 2A/(L1)两者对比如下表窗口长度L简单平均脉冲保留率线性加权脉冲保留率线性加权的相对优势520%33.3%1.67倍1010%18.2%1.82倍205%9.5%1.90倍窗口越长线性加权对峰值保留率的好处越明显。这个峰值保留率在做特征检测时是核心指标因为它直接决定了峰值能不能被后续的阈值或寻峰算法看见。场景二纯白噪声。简单滑动平均把噪声标准差衰减为原来的 1/√L线性加权滑移平均的噪声衰减系数是sqrt(1²2²...L²) / [L(L1)/2]代入L10简单平均约为0.316线性加权约为0.357。两者噪声抑制水平接近但线性加权明显保留了更多信号幅度。相当于同样的噪声风险下信号的收益更高是一个相当划算的权衡。需要说明的是以上是单点脉冲的极限情况。如果局部峰值本身占了好几个采样点简单平均的保留率会随着峰值实际宽度的增加而变高差距会缩小。但方向是确定的让新数据获得更多权重峰值就更容易从平滑结果中露头。2. Matlab实现从movavg到一行卷积2.1 有Financial Toolbox时直接用movavg如果你的环境里安装了Financial Toolbox最快的方式是直接用movavg函数x randn(100, 1); L 10; ma_linear movavg(x, L, 0, linear);这里的第一个参数是原始序列第二和第三个参数是一对Lead/Lag窗口长度第四个参数指定类型为linear。对捕捉当前时刻附近的局部峰值来说Lag设为0就足够了意思是只看当前以及之前的L个数据点。需要提醒的是movavg默认的计算类型可能是指数型所以一定要显式写上linear。不过这个函数通常用在金融时序分析场景中很多人处理通用信号时不会第一时间想到它。另外如果目标环境没有Financial Toolbox这段代码就跑不了所以我更推荐下面这种不依赖工具箱的写法。2.2 不依赖工具箱filter实现与原理用Matlab自带的filter函数可以一行搞定线性加权滑移平均而且完全不需要额外工具箱。先看一下filter做的事y(i) num(1)·x(i) num(2)·x(i-1) ...这个形式恰恰就是线性加权滑移平均要的当前点乘以最大权重、前一个点乘以次大权重L 10; w L:-1:1; % 最新点权重L最旧点权重1 y filter(w, sum(w), x);sum(w)等于L(L1)/2就是归一化分母。filter会把每一项求和再除以这个分母输出长度和x一致实时流式数据也能用。这段代码的可移植性极强任何装了Matlab的机器都能跑是我实际项目中用的主力版本。还有一种离线快速写法y conv(x, L:-1:1, valid) ./ (L*(L1)/2);conv用valid选项时输出长度是N-L1每个输出点对应原始信号第L到第N个点等价于把窗口完整滑过整个序列。如果需要和原数组对齐使用要自己补前L-1个NaN。这个写法适合一次性离线验证快速看效果。2.3 边界对齐三种实现的关键差异这是最容易踩坑的地方我在这里吃过不少亏。三种实现在边界上的处理完全不同movavg专门处理了边界前L-1个点是基于已有数据计算的部分窗口结果没有NaN输出长度和输入一致。做整段分析时最省心。filter默认初始状态为零前L-1个点实际只用了部分输入数据但分母仍然用的是完整权重和。结果是边界处的数值系统性偏小看起来像是开头被压低了一截。conv valid直接丢弃边界输出从第L个点开始长度最短但边界最干净。我的建议是离线分析优先用movavg或conv如果为了可移植性用filter一定要主动把前L-1个点设为NaN或者直接忽略不要用它们做峰值检测。2.4 验证合成信号上的对比测试为了确认效果我构造了一个包含缓慢基线、高斯脉冲和随机噪声的测试信号代码如下rng(0); N 1000; t (0:N-1); x 0.5 * sin(2*pi*0.01*t) ... 3.0 * exp(-((t-500).^2) / 2) ... % 高斯小峰值 0.15 * randn(N, 1); L 10; y_simple movmean(x, L); % 简单滑动平均 w L:-1:1; y_linear filter(w, sum(w), x); y_linear(1:L-1) NaN; % 边界数据不参与比较 fprintf(原始信号峰值: %.3f\n, max(x(480:520))); fprintf(简单滑动平均峰值: %.3f\n, max(y_simple(480:520))); fprintf(线性加权滑移平均峰值: %.3f\n, max(y_linear(480:520)));我实际跑过很多次结果稳定。原始信号里峰值大约在3.0附近经过10点简单滑动平均峰值幅度掉到1.0以下改用线性加权滑移平均后峰值能保留下来到1.3左右。也就是说同样的窗口长度仅仅把权重从平权改成线性递减峰值幅度就能提高三四成。再看噪声区域两者都被压得很平肉眼几乎看不出区别。后续做峰值定位时直接对线性加权后的信号调用findpeaks设置一个略高于噪声底的高度阈值就能稳定把小峰值找出来。这在后续的窗口参数调整里非常有用。3. 参数怎么定窗口长度、峰值宽度与滞后折中3.1 峰值保留率如何随窗口长度变化窗口长度L是整个算法最关键的旋钮。从第1章的公式可以看出单点脉冲的保留率是2/(L1)L越大保留率越低。但L太小平滑效果就弱噪声压不住。这里存在一个典型的折中。如果目标峰值本身有一定宽度比如一个高斯脉冲占了K个采样点经验规律是当窗口长度约为峰值宽度的1.5到3倍时线性加权可以保留到原始幅度的50%到70%左右一旦L超过峰值宽度的3倍峰值幅度衰减会变得非常快继续加大窗口就没有意义了。所以第一条准则就是窗口长度应尽量与目标峰值宽度匹配而不是拍脑袋选个固定值。3.2 时间滞后线性加权比简单平均快多少滑动平均本质上会引入滞后这个滞后可以用权重的一阶矩来估算。简单滑动平均的权重重心在窗口正中间对斜坡信号的滞后约为(L-1)/2个采样点。线性加权滑移平均的权重重心偏向当前时刻滞后约为(L-1)/3个采样点。算一笔账L10时简单平均滞后约4.5个采样点线性加权滞后约3个采样点。也就是说用线性加权后的峰值时刻更接近真实发生时刻。如果监控系统对报警延迟敏感这个差异值得写进技术方案里。更重要的是滞后减小意味着峰值走形也小后续做特征提取时峰宽、峰高这些参数更接近真实值。3.3 一个实用流程用小峰值检测来定窗口我常用的定参流程是这样的先不做平滑直接在原始信号上粗看目标峰值的宽度K确认大概范围把L设置为K的1.5倍运行线性加权滑移平均对比平滑前后的峰值幅度和相邻波谷如果峰值还是太毛把L往大调如果峰被压扁了把L往小调固定一个L之后不要再为了画面好看频繁更改。有条件的还可以做一个一维扫描直观观察窗口长度对峰值幅度的影响Ls 3:2:25; peak_amps zeros(size(Ls)); for i 1:numel(Ls) y_tmp filter(Ls(i):-1:1, sum(Ls(i):-1:1), x); peak_amps(i) max(y_tmp(480:520)); end plot(Ls, peak_amps); xlabel(窗口长度 L); ylabel(峰值幅度);做出来的曲线通常是一条先平缓后陡降的线拐点附近就是合适的L。这个方法比直接猜参数可靠得多尤其适合数据特点不明确、需要快速探索的场景。3.4 不同应用场景的窗口推荐基于我自己的项目经验整理了一张常用窗口范围表可以当作起点应用场景推荐窗口长度L备注振动/加速度脉冲检测5~15建议先中值滤噪再做线性加权心电/脑电特征波QRS波、棘波10~30基线漂移慢窗口可以稍大金融收益率序列动量拐点5~20优先选小窗口避免滞后太大工业传感器趋势跟踪20~50主要用于趋势滤波而非峰检测这张表只是经验值真正的决定因素是采样率、目标峰值形状和噪声水平。团队里如果对参数有分歧最好的办法就是回到3.3的扫描流程用数据说话。4. 实战中的坑与延伸噪声尖峰、边界效应、实时流式4.1 边界数据不可直接用于峰检测第2.3节提过filter输出的前L-1个点相当于部分窗口完整归一化数值系统性偏小。我实际踩过这个坑一段数据的开头有一个起始脉冲用findpeaks怎么都找不到后来打印前几个输出点才发现第一个点的数值只有正常状态的18%左右。处理办法有三种主动把前L-1个点设为NaN后续任何分析直接跳过用movavg函数处理边界对数据做重叠切片每段往前多取L-1个点分析完再丢弃前半段。边界问题的本质是窗口内数据不足L个点时平均结果没有足够的统计意义。不能因为某个现成函数看起来没报错就忽略边界峰检测这种任务里边界错误经常导致漏检。4.2 噪声尖峰反而会被优待建议配合中值滤波线性加权对真实峰值和噪声尖峰的判断并不智能。只要一个突变点落在窗口最新位置它就能拿到最高的权重输出会明显跳一下。换句话说线性加权会放大高频突变的可见性。如果做异常检测噪声尖峰可能变成误报。我的做法是两步走先用medfilt1对原始信号做中值滤波窗口取3~5个点把孤立的噪声尖峰去掉再用线性加权滑移平均做平滑和峰值增强。中值滤波对孤立尖峰特别有效因为只看窗口内排序后的中位数一个异常点不会改变中位数。线性加权再对整理后的信号做加权平均既能保留局部小峰值又不会把噪声尖峰放大。这个组合我在振动数据和心电数据上都验证过效果比单用任何一种都稳定。4.3 实时流式数据的filter状态维护filter版本很适合实时流式处理每个新样本进来只需要更新内部状态。具体写法是把滤波状态变量透传下去L 10; w L:-1:1; den sum(w); [y, zf] filter(w, den, x_new, zi); zi zf; % 每次更新状态第一次调用时不传zi默认初始状态为0。系统刚启动时的前L个点仍然有边界效应建议做一个预热段前L个点只累计不判定等窗口填满了再开始做峰值检测。filter默认使用双精度浮点状态长时间运行精度没有问题。另外如果你在Simulink里做嵌入式代码生成filter也能顺利生成对应的C代码只是要注意输入数据类型尽量统一成形避免自动引入额外的类型转换开销。4.4 其他加权滑移平均的适用场景三角形与指数型线性加权只是越新越重要的一种实现方式。实际工程中还会用到另外两种三角形加权窗口中间权重最大两头小。它更像以窗口中心为参照的加权平均适合数据在一个窗口内相对平稳、你想突出中心状态而不是最新状态的场景比如光谱数据的局部平滑。指数加权权重按指数规律衰减对近期数据最敏感但理论上没有窗口边界所有历史数据都会以衰减形式参与。实时流式处理非常方便因为只需要保存一个累加状态但参数衰减因子的可解释性不如线性权重直观。选型逻辑可以根据需求来定需要明确窗口边界、希望权重可解释选线性做无限记忆的递推平滑选指数想在窗口内突出中心对称特征选三角形。我大部分峰值检测场景都选线性因为它把时间性和可控窗口结合得最好。4.5 扩展多尺度滑移平均或两级套用最后一个思路不要只用一组窗口。如果目标峰值宽度在不同时间尺度上变化比如一会儿是窄脉冲一会儿是宽缓波包可以同时跑几组不同L的线性加权滑移平均然后取它们的逐点最大值或交集。计算量也不算大本质上就是多跑几次filter。我自己的做法是两级方案L15负责宽带信号L215负责窄峰信号然后取两组输出的逐点最大值。这样窄峰不会被L2拖垮宽峰也不会被L1切碎多数情况都能覆盖到。这个做法在局部小峰值捕捉任务里比死磕一个窗口参数省心得多。最后说点个人体会线性加权滑移平均不是什么高深算法但它的核心思想——让新数据比旧数据更有发言权——几乎适用于所有时序平滑任务。在我的项目里从还算平滑变成既平滑又保真只是改了一行权重代码代价几乎为零。如果你也在被平滑后峰值消失困扰建议先别急着把窗口调小试着把简单平均换成线性加权再配合中值滤波大多数情况下会比你想的更管用。