2023年6月29日发(作者：)

Savitzky-Golay平滑滤波的python实现Savitzky-Golay平滑滤波是光谱预处理中常⽤滤波⽅法，它的核⼼思想是对⼀定长度窗⼝内的数据点进⾏k阶多项式拟合，从⽽得到拟合后的结果。对它进⾏离散化处理后后，Ｓ-Ｇ 滤波其实是⼀种移动窗⼝的加权平均算法，但是其加权系数不是简单的常数窗⼝，⽽是通过在滑动窗⼝内对给定⾼阶多项式的最⼩⼆乘拟合得出。

接下来以⼆阶多项式情况下的直线滑动平均法作为特殊例⼦，进⾏详细分析，并给出通⽤的多阶加权系数计算⽅法以及python实现的代码。

⽂章结构：直线滑动平均法对⾃变量x按等距△x作实验观测得数据如下:

令，上述数据变为

⽤下⾯的⽅法修正的值:取定正整数n，使⾄少有⼀个成⽴。当同时成⽴，及窗⼝不在数据边缘时，选取使得

解得：

于是得时的修正公式:

假设当t = 0, 可得修正值：

当只有⼀个成⽴，既窗⼝在数据边缘时，这⾥不做特别讨论，在实现时的思路是将数据按照⾸尾的值进⾏拓展，当然这个思路不够严谨，读者如果希望更为可靠可以⾃⾏研究这⼀部分。由此，可以得到直线滑动平均法在窗⼝⼤⼩为3点，5点，7点…时的结果分别为

3点

5点

通⽤的多阶加权系数计算⽅法通过以上⼆阶多项式拟合的例⼦分析，我们可以对S-G⽅法有个形象的认识。接下来给出通⽤的多阶加权系数计算⽅法。设滤波窗⼝的宽度为，各测量点为采⽤次多项式对窗⼝内的数据点进⾏拟合

于是就有了n个这样的⽅程，扣成了k元线性⽅程组。要使⽅程组有解则n应⼤于等于k，⼀般选择n>k，通过最⼩⼆乘法拟合确定拟合参数A。由此得到

其中， .

由此可得，修正值此处同样不考虑窗⼝在数据边缘时的加权系数计算，实现时处理思路同上。python实现# -*- coding: utf-8 -*-

import numpy as np"""* 创建系数矩阵X* size - 2×size+1 = window_size* rank - 拟合多项式阶次* x - 创建的系数矩阵"""def create_x(size, rank): x = [] for i in range(2 * size + 1): m = i - size row = [m**j for j in range(rank)] (row)

x = (x) return x""" * Savitzky-Golay平滑滤波函数 * data - list格式的1×n纬数据 * window_size - 拟合的窗⼝⼤⼩ * rank - 拟合多项式阶次 * ndata - 修正后的值"""def savgol(data, window_size, rank): m = (window_size - 1) / 2 odata = data[:] # 处理边缘数据，⾸尾增加m个⾸尾项 for i in range(m): (0,odata[0]) (len(odata),odata[len(odata)-1]) # 创建X矩阵 x = create_x(m, rank) # 计算加权系数矩阵B b = (x * (x.T * x).I) * x.T a0 = b[m] a0 = a0.T # 计算平滑修正后的值 ndata = [] for i in range(len(data)): y = [odata[i + j] for j in range(window_size)] y1 = (y) * a0 y1 = float(y1) (y1) return ndata边缘数据处理⽅法不是⼗分严谨，如读者有更⾼的要求可以⾃⾏修改。实验结果对⼀段实际的光谱数据进⾏滑动窗⼝为5，多项式阶数为3的S-G平滑滤波后，得到的结果如下图所⽰：

由于实验数据点数较多，⼤约为3000多个点，此处实验效果不是⾮常明显，但是可以看到的是在400-500nm波段，500-600nm波段内数据明显变得更加光滑。本⽂综合参考了以下⽂献，在这⾥致以衷⼼的感谢：1. 《Savitzky-Golay平滑去噪》：

2. 黄耀欢, 王建华, 江东, 等. 利⽤ SG 滤波进⾏ MODIS-EVI 时间序列数据重构[J]. 武汉⼤学学报信息科学版, 2009, 34(12): 1440-1443.3. 毕志伟, 叶鹰. 数学⼿册:⼤学⽣⽤[M]. ⾼等教育出版社, 2014.