简介这份MATLAB代码实现五点三次平滑滤波基于多项式最小二乘法逼近采样点专门用于波动曲线的去毛刺与趋势分析适合科研数据处理、论文图表绘制及实验信号预处理等场景。与样条插值、平均计算等方法相比该算法既保持滤波效果又兼顾计算简便性在保留数据主要特征的同时有效抑制高频噪声是平滑滤波中较为实用的选择。压缩包体积仅407B包含1个可直接运行的M文件代码由作者编写调试无复杂依赖下载后即可调用。目前已有2713人学习使用属于小而精的实用工具。通过该文件读者能掌握五点三次平滑的核心实现思路并可直接嵌入自己的数据处理流程中快速得到平滑后的曲线辅助判断数据走向与趋势变化。1. 五点三次平滑滤波是去毛刺最稳的滤波之一拿到一段传感器或业务曲线第一眼能看到总体趋势但锯齿状毛刺把局部形态搅得没法看。用移动平均会削峰、用中值滤波会丢坡度信息、用低通滤波要调截止频率简单粗暴都不合适。这时往往会想到一个在教科书里常驻但实际工程里容易忽视的工具——五点三次平滑滤波。它属于 Savitzky-Golay 滤波在窗口宽度 5、多项式阶数 3 条件下的特例本质上是在滑动窗口内做局部最小二乘多项式拟合用拟合值替代原始观测值。它对高频毛刺抑制明显且不会像普通平均那样把波峰波谷压扁适合处理振荡型、带趋势的数据曲线。做数据分析、设备监测、实验数据处理的人用它能同时满足曲线整形和趋势判断两大需求而且原理不复杂、三五十行代码就能跑起来。2. 为什么选“五点”和“三次”滑动窗口内的最小二乘拟合2.1 从局部拟合到加权平均五点三次平滑的核心操作可以这样理解任取原始序列中连续的 5 个点用一个三次多项式对该局部区间做最小二乘拟合取拟合多项式在这个区间中间点的值作为平滑后的输出。然后窗口向右滑动一个点继续重复同样的拟合。这与直接做移动平均的根本差别在于移动平均假设窗口内信号近似为常数而三次拟合假设窗口内信号是一个连续可微的低阶曲线因此能在平滑的同时保留二阶以内的变化特征。把这个过程展开会得到一个非常有用的结论对于等间隔采样的数据中间点的平滑值是这 5 个原始值的固定加权组合权重只与窗口宽度和多项式阶数有关和具体数据无关。对于三次多项式拟合的 5 点窗口中间点的权重系数是对称的越靠近中心点的数据权重越大两端的权重为负值。这个负权重让拟合可以补偿窗口内的曲率是它比普通滑动平均能更好地保留波形峰谷形态的真正原因。2.2 权重系数的来源与归一化设 5 个点对应的横坐标为 -2、-1、0、1、2三次多项式形式为 f(t) a₀ a₁t a₂t² a₃t³。最小二乘求解就是最小化残差平方和S Σ [x(tᵢ) - f(tᵢ)]²由于横坐标对称正规方程中奇偶次项彼此解耦a₀ 和 a₁ 的求解只涉及偶次项列向量。对中心点 t0 来说平滑值就是常数项 a₀通过正规方程展开可以得到a₀ (-3x(-2) 12x(-1) 17x(0) 12x(1) - 3x(2)) / 35于是中心点的 5 个权重就是 -3/35、12/35、17/35、12/35、-3/35权重之和等于 1。对常数信号完全不改变幅值这一性质保证平滑不会引入系统性偏移。2.3 为什么不用高次多项式或更大窗口把窗口增大到 7 点、9 点能获得更强的平滑效果代价是局部细节被抹掉实信号中的转折点会被明显钝化。把阶数升到四次、五次拟合曲线会更贴近原始数据毛刺抑制能力反而下降失去了平滑的本意。五点三次的组合是平滑能力和保形能力之间一个很实用的平衡点。尤其对采样率不高、信号本身又带明显趋势的曲线比如温度爬升曲线、电池放电电压曲线、负荷变化曲线选用这一组参数几乎不需要先验知识就能拿到可用的结果。3. 动手实现Python、NumPy 与嵌入式 C 的最小可运行版本3.1 用纯 NumPy 完成五点三次平滑不依赖 SciPy 的最小实现如下import numpy as np def smooth_5_3(x): 五点三次平滑滤波 Parameters ---------- x : np.ndarray 一维等间隔采样序列 Returns ------- y : np.ndarray 平滑后的序列 n len(x) y np.copy(x).astype(np.float64) if n 5: # 点数不足时不作处理避免边界公式越界 return y # 起点与终点采用专门推导的边界公式 y[0] (3 * x[0] 2 * x[1] x[2] - x[4]) / 5.0 y[1] (4 * x[0] 3 * x[1] 2 * x[2] x[3]) / 10.0 y[n - 2] (4 * x[n - 1] 3 * x[n - 2] 2 * x[n - 3] x[n - 4]) / 10.0 y[n - 1] (3 * x[n - 1] 2 * x[n - 2] x[n - 3] - x[n - 5]) / 5.0 # 中间段统一用五点权重系数 # 对应公式: y[i] (-3*x[i-2] 12*x[i-1] 17*x[i] 12*x[i1] - 3*x[i2]) / 35 w np.array([-3.0, 12.0, 17.0, 12.0, -3.0]) / 35.0 for i in range(2, n - 2): y[i] np.dot(x[i - 2 : i 3], w) return y这段代码里需要特别说明的是边界处理。因为窗口宽度为 5序列开头两个点和结尾两个点无法取满完整窗口如果直接跳过不处理把原始值原样输出边界处会保留毛刺后续做趋势分析时端点的导数也会异常。标准做法是用满足同样最小二乘准则、但在窗口端点的拟合值表达式来填补也就是代码中 y[0]、y[1]、y[n-2]、y[n-1] 这四个式子。它们不是凭空取的近似而是把拟合多项式分别代入 t 2、1、-1、-2 后得到的值。3.2 与 SciPy 的 savgol_filter 对照验证如果你的环境里有 SciPy可以直接用现成函数验证上面的手写实现是否正确。savgol_filter是通用 Savitzky-Golay 滤波器指定窗口大小 5、多项式阶数 3 时结果应与手写版一致from scipy.signal import savgol_filter # 生成带毛刺的测试信号 t np.linspace(0, 4 * np.pi, 200) x np.sin(t) 0.15 * np.sin(30 * t) 0.05 * np.random.randn(len(t)) y_manual smooth_5_3(x) y_scipy savgol_filter(x, window_length5, polyorder3, modeinterp) # 最大差异应该在浮点误差范围内 print(np.max(np.abs(y_manual - y_scipy)))savgol_filter的modeinterp表示边界用多项式拟合值推算默认情况不会改变长度也和手写边界公式思路一致。运行后最大差异通常在 1e-10 量级说明实现没有问题。实际项目里如果数据量很大、窗口需要动态调整直接调用 SciPy 更省事如果只在某个嵌入式模块里用手写版本更合适。3.3 嵌入式场景下的定点化写法单片机上没有浮点协处理器时用浮点数计算五点三次公式会拖慢中断处理流程。常见做法是把权重换算成整数系数再用移位完成除法。以下是一个适配整数采集值的版本#include stdint.h void smooth_5_3_int16(const int16_t *in, int16_t *out, uint32_t n) { if (n 5) { for (uint32_t i 0; i n; i) out[i] in[i]; return; } // 边界: 用 Q12 定点系数, 最终右移 12 位 out[0] (int16_t)((int32_t)(3 * in[0] 2 * in[1] in[2] - in[4]) * 819 / 4096); out[1] (int16_t)((int32_t)(4 * in[0] 3 * in[1] 2 * in[2] in[3]) * 819 / 4096); out[n - 2] (int16_t)((int32_t)(4 * in[n-1] 3 * in[n-2] 2 * in[n-3] in[n-4]) * 819 / 4096); out[n - 1] (int16_t)((int32_t)(3 * in[n-1] 2 * in[n-2] in[n-3] - in[n-5]) * 819 / 4096); for (uint32_t i 2; i n - 2; i) { int32_t acc -11 * in[i-2] 44 * in[i-1] 62 * in[i] 44 * in[i1] - 11 * in[i2]; out[i] (int16_t)(acc / 128); } }中间段把 -3、12、17、12、-3 这组系数乘到了 128 倍变成 -11、44、62、44、-11最后统一除以 128。这里选用 2 的幂次做定点缩放编译器会把除法优化成移位运算速度和代码体积都友好。边界公式里的系数 819 / 4096 是对 1/5 的 Q12 近似虽然只有约 0.02% 的误差但在整数除法下已经足够满足 16 位 ADC 数据的精度要求。4. 参数调优与踩坑窗口宽度、边界效应和迭代平滑的影响4.1 三个直接影响输出形态的操作参数五点三次平滑虽然只有窗口宽度和多项式阶数两个核心参数但在工程实现中还有几个容易被忽视的选择点。参数/操作典型取值影响适用情况窗口宽度5、7、9越大越平滑但峰谷越钝采样密集、毛刺频率高时用 7 或 9多项式阶数2、3、4阶数越高越贴合原数据保留波形细节用 3趋势线提取用 2边界模式interp、mirror、constant决定首尾 2 个点是拟合还是镜像默认 interp数据段短时慎用 zero 填充是否迭代13 次二次平滑可进一步提升去噪强度毛刺严重时迭代 2 次超过 3 次会失真多数情况下用窗口 5、阶数 3 已经足够。但要注意如果信号本身采样率很高相邻点之间差异极小窗口内 5 个点覆盖的时间区间太短平滑效果会非常有限。这时候应该增大窗口宽度至 7 或 9而不是保持 5 不变。反过来采样率不高、一个周期只有十来个点再用 9 点窗口会把整个波峰抹平到无法辨认。4.2 不要盲目迭代迭代次数存在临界点对同一段数据多次应用五点三次平滑等效于多次低通滤波会让频谱能量向低频集中。第一次迭代通常能压掉绝大部分高频毛刺第二次迭代让曲线更加顺滑但到第三次以后波形前缘的上升沿开始出现明显的滞后和预摆看起来像过冲反弹。判断迭代是否过度的简单方法是检查平滑后曲线是否出现原始数据中不存在的“新极值”。如果 1 秒内的抖动曲线在平滑后多出几个小折点说明迭代次数太多或窗口形状与数据不匹配。一个更稳妥的方法是每次迭代前记录平滑前后残差的标准差residual_std np.std(x - y_manual)然后按比例设定停止条件比如残差标准差小于原信号标准差的 10% 就不再迭代。这个阈值不是绝对标准但能在不同量纲的信号之间提供一个可迁移的尺度。4.3 边界失真数据段越短越明显整段数据只有 2030 个点时边界外的 4 个点会占很大比例边界拟合公式的作用会被放大。假设数据本身不是零均值直接套用五点三次公式会让首尾点往中间值方向轻微收缩导致端点出现明显的人为斜率。常见做法是在滤波前对整段数据做线性去趋势把两端平移到接近零的位置平滑之后再按原斜率加回来。这样边界公式处理的是围绕水平线的波动失真会小很多。如果要在滑窗实时流上做每次来一个新点就更新整个数组再全量滤波不仅耗时还会让新旧数据边界反复抖动。更好的方式是用环形缓冲区保留过去 4 个点每次只输出窗口中心点的平滑值。这种做法等价于因果滤波单次延迟为 2 个采样周期适合对实时性要求不高的采集系统。4.4 与移动平均、中值滤波的分工三点五次适合做趋势和低速波形的整形但不应替代所有滤波手段。遇到脉冲型毛刺比如开关切换带来的尖峰五点三次的负权重会让尖峰产生一圈小振荡不如中值滤波干净。遇到随机噪声均匀叠加在方波信号上移动平均的时间常数更好控制。实际流程中我会先用中值滤波剔除明显离群点再对清后的序列做五点三次平滑最后交给趋势分析模块三步各司其职效果优于任何单一滤波器。5. 数据预处理中的实战接上毛刺较多的曲线直接看趋势和转折5.1 用平滑后序列提取走向趋势与转折点拿到一条毛刺较多的曲线最常见的需求是把噪声压掉之后回答三个问题整体在升、降还是震荡趋势何时发生变化极值点到底出现在哪个时刻。五点三次平滑之后趋势可以用差分判断极值位置可以用一阶导过零定位def extract_trend(y, slope_window3): 基于平滑序列计算每点斜率并判断趋势段 dy np.diff(y) # 对差分再做一次平滑, 减少 residual 中的随机噪声干扰 dy_smooth smooth_5_3(np.concatenate([[dy[0]], dy, [dy[-1]]])) trend [] for i in range(len(y)): seg dy_smooth[max(0, i - slope_window) : i slope_window] mean_slope np.mean(seg) if mean_slope 0.005 * np.std(y): trend.append(up) elif mean_slope -0.005 * np.std(y): trend.append(down) else: trend.append(flat) return np.array(trend) # 定位极值点: 平滑后相邻差分变号的位置 cross np.where(np.diff(np.sign(dy_smooth)) ! 0)[0]判断趋势的阈值用0.005 * np.std(y)做相对化避免了不同信号之间需要反复试凑绝对阈值的麻烦。极值点定位要在平滑序列上做因为原始数据的大量毛刺会让np.sign(dy)在每一点都翻转找不到真正有效的转折。5.2 评估平滑效果的两个量化指标只有视觉判断不够还要给出可靠量化。一般用两个指标平滑前后残差的均方根误差描述保真度平滑序列的一阶差分绝对和描述毛刺抑制能力。rmse np.sqrt(np.mean((x - y_manual) ** 2)) roughness_original np.mean(np.abs(np.diff(x))) roughness_smoothed np.mean(np.abs(np.diff(y_manual))) print(f保真度 RMSE{rmse:.4f}, 粗糙度 {roughness_original:.4f} - {roughness_smoothed:.4f})粗糙度下降越大说明滤波越强而 RMSE 过大说明失真严重。用这两组数据配合可以快速判断当前参数是该加大窗口还是减少迭代远好过对着图形猜。检验滤波是否破坏原始趋势的另一个技巧对平滑序列和原始序列分别计算 Spearman 等级相关系数若相关系数仍高于 0.98基本可以断定这段数据的整体走向被完整保留了。本文还有配套的精品资源点击获取