Python开发之实现SG滤波

前言：主要介绍SG滤波的Python实现，顺带介绍SG滤波的实现原理。

---

1 SG滤波

• Savitzky-Golay滤波器（通常简称为S-G滤波器）最初由Savitzky和Golay于1964年提出，发表于Analytical Chemistry 杂志。之后被广泛地运用于数据流平滑除噪，是一种在时域内基于局域多项式最小二乘法拟合的滤波方法。这种滤波器最大的特点在于在滤除噪声的同时可以确保信号的形状、宽度不变。
• 平滑滤波是光谱分析中常用的预处理方法之一。用 Savitzky. Golay 方法进行平滑滤波，可以提高光谱的平滑性，并降低噪音的干扰。S-G 平滑滤波的效果，随着选取窗宽不同而不同，可以满足不同场合的需求。
• Savitzy-Golay 卷积平滑算法是移动平滑算法的改进。关键在于矩阵算子的求解。
  更多理论介绍请查看：
  《【UWB】Savitzky Golay filter SG滤波器原理讲解》
  《Savitzky-Golay 滤波器》
  《Savitzky-Golay平滑去噪》
  《SG平滑算法（又称多项式平滑算法）》
  《最小二乘法》

2 借助Python中的scipy.signal库实现SG滤波

from scipy.signal import savgol_filter
import matplotlib.pyplot as plt

def SG(data,window_size,polyorder):
    smoothed_data = savgol_filter(data, window_size, polyorder).tolist()
    return smoothed_data

if __name__ == '__main__':
    data = [0.3962, 0.4097, 0.2956, 0.4191, 0.3456, 0.3056, 0.6346, 0.7025, 0.6568, 0.4719, 0.5645, 0.6514, 0.5717,
            0.6072, 0.7076, 0.7062, 0.7086, 0.677, 0.8141, 0.7985, 0.7037, 0.7961, 0.6805, 0.5463, 0.2766]

    smoothed_data = SG(data,5,2)
    smoothed_data = [round(i, 4) for i in smoothed_data]#保留四位小数
    smoothed_data2 = SG(data,9,3)
    smoothed_data2 = [round(i, 4) for i in smoothed_data2]#保留四位小数
    print("data：", data)
    print("smoothed_data：", smoothed_data)
    print("smoothed_data2：", smoothed_data2)
    plt.plot(data, label='data')
    plt.plot(smoothed_data, label='smoothed_data')
    plt.plot(smoothed_data2, label='smoothed_data2')
    plt.xlabel('Time')
    plt.ylabel('Value')
    plt.title('Line Plot')
    plt.legend()
    plt.show()

data： [0.3962, 0.4097, 0.2956, 0.4191, 0.3456, 0.3056, 0.6346, 0.7025, 0.6568, 0.4719, 0.5645, 0.6514, 0.5717, 0.6072, 0.7076, 0.7062, 0.7086, 0.677, 0.8141, 0.7985, 0.7037, 0.7961, 0.6805, 0.5463, 0.2766]
smoothed_data： [0.4007, 0.3779, 0.3642, 0.3621, 0.3366, 0.3884, 0.5679, 0.7173, 0.6189, 0.5319, 0.554, 0.6135, 0.6002, 0.6172, 0.6843, 0.7185, 0.688, 0.7219, 0.7803, 0.782, 0.7604, 0.746, 0.7068, 0.5433, 0.2737]
smoothed_data2： [0.4295, 0.3527, 0.326, 0.3402, 0.3864, 0.4741, 0.5428, 0.5704, 0.6193, 0.6247, 0.5715, 0.5688, 0.6023, 0.6561, 0.6591, 0.6804, 0.7305, 0.7409, 0.7562, 0.7778, 0.7894, 0.7619, 0.6805, 0.5283, 0.2882]

-注意：window_size,polyorder两个参数决定着平滑的效果，要根据自己的数据，调整合适的参数，注意不要欠拟合和过拟合。

3 手动代码实现SG滤波

import matplotlib.pyplot as plt
import numpy as np

