简介

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

平滑滤波是光谱分析中常用的预处理方法之一。用 Savitzky. Golay 方法进行平滑滤波，可以提高光谱的平滑性，并降低噪音的干扰。S-G 平滑滤波的效果，随着选取窗宽不同而不同，可以满足不同场合的需求。

Savitzy-Golay 卷积平滑算法是移动平滑算法的改进。关键在于矩阵算子的求解。

推导

设滤波窗口的宽度为 n=2m+1n=2m+1n=2m+1，各测量点为 x=(−m,−m+1,⋯,0,1,⋯,m−1,m)x=(-m, -m+1, \cdots, 0, 1, \cdots, m-1, m )x=(−m,−m+1,⋯,0,1,⋯,m−1,m)，采用 k−1k-1k−1 次多项式对窗口内的数据点进行拟合
y=a0+a1x+a2x2+⋯+ak−1xk−1y = a_0 + a_1 x + a_2 x^2 + \cdots + a_{k-1} x^{k-1}y=a0+a1x+a2x2+⋯+ak−1xk−1

于是就有了 nnn 个这样的方程，构成了 kkk 元线性方程组。要使方程组有解则 nnn 应大于等于 kkk，一般选择 n>kn > kn>k，通过最小二乘法拟合确定拟合参数 AAA。由此可得到
[y−my−m+1⋮y−1y0y1⋮ym−1ym]=[1(x−m)1⋯(x−m)k−12−1(x−m)k−12(x−m)k−12+1⋯(x−m)k−2(x−m)k−11(x−m+1)1⋯(x−m+1)k−12−1(x−m+1)k−12(x−m+1)k−12+1⋯(x−m+1)k−2(x−m+1)k−1⋮1(x−1)1⋯(x−1)k−12−1(x−1)k−12(x−1)k−12+1⋯(x−1)k−2(x−1)k−11(x0)1⋯(x0)k−12−1(x0)k−12(x0)k−12+1⋯(x0)k−2(x0)k−11(x1)1⋯(x1)k−12−1(x1)k−12(x1)k−12+1⋯(x1)k−2(x1)k−1⋮1(xm−1)1⋯(xm−1)k−12−1(xm−1)k−12(xm−1)k−12+1⋯(xm−1)k−2(xm−1)k−11(xm)1⋯(xm)k−12−1(xm)k−12(xm)k−12+1⋯(xm)k−2(xm)k−1][a0a1⋮ak−12−1ak−12ak−12+1⋮ak−2ak−1]+[e−me−m+1⋮e−1e0e1⋮em−1em]\left[\begin{matrix} y_{-m} \\ y_{-m+1} \\ \vdots \\ y_{-1} \\ y_{0} \\ y_{1} \\ \vdots \\ y_{m-1} \\ y_{m} \\ \end{matrix}\right]= \left[\begin{matrix} 1 & (x_{-m})^{1} & \cdots & (x_{-m})^{\frac{k-1}{2}-1} & (x_{-m})^{\frac{k-1}{2}} & (x_{-m})^{\frac{k-1}{2}+1} & \cdots & (x_{-m})^{k-2} & (x_{-m})^{k-1} \\ 1 & (x_{-m+1})^{1} & \cdots & (x_{-m+1})^{\frac{k-1}{2}-1} & (x_{-m+1})^{\frac{k-1}{2}} & (x_{-m+1})^{\frac{k-1}{2}+1} & \cdots & (x_{-m+1})^{k-2} & (x_{-m+1})^{k-1} \\ &&&& \vdots \\ 1 & (x_{-1})^{1} & \cdots & (x_{-1})^{\frac{k-1}{2}-1} & (x_{-1})^{\frac{k-1}{2}} & (x_{-1})^{\frac{k-1}{2}+1} & \cdots & (x_{-1})^{k-2} & (x_{-1})^{k-1} \\ 1 & (x_{0})^{1} & \cdots & (x_{0})^{\frac{k-1}{2}-1} & (x_{0})^{\frac{k-1}{2}} & (x_{0})^{\frac{k-1}{2}+1} & \cdots & (x_{0})^{k-2} & (x_{0})^{k-1} \\ 1 & (x_{1})^{1} & \cdots & (x_{1})^{\frac{k-1}{2}-1} & (x_{1})^{\frac{k-1}{2}} & (x_{1})^{\frac{k-1}{2}+1} & \cdots & (x_{1})^{k-2} & (x_{1})^{k-1} \\ &&&& \vdots \\ 1 & (x_{m-1})^{1} & \cdots & (x_{m-1})^{\frac{k-1}{2}-1} & (x_{m-1})^{\frac{k-1}{2}} & (x_{m-1})^{\frac{k-1}{2}+1} & \cdots & (x_{m-1})^{k-2} & (x_{m-1})^{k-1} \\ 1 & (x_{m})^{1} & \cdots & (x_{m})^{\frac{k-1}{2}-1} & (x_{m})^{\frac{k-1}{2}} & (x_{m})^{\frac{k-1}{2}+1} & \cdots & (x_{m})^{k-2} & (x_{m})^{k-1} \\ \end{matrix}\right] \left[\begin{matrix} a_{0} \\ a_{1} \\ \vdots \\ a_{\frac{k-1}{2}-1} \\ a_{\frac{k-1}{2}} \\ a_{\frac{k-1}{2}+1} \\ \vdots \\ a_{k-2} \\ a_{k-1} \\ \end{matrix}\right] + \left[\begin{matrix} e_{-m} \\ e_{-m+1} \\ \vdots \\ e_{-1} \\ e_{0} \\ e_{1} \\ \vdots \\ e_{m-1} \\ e_{m} \\ \end{matrix}\right] ⎣⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎡y−my−m+1⋮y−1y0y1⋮ym−1ym⎦⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎤=⎣⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎡1111111(x−m)1(x−m+1)1(x−1)1(x0)1(x1)1(xm−1)1(xm)1⋯⋯⋯⋯⋯⋯⋯(x−m)2k−1−1(x−m+1)2k−1−1(x−1)2k−1−1(x0)2k−1−1(x1)2k−1−1(xm−1)2k−1−1(xm)2k−1−1(x−m)2k−1(x−m+1)2k−1⋮(x−1)2k−1(x0)2k−1(x1)2k−1⋮(xm−1)2k−1(xm)2k−1(x−m)2k−1+1(x−m+1)2k−1+1(x−1)2k−1+1(x0)2k−1+1(x1)2k−1+1(xm−1)2k−1+1(xm)2k−1+1⋯⋯⋯⋯⋯⋯⋯(x−m)k−2(x−m+1)k−2(x−1)k−2(x0)k−2(x1)k−2(xm−1)k−2(xm)k−2(x−m)k−1(x−m+1)k−1(x−1)k−1(x0)k−1(x1)k−1(xm−1)k−1(xm)k−1⎦⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎤⎣⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎡a0a1⋮a2k−1−1a2k−1a2k−1+1⋮ak−2ak−1⎦⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎤+⎣⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎡e−me−m+1⋮e−1e0e1⋮em−1em⎦⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎤

