高精度保形滤波器Savitzky-Golay的数学原理、Python实现与工程应用

面向信号处理的特征保持平滑技术

在数据分析领域，信号处理中的噪声问题始终是一个重要议题。无论是实验数据、金融时间序列还是其他形式的信号处理，噪声都会干扰目标模式和趋势的识别。尽管存在多种降噪方法，但在处理短时信号时，算法的性能往往比执行效率更为重要。在众多方法中Savitzky-Golay滤波器因其独特的特征保持能力而脱颖而出。

Savitzky-Golay滤波器由Abraham Savitzky和Marcel J. E. Golay于1964年提出，是一种应用广泛的数字滤波器，可用于数据平滑和微分运算。与传统的中值滤波或均值滤波等容易造成信号特征损失的方法相比，Savitzky-Golay滤波器能够在实现信号平滑的同时保持原始信号的关键特征。这一特性使其在信号形状和特征保持要求较高的应用场景中具有显著优势。

本文将系统地介绍Savitzky-Golay滤波器的原理、实现和应用。我们将从基本原理出发，通过数学推导和直观解释，深入理解该滤波器的工作机制。同时将结合Python实现，展示其在实际应用中的效果。

Savitzky-Golay滤波器原理

Savitzky-Golay滤波器是一种基于局部多项式回归的数字滤波器，其核心是通过线性最小二乘法将低阶多项式拟合到相邻数据点的滑动窗口中。该方法的主要优势在于能够在降低噪声的同时保持信号的高阶矩，这意味着信号的峰值、谷值等特征可以得到较好的保持。

滤波器的工作过程可以概括为：在信号序列上滑动固定大小的窗口，对窗口内的数据点进行多项式拟合。窗口大小和多项式阶数是该算法的两个关键参数。算法在每个窗口位置计算多项式在中心点处的值，将其作为该点的滤波输出。通过对每个数据点重复此过程，最终得到完整的滤波信号。

数学原理

多项式拟合

Savitzky-Golay滤波器的核心是局部多项式拟合。设数据序列为(xi, yi)，其中i∈[1, N]，目标是用p阶多项式对局部数据进行拟合。

多项式表达式为：

对于中心位于x_k的窗口，需要确定系数向量[a0, a1, …, ap]，使得多项式能最佳拟合窗口内的数据点。这个优化问题可以通过最小化均方误差来解决：

其中，2m+1表示窗口大小，窗口中心为点x_k。

拟合实例

为了说明算法的具体实现过程，我们考虑一个简单的例子：窗口大小为5（即m=2）的2阶多项式拟合。

假设窗口内的数据点为：

采用2阶多项式进行拟合：

最小化误差函数：

求解得到系数后，滤波后的值yhat_k由多项式在中心点x_k处的值给出：

这个过程体现了Savitzky-Golay滤波器的本质：通过局部多项式拟合来实现数据平滑，同时保持信号的高阶特征。

Python实现与应用示例

以下通过一个完整的示例演示Savitzky-Golay滤波器的应用过程。首先生成含噪声的测试信号：

 importnumpyasnp  
 importmatplotlib.pyplotasplt  
 fromscipy.signalimportsavgol_filter  
       
 np.random.seed(0)  
 x=np.linspace(0, 2*np.pi, 100)  
 y=np.sin(x) +np.random.normal(0, 0.1, x.size)  
       
 plt.plot(x, y, label='Noisy Signal')  # 原始含噪信号
 plt.grid(lw=2,ls=':')  
 plt.xlabel('Time Step')  # 时间步长
 plt.ylabel("Value")      # 信号值
 plt.legend()  
 plt.show()

上图展示了添加高斯噪声后的正弦信号。

使用scipy.signal模块中的savgol_filter函数实现滤波。选择窗口大小为11，多项式阶数为3：

 window_size=11  
 poly_order=3  
 y_smooth=savgol_filter(y, window_size, poly_order)  
       
 plt.plot(x, y, label='Noisy Signal')        # 原始含噪信号
 plt.plot(x, y_smooth, label='Smoothed Signal', color='red')  # 滤波后信号
 plt.grid(lw=2,ls=':')  
 plt.xlabel('Time Step')   # 时间步长
 plt.ylabel("Value")       # 信号值
 plt.legend()  
 plt.show()

滤波结果显示，算法成功地去除了噪声同时保持了信号的基本形状。

上述动画展示了滤波过程中窗口滑动和局部拟合的过程。

参数影响分析

以下代码比较了不同窗口大小和多项式阶数对滤波效果的影响：

 fig, axs=plt.subplots(2, 2, figsize=(20, 12))  
       
 # 配置1：小窗口，低阶多项式
 y_smooth_1=savgol_filter(y_complex, 5, 2)  
 axs[0, 0].plot(x, y_complex, label='Noisy Signal')  
 axs[0, 0].plot(x, y_smooth_1, label='Smoothed Signal (5, 2)', color='red')  
 axs[0, 0].legend()  
 axs[0, 0].set_title('Window Size: 5, Poly Degree: 2')  
 plt.xlabel('Time Step')   # 时间步长
 plt.ylabel("Value")       # 信号值
 plt.legend()  
 
 # 配置2：小窗口，高阶多项式
 y_smooth_2=savgol_filter(y_complex, 5, 4)  
 axs[0, 1].plot(x, y_complex, label='Noisy Signal')  
 axs[0, 1].plot(x, y_smooth_2, label='Smoothed Signal (5, 4)', color='red')  
 axs[0, 1].legend()  
 axs[0, 1].set_title('Window Size: 5, Poly Degree: 4')  
       
 # 配置3：大窗口，低阶多项式
 y_smooth_3=savgol_filter(y_complex, 21, 2)  
 axs[1, 0].plot(x, y_complex, label='Noisy Signal')  
 axs[1, 0].plot(x, y_smooth_3, label='Smoothed Signal (21, 2)', color='red')  
 axs[1, 0].legend()  
 axs[1, 0].set_title('Window Size: 21, Poly Degree: 2')  
       
 # 配置4：大窗口，高阶多项式
 y_smooth_4=savgol_filter(y_complex, 21, 4)  
 axs[1, 1].plot(x, y_complex, label='Noisy Signal')  
 axs[1, 1].plot(x, y_smooth_4, label='Smoothed Signal (21, 4)', color='red')  
 axs[1, 1].legend()  
 axs[1, 1].set_title('Window Size: 21, Poly Degree: 4')  
       
 plt.tight_layout()  
 plt.show()