def SG01(data,window_size):
    # 前后各m个数据，共2m+1个数据，作为滑动窗口内滤波的值
    m = int((window_size - 1) / 2)  # (59-1)  /2 = 29
    # 计算 矩阵X 的值 ，就是将自变量x带进去的值算 0次方，1次方，2次方.....k-1次方,一共window_size行，k列
    # 大小为（2m+1,k）
    X_array = []
    for i in range(window_size):  #
        arr = []
        for j in range(3):
            X0 = np.power(-m + i, j)
            arr.append(X0)
        X_array.append(arr)
    X_array = np.mat(X_array)
    # B = X*（X.T*X）^-1*X.T
    B = X_array * (X_array.T * X_array).I * X_array.T
    data = np.insert(data, 0, [data[0] for i in range(m)])  # 首位插入m-1个data[0]
    data = np.append(data, [data[-1] for i in range(m)])  # 末尾插入m-1个data[-1]
    # 取B中的第m行 进行拟合  因为是对滑动窗口中的最中间那个值进行滤波，所以只要获取那个值对应的参数就行， 固定不变
    B_m = B[m]
    # 存储滤波值
    y_array = []
    # 对扩充的data 从第m个数据开始遍历一直到（data.shape[0] - m）  ：（第m个数据就是原始data的第1个，（data.shape[0] - m）为原始数据的最后一个
    for n in range(m, data.shape[0] - m):
        y_true = data[n - m: n + m + 1]  # 取出真实y值的前后各m个，一共2m+1个就是滑动窗口的大小
        y_filter = np.dot(B_m, y_true)  # 根据公式 y_filter = B * X 算的  X就是y_true
        y_array.append(float(y_filter))  # float(y_filter) 从矩阵转为数值型
    return y_array

if __name__ == '__main__':
    data = [0.3962, 0.4097, 0.2956, 0.4191, 0.3456, 0.3056, 0.6346, 0.7025, 0.6568, 0.4719, 0.5645, 0.6514, 0.5717,
            0.6072, 0.7076, 0.7062, 0.7086, 0.677, 0.8141, 0.7985, 0.7037, 0.7961, 0.6805, 0.5463, 0.2766]
    smoothed_data = SG01(data,5)
    smoothed_data = [round(i, 4) for i in smoothed_data]
    smoothed_data2 = SG01(data,11)
    smoothed_data2 = [round(i, 4) for i in smoothed_data2]
    print("data：", data)
    print("smoothed_data：", smoothed_data)
    print("smoothed_data2：", smoothed_data2)
    plt.plot(data, label='data')
    plt.plot(smoothed_data, label='smoothed_data')
    plt.plot(smoothed_data2, label='smoothed_data2')
    plt.xlabel('Time')
    plt.ylabel('Value')
    plt.title('Line Plot')
    plt.legend()
    plt.show()

data： [0.3962, 0.4097, 0.2956, 0.4191, 0.3456, 0.3056, 0.6346, 0.7025, 0.6568, 0.4719, 0.5645, 0.6514, 0.5717, 0.6072, 0.7076, 0.7062, 0.7086, 0.677, 0.8141, 0.7985, 0.7037, 0.7961, 0.6805, 0.5463, 0.2766]
smoothed_data： [0.4095, 0.3663, 0.3642, 0.3621, 0.3366, 0.3884, 0.5679, 0.7173, 0.6189, 0.5319, 0.554, 0.6135, 0.6002, 0.6172, 0.6843, 0.7185, 0.688, 0.7219, 0.7803, 0.782, 0.7604, 0.746, 0.7068, 0.5015, 0.3344]
smoothed_data2： [0.3916, 0.3559, 0.3443, 0.368, 0.4265, 0.4713, 0.5107, 0.5667, 0.5903, 0.6093, 0.6122, 0.5879, 0.6034, 0.6261, 0.6685, 0.6965, 0.7121, 0.7497, 0.7714, 0.7859, 0.7623, 0.7017, 0.6148, 0.5048, 0.4128]

• 注意：效果不佳，不推荐使用。