Python求离散序列导数的示例
作者:CaspianR 发布时间:2022-10-04 11:18:52
有一组4096长度的数据,需要找到一阶导数从正到负的点,和三阶导数从负到正的点,截取了一小段。
394.0
388.0
389.0
388.0
388.0
392.0
393.0
395.0
395.0
394.0
394.0
390.0
392.0
按照之前所了解的,对离散值求导其实就是求差分,例如第i点的导数(差分)为:
即在一个宽度为2m+1的窗口内通过计算前后m个值加权后的和得到。但是在实际使用过程中效果不是很好。于是想到了同样在一个宽度为2k+1的窗口内,将这2k+1个点拟合成一个函数,然后求导就可以得到任意阶数的导数值。
首先是函数拟合,使用from scipy.optimize import leastsq即最小二乘拟合
from scipy.optimize import leastsq
class search(object):
def __init__(self, filename):
self.filename = filename
def func(self, x, p):
f = np.poly1d(p)
return f(x)
def residuals(self, p, x, y, reg):
regularization = 0.1 # 正则化系数lambda
ret = y - self.func(x, p)
if reg == 1:
ret = np.append(ret, np.sqrt(regularization) * p)
return ret
def LeastSquare(self, data, k=100, order=4, reg=1, show=1): # k为求导窗口宽度,order为多项式阶数,reg为是否正则化
l = self.len
step = 2 * k + 1
p = [1] * order
for i in range(0, l, step):
if i + step < l:
y = data[i:i + step]
x = np.arange(i, i + step)
else:
y = data[i:]
x = np.arange(i, l)
try:
r = leastsq(self.residuals, p, args=(x, y, reg))
except:
print("Error - curve_fit failed")
fun = np.poly1d(r[0]) # 返回拟合方程系数
df_1 = np.poly1d.deriv(fun) # 求得导函数
df_2 = np.poly1d.deriv(df_1)
df_3 = np.poly1d.deriv(df_2)
df_value = df_1(x)
df3_value = df_3(x)
fun = np.poly1d(r[0]),fun返回的是一个 polynomial class,具体使用可以见官方文档numpy.poly1d
polynomial对象可以使用deriv方法求导数,求得的依然是 polynomial对象。 df_value = df_1(x)所得到的就是x这个几个点求得的导数值。
看似大功告成,但是求导的结果并不是很好,如下图,实际最高点在100左右,但是拟合出来的曲线最高点在120左右,而原因在于使用多项式拟合很难准确拟合曲线。
于是想用高斯函数来实现对曲线的拟合,在matlab中试了下,三阶高斯拟合可以很好的拟合曲线,
但是numpy以及sicpy中没有找到类似poly1d这种对象,虽然可以自己定义高斯函数,如下
def gaussian(self, x, *param):
fun = param[0]*np.exp(-np.power(x - param[2], 2.) / (2 * np.power(param[4], 2.)))+param[1]*np.exp(-np.power(x - param[3], 2.) / (2 * np.power(param[5], 2.)))
return fun
但是,在通过最小二乘拟合得到函数参数后只能得到拟合后的点,无法直接求导数..所以并不适合。
所以还是只能回到多项式拟合,如果4阶多项式不能表征的话,更高阶的呢
总体来说,效果还是可以接受的。
如果下阶段找到好的高斯函数拟合方法,会继续更新。
来源:https://blog.csdn.net/renjunsong0/article/details/70176119


猜你喜欢
- 在 Go语言中通过调用 reflect.TypeOf 函数,我们可以从一个任何非接口类型的值创建一个 reflect.Type 值。refl
- 背景说明服务部署在阿里云的K8s上,配置了基于Prometheus的Grafana监控。原本用的是自定义的Metrics接口统计,上报一些字
- 将np图片(imread后的图片)转码为base64格式def image_to_base64(image_np):image = cv2.
- 本文实例讲述了Django框架静态文件处理、中间件、上传文件操作。分享给大家供大家参考,具体如下:Django静态文件处理、中间件、上传文件
- 基本思路是使用opencv来把随机生成的字符,和随机生成的线段,放到一个随机生成的图像中去。虽然没有加复杂的形态学处理,但是目前看起来效果还
- 本文实例讲述了Go语言Web编程实现Get和Post请求发送与解析的方法。分享给大家供大家参考,具体如下:这是一篇入门文章,通过一个简单的例
- 首先是安装好PHP之后配置环境变量然后在cmd中输入php -v 能看到版本号即为配置好了之后在sublime中新建编译系统,输入代码{&q
- 本文将介绍如何使用 Vue 和第三方组件库 Element UI 实现轮播图功能。我们将从以下几个方面进行讲解:安装 Element UI创
- 1.变量的赋值操作只是多生成了一个变量,实际上还是指向同一个对象# -*- coding: utf-8 -*-class CPU: &nbs
- --创建测试数据库 CREATE DATABASE Db GO --对数据库进行备份 BACKUP DATABASE Db TO DISK=
- 在前几章节中,我们已经学会了如果在一张表中读取数据,这是相对简单的,但是在真正的应用中经常需要从多个数据表中读取数据。本章节我们将向大家介绍
- js删除数组中某一项或几项的几种方法一、删除第一个元素1、shift() 方法用于把数组的第一个元素从其中删除,并返回第一个元素的值。注意:
- var obj = document.getElementById("testSelect"); //定位idvar i
- 一、概述及完整代码对MNIST(MixedNational Institute of Standard and Technology dat
- 1.在Qt Designer中设计两个简单窗口2.将.ui文件转换成.py文件3.新建**.py文件#-*- coding:utf-8 -*
- 本文实例讲述了JS小游戏的象棋暗棋源码,分享给大家供大家参考。具体如下:游戏运行后如下图所示:Javascript 部分:/** chine
- 其实,这是一个非常容易解决掉的问题。在我看来,似曾相识,呵呵,最近学JavaScript可是学会了使用var声明变量。其实,在PHP中根本不
- 高阶函数除了可以接受函数作为参数外,还可以把函数作为结果值返回。看代码:# -*- coding: utf-8 -*-# @File &nb
- 树,因其清晰明了的展现形式而被广泛的使用日常的开发过程中我们需要经常与“树”打交道,例如公司的组织架构树、服务器的项目归属树,管理后台侧边树
- 1.静态传值(在父组件中赋值好props中属性的值传递给子组件)父组件<template> <div>