嵌入式信号处理:Savitzky-Golay滤波在STM32上的原理与实现
2026/9/12 11:37:36 网站建设 项目流程

简介: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 滤波移植中最常见的三类错误,值得写进你的固件自检流程。

参考文献

  1. Savitzky, A., Golay, M. J. E. (1964). Smoothing and Differentiation of Data by Simplified Least Squares Procedures. Analytical Chemistry, 36(8), 1627–1639.
  2. 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.
  3. Schafer, R. W. (2011). What Is a Savitzky-Golay Filter? IEEE Signal Processing Magazine, 28(4), 111–117.
  4. STMicroelectronics. (2023). STM32F103x8/xB Datasheet. Available at: ST.com
  5. Virtanen, P., et al. (2020). SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python. Nature Methods, 17, 261–272.

本文还有配套的精品资源,点击获取

需要专业的网站建设服务?

联系我们获取免费的网站建设咨询和方案报价,让我们帮助您实现业务目标

立即咨询