文章目录

  • 1 插值法对曲线平滑处理
    • 1.1 插值法的常见实现方法
    • 1.2 拟合和插值的区别
    • 1.3 代码实例
  • 2 Savitzky-Golay 滤波器实现曲线平滑
    • 2.1 问题描述
    • 2.2 Savitzky-Golay 滤波器–调用讲解
    • 2.3 Savitzky-Golay 曲线平滑处理 示例
    • 2.4 Savitzky-Golay原理剖析
  • 3 基于Numpy.convolve实现滑动平均滤波
    • 3.1 滑动平均概念
    • 3.2 滑动平均的数学原理
    • 3.3 语法
    • 3.4 滑动平均滤波示例

有时我们得到曲线震荡或者噪声比较多，不利于观察曲线的趋势走向，需要对其平滑处理，

本文结介绍Savitzky-Golay 滤波器、make_interp_spline插值法和convolve滑动平均滤波，三种平滑处理方法

1 插值法对曲线平滑处理

实现所需的库： numpy、scipy、matplotlib

1.1 插值法的常见实现方法

nearest：最邻近插值法zero：阶梯插值slinear：线性插值quadratic、cubic：2、3阶B样条曲线插值

1.2 拟合和插值的区别

1、插值：简单来说，插值就是根据原有数据进行填充，最后生成的曲线一定过原有点。

2拟合：拟合是通过原有数据，调整曲线系数，使得曲线与已知点集的差别（最小二乘）最小，最后生成的曲线不一定经过原有点。

1.3 代码实例

代码语法：

通过执行from scipy.interpolate import make_interp_spline，导入make_interp_spline模块，之后调用make_interp_spline(x, y)(x_smooth)函数实现。

官方帮助文档: scipy.interpolate.make_interp_spline

import numpy as npfrom matplotlib import pyplot as pltfrom scipy.interpolate import make_interp_splineSize = 30x = np.arange(Size)y = np.random.randint(1, Size, Size)#平滑前plt.plot(x, y,'r')plt.show()#平滑处理后x_smooth = np.linspace(x.min(), x.max(), 500)# np.linspace 等差数列,从x.min()到x.max()生成300个数，便于后续插值y_smooth = make_interp_spline(x, y)(x_smooth)plt.plot(x_smooth, y_smooth,'b')plt.show()

插值法 平滑处理 前后对比

2 Savitzky-Golay 滤波器实现曲线平滑

2.1 问题描述

在寻找曲线的波峰、波谷时，由于数据帧数多的原因，导致生成的曲线图噪声很大，不易寻找规律。如下图：

由于高频某些点的波动导致高频曲线非常难看，为了降低噪声干扰，需要对曲线做平滑处理，让曲线过渡更平滑。常见的对曲线进行平滑处理的方法包括： Savitzky-Golay 滤波器、插值法等。

2.2 Savitzky-Golay 滤波器–调用讲解

对曲线进行平滑处理，通过Savitzky-Golay 滤波器，可以在scipy库里直接调用，不需要再定义函数。

python中Savitzky-Golay滤波器调用如下：

y_smooth = scipy.signal.savgol_filter(y,53,3)# 亦或y_smooth2 = savgol_filter(y, 99, 1, mode= 'nearest')# 备注：y：代表曲线点坐标（x,y）中的y值数组window_length：窗口长度，该值需为正奇整数。例如：此处取值53k值：polyorder为对窗口内的数据点进行k阶多项式拟合，k的值需要小于window_length。例如：此处取值3mode：确定了要应用滤波器的填充信号的扩展类型。（This determines the type of extension to use for the padded signal to which the filter is applied. ）

调参规律：

现在看一下window_length和k这两个值对曲线的影响。

1）window_length对曲线的平滑作用：
（ window_length的值越小，曲线越贴近真实曲线；window_length值越大，平滑效果越厉害（备注：该值必须为正奇整数）。

1）2）k值对曲线的平滑作用：
（ k值越大，曲线越贴近真实曲线；k值越小，曲线平滑越厉害。另外，当k值较大时，受窗口长度限制，拟合会出现问题，高频曲线会变成直线。

2.3 Savitzky-Golay 曲线平滑处理 示例

# 用于生成问题描述中示例曲线的代码如下：import numpy as npfrom matplotlib import pyplot as pltSize = 100x = np.linspace(1, Size,Size)#生成随机矩阵data = np.random.randint(1, Size, Size)print(data)# 可视化图线plt.plot(x, data,'r')# 使用Savitzky-Golay 滤波器后得到平滑图线from scipy.signal import savgol_filtery = savgol_filter(data, 15, 2, mode= 'nearest')# 可视化图线plt.plot(x, y, 'b', label = 'savgol')#显示曲线plt.show()

#生成的随机矩阵

>>>[61 36 90 88 89 29 36 39 92 62 89 108 66 37 92 14 45 97 35 941 10 15 14 65 55 55 108 57 39 28 62 20 19 30 75 82 71 54 24 40 48 64 65 22 97 61 13 14 69 35 58 612 42 93 43 62 75 39 63 75 82 53 32 86 17 95 89 25 73 47 22 57 85 27 49 47 63 54 616 99 84 78 41 882 41 63 32 43 81 70 75 86 13 57]

Savitzky-Golay 平滑曲线 效果

2.4 Savitzky-Golay原理剖析

在scipy官方帮助文档里可以看到关于Savitzky-Golay 滤波器中关于 savgol_filter()函数 的详细说明。

以下是关于 Savitzky-Golay平滑滤波 的简单介绍（参考Python 生成曲线进行快速平滑处理）：

