12.2 Savitzky-Golay平滑滤波

Savitzky-Golay平滑滤波命名自发明者Abraham Savitzky和Marcel J. E. Golay,简称为S-G平滑滤波。根据原理亦称为移动窗口最小二乘多项式平滑,是一种特殊的低通滤波器。该方法结合了卷积和多项式回归的方法来实现平滑滤波,其在众多平滑滤波中是最广泛应用于分析化学数据处理的平滑算子。该算法的另一个优点是亦可用于数值平滑求导。

原理部分

假设数据矢量为x,选定的平滑窗口大小为奇数m,多项式次数为n,当前待平滑的点为x0,前后m个点分别记为:… xm-3,xm-2,xm-1,x0,xm+1,xm+2,xm+3…。

Savitzky-Golay平滑利用中心点以及其前后m-1个点进行多项式的最小二乘拟和。每一个点可以表示为不同的多项式的结果,从而m个点xi可以表示成为含有n+1个未知数b,m个方程的方程组。以下以m=7,n=3为例进行介绍:

\left\{ \begin{gathered} {x_{ - 3}} = {b_0} + {b_1}*( - 3) + {b_2}*{( - 3)^2} + {b_3}*{( - 3)^3} = {b_0} - 3{b_1} + 9{b_2} - 27{b_3} \hfill \\ {x_{ - 2}} = {b_0} + {b_1}*( - 2) + {b_2}*{( - 2)^2} + {b_3}*{( - 2)^3} = {b_0} - 2{b_1} + 4{b_2} - 8{b_3} \hfill \\ {x_{ - 1}} = {b_0} + {b_1}*( - 1) + {b_2}*{( - 1)^2} + {b_3}*{( - 1)^3} = {b_0} - {b_1} + {b_2} - {b_3} \hfill \\ {x_0} = {b_0} + {b_1}*( - 0) + {b_2}*{( - 0)^2} + {b_3}*{( - 0)^3} = {b_0} \hfill \\ {x_1} = {b_0} + {b_1}*(1) + {b_2}*{(1)^2} + {b_3}*{(1)^3} = {b_0} + {b_1} + {b_2} + {b_3} \hfill \\ {x_2} = {b_0} + {b_1}*(2) + {b_2}*{(2)^2} + {b_3}*{(2)^3} = {b_0} + 2{b_1} + 4{b_2} + 8{b_3} \hfill \\ {x_3} = {b_0} + {b_1}*(3) + {b_2}*{(3)^2} + {b_3}*{(3)^3} = {b_0} + 3{b_1} + 9{b_2} + 27{b_3} \hfill \\ \end{gathered} \right.

对于上述线性方程,可以表示为以下矩阵形式:

\left( {\begin{array}{*{20}{c}} 1&{ - 3}&9&{ - 27} \\ 1&{ - 2}&4&{ - 8} \\ 1&{ - 1}&1&{ - 1} \\ 1&0&0&0 \\ 1&1&1&1 \\ 1&2&4&8 \\ 1&3&9&{27} \end{array}} \right){\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} *{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \left( \begin{gathered} {b_0} \hfill \\ {b_1} \hfill \\ {b_2} \hfill \\ {b_3} \hfill \\ \end{gathered} \right){\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} = {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \left( \begin{gathered} {x_{ - 3}} \hfill \\ {x_{ - 2}} \hfill \\ {x_{ - 1}} \hfill \\ {x_0} \hfill \\ {x_1} \hfill \\ {x_2} \hfill \\ {x_3} \hfill \\ \end{gathered} \right)

即:

x=Mb

进而得到方程组的最小二乘解估计值\hat x

\hat x=M(M^T M)M^T x

或其一般形式:

\left\{ \begin{gathered} {x_{-3} = \frac{1}{42}*({39}{x}_{-3}{ + 8}{x}_{-2}{-4}{x}_{-1}{-4}{x}_{0}{ + }{x}_{1}{ + 4}{x}_{2}{-2}{x}_{3}) \hfill \\ {x_{-2} = \frac{1}{42}*({8}{x}_{-3}{ + 19}{x}_{-2}{ + 16}{x}_{-1}{ + 6}{x}_{0}{-4}{x}_{1}{-7}{x}_{2}{ + 4}{x}_{3}) \hfill \\ {x_{-1} = \frac{1}{42}*({-4}{x}_{-3}{ + 16}{x}_{-2}{ + 19}{x}_{-1}{ + 12}{x}_{0}{ + 2}{x}_{1}{-4}{x}_{2}{ + }{x}_{3}) \hfill \\ {x_0} = \frac{1}{42}*({-4}{x}_{-3}{ + 6}{x}_{-2}{ + 12}{x}_{-1}{ + 14}{x}_{0}{ + 12}{x}_{1}{ + 6}{x}_{2}{-4}{x}_{3}) \hfill \\ {x_1} = \frac{1}{42}*({x}_{-3}{-4}{x}_{-2}{ + 2}{x}_{-1}{ + 12}{x}_{0}{ + 19}{x}_{1}{ + 16}{x}_{2}{-4}{x}_{3}) \hfill \\ {x_2} = \frac{1}{42}*({4}{x}_{-3}{-7}{x}_{-2}{-4}{x}_{-1}{ + 6}{x}_{0}{ + 16}{x}_{1}{ + 19}{x}_{2}{ + 8}{x}_{3}) \hfill \\ {x_3} = \frac{1}{42}*({-2}{x}_{-3}{ + 4}{x}_{-2}{ + }{x}_{-1}{-4}{x}_{0}{-4}{x}_{1}{ + 8}{x}_{2}{ + 39}{x}_{3}) \hfill \\ \end{gathered} \right.

包括要计算的点x0在内,7个方程所代表的7个点都是窗口内部各个点的不同加权线性组合。因此从本质上说,移动窗口多项式平滑就是利用窗口内部各个点之间的加权来计算平滑后的新值。