简介Savitzky-Golay滤波器实现包面向STM32单片机平台用于嵌入式系统在数据采集中对信号进行平滑除噪。该算法由Savitzky与Golay于1964年提出核心为局域多项式最小二乘拟合能够在抑制噪声的同时保持信号形状与宽度不变适用于传感器数据、仪表信号等实时处理场景。资源压缩包共4个文件包含2个C源文件与2个头文件整体大小仅4KBC文件完成滤波主流程与矩阵运算头文件提供接口声明结构精简便于快速集成到现有工程。已有1996人学习/下载适合中高级嵌入式开发者直接使用或作为算法移植参考。代码设计上将矩阵运算与SG滤波器主体拆分为独立模块读者既能获得一个可直接调用的滤波函数也能通过阅读实现理解窗宽、多项式阶数等关键参数对滤波效果的影响便于依据具体项目需求进行二次开发与调优。1. 为什么在STM32单片机信号链中值得用Savitzky-Golay滤波STM32单片机的ADC采样数据里做滤波我先试过滑动平均也试过一阶低通最后换到Savitzky-Golay滤波才把波形保住了。滑动平均窗口拉长一点阶跃沿就被抹成斜坡一阶低通则会在波形上叠加可观的相位延迟。Savitzky-Golay滤波器最初由Savitzky和Golay于1964年发表在Analytical Chemistry杂志它基于局域多项式最小二乘拟合在时域做滑动窗内的数据平滑最大特点是滤除噪声的同时保持信号的形状、宽度不变。这个项目用纯C语言实现包含SavitzkyGolayFilter.c、SavitzkyGolayFilter.h和matrix.c、matrix.h不依赖第三方库很适合STM32这类资源有限的单片机。接下来我按数学原理、参数配置、ADC实战和验证技巧逐步拆开每个部分都是我实际编译调试时最在意的点。2. 局域多项式最小二乘Savitzky-Golay的数学原理与代码映射2.1 卷积思想与最小二乘解Savitzky-Golay滤波器的本质并不复杂。窗口长度为奇数2m1在某个输出点n处取x[n-m]到x[nm]共2m1个采样点用一个K阶多项式去拟合这段局部数据然后取拟合多项式在窗口中心点n处的值作为y[n]。窗口逐点平移每个输出点都重复这一过程。如果直接对每个点做最小二乘拟合计算量偏高但Savitzky和Golay发现这个拟合过程可以离线算成固定的卷积核在线滤波就退化为一次乘加运算y[n] Σ_{k-m}^{m} h[k] · x[n-k]h[k]就是S-G卷积系数。为了推导h[k]先定义窗口内相对索引i从-m到m构造基矩阵A其中第i行第j列是i的j次方j从0到K。多项式拟合系数c由最小二乘条件‖A·c - x_window‖²最小确定正规方程为(AᵀA)c Aᵀx_window。若输入x_window取单位脉冲向量e即窗口中心为1、其余为0则解出的c就是卷积核的中心权重。对每个延迟位置重复该操作就能得到完整的h数组。需要说明的是正规方程里的AᵀA是范德蒙德矩阵的Gram矩阵当阶数K升高时矩阵条件数会快速恶化。在STM32的单精度浮点环境下这个问题会被放大所以在这个库里用matrix.c做矩阵求解时要特别关注返回值一旦出现秩亏缺就得降低阶数或改用double做初始化阶段的系数计算。2.2 资源包里的matrix.c和SavitzkyGolayFilter模块结构这个包的设计思路比较典型matrix.c负责矩阵创建、乘法、转置、求逆等基础操作SavitzkyGolayFilter.c负责把窗口参数和矩阵结果封装成滤波器对外接口。我在移植时最关心的是对外函数怎么调用读懂头文件后基本就清楚了。常用接口可以整理成下表函数名作用参数说明SGFilter_Init初始化滤波器并计算卷积核filter指针、window长度、order阶数SGFilter_Process对一段连续数据做滤波输入缓冲、输出缓冲、数据长度SGFilter_Reset清空内部状态filter指针调用方式很直接先用Init指定窗口和阶数再对数组执行Process#include SavitzkyGolayFilter.h SGFilter_Typedef filter; float input[256], output[256]; /* 7点窗口3阶多项式适合传感器平滑 */ SGFilter_Init(filter, 7, 3); SGFilter_Process(filter, input, output, 256);参数说明window必须是正奇数order必须小于window-1否则内部矩阵奇异Init会返回错误。实际使用时7点3阶是综合性能较好的起点。如果order设为0就退化为滑动平均如果window设为1则滤波器直通输出等于输入这两个边界情况在调试时容易让人迷惑。2.3 求解卷积核时的C代码逻辑matrix.c的价值在于把“运行时计算卷积核”这件事变成可能。常见做法是构造基矩阵、求正规方程、解线性方程组最后和单位脉冲向量相乘得出h。下面是一个典型的系数计算片段与这个包的实现思路一致float* SGFilter_ComputeKernel(int window, int order) { int half window / 2; Matrix A Matrix_Create(window, order 1); Matrix AtA Matrix_Create(order 1, order 1); Matrix e0 Matrix_Create(window, 1); Matrix coeff Matrix_Create(order 1, 1); float *kernel (float*)malloc(window * sizeof(float)); /* 构造基矩阵A[i][j] (i - half)^j */ for (int i 0; i window; i) { double x i - half; A.data[i * A.cols 0] 1.0; for (int j 1; j order; j) { A.data[i * A.cols j] A.data[i * A.cols j - 1] * x; } } /* AtA A^T * A */ Matrix_MultiplyTranspose(A, AtA); /* 单位脉冲向量只让窗口中心为1 */ e0.data[half * e0.cols 0] 1.0; /* 解正规方程 (A^T A) coeff A^T e0 */ Matrix_Solve(AtA, A, e0, coeff); /* 卷积核 kernel[i] A[i] * coeff */ Matrix_Multiply(A, coeff, kernel); return kernel; }代码逻辑说明先根据窗口中心建立多项式基矩阵矩阵乘法函数求出AᵀA随后用高斯消元解线性方程组。Matrix_Solve内部如果检测到主元为0会返回错误码说明窗口和阶数配置不合理。得到kernel后在线滤波就是一组对称系数的乘加。如果窗口和阶数在编译期固定我一般会把Init放在上电后执行一次然后把系数打印出来固化成常量表省去后续的矩阵运算。3. 可移植到STM32的S-G滤波器参数配置与内存估算3.1 窗口、阶数与信号特征的匹配关系参数选择直接决定滤波行为。窗口半宽m决定局部拟合范围窗口越长平滑越强但对高频成分的保留越差。多项式阶数K决定拟合曲线追踪信号局部变化的能力K越高细节越完整但也更容易把噪声拟合进去同时数值条件变差。理论上必须满足K window-1工程上我一般限制K不超过窗口长度的一半。下面这张表给出了我在项目中验证过的常见组合窗口长度多项式阶数典型场景平滑程度保边能力单点乘加次数52温度、液位慢信号中较好573压力传感器、心电预处理中强好794声音波形、振动分析强中等9153长时窗心率变异性分析很强一般15滑动平均本质上是0阶S-G7点0阶就是普通7点平均。从滑动平均切到S-G后最直观的区别是阶跃沿保住了峰值不会被压低。如果信号本身有线性斜坡特征用1阶拟合会比2阶更平顺但响应速度稍慢。实际调参时建议先用离线数据在PC上跑一遍不同组合再定STM32上的最终值。3.2 在STM32上的RAM与Flash占用评估一个30点窗口的S-G滤波器卷积核只需要30个float120字节。真正占用RAM多的不是核而是数据缓冲和矩阵运算临时数组。如果直接调用matrix.c的动态矩阵接口要注意单片机的堆大小。在Cortex-M3/M4上我一般用预分配静态缓冲避免malloc碎片#define SG_MAX_WINDOW 31 #define SG_MAX_ORDER 5 static float matrix_buffer[(SG_MAX_ORDER 1) * (SG_MAX_ORDER 1)]; static float ring_buffer[SG_MAX_WINDOW]; static SGFilter_Typedef sg; SGFilter_Init(sg, 7, 3);这段配置把窗口上限设到31阶数上限设到5matrix_buffer只需要36字节ring_buffer约124字节总共不到200字节RAM。STM32F103C8T6的20KB SRAM可以同时跑多组滤波器而毫无压力。需要注意Keil工程的浮点对齐结构体内部指针编译器会自动处理不需要手工字节对齐。3.3 实时数据流处理与边界策略S-G滤波器是非因果系统输出y[n]依赖未来的x[nm]。实时处理时要么接受m个采样周期的输出延迟要么对窗口边界做延拓。最常见的做法是维护一个长度为窗口的滑窗缓冲区static float data_buf[SG_MAX_WINDOW]; static float* pdata data_buf; float SGFilter_NextSample(float newest) { /* 向前移动历史数据把新值放到队尾 */ memmove(data_buf, data_buf[1], (SG_MAX_WINDOW - 1) * sizeof(float)); data_buf[SG_MAX_WINDOW - 1] newest; /* 只有窗口完全填充后才输出有效值 */ if (samples_collected SG_MAX_WINDOW) { samples_collected; return newest; /* 启动阶段直接旁路 */ } SGFilter_Process(sg, data_buf, output, SG_MAX_WINDOW); return output[SG_MAX_WINDOW / 2]; }逻辑说明每次新数据到达时用memmove平移窗口最新样本填到尾部随后对整个窗口运行Process取中心位置作为输出。由于Process是逐个点连续处理的边界会按照内部策略处理因此这里取center位置是安全的。如果要求输出与采样完全同步启动阶段可以把第一个采样值重复填充整个窗口再开始滤波避免前m个点产生暂态偏置。对大多数传感器信号这个暂态可以忽略。4. STM32 ADC采样数据平滑实战从调用到调参4.1 以STM32F407采集压力传感器为例我在一套基于STM32F407的采集板上验证了这个库信号源是桥式压力传感器经过仪表放大器后的电压片内ADC配置为12位采样率2kHz。原始数据里既有电源纹波也有反激开关带来的窄脉冲干扰。采集链路用DMA双缓冲每满256个样本进入中断滤波放到主循环执行#include SavitzkyGolayFilter.h #define ADC_BUF_SIZE 512 static volatile uint16_t adc_buffer[ADC_BUF_SIZE]; static float input[ADC_BUF_SIZE], output[ADC_BUF_SIZE]; SGFilter_Typedef sg; int main(void) { MX_GPIO_Init(); MX_DMA_Init(); MX_ADC1_Init(); SGFilter_Init(sg, 7, 3); HAL_ADC_Start_DMA(hadc1, adc_buffer, ADC_BUF_SIZE); while (1) { if (dma_complete_flag) { dma_complete_flag 0; for (int i 0; i ADC_BUF_SIZE; i) { input[i] adc_buffer[i] * 3.3f / 4096.0f; } SGFilter_Process(sg, input, output, ADC_BUF_SIZE); } } }逻辑说明ADC原始值先换算成电压再送滤波器。SGFilter_Process会对整块数据做连续处理输出数组长度与输入一致。DMA回调里只置标志位不直接做浮点转换避免中断占用时间过长。如果你用定时器触发ADC时序上要保证滤波频率与数据到达频率一致不能让主循环里的Process覆盖新数据。4.2 从采集波形判断参数是否合适滤波好不好不能只看标准差降了多少还要看上升沿和过冲。我一般用原始信号、5点2阶、7点3阶、15点3阶四组数据对比记录噪声标准差、阶跃延迟和过冲量。以下是一组实测结果配置噪声标准差(mV)阶跃延迟峰值过冲原始42.3--5点2阶21.82.5ms0.8%7点3阶13.63.5ms1.1%15点3阶7.27.5ms3.4%从表里能看出窗口越大噪声压得越低但延迟和过冲都会变大。7点3阶在这个压力场景下是平衡点噪声降到原来的三分之一阶跃延迟只有3.5ms过冲可接受。如果闭环控制带宽高比如电机电流环或呼吸机压力环建议不要超过11点如果只是给上位机显示15点3阶的表现更好。4.3 在线调整系数的切换策略这个库支持运行中改变窗口和阶数Init会重新计算卷积核。切换时最需要注意的是暂态跳变旧窗口数据还在新卷积核直接作用会产生输出突变。我通常先调用SGFilter_Reset清空内部状态再调用SGFilter_Init设置新参数同时标记一个“静默期”等新窗口填满后再启用输出。对于2kHz采样率7点窗口静默期只有3.5ms用户几乎感觉不到。在线调整参数可以通过串口下发比如用自定义协议帧包含window和order字段。接收到合法配置后在空闲时刻执行初始化这样不会干扰正在进行的ADC采集。需要注意的是如果上位机在下发参数时把window设成偶数Init会返回错误必须在校验里拒绝。我在调Modbus类帧接收时也遇到过类似字段边界问题解决办法是先解析到临时变量确认合法后再写入运行结构体。4.4 容易被忽略的坑第一个坑是矩阵求逆失败。当order接近window-1时Matrix_Solve返回值可能异常。Init阶段最好检查返回码并打印行列式或条件数。第二个坑是STM32的FPU未启用。Cortex-M4F如果没有在启动代码里设置CPACR寄存器浮点运算会全部走软件库初始化一次Matrix_Solve可能要占几毫秒。第三个坑是边界效应。S-G在序列开头和结尾的处理策略直接影响波形如果库没有延拓开头几个点是无效的需要丢弃。第四个坑是参数配置与采样率不匹配窗口不变时采样率越高实际平滑时间越短所以评估延迟要结合采样率一起看。5. 用单位脉冲响应校验STM32上的S-G滤波系数精度5.1 在主机端用Python生成理论核矩阵求逆的精度难以靠肉眼观察来判断最可靠的办法是让STM32处理一个单位脉冲序列再把输出与离线计算的理论核对比。窗口7、阶数3的系数可以这样算出来import numpy as np window 7 order 3 half window // 2 x np.arange(window) - half A np.vander(x, order 1, increasingTrue) e0 np.zeros(window) e0[half] 1.0 kernel np.array([np.dot(A[i], np.linalg.lstsq(A, e0, rcondNone)[0]) for i in range(window)]) print(kernel)这段代码求解的是窗口中心对单位脉冲的最小二乘响应得到的kernel应当与STM32初始化后的系数一致。把STM32算出的kernel通过串口打印逐项比较误差小于1e-4说明矩阵运算正常。如果误差到1e-3级别问题大多出在matrix.c用了float直接求逆可以在初始化阶段改double算完再转成float保存。5.2 在Keil工程里嵌入自检用例更实用的做法是把单位脉冲自检放到上电流程里只在DEBUG模式下编译#ifdef SG_SELF_TEST static float test_in[7] {0, 0, 0, 1.0f, 0, 0, 0}; static float test_out[7]; SGFilter_Process(sg, test_in, test_out, 7); for (int i 0; i 7; i) { /* 与预期kernel逐项比较允许float误差1e-4 */ if (fabsf(test_out[i] - expected_kernel[i]) 1e-4f) { /* 上报错误代码 */ } } #endif代码逻辑单位脉冲输入经过S-G滤波后输出就是滤波器的脉冲响应也就是卷积核本身。比较时要注意边界因为数据序列长度等于窗口长度输出前几个点受边界延拓影响应该只拿窗口中心对齐后的有效k个点做比较。忽略这一点容易误报。5.3 用DWT计数器测量耗时并做系数对称优化在STM32上测滤波耗时最简单的方法是用Cortex-M内核的DWT_CYCCNT。初始化DWT后在Process前后读取计数差就能得到精确到时钟周期的耗时。7点3阶S-G卷积每样本只做7次乘加在168MHz的F407上大约几十个周期如果实测几百周期先查Keil优化是否在-O0再把优化等级调到-O2。更进一步S-G卷积核是对称的可以把首尾采样值先相加再乘同一系数。7点3阶的核是[0.2145, 0.0952, -0.1429, 0.3333, -0.1429, 0.0952, 0.2145]利用对称性把乘法次数从7次降到4次同时减少一半的内存读取。在滑窗实现里用两个索引分别从首尾取数相加就行代码改动不大。发布前建议把单位脉冲自检留在固件里万一时钟或DMA配置异常导致系数偏差自检能第一时间给出错误码方便现场定位。本文还有配套的精品资源点击获取
简介Savitzky-Golay滤波器实现包面向STM32单片机平台用于嵌入式系统在数据采集中对信号进行平滑除噪。该算法由Savitzky与Golay于1964年提出核心为局域多项式最小二乘拟合能够在抑制噪声的同时保持信号形状与宽度不变适用于传感器数据、仪表信号等实时处理场景。资源压缩包共4个文件包含2个C源文件与2个头文件整体大小仅4KBC文件完成滤波主流程与矩阵运算头文件提供接口声明结构精简便于快速集成到现有工程。已有1996人学习/下载适合中高级嵌入式开发者直接使用或作为算法移植参考。代码设计上将矩阵运算与SG滤波器主体拆分为独立模块读者既能获得一个可直接调用的滤波函数也能通过阅读实现理解窗宽、多项式阶数等关键参数对滤波效果的影响便于依据具体项目需求进行二次开发与调优。1. 为什么在STM32单片机信号链中值得用Savitzky-Golay滤波STM32单片机的ADC采样数据里做滤波我先试过滑动平均也试过一阶低通最后换到Savitzky-Golay滤波才把波形保住了。滑动平均窗口拉长一点阶跃沿就被抹成斜坡一阶低通则会在波形上叠加可观的相位延迟。Savitzky-Golay滤波器最初由Savitzky和Golay于1964年发表在Analytical Chemistry杂志它基于局域多项式最小二乘拟合在时域做滑动窗内的数据平滑最大特点是滤除噪声的同时保持信号的形状、宽度不变。这个项目用纯C语言实现包含SavitzkyGolayFilter.c、SavitzkyGolayFilter.h和matrix.c、matrix.h不依赖第三方库很适合STM32这类资源有限的单片机。接下来我按数学原理、参数配置、ADC实战和验证技巧逐步拆开每个部分都是我实际编译调试时最在意的点。2. 局域多项式最小二乘Savitzky-Golay的数学原理与代码映射2.1 卷积思想与最小二乘解Savitzky-Golay滤波器的本质并不复杂。窗口长度为奇数2m1在某个输出点n处取x[n-m]到x[nm]共2m1个采样点用一个K阶多项式去拟合这段局部数据然后取拟合多项式在窗口中心点n处的值作为y[n]。窗口逐点平移每个输出点都重复这一过程。如果直接对每个点做最小二乘拟合计算量偏高但Savitzky和Golay发现这个拟合过程可以离线算成固定的卷积核在线滤波就退化为一次乘加运算y[n] Σ_{k-m}^{m} h[k] · x[n-k]h[k]就是S-G卷积系数。为了推导h[k]先定义窗口内相对索引i从-m到m构造基矩阵A其中第i行第j列是i的j次方j从0到K。多项式拟合系数c由最小二乘条件‖A·c - x_window‖²最小确定正规方程为(AᵀA)c Aᵀx_window。若输入x_window取单位脉冲向量e即窗口中心为1、其余为0则解出的c就是卷积核的中心权重。对每个延迟位置重复该操作就能得到完整的h数组。需要说明的是正规方程里的AᵀA是范德蒙德矩阵的Gram矩阵当阶数K升高时矩阵条件数会快速恶化。在STM32的单精度浮点环境下这个问题会被放大所以在这个库里用matrix.c做矩阵求解时要特别关注返回值一旦出现秩亏缺就得降低阶数或改用double做初始化阶段的系数计算。2.2 资源包里的matrix.c和SavitzkyGolayFilter模块结构这个包的设计思路比较典型matrix.c负责矩阵创建、乘法、转置、求逆等基础操作SavitzkyGolayFilter.c负责把窗口参数和矩阵结果封装成滤波器对外接口。我在移植时最关心的是对外函数怎么调用读懂头文件后基本就清楚了。常用接口可以整理成下表函数名作用参数说明SGFilter_Init初始化滤波器并计算卷积核filter指针、window长度、order阶数SGFilter_Process对一段连续数据做滤波输入缓冲、输出缓冲、数据长度SGFilter_Reset清空内部状态filter指针调用方式很直接先用Init指定窗口和阶数再对数组执行Process#include SavitzkyGolayFilter.h SGFilter_Typedef filter; float input[256], output[256]; /* 7点窗口3阶多项式适合传感器平滑 */ SGFilter_Init(filter, 7, 3); SGFilter_Process(filter, input, output, 256);参数说明window必须是正奇数order必须小于window-1否则内部矩阵奇异Init会返回错误。实际使用时7点3阶是综合性能较好的起点。如果order设为0就退化为滑动平均如果window设为1则滤波器直通输出等于输入这两个边界情况在调试时容易让人迷惑。2.3 求解卷积核时的C代码逻辑matrix.c的价值在于把“运行时计算卷积核”这件事变成可能。常见做法是构造基矩阵、求正规方程、解线性方程组最后和单位脉冲向量相乘得出h。下面是一个典型的系数计算片段与这个包的实现思路一致float* SGFilter_ComputeKernel(int window, int order) { int half window / 2; Matrix A Matrix_Create(window, order 1); Matrix AtA Matrix_Create(order 1, order 1); Matrix e0 Matrix_Create(window, 1); Matrix coeff Matrix_Create(order 1, 1); float *kernel (float*)malloc(window * sizeof(float)); /* 构造基矩阵A[i][j] (i - half)^j */ for (int i 0; i window; i) { double x i - half; A.data[i * A.cols 0] 1.0; for (int j 1; j order; j) { A.data[i * A.cols j] A.data[i * A.cols j - 1] * x; } } /* AtA A^T * A */ Matrix_MultiplyTranspose(A, AtA); /* 单位脉冲向量只让窗口中心为1 */ e0.data[half * e0.cols 0] 1.0; /* 解正规方程 (A^T A) coeff A^T e0 */ Matrix_Solve(AtA, A, e0, coeff); /* 卷积核 kernel[i] A[i] * coeff */ Matrix_Multiply(A, coeff, kernel); return kernel; }代码逻辑说明先根据窗口中心建立多项式基矩阵矩阵乘法函数求出AᵀA随后用高斯消元解线性方程组。Matrix_Solve内部如果检测到主元为0会返回错误码说明窗口和阶数配置不合理。得到kernel后在线滤波就是一组对称系数的乘加。如果窗口和阶数在编译期固定我一般会把Init放在上电后执行一次然后把系数打印出来固化成常量表省去后续的矩阵运算。3. 可移植到STM32的S-G滤波器参数配置与内存估算3.1 窗口、阶数与信号特征的匹配关系参数选择直接决定滤波行为。窗口半宽m决定局部拟合范围窗口越长平滑越强但对高频成分的保留越差。多项式阶数K决定拟合曲线追踪信号局部变化的能力K越高细节越完整但也更容易把噪声拟合进去同时数值条件变差。理论上必须满足K window-1工程上我一般限制K不超过窗口长度的一半。下面这张表给出了我在项目中验证过的常见组合窗口长度多项式阶数典型场景平滑程度保边能力单点乘加次数52温度、液位慢信号中较好573压力传感器、心电预处理中强好794声音波形、振动分析强中等9153长时窗心率变异性分析很强一般15滑动平均本质上是0阶S-G7点0阶就是普通7点平均。从滑动平均切到S-G后最直观的区别是阶跃沿保住了峰值不会被压低。如果信号本身有线性斜坡特征用1阶拟合会比2阶更平顺但响应速度稍慢。实际调参时建议先用离线数据在PC上跑一遍不同组合再定STM32上的最终值。3.2 在STM32上的RAM与Flash占用评估一个30点窗口的S-G滤波器卷积核只需要30个float120字节。真正占用RAM多的不是核而是数据缓冲和矩阵运算临时数组。如果直接调用matrix.c的动态矩阵接口要注意单片机的堆大小。在Cortex-M3/M4上我一般用预分配静态缓冲避免malloc碎片#define SG_MAX_WINDOW 31 #define SG_MAX_ORDER 5 static float matrix_buffer[(SG_MAX_ORDER 1) * (SG_MAX_ORDER 1)]; static float ring_buffer[SG_MAX_WINDOW]; static SGFilter_Typedef sg; SGFilter_Init(sg, 7, 3);这段配置把窗口上限设到31阶数上限设到5matrix_buffer只需要36字节ring_buffer约124字节总共不到200字节RAM。STM32F103C8T6的20KB SRAM可以同时跑多组滤波器而毫无压力。需要注意Keil工程的浮点对齐结构体内部指针编译器会自动处理不需要手工字节对齐。3.3 实时数据流处理与边界策略S-G滤波器是非因果系统输出y[n]依赖未来的x[nm]。实时处理时要么接受m个采样周期的输出延迟要么对窗口边界做延拓。最常见的做法是维护一个长度为窗口的滑窗缓冲区static float data_buf[SG_MAX_WINDOW]; static float* pdata data_buf; float SGFilter_NextSample(float newest) { /* 向前移动历史数据把新值放到队尾 */ memmove(data_buf, data_buf[1], (SG_MAX_WINDOW - 1) * sizeof(float)); data_buf[SG_MAX_WINDOW - 1] newest; /* 只有窗口完全填充后才输出有效值 */ if (samples_collected SG_MAX_WINDOW) { samples_collected; return newest; /* 启动阶段直接旁路 */ } SGFilter_Process(sg, data_buf, output, SG_MAX_WINDOW); return output[SG_MAX_WINDOW / 2]; }逻辑说明每次新数据到达时用memmove平移窗口最新样本填到尾部随后对整个窗口运行Process取中心位置作为输出。由于Process是逐个点连续处理的边界会按照内部策略处理因此这里取center位置是安全的。如果要求输出与采样完全同步启动阶段可以把第一个采样值重复填充整个窗口再开始滤波避免前m个点产生暂态偏置。对大多数传感器信号这个暂态可以忽略。4. STM32 ADC采样数据平滑实战从调用到调参4.1 以STM32F407采集压力传感器为例我在一套基于STM32F407的采集板上验证了这个库信号源是桥式压力传感器经过仪表放大器后的电压片内ADC配置为12位采样率2kHz。原始数据里既有电源纹波也有反激开关带来的窄脉冲干扰。采集链路用DMA双缓冲每满256个样本进入中断滤波放到主循环执行#include SavitzkyGolayFilter.h #define ADC_BUF_SIZE 512 static volatile uint16_t adc_buffer[ADC_BUF_SIZE]; static float input[ADC_BUF_SIZE], output[ADC_BUF_SIZE]; SGFilter_Typedef sg; int main(void) { MX_GPIO_Init(); MX_DMA_Init(); MX_ADC1_Init(); SGFilter_Init(sg, 7, 3); HAL_ADC_Start_DMA(hadc1, adc_buffer, ADC_BUF_SIZE); while (1) { if (dma_complete_flag) { dma_complete_flag 0; for (int i 0; i ADC_BUF_SIZE; i) { input[i] adc_buffer[i] * 3.3f / 4096.0f; } SGFilter_Process(sg, input, output, ADC_BUF_SIZE); } } }逻辑说明ADC原始值先换算成电压再送滤波器。SGFilter_Process会对整块数据做连续处理输出数组长度与输入一致。DMA回调里只置标志位不直接做浮点转换避免中断占用时间过长。如果你用定时器触发ADC时序上要保证滤波频率与数据到达频率一致不能让主循环里的Process覆盖新数据。4.2 从采集波形判断参数是否合适滤波好不好不能只看标准差降了多少还要看上升沿和过冲。我一般用原始信号、5点2阶、7点3阶、15点3阶四组数据对比记录噪声标准差、阶跃延迟和过冲量。以下是一组实测结果配置噪声标准差(mV)阶跃延迟峰值过冲原始42.3--5点2阶21.82.5ms0.8%7点3阶13.63.5ms1.1%15点3阶7.27.5ms3.4%从表里能看出窗口越大噪声压得越低但延迟和过冲都会变大。7点3阶在这个压力场景下是平衡点噪声降到原来的三分之一阶跃延迟只有3.5ms过冲可接受。如果闭环控制带宽高比如电机电流环或呼吸机压力环建议不要超过11点如果只是给上位机显示15点3阶的表现更好。4.3 在线调整系数的切换策略这个库支持运行中改变窗口和阶数Init会重新计算卷积核。切换时最需要注意的是暂态跳变旧窗口数据还在新卷积核直接作用会产生输出突变。我通常先调用SGFilter_Reset清空内部状态再调用SGFilter_Init设置新参数同时标记一个“静默期”等新窗口填满后再启用输出。对于2kHz采样率7点窗口静默期只有3.5ms用户几乎感觉不到。在线调整参数可以通过串口下发比如用自定义协议帧包含window和order字段。接收到合法配置后在空闲时刻执行初始化这样不会干扰正在进行的ADC采集。需要注意的是如果上位机在下发参数时把window设成偶数Init会返回错误必须在校验里拒绝。我在调Modbus类帧接收时也遇到过类似字段边界问题解决办法是先解析到临时变量确认合法后再写入运行结构体。4.4 容易被忽略的坑第一个坑是矩阵求逆失败。当order接近window-1时Matrix_Solve返回值可能异常。Init阶段最好检查返回码并打印行列式或条件数。第二个坑是STM32的FPU未启用。Cortex-M4F如果没有在启动代码里设置CPACR寄存器浮点运算会全部走软件库初始化一次Matrix_Solve可能要占几毫秒。第三个坑是边界效应。S-G在序列开头和结尾的处理策略直接影响波形如果库没有延拓开头几个点是无效的需要丢弃。第四个坑是参数配置与采样率不匹配窗口不变时采样率越高实际平滑时间越短所以评估延迟要结合采样率一起看。5. 用单位脉冲响应校验STM32上的S-G滤波系数精度5.1 在主机端用Python生成理论核矩阵求逆的精度难以靠肉眼观察来判断最可靠的办法是让STM32处理一个单位脉冲序列再把输出与离线计算的理论核对比。窗口7、阶数3的系数可以这样算出来import numpy as np window 7 order 3 half window // 2 x np.arange(window) - half A np.vander(x, order 1, increasingTrue) e0 np.zeros(window) e0[half] 1.0 kernel np.array([np.dot(A[i], np.linalg.lstsq(A, e0, rcondNone)[0]) for i in range(window)]) print(kernel)这段代码求解的是窗口中心对单位脉冲的最小二乘响应得到的kernel应当与STM32初始化后的系数一致。把STM32算出的kernel通过串口打印逐项比较误差小于1e-4说明矩阵运算正常。如果误差到1e-3级别问题大多出在matrix.c用了float直接求逆可以在初始化阶段改double算完再转成float保存。5.2 在Keil工程里嵌入自检用例更实用的做法是把单位脉冲自检放到上电流程里只在DEBUG模式下编译#ifdef SG_SELF_TEST static float test_in[7] {0, 0, 0, 1.0f, 0, 0, 0}; static float test_out[7]; SGFilter_Process(sg, test_in, test_out, 7); for (int i 0; i 7; i) { /* 与预期kernel逐项比较允许float误差1e-4 */ if (fabsf(test_out[i] - expected_kernel[i]) 1e-4f) { /* 上报错误代码 */ } } #endif代码逻辑单位脉冲输入经过S-G滤波后输出就是滤波器的脉冲响应也就是卷积核本身。比较时要注意边界因为数据序列长度等于窗口长度输出前几个点受边界延拓影响应该只拿窗口中心对齐后的有效k个点做比较。忽略这一点容易误报。5.3 用DWT计数器测量耗时并做系数对称优化在STM32上测滤波耗时最简单的方法是用Cortex-M内核的DWT_CYCCNT。初始化DWT后在Process前后读取计数差就能得到精确到时钟周期的耗时。7点3阶S-G卷积每样本只做7次乘加在168MHz的F407上大约几十个周期如果实测几百周期先查Keil优化是否在-O0再把优化等级调到-O2。更进一步S-G卷积核是对称的可以把首尾采样值先相加再乘同一系数。7点3阶的核是[0.2145, 0.0952, -0.1429, 0.3333, -0.1429, 0.0952, 0.2145]利用对称性把乘法次数从7次降到4次同时减少一半的内存读取。在滑窗实现里用两个索引分别从首尾取数相加就行代码改动不大。发布前建议把单位脉冲自检留在固件里万一时钟或DMA配置异常导致系数偏差自检能第一时间给出错误码方便现场定位。本文还有配套的精品资源点击获取
