Savitsky-Golay 平滑算法

本文首发于算法社区,转载请注明出处，谢谢。

名称

Savitsky-Golay平滑算法

作用

信号平滑处理可以除去高频噪音对数据信号的干扰，是消除噪音最常用的一种方法

算法原理

通过多项式对移动窗口内的数据进行多项式最小二乘拟合，算出窗口内中心点关于其周围点的加权平均和。
计算公式：

\[X_i^*=\frac{\sum\limits_{j=-r}^r X_{i+j}W_j}{\sum\limits_{j=-r}^r W_j} (3-1) \]

式（3-1）中的 \(X_{i}\) 和 \(X_i^*\) 为平滑前、后的数据，\(W_j\) 是移动窗口平滑中的权重因子（窗口长度 2r + 1）。
窗口移动多项式拟合的基本思想是：利用多项式对窗口内 N = 2r + 1 的波谱点进行多项式最小二乘拟合，将窗口内等距离的数据拟合成 k 次多项式：

\[X_{i,j}^{\star}=a_0+a_1\lambda_j+...+a_k\lambda_j^k\quad i=1,2,...,n,j=j-r,...j,...,j+r \]

其中 \(a_0,a_1,...,a_k\) 为待定系数或权重。

实例

以窗口为5的二次多项式为例进行拟合。
设有5个数据对（等间隔）：\((x_1,y_1),(x_2,y_2),(x_3,y_3),(x_4,y_4),(x_5,y_5)\)，进行二次多项式拟合：\(\Big|Y=a_0+a_1X+a_2X^2\)
应用最小二乘法，使误差 \(e=\sum_{i=1}^5\left(a_0+a_1x_i+a_2x_i^2-y_i\right)^2;\) 达到最小，所以分别对 \(a_0,a_1,a_2\) 进行求导 \(\frac{\partial e}{\partial a_i} =0\)，得：

\(\begin{cases}a_0\sum_{i=1}^{5}1+a_1\sum_{i=1}^{5}x_i+a_2\sum_{i=1}^{5}x_i^2=\sum_{i=1}^{5}y_i\\ a_0\sum_{i=1}^{5}x_i+a_1\sum_{i=1}^{5}x_i^2+a_2\sum_{i=1}^{5}x_i^3=\sum_{i=1}^{5}y_i x_i\\ a_0\sum_{i=1}^{5}x_i^2+a_1\sum_{i=1}^{5}x_i^3+a_2\sum_{i=1}^{5}x_i^4=\sum_{i=1}^{5}y_i x_i^2\\ \end{cases}\)（4-1）

由于等间隔，所以 \(x_i=i*\Delta x,i=-2,-1,0,1,2\)，式（4-1）可化为

\(\begin{cases}a_0\sum_{i=-2}^{2}1+a_1\sum_{i=-2}^{2}x_i+a_2\sum_{i=-2}^{2}x_i^2=\sum_{i=-2}^{2}y_i\\ a_0\sum_{i=-2}^{2}x_i+a_1\sum_{i=-2}^{2}x_i^2+a_2\sum_{i=-2}^{2}x_i^3=\sum_{i=-2}^{2}y_i x_i\\ a_0\sum_{i=-2}^{2}x_i^2+a_1\sum_{i=-2}^{3}x_i^3+a_2\sum_{i=-2}^{2}x_i^4=\sum_{i=-2}^{2}y_i x_i^2\\ \end{cases}\)（4-2）

注意式（4-2）下标的改变

而

\(\begin{aligned} & \\ &\sum_{i=-2}^{2}x_{i}=\sum_{i=-2}^{2}i\Delta x=\Delta x\sum_{i=-2}^{2}i=0, \\ &\sum_{i=-2}^{2}x_{i}^{2}= \Delta x^{2}\sum_{i=-2}^{2}i^{2}=10\Delta x^{2}, \\ &\sum_{i=-2}^{2}x_{i}^{3}= \Delta x^{3}\sum_{i=-2}^{2}i^{3}=0, \\ &\sum_{i=-2}^{2}x_{i}^{4}= \Delta x^{4}\sum_{i=2}^{2}i^{4}=34\Delta x^{4} \end{aligned}\)（4-3）

由式（4-2）和式（4-3）即可求得：
\(\begin{aligned} &a_0 =(-3y_{-2}+12y_{-1}+17y_0+12y_1-3y_2)/35 \\ & a_1=(-2y_{-2}-y_{-1}+y_1+2y_2)/(10\Delta x) \\ &a_2 =(2y_{-2}-y_{-1}-2y_0-y_1+2y_2)/\left(14\Delta x^2\right) \end{aligned}\) （4-4）

对于中心点，也即所求的平滑点，其值为：

\(y=a_{0}+a_{1}*x_{0}+a_{2}*{x_{0}}^{2}=a_{0},x_0=0*\bigtriangleup x\)

由式（4-4）知

\(a_0=(-3y_{-2}+12y_{-1}+17y_0+12y_1-3y_2)/35\),

故权重为

\(W=[\frac{-3}{35},\frac{12}{35},\frac{17}{35},\frac{12}{35},\frac{-3}{35}]\), 即可利用权重对中心点进行平滑