Savitzky-Golay滤波器原理-01

0 前言

最近在看文章时，发现有人使用Savitzky-Golay滤波器对数据进行预处理，提取背景值，尝试之后发现效果显著。因此，简单研究下Savitzky-Golay滤波器的原理。

1 Savitzky-Golay滤波器原理

1.1 简介

Savitzky-Golay滤波器（通常简称为S-G滤波器）最初由Savitzky和Golay于1964年提出，发表于Analytical Chemistry 杂志。之后被广泛地运用于数据流平滑除噪，是一种在时域内基于局域多项式最小二乘法拟合的滤波方法。这种滤波器最大的特点在于在滤除噪声的同时可以确保信号的形状、宽度不变${}^{[1]}$。

Savitzy-Golay 卷积平滑算法是移动平滑算法的改进。用 Savitzky-Golay 方法进行平滑滤波，可以提高光谱的平滑性，并降低噪音的干扰。Savitzy-Golay 平滑滤波的效果，随着选取窗宽不同而不同，可以满足不同场合的需求${}^{[1]}$。

1.2 推导

思想：使用一定长度窗口的数据对窗口中心的数据进行平滑处理。

设窗口的宽度为$2m+1$，窗口内的数据值为$y_i$， $i=-m,-m+1,\dots ,-1,0,1,\dots ,m-1,m$，使用$n$阶多项式对其进行拟合，有
$$f_i=f(i)=\sum _{k=0}^n{b}_k{i}^k$$
那么拟合数据与原始数据的残差平方和为
$$R=\sum _{i=-m}^m(f_i-y_i{)}^2=\sum _{i=-m}^m{(\sum _{k=0}^n{b}_k{i}^k-y_i)}^2$$
使用最小二乘的思想，为了使拟合效果最好，需要使残差的平方和最小，则需要使得所有拟合系数的一阶偏导为零。对于第$l$个系数$b_l$有
$$\frac{\partial R}{\partial b_l}=2\sum _{i=-m}^m(\sum _{k=0}^n{b}_k{i}^k-y_i)i^l=0$$
整理后可得
$$\sum _{i=-m}^m{y}_i{i}^l=\sum _{k=0}^n{b}_k\sum _{i=-m}^m{i}^{k+l}$$
对于$l\in [0,n]$可以得到$n+1$个方程，联立后可以解得所有系数$b_l$。写成矩阵形式，有
$$\left[\begin{array}{ccccc}(-m{)}^0& \dots & 0^0& \dots & m^0\\ ⋮& & ⋮& & ⋮\\ (-m{)}^l& \dots & 0^l& \dots & m^l\\ ⋮& & ⋮& & ⋮\\ (-m{)}^n& \dots & 0^n& \dots & m^n\end{array}\right]\left[\begin{array}{ccccc}(-m{)}^0& \dots & (-m{)}^k& \dots & (-m{)}^n\\ ⋮& & ⋮& & ⋮\\ 0^0& \dots & 0^k& \dots & m^n\\ ⋮& & ⋮& & ⋮\\ m^0& \dots & m^k& \dots & m^n\end{array}\right]\left[\begin{array}{c}{b}_0\\ ⋮\\ b_k\\ ⋮\\ b_n\end{array}\right]=\left[\begin{array}{ccccc}(-m{)}^0& \dots & 0^0& \dots & m^0\\ ⋮& & ⋮& & ⋮\\ (-m{)}^l& \dots & 0^l& \dots & m^l\\ ⋮& & ⋮& & ⋮\\ (-m{)}^n& \dots & 0^n& \dots & m^n\end{array}\right]\left[\begin{array}{c}{y}_{-m}\\ ⋮\\ y_0\\ ⋮\\ y_m\end{array}\right]$$
得到所有系数后，将其代回拟合方程，并计算中心点处的拟合值
$$f_0=f(0)=b_0$$
此时，我们就得到了Savitzky-Golay滤波器平滑后的中心点的数值大小。

1.3 一阶形式

当使用一阶多项式进行拟合时，代入上一节矩阵形式，很容易可以得到
$$\left[\begin{array}{ccccc}1& \dots & 1& \dots & 1\\ -m& \dots & 0& \dots & m\end{array}\right]\left[\begin{array}{cc}1& -m\\ ⋮& ⋮\\ 1& 0\\ ⋮& ⋮\\ 1& m\end{array}\right]\left[\begin{array}{c}{b}_0\\ b_1\end{array}\right]=\left[\begin{array}{ccccc}1& \dots & 1& \dots & 1\\ -m& \dots & 0& \dots & m\end{array}\right]\left[\begin{array}{c}{y}_{-m}\\ ⋮\\ y_0\\ ⋮\\ y_m\end{array}\right]$$
整理后可得
$$\{\begin{array}{rl}& b_0=\sum _{i=-m}^m{y}_i/\sum _{i=-m}^{m}1\\ & b_1=\sum _{i=-m}^m{y}_{i}i/\sum _{i=-m}^m{i}^2\end{array}$$
从而可以得到中心点处的拟合值为
$$f_0=b_0=\sum _{i=-m}^m{y}_i/(2m+1)$$
此即为滑窗平均的表达式。

参考：

1. 【UWB】Savitzky Golay filter SG滤波器原理讲解
2. Savitzky-Golay滤波器原理阐述