上述公式可能过于复杂不方便理解，我们假设一个五点三次平滑公式，即 m=2,n=2∗2+1=5,k=3m=2, n= 2*2+1 = 5, k=3m=2,n=2∗2+1=5,k=3 代入公式

[y−2y−1y0y1y2]=[1(x−2)1(x−2)21(x−1)1(x−1)21(x0)1(x0)21(x1)1(x1)21(x2)1(x2)2][a0a1a2]+[e−2e−1e0e1e2]\begin{aligned} \left[\begin{matrix} y_{-2} \\ y_{-1} \\ y_{0} \\ y_{1} \\ y_{2} \\ \end{matrix}\right] &= \left[\begin{matrix} 1 & (x_{-2})^{1} & (x_{-2})^{2} \\ 1 & (x_{-1})^{1} & (x_{-1})^{2} \\ 1 & (x_{0})^{1} & (x_{0})^{2} \\ 1 & (x_{1})^{1} & (x_{1})^{2} \\ 1 & (x_{2})^{1} & (x_{2})^{2} \\ \end{matrix}\right] \left[\begin{matrix} a_{0} \\ a_{1} \\ a_{2} \\ \end{matrix}\right] + \left[\begin{matrix} e_{-2} \\ e_{-1} \\ e_{0} \\ e_{1} \\ e_{2} \\ \end{matrix}\right] \\ \end{aligned}⎣⎢⎢⎢⎢⎡y−2y−1y0y1y2⎦⎥⎥⎥⎥⎤=⎣⎢⎢⎢⎢⎡11111(x−2)1(x−1)1(x0)1(x1)1(x2)1(x−2)2(x−1)2(x0)2(x1)2(x2)2⎦⎥⎥⎥⎥⎤⎣⎡a0a1a2⎦⎤+⎣⎢⎢⎢⎢⎡e−2e−1e0e1e2⎦⎥⎥⎥⎥⎤