Savitzky-Golay平滑滤波是光谱预处理中的常用滤波方法，其 核心思想：是对一定长度窗口内的数据点进行k阶多项式拟合，从而得到拟合后的结果。 对它进行离散化处理后，Ｓ-Ｇ 滤波其实是一种移动窗口的加权平均算法，但是其加权系数不是简单的常数窗口，而是通过在滑动窗口内对给定高阶多项式的最小二乘拟合得出。

Savitzky-Golay平滑滤波被广泛地运用于数据流平滑除噪，是一种在时域内基于局域多项式最小二乘法拟合的滤波方法。这种滤波器的 最大特点：在滤除噪声的同时可以确保信号的形状、宽度不变。

使用平滑滤波器对信号滤波时，实际上是拟合了信号中的低频成分，而将高频成分平滑出去了。 如果噪声在高频端，那么滤波的结果就是去除了噪声，反之，若噪声在低频段，那么滤波的结果就是留下了噪声。

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

参考链接：Savitzky-Golay平滑滤波的python实现

3 基于Numpy.convolve实现滑动平均滤波

3.1 滑动平均概念

滑动平均滤波法 （又称：递推平均滤波法），它把连续取N个采样值看成一个队列 ，队列的长度固定为N ，每次采样到一个新数据放入队尾，并扔掉原来队首的一次数据(先进先出原则) 。把队列中的N个数据进行算术平均运算，就可获得新的滤波结果。

N值的选取：流量，N=12；压力：N=4；液面，N=4 ~ 12；温度，N=1~4

滑动平均的优缺点：

优点： 对周期性干扰有良好的抑制作用，平滑度高，适用于高频振荡的系统。
缺点： 灵敏度低，对偶然出现的脉冲性干扰的抑制作用较差，不易消除由于脉冲干扰所引起的采样值偏差，不适用于脉冲干扰比较严重的场合，比较浪费RAM 。

3.2 滑动平均的数学原理

滑动平均滤波法计算类似一维卷积的工作原理，滑动平均的N就对应一维卷积核大小（长度）。 区别在于：

（1）步长会有些区别，滑动平均滤波法滑动步长为1，而一维卷积步长可以自定义；
（2）一维卷积的核参数是需要更新迭代的，而滑动平均滤波法核参数都是1。

我们应该怎么利用这个相似性呢？其实也很简单，只需要把一维卷积核大小（长度）和N相等，步长设置为1，核参数都初始为1就可以了。由于一维卷积计算速度快，因此我们可以使用一维卷积来快速高效地实现这个功能。

滑动平均值是卷积数学运算的一个例子。对于滑动平均值，沿着输入滑动窗口并计算窗口内容的平均值。对于离散的1D信号，卷积是相同的，除了代替计算任意线性组合的平均值，即将每个元素乘以相应的系数并将结果相加。那些系数，一个用于窗口中的每个位置，有时称为卷积核。现在，N值的算术平均值是( x1 + x2 + . . . + xN ) / N ，所以相应的内核是( 1/N , 1/N , . . . , 1 / N ) ，这正是我们通过使用得到的np.ones((N,))/N。

3.3 语法

通过Numpy库中的convolve()函数可以实现这些功能。

def np_move_avg(a,n,mode="same"):return(np.convolve(a, np.ones((n,))/n, mode=mode))

Numpy.convolve函数：（numpy.convolve函数官方文档）

参数说明：

  • a:(N,)输入的第一个一维数组
  • v:(M,)输入的第二个一维数组
  • mode:{‘full’, ‘valid’, ‘same’}参数可选，该参数指定np.convolve函数如何处理边缘。
mode可能的三种取值情况：full’　默认值，返回每一个卷积值，长度是N+M-1,在卷积的边缘处，信号不重叠，存在边际效应。‘same’　返回的数组长度为max(M, N),边际效应依旧存在。‘valid’ 　返回的数组长度为max(M,N)-min(M,N)+1,此时返回的是完全重叠的点。边缘的点无效。

和一维卷积参数类似，a就是被卷积数据，v是卷积核大小。

3.4 滑动平均滤波示例

np.convolve函数中通过mode参数指定如何处理边缘。

下面是一个说明模式不同取值之间差异的图：

import numpy as npimport matplotlib.pyplot as plt def np_move_avg(a,n,mode="same"):return(np.convolve(a, np.ones((n,))/n, mode=mode)) modes = ['full', 'same', 'valid']for m in modes:plt.plot(np_move_avg(np.ones((200,)), 50, mode=m)) plt.axis([-10, 251, -.1, 1.1]) plt.legend(modes, loc='lower center') plt.show()


参考链接：
[开发技巧]·Python极简实现滑动平均滤波（基于Numpy.convolve）

numpy中的convolve的理解

典型范例：

# 实现数据可视化中的数据平滑import numpy as npimport matplotlib.pylab as plt '''其它的一些知识点：raise：当程序发生错误，python将自动引发异常，也可以通过raise显示的引发异常一旦执行了raise语句，raise语句后面的语句将不能执行''' def moving_average(interval, windowsize):window = np.ones(int(windowsize)) / float(windowsize)re = np.convolve(interval, window, 'same')return re def LabberRing():t = np.linspace(-4, 4, 100) # np.linspace 等差数列,从-4到4生成100个数print('t=', t) # np.random.randn 标准正态分布的随机数，np.random.rand 随机样本数值y = np.sin(t) + np.random.randn(len(t)) * 0.1 # 标准正态分布中返回1个，或者多个样本值print('y=', y)plt.plot(t, y, 'k') # plot(横坐标，纵坐标， 颜色)y_av = moving_average(y, 10)plt.plot(t, y_av, 'b')plt.xlabel('Time')plt.ylabel('Value')# plt.grid()网格线设置plt.grid(True)plt.show()return LabberRing()# 调用函数