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

序号	内容
1	【小项目关键技术六】控制北斗 GPS 定位 / UWB 室内定位
2	【UWB】Savitzky Golay filter SG滤波器原理讲解
3	【UWB】Savitzky Golay filter SG滤波器快速入门并上手使用
4	【UWB】MSE 均方误差、RMSE 均方根误差
5	【UWB】Kalman filter, KF卡尔曼滤波, EKF 扩展卡尔曼滤波
6	【UWB】公式推导计算坐标值
7	【UWB】ELM 极限学习机原理及公式推导
8	【UWB】ELM，Extreme Learning Machine 极限学习机
9	【UWB】数学建模 E 题目个人解题答案 - 2021年第十八届华为杯
10	【UWB】使用 python 操作 jetson 提取 UWB 定位模块的数据，处理成坐标的格式，模块厂家为维特智能

关于 Matlab 程序的操作请参考：【UWB】Savitzky Golay filter SG滤波器快速入门并上手使用

简介

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

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

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

推导

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

其中
$x$ 为待拟合的数据
$y$ 为拟合后的输出数据
$a$ 为待求解的参数。

其本质就是通过求解一个 $k-1$ 次方程的到一个输出值 $y$。
因变量是 $y$，自变量是 $x$，参数是 $a$，后边还会加上参数 $b$。

---

于是就有了 $n=2m+1$ 个这样的方程，构成了 $k$ 元线性方程组。要使方程组有解则 $n$ 应大于等于 $k$，一般选择 $n>k$，通过最小二乘法拟合确定拟合参数 $A$。由此可得到
KaTeX parse error: Got function '\textcolor' with no arguments as subscript at position 1: \̲t̲e̲x̲t̲c̲o̲l̲o̲r̲{##df0030}{#1}

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

$$\begin{array}{rl}\left[\begin{array}{c}{y}_{-2}\\ y_{-1}\\ y_0\\ y_1\\ y_2\end{array}\right]& =\left[\begin{array}{ccc}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{array}\right]\left[\begin{array}{c}{a}_0\\ a_1\\ a_2\end{array}\right]+\left[\begin{array}{c}{b}_{-2}\\ b_{-1}\\ b_0\\ b_1\\ b_2\end{array}\right]\end{array}$$

---

转换成矩阵表示的形式为
$$Y_{(2m+1)\times 1}=X_{(2m+1)\times k}\cdot A_{k\times 1}+B_{(2m+1)\times 1}$$

规整下符号：$Y$ 表示待求解的输出值，$X$ 为观测值

$A$ 的最小二乘解 $\hat{A}$ 为
$$\hat{A}=(X^{\text{T}}\cdot X{)}^{-1}\cdot X^{\text{T}}\cdot Y$$

$Y$ 的模型预测值或滤波值 $\hat{Y}$ 为
$$\hat{Y}=X\cdot A=X\cdot (X^{\text{T}}\cdot X{)}^{-1}\cdot X^{\text{T}}\cdot Y=E\cdot Y$$$$E=X\cdot (X^{\text{T}}\cdot X{)}^{-1}\cdot X^{\text{T}}$$

举例

为了进一步理解上述过程，我们假设信号为 $x_{-2}=-2,x_{-1}=-1,x_0=0,x_1=1,x_2=2$，那么就有

$$\begin{array}{rl}\left[\begin{array}{c}{y}_{-2}\\ y_{-1}\\ y_0\\ y_1\\ y_2\end{array}\right]& =\left[\begin{array}{ccc}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{array}\right]\left[\begin{array}{c}{a}_0\\ a_1\\ a_2\end{array}\right]+\left[\begin{array}{c}{b}_{-2}\\ b_{-1}\\ b_0\\ b_1\\ b_2\end{array}\right]\\ & =\left[\begin{array}{ccc}1& -2& 4\\ 1& -1& 1\\ 1& 0& 0\\ 1& 1& 1\\ 1& 2& 4\end{array}\right]\left[\begin{array}{c}{a}_0\\ a_1\\ a_2\end{array}\right]+\left[\begin{array}{c}{b}_{-2}\\ b_{-1}\\ b_0\\ b_1\\ b_2\end{array}\right]\end{array}$$

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

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

可以得到 $E$ 矩阵为
$$E=\frac{1}{35}\left[\begin{array}{ccccc}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{array}\right]$$

$$\hat{Y}=\frac{1}{35}\left[\begin{array}{ccccc}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{array}\right]\cdot Y$$

展开就可以得到根据 5 个实际值通过多项式计算得到对应的预测值的计算过程
$$\begin{array}{rlr}& {\hat{y}}_{-2}=\frac{1}{35}(31y_{-2}+9y_{-1}-3y_0-5y_1+3y_2)& =-2\\ & {\hat{y}}_{-1}=\frac{1}{35}(9y_{-2}+13y_{-1}+12y_0+6y_1-5y_2)& =-1\\ & {\hat{y}}_{-0}=\frac{1}{35}(-3y_{-2}+12y_{-1}+17y_0+12y_1-3y_2)& =0\\ & {\hat{y}}_{-1}=\frac{1}{35}(-5y_{-2}+6y_{-1}+12y_0+13y_1+9y_2)& =1\\ & {\hat{y}}_2=\frac{1}{35}(3y_{-2}-5y_{-1}-3y_0+9y_1+31y_2)& =2\end{array}$$

最小二乘法原理

---

针对矩阵形式，其控制目标是使误差最小，即
$$\underset{A}{\mathrm{argmin}}=|Y-(AX+B){|}^2$$

将上述等式右侧展开得
$$\begin{array}{rl}|Y-(AX+B){|}^2& =(Y-(AX+B){)}^{\text{T}}(Y-(AX+B))\\ & =(Y^{\text{T}}-X^{\text{T}}{A}^{\text{T}}-B^{\text{T}})(Y-AX-B)\\ & =Y^{\text{T}}Y-Y^{\text{T}}AX-Y^{\text{T}}B-X^{\text{T}}{A}^{\text{T}}Y+X^{\text{T}}{A}^{\text{T}}AX+X^{\text{T}}{A}^{\text{T}}B-B^{\text{T}}Y+B^{\text{T}}AX+B^{\text{T}}B\\ & =\end{array}$$

Ref：

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