Savitzky Golay filter SG滤波器原理讲解

https://zhaojichao.blog.csdn.net/article/details/121238628

Y= X*A

当窗口长度确定，拟合阶数固定，可认为X矩阵的值也固定了
则每一次滤波就是根据已知的测量值Y和X矩阵计算系数A
再通过系数A计算拟合的y值

如果已经通过实验确定了窗口长度和拟合阶数，那么直接取出X矩阵第M行的系数
将该系数与测量值相乘计算量会更小

import numpy as np
import math
import matplotlib.pyplot as plt
from scipy.signal import savgol_filter
import random

def savgol(data: list, window_size: int, rank: int):
    m = int((window_size - 1) / 2)
    odata = data[:]
    # 处理边缘数据，首尾增加m个首尾项
    for i in range(m):
        odata.insert(0, odata[0])
        odata.insert(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 = np.mat(y) * a0
        y1 = float(y1)
        ndata.append(y1)
    return ndata

def create_x(size, rank):
    x = []
    for i in range(2 * size + 1):
        m = i - size
        row = [m**j for j in range(rank)]
        x.append(row)
    x = np.mat(x)
    return x



if __name__=='__main__':
    n_size = 500
    # generate orin data
    orin = [2*math.sin(i/20) for i in range(n_size)]
    # generate noise
    noise = np.random.normal(0,0.2,n_size)
    # generate data
    data = orin + noise

    # method1 use the Scipy Library Functions
    newans1 = savgol_filter(data, 5, 3, mode= 'nearest')

    # method2 use the Written Functions
    newans2 =  savgol(list(data), 5, 3)

    # method3 use the coefficient
    a_coff = 1/35 * np.array( [-3, 12, 17, 12 ,-3])
    newans3 = np.convolve(data, a_coff , mode='valid')
    newans3 = np.insert(newans3, 0, [0,0]) # Add two data at the beginning



    # plot filter1 filter2 filter3
    # the result is the same
    plt.plot(orin, label = "orin")
    plt.plot(data, label = "noise")
    plt.plot(newans1, label = "filter1")
    # plt.plot(newans2, label = "filter2",linestyle="--")
    plt.plot(newans3, label = "filter3",linestyle="--")
    plt.legend()
    plt.show()