Python处理离散序列导数：从多项式拟合到高斯函数的探索

Miles

发布时间：2026-10-13 | 5人浏览过

原创

本文介绍在Python中处理离散数据求导的方法。通过对比多项式拟合与高斯函数拟合，探讨如何利用scipy.optimize.leastq进行曲线拟合并计算各阶导数。文中指出高阶多项式在局部峰值拟合上的局限性，并尝试引入高斯函数，但发现其难以直接求导，最终回归多项式方案并评估其效果。

离散序列导数求解目标与数据背景

面对长度为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])返回一个多项式类对象，其具体用法可参考官方文档numpy.poly1d。该多项式对象提供了deriv方法用于求导，求导结果仍为多项式对象。通过df_value = df_1(x)即可计算出特定x点对应的导数值。

虽然看似实现顺利，但求导结果并不理想。如图所示，实际数据的最高峰值出现在100左右，而拟合曲线显示的峰值却高达120。造成这一偏差的主要原因是多项式拟合难以精确捕捉曲线的局部细节。

高斯函数拟合的尝试与局限

鉴于此，尝试改用高斯函数进行曲线拟合。在Matlab环境中测试发现，三阶高斯拟合能够较好地还原曲线形态。

然而，在NumPy及SciPy库中并未找到类似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

但在通过最小二乘法确定函数参数后，仅能获得拟合后的离散点数据，无法直接对连续函数形式进行解析求导，因此该方案并不适用。