Savitzky-Golay平滑滤波的python

直接调用库： 

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

data = np.exp(np.linspace(0,5,100)) * np.random.random(100)
zs=savgol_filter(data, 5, 3) # window size 5, polynomial order 3
plt.figure()
plt.plot(zs)
plt.plot(data)
plt.show()

自写：

# -*- coding: utf-8 -*-
import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import savgol_filter

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

"""
 * Savitzky-Golay平滑滤波函数
 * data - list格式的1×n纬数据
 * window_size - 拟合的窗口大小
 * rank - 拟合多项式阶次
 * ndata - 修正后的值
"""
def savgol(data, window_size, rank):
    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

if __name__ == '__main__':
    data = np.exp(np.linspace(0, 5, 100)) * np.random.random(100)

    fig, ax = plt.subplots(figsize=(4, 3))
    ax.plot(range(len(data)), data)

    zs = savgol_filter(data, 5, 3)
    ax.plot(range(len(zs)), zs)

    datenew = savgol(list(data),5,3)
    ax.plot(range(len(datenew)), datenew)
    plt.show()