一种快速做到这一点的方法是通过与高斯内核的导数卷积。简单的情况是你的数组与[-1,1]的卷积,这就是简单的有限差分公式。除此之外,(f * g)’= f’* g = f * g’其中*是卷积,所以你最终得到用高斯表达式进行的衍生,所以当然这样会使你的数据平滑一些可以通过选择最小的合理内核来最小化。
import numpy as np
from scipy import ndimage
import matplotlib.pyplot as plt
#Data:
x = np.linspace(0,2*np.pi,100)
f = np.sin(x) + .02*(np.random.rand(100)-.5)
#Normalization:
dx = x[1] - x[0] # use np.diff(x) if x is not uniform
dxdx = dx**2
#First derivatives:
df = np.diff(f) / dx
cf = np.convolve(f, [1,-1]) / dx
gf = ndimage.gaussian_filter1d(f, sigma=1, order=1, mode='wrap') / dx
#Second derivatives:
ddf = np.diff(f, 2) / dxdx
ccf = np.convolve(f, [1, -2, 1]) / dxdx
ggf = ndimage.gaussian_filter1d(f, sigma=1, order=2, mode='wrap') / dxdx
#Plotting:
plt.figure()
plt.plot(x, f, 'k', lw=2, label='original')
plt.plot(x[:-1], df, 'r.', label='np.diff, 1')
plt.plot(x, cf[:-1], 'r--', label='np.convolve, [1,-1]')
plt.plot(x, gf, 'r', label='gaussian, 1')
plt.plot(x[:-2], ddf, 'g.', label='np.diff, 2')
plt.plot(x, ccf[:-2], 'g--', label='np.convolve, [1,-2,1]')
plt.plot(x, ggf, 'g', label='gaussian, 2')
因为你提到np.gradient我假定你至少有二维数组,所以以下是适用的:如果你想为ndarrays做,这是内置在scipy.ndimage包中的。不过要谨慎,因为这当然不会给你全部的渐变,但我相信各方面的产物。有更好的专业知识的人希望能够说出来。
以下是一个例子:
from scipy import ndimage
x = np.linspace(0,2*np.pi,100)
sine = np.sin(x)
im = sine * sine[...,None]
d1 = ndimage.gaussian_filter(im, sigma=5, order=1, mode='wrap')
d2 = ndimage.gaussian_filter(im, sigma=5, order=2, mode='wrap')
plt.figure()
plt.subplot(131)
plt.imshow(im)
plt.title('original')
plt.subplot(132)
plt.imshow(d1)
plt.title('first derivative')
plt.subplot(133)
plt.imshow(d2)
plt.title('second derivative')
使用gaussian_filter1d可以让您沿着某个轴采取方向导数:
imx = im * x
d2_0 = ndimage.gaussian_filter1d(imx, axis=0, sigma=5, order=2, mode='wrap')
d2_1 = ndimage.gaussian_filter1d(imx, axis=1, sigma=5, order=2, mode='wrap')
plt.figure()
plt.subplot(131)
plt.imshow(imx)
plt.title('original')
plt.subplot(132)
plt.imshow(d2_0)
plt.title('derivative along axis 0')
plt.subplot(133)
plt.imshow(d2_1)
plt.title('along axis 1')
上面的第一组结果对我来说有点混乱(当曲率应该下降时,原始曲线的峰值显示为二阶导数的峰值)。没有进一步了解2d版本如何工作,我只能真正推荐1d版本。如果你想要的大小,只需做一些像:
d2_mag = np.sqrt(d2_0**2 + d2_1**2)