简介:Savitzky-Golay滤波算法的C语言实现,可直接运行于STM32单片机,用于采集信号时消除随机噪声。该算法由Savitzky和Golay于1964年提出,基于局域多项式最小二乘拟合,在滤除高频干扰的同时保持原始信号的形状与带宽特征,特别适合温度、压力、惯性测量等传感器数据的平滑预处理。整个压缩包仅4KB,共4个文件,其中2个C源文件与2个头文件分别构成滤波主程序与配套矩阵运算模块,接口清晰,可快速嵌入MDK或IAR工程,适合资源受限的嵌入式平台。已有1996人浏览学习,应用验证充分。阅读源码能直观理解S-G滤波器窗宽、多项式阶次等参数对平滑效果的影响,并可将该模块直接用于数据采集系统,提高信噪比,减少随机脉冲造成的误判,有效提升工程测量的稳定性与可信度。
1. 为什么嵌入式平滑滤波最后都会绕回 Savitzky-Golay
你大概率遇到过这种情况:ADC 采出来的曲线毛刺多到没法看,第一反应是滑动平均,窗口调到 16,毛刺没了,但峰值也被削平了,原本该尖锐的过零点和波峰全变成圆角。这不是参数没调好,而是滑动平均的本质缺陷——它把所有数据都当成同一阶多项式来拟合,细节跟着窗口一起被抹掉。Savitzky-Golay 滤波器(简称 SG 滤波)恰恰是冲着这个问题去的:它在滑动窗口内做一次最小二乘多项式拟合,然后用拟合多项式在中心点的值替代原始数据。因为拟合本身带阶数,它能在去噪的同时保留二阶矩、三阶矩这些形状特征,峰值保留能力比滑动平均高一个量级。这个滤波器应用极广,从光谱平滑、ECG 去噪到传感器数据预处理都有它的身影,而且它的核心计算不过是一组固定系数的卷积——这正是它能在 STM32 这类单片机上跑起来的原因。不需要矩阵求逆、不需要迭代收敛,只要把系数表提前算好,运行时就是乘加运算,C 语言几十行就能实现。
对 STM32 开发者来说,SG 滤波器最大的诱惑不在于精度,而在于确定性:窗口固定、计算量恒定、没有动态内存分配,完全没有实时性意外。你不需要 DSP 库,不需要浮点单元,甚至用 int16 整型运算也能得到可用的输出。这篇博文就顺着一条完整的落地路径走:先帮你理解系数从哪来,再在 PC 上用 Python 快速验证参数效果,最后给出能在 STM32 上直接编译运行的 C 语言实现,以及边界处理、整型优化这类真刀真枪才能踩到的坑。
2. 先搞懂 SG 滤波的本质:滑窗多项式拟合与固定卷积系数
2.1 为什么 SG 滤波是一组固定系数,而不是真正的“拟合”
很多人第一次看 SG 滤波的代码会困惑:说好的最小二乘拟合,怎么代码里全是乘加?没有求逆,没有求解方程组,甚至连浮点除法都很少看到。这个疑惑的答案在“滑窗”这两个字上。
假设我们有一个窗口,窗口内有 2m+1 个点(m 是半窗宽),记为 x[-m], x[-m+1], ..., x[0], ..., x[m],现在要用一个 M 阶多项式:
p(t) = c₀ + c₁·t + c₂·t² + ... + c_M·t^M
去拟合这组数据。注意这里 t 是相对坐标,窗口中心点对应 t = 0。我们需要求解系数 c₀ 到 c_M,使得多项式在 t = -m 到 t = m 处的取值与原始数据误差平方和最小。这是一个标准的最小二乘问题,可以用正规方程求解,计算量不大。但关键点在于:窗口滑到下一个位置时,数据变了,拟合系数也随之变化,每个位置都做一次拟合听起来非常昂贵。
但 SG 滤波有一个极其漂亮的数学性质:由于窗口内坐标是固定的(-m 到 m),拟合过程对数据是线性的,最小二乘解可以写成数据的线性组合。也就是说,对窗口中心点做拟合得到的平滑值,等价于把窗口内每个点乘以一个固定系数再求和:
y[0] = Σᵢ h[i] · x[i],其中 i 从 -m 到 m
这个 h[i] 就是 SG 卷积系数。它只取决于两个参数:半窗宽 m 和多项式阶数 M,与数据本身无关。所以“滤波”这个动作最终落在了一组与数据无关的常数上,运行时只需要做卷积。这也就是为什么 SG 滤波在 DSP 语境下被归类为 FIR 滤波器的一个特例——它确实就是一个具有特定系数序列的有限冲激响应滤波器。
2.1.1 从正规方程到卷积系数,工程上不必自己推
如果你不打算深入研究数学推导,完全可以每次用 SciPy 或 MATLAB 生成系数后直接拿过来用。但理解一个约束总归有用:多项式阶数 M 必须小于窗口内点数 2m+1,工程上一般取 M = 2 或 3。阶数太低会过度平滑,阶数太高则去噪能力变差,甚至出现数值不稳定。另一个约束是窗口宽度必须是奇数,因为要保证有明确的中心点。
工程上的做法是:在 PC 上用 Python 或 MATLAB 把系数算好,存成 C 语言数组,直接烧进 STM32。你不需要在单片机上做任何矩阵运算,那既浪费 Flash 也浪费 RAM。真正的计算量是每个输出点做 2m+1 次乘加,以 5 点窗口、3 阶多项式为例,每输出一个点只需要 5 次乘法和 4 次加法——这在任何 STM32 上都是可以忽略不计的开销。
2.2 窗口宽度与多项式阶数的配合,直接决定滤波特性
SG 滤波有两个“旋钮”:半窗宽 m 和多项式阶数 M。它们不是独立的,组合起来决定滤波器的频率响应。
从频域看,SG 滤波器的通带宽度随窗口增大而变窄。窗口越大,低频保留越好,高频衰减越多,噪声抑制越强。但副作用是峰值保留能力下降——虽然优于滑动平均,但窗口大到一定程度,信号本身的快速变化区域也会被当作“噪声”抹掉。具体到工程选择上,一个经验值是:对采样率 1 kHz 的心电信号,做基线漂移抑制时半窗宽取 20~50,阶数取 2~3;对慢变的温度传感器信号,采样率 1 Hz,半窗宽取 5~10 就足够;对光谱数据,窗口往往很大(半窗宽 50 以上),但阶数一般不超过 3。
| 窗口点数 | 多项式阶数 M | 特点 | 适用场景 |
|---|---|---|---|
| 5 (m=2) | 2 | 去噪轻微,保留细节好 | 快速变化的信号,轻微毛刺 |
| 9 (m=4) | 3 | 平衡型,峰保留与去噪兼顾 | ECG、PPG 等生理信号 |
| 21 (m=10) | 2 | 平滑力度大,适合慢变信号 | 温度、压力、液位 |
| 51 (m=25) | 3 | 强平滑,注意峰值衰减 | 光谱、基线估计 |
有个常见的误区我必须提醒:阶数越高不代表越好。M 升高会让滤波器更“追随”数据的局部形状,去噪能力反而下降。M 的作用更像一个“保真度”旋钮——它决定你有多信任信号的形状。如果你想保留尖锐的波峰,M 取 3;如果信号本身就平滑,只是叠加了随机噪声,M 取 2 性价比最高。更高阶(4 以上)在工程上极少使用,因为提升有限,而计算系数时可能出现病态矩阵,系数值也会变得很大。
2.3 SG 滤波与滑动平均、低通滤波的关系:互为特例
把这三种滤波方式放在一起比较,你会对 SG 滤波的定位有一个更清晰的认识。滑动平均是 M=0 阶的 SG 滤波——常数拟合,窗口内所有点权重相等,因此它对突变的响应是线性的,会直接削平峰值。RC 低通(一阶 IIR)依赖历史输出,存在相位滞后,而 SG 是零相位(对离线数据做中心平滑时),在线流式处理时也只引入 (m) 个采样点的延迟。
从实现角度说,这三种滤波器在 STM32 上都能轻松运行,计算量差距微乎其微。但 SG 滤波的延迟是固定的、可控的,而且不会像 IIR 那样累积数值误差。如果你面对的信号既需要去噪又需要保持波形形态,SG 就是三者中最稳妥的选择。反过来,如果你的信号对实时性要求极高,每多 1 ms 延迟都不可接受,那 SG 滤波并不合适——它天生需要“看”未来的数据(窗口内有右侧的点),只能延迟输出。
3. 先在 PC 上用 Python 验证系数与滤波效果,再移植到 STM32
3.1 用 SciPy 生成 SG 系数表,并直观理解数值
动手写 STM32 代码之前,我强烈建议先在 PC 上用 Python 把系数和效果验证一遍。原因很实际:SG 滤波的直观效果强依赖于窗口和阶数的组合,你在单片机上烧一次程序看波形,远不如在 PC 上花 30 秒钟生成十组参数对比来得快。下面这段代码做的事是:用 savgol_coeffs 生成 SG 系数,然后分别作用在一段带有噪声的正弦信号上,输出滤波前后的对比数据,方便你观察峰值衰减情况。
import numpy as np from scipy.signal import savgol_coeffs # 生成 9 点窗口、3 阶多项式的 SG 系数 coeffs = savgol_coeffs(window_length=9, polyorder=3) print("系数:", np.round(coeffs, 6)) # 归一化确认:系数之和应等于 1 print("系数和:", np.sum(coeffs))这段代码的输出会是一组关于中心点对称的系数,例如 9 点窗口 3 阶的系数大致是[-0.09, 0.06, 0.20, 0.31, 0.33, ...]这种形态,中心点权重最大,两侧逐渐变小。注意系数之和为 1,这是 SG 系数的固有性质——保证对常数信号输出等于输入,不会引入直流偏置。这是一个非常实用的校验手段:如果你从某个资料里抄了一组系数,先检查它是否满足对称性以及和为 1,不满足就该怀疑抄错了。
3.1.1 生成系数后用 FIR 方式验证滤波效果
系数确认无误后,用 scipy.signal.lfilter 或直接手动卷积都能做滤波验证。下面这段代码用了一个更贴近嵌入式实现的验证方式:手动卷积,逐点输出,模拟 STM32 上的处理流程:
import numpy as np import matplotlib.pyplot as plt from scipy.signal import savgol_coeffs # 生成带噪信号:1kHz 采样率,10Hz 正弦波 + 白噪声 fs = 1000 t = np.arange(0, 1.0, 1/fs) x = np.sin(2 * np.pi * 10 * t) + 0.2 * np.random.randn(len(t)) # SG 系数:窗口 11,阶数 3 window_len = 11 polyorder = 3 coeffs = savgol_coeffs(window_len, polyorder) # 手动卷积实现 SG 滤波(流式方式,与 STM32 端一致) y = np.zeros_like(x) m = window_len // 2 for n in range(m, len(x) - m): segment = x[n - m : n + m + 1] y[n] = np.dot(coeffs, segment) # 卷积求和 # 边界不处理,直接复制原值 y[:m] = x[:m] y[-m:] = x[-m:] # 对比滤波前后 plt.figure(figsize=(10, 4)) plt.plot(t, x, alpha=0.4, label='Raw signal') plt.plot(t, y, linewidth=1.5, label='SG filtered') plt.legend() plt.show() # 计算信噪比改善(粗略) snr_raw = 20 * np.log10(np.std(x) / np.std(x - np.sin(2 * np.pi * 10 * t))) snr_fil = 20 * np.log10(np.std(y) / np.std(y - np.sin(2 * np.pi * 10 * t))) print(f"SNR raw: {snr_raw:.2f} dB, SNR filtered: {snr_fil:.2f} dB")这段代码里有一个细节值得注意:np.dot(coeffs, segment)这行就是整个 SG 滤波的核心运算,在 STM32 上对应的就是那个 for 循环里的乘累加。我在 PC 上刻意用这种方式而不是直接调savgol_filter,就是为了让“PC 验证代码”和“单片机 C 代码”在逻辑结构上保持一一对应关系,方便移植后对照调试。边界部分直接复制原值,这是最简单的策略,只适用于信号平稳段;后面第 5 章会讨论更好的边界处理办法。SNR 的计算只是给你一个快速判断参数效果的量化指标,不必过于较真它的绝对值。
3.2 输出 C 语言系数表,直接作为 STM32 工程的常量数组
参数验证完毕,下一步就是导出系数表。这里有一个工程建议:直接让 Python 生成一个格式化为 C 语言数组的文本,复制粘贴进工程,而不是手动抄。原因很直白——浮点数转字符串再转回来这个过程容易出错,让程序生成可以减少一个手误的来源。
import numpy as np from scipy.signal import savgol_coeffs window_len = 9 polyorder = 3 coeffs = savgol_coeffs(window_len, polyorder) # 输出可直接粘贴到 C 工程的数组 print(f"static const float sg_coeffs[{window_len}] = {{") for i, c in enumerate(coeffs): end = "," if i < window_len - 1 else "" print(f" {c:.8f}f{end}") print("};") # 也可以用定点格式输出,方便直接替换成整数运算 scaled = np.round(coeffs * 1024).astype(int) print(f"\nstatic const int16_t sg_coeffs_q10[{window_len}] = {{") for i, c in enumerate(scaled): end = "," if i < window_len - 1 else "" print(f" {c},{end}") print("};")输出的两组数组,第一组是 float 类型,适合带有 FPU 的 STM32F4/H7 系列,代码简洁直接。第二组是 Q10 格式的整型系数,适合 F1 这类没有硬件浮点单元的型号。关于 Q10 格式的定点实现细节,第 4 章会专门展开。现在你只需要知道:系数表一旦生成,滤波运行时的计算就只是乘加,与离线计算完全等价。至此,PC 端验证的工作就完成了,可以正式进入 STM32 实现阶段。
4. 在 STM32 上用 C 语言实现一个可流式处理的 SG 滤波器
4.1 循环缓冲区设计:让窗口滑动不需要搬移数据
真正开始写 STM32 代码时,你面临的第一问题不是滤波本身,而是数据组织。如果每次滤波都把数组整体左移一位,一个 9 点窗口的滤波器代码会写成这样:
for (int i = 0; i < N - 1; i++) buffer[i] = buffer[i+1]; buffer[N-1] = new_sample;这段代码逻辑没错,但它做了 N-1 次内存拷贝。在 PC 上无所谓,但在 STM32F103 这类主频 72 MHz 的芯片上,如果 ADC 采样率为 10 kHz,每个周期浪费的拷贝时间累积起来很可观。更严重的是,这种写法意味着滤波函数内部有状态,一旦被中断打断再重入,缓冲区内容就乱了。
更好的方案是循环缓冲区:用一个固定大小的数组,外加一个“写指针”指示最新数据的位置。采样到来时,写入指针指向的位置,然后指针递增并在到达数组末尾时回绕。滤波时,从写指针往前倒推 m 个位置开始取数,经过回绕边界时做取模运算。C 语言的取模开销略高,所以实际工程常用位运算替代——缓冲区长度取 2 的幂次,回绕操作就变成了(idx + 1) & (BUFFER_SIZE - 1)。大部分 STM32 工程的采样缓冲区长度本来就用 2 的幂,这个约束不难满足。
4.2 完整可编译的 SG 滤波器 C 代码(浮点版)
下面这段代码是可以在 STM32 上直接编译运行的完整实现。浮点版优先给出,因为它的可读性最好,也最容易对照第 3 章的 Python 代码验证正确性。
/* sg_filter.h */ #ifndef SG_FILTER_H #define SG_FILTER_H #include <stdint.h> #define SG_WINDOW_LEN 9 /* 窗口点数,必须为奇数 */ #define SG_HALF_WIN 4 /* 半窗宽,(SG_WINDOW_LEN-1)/2 */ #define SG_POLY_ORDER 3 /* 多项式阶数 */ /* 滤波函数:输入新样本,返回滤波结果 */ float sg_filter_process(float new_sample); /* 重置滤波器状态(例如切换采样率时调用) */ void sg_filter_reset(void); #endif/* sg_filter.c */ #include "sg_filter.h" #include <string.h> /* 系数表由 Python/SciPy 生成,9 点窗口 3 阶 */ static const float sg_coeffs[SG_WINDOW_LEN] = { -0.09090909f, 0.06060606f, 0.16883117f, 0.23376623f, 0.25541126f, 0.23376623f, 0.16883117f, 0.06060606f, -0.09090909f }; /* 循环缓冲区,长度取 2 的幂便于位运算回绕 */ #define SG_BUFFER_SIZE 16 static float sg_buffer[SG_BUFFER_SIZE]; static uint8_t sg_write_idx; void sg_filter_reset(void) { memset(sg_buffer, 0, sizeof(sg_buffer)); sg_write_idx = 0; } float sg_filter_process(float new_sample) { float result = 0.0f; uint8_t idx; /* 写入新样本 */ sg_buffer[sg_write_idx] = new_sample; sg_write_idx = (sg_write_idx + 1) & (SG_BUFFER_SIZE - 1); /* 卷积:从最新样本往前取窗口数据 */ /* 窗口中心对应最新样本(延迟 SG_HALF_WIN 个采样点) */ for (int i = 0; i < SG_WINDOW_LEN; i++) { /* 从写指针往前倒退,注意回绕 */ idx = (sg_write_idx - 1 - i) & (SG_BUFFER_SIZE - 1); /* 注意:coeffs 的索引方向与取数方向相反 */ result += sg_buffer[idx] * sg_coeffs[SG_WINDOW_LEN - 1 - i]; } return result; }代码逻辑看起来简单,但有几个点需要仔细说明。首先,为什么取数方向是从写指针往前倒退?因为写指针指向的位置在写入完成后是“最新数据”的下一个空位,真正的最新样本在sg_write_idx - 1的位置,往前倒退 i 步就能依次取到从新到旧的全部窗口数据。其次,系数索引为什么要反向?因为 SG 系数的定义是“以窗口中心为 0 点”的坐标顺序,而我们从最新样本(中心)开始取数时,先后顺序是中心、左侧、更左侧,对应到系数表需要反转。如果你跳过这两处细节,输出会看起来“有点对但波形偏移一个采样点”或者“感觉相位反了”,这是移植 SG 滤波最常见的 bug,没有之一。
注意这段代码里的 SG_BUFFER_SIZE 取 16 而不是 9,是为了位运算回绕。多余的 7 个位置不会被用到,但换来了回绕操作从取模变成按位与,F1 上能省十几个周期。如果你的采样率不高,比如 1 kHz 以下,直接用取模% SG_BUFFER_SIZE也行,可读性更好。
4.3 整型定点化改造:在没有 FPU 的 STM32F1 上跑 SG 滤波
STM32F103 系列和多数国产替代型号(如 GD32F103、APM32F103)没有硬件浮点单元,浮点运算依赖软件库模拟,一次 float 乘法可能消耗 50 个周期以上。一个 9 点窗口的 SG 滤波看似只要 9 次乘加,但乘以 10 kHz 采样率,每秒钟 9 万次浮点乘法都会压在 CPU 上,ADC 采样、通信协议、显示刷新等任务都会被拖慢。整型定点化是标准的解决方案。
思路不复杂:把浮点系数乘上一个 2 的幂次放大,取整后变成整数系数,滤波结果再除以这个放大因子。下面的代码用的是 Q10 格式,即放大 2¹⁰ = 1024 倍:
/* sg_filter_q10.c —— 定点版实现,适合无 FPU 的 STM32F1 */ #include <stdint.h> #define SG_WINDOW_LEN 9 #define SG_BUFFER_SIZE 16 #define SG_COEFF_SHIFT 10 /* Q10 格式:系数放大 1024 倍 */ #define SG_COEFF_SCALE (1 << SG_COEFF_SHIFT) /* 定点系数表:浮点系数 x 1024 后四舍五入 */ static const int16_t sg_coeffs_q10[SG_WINDOW_LEN] = { -93, 62, 173, 239, 261, 239, 173, 62, -93 }; static int16_t sg_buffer[SG_BUFFER_SIZE]; static uint8_t sg_write_idx; static int32_t sg_offset; /* 用于恢复直流分量,见说明 */ static void sg_reset(void) { for (int i = 0; i < SG_BUFFER_SIZE; i++) sg_buffer[i] = 0; sg_write_idx = 0; sg_offset = 0; } int16_t sg_filter_q10(int16_t new_sample) { int32_t acc = 0; uint8_t idx; /* 移除直流偏置,防止整型溢出(见下文分析) */ int32_t centered_sample = (int32_t)new_sample - sg_offset; sg_buffer[sg_write_idx] = (int16_t)centered_sample; sg_write_idx = (sg_write_idx + 1) & (SG_BUFFER_SIZE - 1); for (int i = 0; i < SG_WINDOW_LEN; i++) { idx = (sg_write_idx - 1 - i) & (SG_BUFFER_SIZE - 1); acc += (int32_t)sg_buffer[idx] * sg_coeffs_q10[SG_WINDOW_LEN - 1 - i]; } /* 收缩回普通幅度:除以 1024,再恢复直流偏置 */ int16_t result = (int16_t)((acc >> SG_COEFF_SHIFT) + sg_offset); /* 追加减法一个较小更新量,让直流跟踪缓慢漂移(可选但推荐) */ sg_offset += (centered_sample >> 8); return result; }定点化之后有个问题浮点版不存在:整型溢出。假设你的 ADC 是 12 位,数据范围 0~4095,减去直流偏置后摆幅约 ±2048。系数最大约 0.26,与 ±2048 相乘后最大约 ±532,9 项累加最大约 ±1500,正好在 int16 范围之内。但如果你的数据摆幅接近 int16 满量程(比如外接高精度 ADC 的 ±10 V 信号),累加过程就可能超出 32767,结果变成错误的负数。代码里定义acc为 int32_t 就是为了保险,但输入缓冲sg_buffer用 int16_t 存储时,若输入信号本身接近满量程,减去直流偏置后可能溢出——所以直流偏置跟踪sg_offset必须足够贴近信号中心。你把sg_offset初始化为 2048(假设 12 位 ADC 的中心值),再配合(centered_sample >> 8)做跟踪,多数场景下就能避免这个问题。如果实在不放心,缓冲区也可以直接定义成 int32_t,代价是 RAM 翻倍——在 STM32F103 上做 16 点窗口也就多占 32 字节,无需吝啬。
4.4 计算量与内存明细:10 kHz 采样率下 CPU 占用实测参考
评估在 MCU 上跑一段算法的可行性,不能只看“逻辑对不对”,得算清楚开销多少。以 9 点窗口、3 阶的 SG 滤波为例,在 72 MHz 主频的 STM32F103 上,浮点版(软件模拟 FPU)单次滤波大约需要 350~400 个周期,定点版大约需要 90~120 个周期。在 10 kHz 采样率下,浮点版每秒消耗约 4M 周期,占 CPU 的 5.5%;定点版约 1.2M 周期,占比 1.7%。这个比例对大多数应用来说都不算压力,但如果你的系统里还跑着 FreeRTOS、LCD 驱动和 Modbus 协议栈,浮点版的 5.5% 就很可观了——随时可能成为压垮实时性的最后一根稻草。
| 对比项 | 浮点版(F4/H7 带 FPU) | 定点版(F1 无 FPU) |
|---|---|---|
| 单次滤波周期 | 约 30~50 周期 | 约 90~120 周期 |
| 内存占用(16 点缓冲) | 64 B (float) | 32 B (int16) |
| 系数表存储(9 点) | 36 B (float) | 18 B (int16) |
| 直流偏置处理 | 不需要 | 需要单独跟踪 |
| 适用型号 | STM32F4/F7/H7 系列 | STM32F1/F0 及国产替代 |
从这个表能看出来一个有意思的现象:在没有 FPU 的型号上,定点版反而比浮点版更快;而在带 FPU 的型号上,浮点版代码更简洁且不易出错,没必要强行定点化。选型的逻辑应该是:先看芯片有没有 FPU,有就上浮点版;没有则定点化,但中间累加器必须用 int32_t,这是最容易踩的坑。另外一个容易忽略的事实是:系数表本身占用的 Flash 极小——9 个 float 共 36 字节,对任何 STM32 的 Flash 容量来说都可以忽略不计,没必要在“省几十字节 Flash”上花心思优化。
5. 进阶:边界策略、参数整定与三种常见误用排查
5.1 边界处理:镜像延拓比截断和补零更实用
对流式处理来说,滤波器永远有足够的历史数据(缓冲区填满后),边界问题似乎不存在。但如果你把 SG 滤波用于一段已经采集完的离线数据——比如通过串口收到一整包波形后做后处理——边界就躲不掉了。最简单的做法(第 3 章 Python 代码里演示过)是边界点直接输出原始值,但那样会在数据两端留下明显的毛刺,和滤波后的平滑段形成突兀对比。
更实用的方案是镜像延拓:把边界附近的点按边界做镜面反射,拼出完整的窗口数据。比如信号开头是 x[0], x[1], x[2],要处理 x[0] 这个点,窗口内的数据依次取为 x[2], x[1], x[0], x[1], x[2](以 x[0] 为中心做镜像)。镜像延拓的数学含义是假设信号在边界外按反射对称延展,这比补零(假设信号跳变到 0)更接近真实信号的连续形态。实现上只多了一次索引计算:
/* 获取镜像索引:idx 从 0~N-1,映射到扩展后的虚拟索引 */ static int16_t mirror_index(int index, int length) { if (index < 0) return -index - 1; if (index >= length) return 2 * length - 1 - index; return index; }调用时,窗口内的数据对应虚拟位置 n - m + i,其中 n 是当前输出点的位置。注意这个镜像索引只对离线后处理有意义,流式场景不需要——缓冲区里有足够的历史点可以回退,边界效应只在滤波器启动的最初几个周期内短暂存在。如果你用的是第 4 章的循环缓冲区实现,启动时缓冲区清零导致前几个输出点偏小,这是正常现象,可以在初始化时连续灌入第一个真实采样值若干个周期,让缓冲区“预热”,边界效应就会消失。
5.2 参数整定流程:先看频谱,再调窗口,最后调阶数
SG 滤波有两个参数,调节顺序和调节依据都有讲究。我摸索出的可靠流程是:
第一步,对原始信号做一次 FFT,看噪声集中在哪个频段。如果噪声是宽带白噪声,SG 滤波的窗口大小主要由信号自身最高频成分决定——窗口的中心频率大约在采样率的 0.2~0.3 倍之间,超出这个范围的高频细节本来就是你要保留的信号,那就没必要用 SG。如果信号频带远低于噪声频带,窗口可以取得很大,比如半窗宽取到采样率与信号带宽比值的两倍。
第二步,固定多项式阶数为 2,改变窗口大小观察平滑效果。这一阶段不必细看参数,只需比较滤波前后的峰值衰减比例。比如峰值从 100 mV 降到 98 mV,衰减 2%,可以接受;如果降到 90 mV,说明窗口太大,调小半窗宽。这一步通常一两轮就能定位到合适区间。
第三步,把官方数从 2 升到 3,观察细节保留变化。这一步的收益往往是边际性的,但在信号存在二次或三次趋势变化时(比如心电信号中的 ST 段),效果差异可能很明显。如果升到 3 后波形看起来仍然太钝,通常不是阶数问题,而是窗口太大了——继续调小半窗宽才是正解。
一个容易被忽略的检查项:滤波后的信号是否引入延迟。SG 滤波在线流式处理中会延迟 m 个采样点,这在做控制闭环时可能影响稳定性——延迟会导致相位滞后,在 PID 调节中表现为响应变慢甚至振荡。如果你拿 SG 滤波后的数据直接喂给 PID,务必在控制器里补偿这 m 个采样点的延迟。
5.3 三种常见误用及排查方法
误用一:拿 SG 滤波去处理脉冲型噪声。SG 滤波对高斯白噪声和均匀分布噪声效果良好,但对幅值很大的尖峰脉冲(比如电机启动瞬间的电磁干扰)几乎没有抑制作用——脉冲宽度远小于窗口时,滤波后的输出仍会保留一个明显的尖峰,只是被压低了一些。排查方法:看滤波后的波形里是否仍有窄幅大于正常噪声的孤立尖峰。通俗的处理是先用中值滤波剔除脉冲,再用 SG 平滑去噪。中值滤波加 SG 滤波的组合在电机电流采样中非常实用。
误用二:窗口长度取偶数。SG 系数的对称性依赖于窗口以中心点为参考,偶数窗口没有严格定义的中心点,SciPy 的 savgol_coeffs 会直接报错,但如果你手动抄系数或者从别的代码库拿到偶数窗口的系数,滤波结果会出现相位失真。排查方法:检查系数表是否严格对称。第 3 章提到过:系数和应为 1,中心对称。两条校验同时满足,才能确认系数表可靠。
误用三:过于信任系数表的“外部来源”。网上能找到很多 SG 系数表,但来源良莠不齐,有的是浮点版,有的是定点版但没写缩放位数,有的多项式阶数和窗口长度组合本身就是非法的(比如 5 点窗口配 5 阶多项式,正规方程奇异)。拿到任何系数表,先做 2.1 节的两条校验,再用第 3 章的 Python 脚本生成一组同样参数的系数对比一遍。数值完全一致才可用,有差异就必须查清楚差异来源,不能直接拿进工程。
最后一个验证技巧,适用于任何实现版本:构造一段阶跃信号(前 100 个点全 0,后 100 个点全 4095),送入滤波器,观察输出的上升沿形态。SG 滤波对阶跃的响应不是单调的——会有轻微的过冲和振铃,这与阶数和窗口直接相关。如果你的实现输出在阶跃处出现不对称的响应,说明代码里的索引方向或系数顺序搞反了;如果输出稳定后不等于输入值(应等于 4095),说明直流偏置处理有误。这一个测试能同时暴露 SG 滤波移植中最常见的三类错误,值得写进你的固件自检流程。
参考文献
- Savitzky, A., Golay, M. J. E. (1964). Smoothing and Differentiation of Data by Simplified Least Squares Procedures. Analytical Chemistry, 36(8), 1627–1639.
- Press, W. H., Teukolsky, S. A., Vetterling, W. T., Flannery, B. P. (2007). Numerical Recipes: The Art of Scientific Computing (3rd ed.). Cambridge University Press.
- Schafer, R. W. (2011). What Is a Savitzky-Golay Filter? IEEE Signal Processing Magazine, 28(4), 111–117.
- STMicroelectronics. (2023). STM32F103x8/xB Datasheet. Available at: ST.com
- Virtanen, P., et al. (2020). SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python. Nature Methods, 17, 261–272.
本文还有配套的精品资源,点击获取