转换成矩阵表示的形式为
Y(2m+1)×1=X(2m+1)×k⋅Ak×1+E(2m+1)×1Y_{(2m+1) \times 1} = X_{(2m+1) \times k} \cdot A_{k \times 1} + E_{(2m+1) \times 1}Y(2m+1)×1=X(2m+1)×k⋅Ak×1+E(2m+1)×1

AAA 的最小二乘解 A^\hat{A}A^ 为
A^=(XT⋅X)−1⋅XT⋅Y\hat{A} = (X^\text{T} \cdot X)^{-1} \cdot X^\text{T} \cdot YA^=(XT⋅X)−1⋅XT⋅Y

YYY 的模型预测值或滤波值 Y^\hat{Y}Y^ 为
Y^=X⋅A=X⋅(XT⋅X)−1⋅XT⋅Y=B⋅Y\hat{Y} = X \cdot A = X \cdot (X^\text{T} \cdot X)^{-1} \cdot X^\text{T} \cdot Y = B \cdot YY^=X⋅A=X⋅(XT⋅X)−1⋅XT⋅Y=B⋅YB=X⋅(XT⋅X)−1⋅XTB = X \cdot (X^\text{T} \cdot X)^{-1} \cdot X^\text{T}B=X⋅(XT⋅X)−1⋅XT

举例

为了进一步理解上述过程，我们假设信号为 x−2=−2,x−1=−1,x0=0,x1=1,x2=2x_{-2} = -2, x_{-1} = -1, x_{0} = 0, x_{1} = 1, x_{2} = 2x−2=−2,x−1=−1,x0=0,x1=1,x2=2，那么就有

[y−2y−1y0y1y2]=[1(x−2)1(x−2)21(x−1)1(x−1)21(x0)1(x0)21(x1)1(x1)21(x2)1(x2)2][a0a1a2]+[e−2e−1e0e1e2]=[1−241−11100111124][a0a1a2]+[e−2e−1e0e1e2]\begin{aligned} \left[\begin{matrix} y_{-2} \\ y_{-1} \\ y_{0} \\ y_{1} \\ y_{2} \\ \end{matrix}\right] &= \left[\begin{matrix} 1 & (x_{-2})^{1} & (x_{-2})^{2} \\ 1 & (x_{-1})^{1} & (x_{-1})^{2} \\ 1 & (x_{0})^{1} & (x_{0})^{2} \\ 1 & (x_{1})^{1} & (x_{1})^{2} \\ 1 & (x_{2})^{1} & (x_{2})^{2} \\ \end{matrix}\right] \left[\begin{matrix} a_{0} \\ a_{1} \\ a_{2} \\ \end{matrix}\right] + \left[\begin{matrix} e_{-2} \\ e_{-1} \\ e_{0} \\ e_{1} \\ e_{2} \\ \end{matrix}\right] \\ &= \left[\begin{matrix} 1 & -2 & 4 \\ 1 & -1 & 1 \\ 1 & 0 & 0 \\ 1 & 1 & 1 \\ 1 & 2 & 4 \\ \end{matrix}\right] \left[\begin{matrix} a_{0} \\ a_{1} \\ a_{2} \\ \end{matrix}\right] + \left[\begin{matrix} e_{-2} \\ e_{-1} \\ e_{0} \\ e_{1} \\ e_{2} \\ \end{matrix}\right] \end{aligned}⎣⎢⎢⎢⎢⎡y−2y−1y0y1y2⎦⎥⎥⎥⎥⎤=⎣⎢⎢⎢⎢⎡11111(x−2)1(x−1)1(x0)1(x1)1(x2)1(x−2)2(x−1)2(x0)2(x1)2(x2)2⎦⎥⎥⎥⎥⎤⎣⎡a0a1a2⎦⎤+⎣⎢⎢⎢⎢⎡e−2e−1e0e1e2⎦⎥⎥⎥⎥⎤=⎣⎢⎢⎢⎢⎡11111−2−101241014⎦⎥⎥⎥⎥⎤⎣⎡a0a1a2⎦⎤+⎣⎢⎢⎢⎢⎡e−2e−1e0e1e2⎦⎥⎥⎥⎥⎤

