文章详情

短信预约-IT技能 免费直播动态提醒

请输入下面的图形验证码

提交验证

短信预约提醒成功

python 曲线平滑处理——方法总结(Savitzky-Golay 滤波器、make_interp_spline插值法和convolve滑动平均滤波)

2023-09-10 15:42

关注

文章目录

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

本文结介绍Savitzky-Golay 滤波器、make_interp_spline插值法和convolve滑动平均滤波,三种平滑处理方法
在这里插入图片描述

1 插值法对曲线平滑处理

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

1.1 插值法的常见实现方法

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

1.2 拟合和插值的区别

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

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

1.3 代码实例

代码语法:

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

官方帮助文档: https://blog.csdn.net/m0_51233386/article/details/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 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阶多项式拟合,从而得到拟合后的结果。 对它进行离散化处理后,S-G 滤波其实是一种移动窗口的加权平均算法,但是其加权系数不是简单的常数窗口,而是通过在滑动窗口内对给定高阶多项式的最小二乘拟合得出。

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函数官方文档

参数说明:

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 等差数列,-44生成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()  # 调用函数

在这里插入图片描述

来源地址:https://blog.csdn.net/m0_51233386/article/details/129878586

阅读原文内容投诉

免责声明:

① 本站未注明“稿件来源”的信息均来自网络整理。其文字、图片和音视频稿件的所属权归原作者所有。本站收集整理出于非商业性的教育和科研之目的,并不意味着本站赞同其观点或证实其内容的真实性。仅作为临时的测试数据,供内部测试之用。本站并未授权任何人以任何方式主动获取本站任何信息。

② 本站未注明“稿件来源”的临时测试数据将在测试完成后最终做删除处理。有问题或投稿请发送至: 邮箱/279061341@qq.com QQ/279061341

软考中级精品资料免费领

  • 历年真题答案解析
  • 备考技巧名师总结
  • 高频考点精准押题
  • 2024年上半年信息系统项目管理师第二批次真题及答案解析(完整版)

    难度     807人已做
    查看
  • 【考后总结】2024年5月26日信息系统项目管理师第2批次考情分析

    难度     351人已做
    查看
  • 【考后总结】2024年5月25日信息系统项目管理师第1批次考情分析

    难度     314人已做
    查看
  • 2024年上半年软考高项第一、二批次真题考点汇总(完整版)

    难度     433人已做
    查看
  • 2024年上半年系统架构设计师考试综合知识真题

    难度     221人已做
    查看

相关文章

发现更多好内容

猜你喜欢

AI推送时光机
位置:首页-资讯-后端开发
咦!没有更多了?去看看其它编程学习网 内容吧
首页课程
资料下载
问答资讯