文章目录
- 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 np
from matplotlib import pyplot as plt
from scipy.interpolate import make_interp_spline
Size = 30
x = 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：窗口长度，该值需为正奇整数。例如：此处取值53
k值：polyorder为对窗口内的数据点进行k阶多项式拟合，k的值需要小于window_length。例如：此处取值3
mode：确定了要应用滤波器的填充信号的扩展类型。（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 np
from matplotlib import pyplot as plt
Size = 100
x = 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_filter
y = 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 10 8 66 37 92 14 45 97 35 94 1 10 15
14 65 55 55 10 8 57 39 28 62 20 19 30 75 82 71 54 24 40 48 64 65 22 97
61 13 14 69 35 58 61 2 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 61 6 99 84 78 41 88 2 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 np
import 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 np
import 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() # 调用函数
文章出处登录后可见！