再由公式计算 BBB 矩阵。此处借用 Matlab 进行计算，m 代码如下

X = [1 -2  4;1 -1  1;1  0  0;1  1  1;1  2  4;];
B = X * inv(X' * X) * X'

可以得到 BBB 矩阵为
B=135[319−3−53913126−5−3121712−3−56121393−5−3931]B = \frac{1}{35} \left[\begin{matrix} 31 & 9 & -3 & -5 & 3 \\ 9 & 13 & 12 & 6 & -5 \\ -3 & 12 & 17 & 12 & -3 \\ -5 & 6 & 12 & 13 & 9 \\ 3 & -5 & -3 & 9 & 31 \\ \end{matrix}\right]B=351⎣⎢⎢⎢⎢⎡319−3−53913126−5−3121712−3−56121393−5−3931⎦⎥⎥⎥⎥⎤

Y^=135[319−3−53913126−5−3121712−3−56121393−5−3931]⋅Y\hat{Y} = \frac{1}{35} \left[\begin{matrix} 31 & 9 & -3 & -5 & 3 \\ 9 & 13 & 12 & 6 & -5 \\ -3 & 12 & 17 & 12 & -3 \\ -5 & 6 & 12 & 13 & 9 \\ 3 & -5 & -3 & 9 & 31 \\ \end{matrix}\right] \cdot YY^=351⎣⎢⎢⎢⎢⎡319−3−53913126−5−3121712−3−56121393−5−3931⎦⎥⎥⎥⎥⎤⋅Y

展开就可以得到根据 5 个实际值通过多项式计算得到对应的预测值的计算过程
y^−2=135(31y−2+9y−1−3y0−5y1+3y2)=−2y^−1=135(9y−2+13y−1+12y0+6y1−5y2)=−1y^−0=135(−3y−2+12y−1+17y0+12y1−3y2)=0y^−1=135(−5y−2+6y−1+12y0+13y1+9y2)=1y^2=135(3y−2−5y−1−3y0+9y1+31y2)=2\begin{aligned} &\hat{y}_{-2} = \frac{1}{35} (31 y_{-2} + 9 y_{-1} - 3 y_{0} - 5 y_{1} + 3 y_{2}) & = -2 \\ &\hat{y}_{-1} = \frac{1}{35} (9 y_{-2} + 13 y_{-1} + 12 y_{0} + 6 y_{1} - 5 y_{2}) & = -1 \\ &\hat{y}_{-0} = \frac{1}{35} (-3 y_{-2} + 12 y_{-1} + 17 y_{0} + 12 y_{1} - 3 y_{2}) & = 0 \\ &\hat{y}_{-1} = \frac{1}{35} (-5 y_{-2} + 6 y_{-1} + 12 y_{0} + 13 y_{1} + 9 y_{2}) & = 1 \\ &\hat{y}_{2} = \frac{1}{35} (3 y_{-2} - 5 y_{-1} - 3 y_{0} + 9 y_{1} + 31 y_{2}) & = 2 \end{aligned}y^−2=351(31y−2+9y−1−3y0−5y1+3y2)y^−1=351(9y−2+13y−1+12y0+6y1−5y2)y^−0=351(−3y−2+12y−1+17y0+12y1−3y2)y^−1=351(−5y−2+6y−1+12y0+13y1+9y2)y^2=351(3y−2−5y−1−3y0+9y1+31y2)=−2=−1=0=1=2

最小二乘法原理

Ref：

Savitzky-Golay 滤波器
Savitzky-Golay平滑去噪
SG平滑算法（又称多项式平滑算法）
最小二乘法 - WikiPedia

【UWB】Savitzky Golay filter SG滤波器原理讲解相关推荐

【UWB】Savitzky Golay filter SG滤波器快速入门并上手使用
文章目录函数介绍示例程序 Ref: 关于 S-G 滤波器原理的讲解请参考: [UWB]Savitzky Golay filter SG滤波器原理讲解函数介绍 Savitzky Golay fil ...
卡尔曼滤波器原理讲解及其matlab实现
目录一:卡尔曼滤波器的信号模型[1-2] 二:其他方程及变量介绍三:卡尔曼滤波器递推公式四:matlab仿真[3] 参考文献: 引言:在进行一些信号处理的过程中,我们通常会采集到一些数据,但是实 ...
低通滤波器转带通滤波器公式由来_高手讲解滤波器原理(二)，轻松搞定LC滤波器原理...
很多朋友看来,滤波器原理属于难以掌控的内容.但事实上,只要耐心学习,滤波器原理可被轻松掌握.本文中,将为大家讲解LC滤波器原理,并附带LC滤波电路实例.希望通过本文,大家能对LC滤波器原理有更深的理解 ...
matlab中SG滤波器的说明
matlab中SG滤波器的说明项目需要最近在进行光谱预处理方面研究,看了很多论文,发现其中涉及到的方法也就那么几种,比如SG filter.MSC.SNV等等.今天对SG filter 进行了学习, ...
bilateral filter双边滤波器的通俗理解
bilateral filter双边滤波器的通俗理解图像去噪的方法很多,如中值滤波,高斯滤波,维纳滤波等等.但这些降噪方法容易模糊图片的边缘细节,对于高频细节的保护效果并不明显.相比较而言,bila ...
【Matlab】扩展卡尔曼滤波器原理及仿真（初学者入门专用）
文章目录 0.引言及友情链接 1.场景预设 2.扩展卡尔曼滤波器 3.仿真及效果 0.引言及友情链接 \qquad卡尔曼滤波器(Kalman Filter, KF)是传感器融合(Sensor Fusi ...
带通带阻滤波器幅频响应_滤波器原理，各式尽在掌握
滤波器作为电子系统中十分常见的工具,在信号处理中占有重要地位.本文将对两种典型的滤波器原理进行分析,并借以理解其他各式滤波器. 滤波器是一种选频装置,可以使信号中特定的频率成分通过,而极大地衰减其它频 ...
c语言双边滤波算法,浅析bilateral filter双边滤波器的理解
图像去噪的方法很多,如中值滤波,高斯滤波,维纳滤波等等.但这些降噪方法容易模糊图片的边缘细节,对于高频细节的保护效果并不明显.相比较而言,bilateral filter双边滤波器可以很好的边缘保护, ...
AMCL算法原理讲解
ROS进阶教程(二)AMCL算法原理讲解 AMCL算法理解蒙特卡洛定位算法蒙特卡洛定位算法自适应变种里程计运动模型测距仪模型波束模型似然域模型 AMCL算法理解 AMCL(adaptive ...

【UWB】Savitzky Golay filter SG滤波器原理讲解

文章目录

简介

推导

举例

最小二乘法原理

Ref：

【UWB】Savitzky Golay filter SG滤波器原理讲解相关推荐

最新文章

热门文章