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} \left( \sum_{k=0}^{n} b_k i^k - y_i \right)^2
$$
使用最小二乘的思想，为了使拟合效果最好，需要使残差的平方和最小，则需要使得所有拟合系数的一阶偏导为零。对于第$l$个系数$b_l$有
$$
\frac{\partial R}{\partial b_l} = 2 \sum_{i=-m}^{m} \left( \sum_{k=0}^{n} b_k i^k - y_i \right)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$。写成矩阵形式，有
$$
\begin{bmatrix}
(-m)^0 & \dots & 0^0 & \dots & m^0 \\
\vdots & & \vdots & & \vdots \\
(-m)^l & \dots & 0^l & \dots & m^l \\
\vdots & & \vdots & & \vdots \\
(-m)^n & \dots & 0^n & \dots & m^n \\
\end{bmatrix}
\begin{bmatrix}
(-m)^0 & \dots & (-m)^k & \dots & (-m)^n \\
\vdots & & \vdots & & \vdots \\
0^0 & \dots & 0^k & \dots & m^n \\
\vdots & & \vdots & & \vdots \\
m^0 & \dots & m^k & \dots & m^n \\
\end{bmatrix}
\begin{bmatrix}
b_0 \\
\vdots \\
b_k \\
\vdots \\
b_n \\
\end{bmatrix}
= \begin{bmatrix}
(-m)^0 & \dots & 0^0 & \dots & m^0 \\
\vdots & & \vdots & & \vdots \\
(-m)^l & \dots & 0^l & \dots & m^l \\
\vdots & & \vdots & & \vdots \\
(-m)^n & \dots & 0^n & \dots & m^n \\
\end{bmatrix}
\begin{bmatrix}
y_{-m} \\
\vdots \\
y_0 \\
\vdots \\
y_m \\
\end{bmatrix}
$$
得到所有系数后，将其代回拟合方程，并计算中心点处的拟合值
$$
f_0 = f(0) = b_0
$$
此时，我们就得到了Savitzky-Golay滤波器平滑后的中心点的数值大小。

1.3 一阶形式

当使用一阶多项式进行拟合时，代入上一节矩阵形式，很容易可以得到
$$
\begin{bmatrix}
1 & \dots & 1 & \dots & 1 \\
-m & \dots & 0 & \dots & m \\
\end{bmatrix}
\begin{bmatrix}
1 & -m \\
\vdots & \vdots \\
1 & 0 \\
\vdots & \vdots \\
1 & m \\
\end{bmatrix}
\begin{bmatrix}
b_0 \\
b_1 \\
\end{bmatrix}
= \begin{bmatrix}
1 & \dots & 1 & \dots & 1 \\
-m & \dots & 0 & \dots & m \\
\end{bmatrix}
\begin{bmatrix}
y_{-m} \\
\vdots \\
y_0 \\
\vdots \\
y_m \\
\end{bmatrix}
$$
整理后可得
$$
\left\{
\begin{split}
&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{split}
\right.
$$
从而可以得到中心点处的拟合值为
$$
f_0 = b_0 = \sum_{i=-m}^m y_i / (2m+1)
$$
此即为滑窗平均的表达式。

参考：

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