一维插值#

分段线性插值#

如果你只需要线性(即折线)插值,你可以使用 numpy.interp 例程。它需要两个用于插值的数据数组,xy,以及一个用于评估插值的点数组 xnew

>>> import numpy as np
>>> x = np.linspace(0, 10, num=11)
>>> y = np.cos(-x**2 / 9.0)

构建插值

>>> xnew = np.linspace(0, 10, num=1001)
>>> ynew = np.interp(xnew, x, y)

并绘制它

>>> import matplotlib.pyplot as plt
>>> plt.plot(xnew, ynew, '-', label='linear interp')
>>> plt.plot(x, y, 'o', label='data')
>>> plt.legend(loc='best')
>>> plt.show()
../../_images/1D-1.png

三次样条插值#

当然,分段线性插值在数据点处会产生拐角,其中线性片段连接。为了生成更平滑的曲线,你可以使用三次样条插值,其中插值曲线由具有匹配的一阶和二阶导数的三次片段组成。在代码中,这些对象通过 CubicSpline 类的实例表示。实例通过数据的 xy 数组构造,然后可以使用目标 xnew 值进行评估:

>>> from scipy.interpolate import CubicSpline
>>> spl = CubicSpline([1, 2, 3, 4, 5, 6], [1, 4, 8, 16, 25, 36])
>>> spl(2.5)
5.57

CubicSpline 对象的 __call__ 方法接受标量值和数组。它还接受第二个参数 nu,用于评估 nu 阶导数。例如,我们绘制样条的导数:

>>> from scipy.interpolate import CubicSpline
>>> x = np.linspace(0, 10, num=11)
>>> y = np.cos(-x**2 / 9.)
>>> spl = CubicSpline(x, y)
>>> import matplotlib.pyplot as plt
>>> fig, ax = plt.subplots(4, 1, figsize=(5, 7))
>>> xnew = np.linspace(0, 10, num=1001)
>>> ax[0].plot(xnew, spl(xnew))
>>> ax[0].plot(x, y, 'o', label='数据')
>>> ax[1].plot(xnew, spl(xnew, nu=1), '--', label='一阶导数')
>>> ax[2].plot(xnew, spl(xnew, nu=2), '--', label='二阶导数')
>>> ax[3].plot(xnew, spl(xnew, nu=3), '--', label='三阶导数')
>>> for j in range(4):
...     ax[j].legend(loc='best')
>>> plt.tight_layout()
>>> plt.show()

注意,一阶和二阶导数在构造上是连续的,而三阶导数在数据点处跳跃。

单调插值#

三次样条在构造上是两次连续可微的。这可能导致样条函数在数据点之间振荡和“过冲”。在这种情况下,一种替代方法是使用所谓的*单调*三次插值:这些插值被构造为仅一次连续可微,并试图保留数据所隐含的局部形状。scipy.interpolate 提供了两种此类对象:PchipInterpolatorAkima1DInterpolator。为了说明,让我们考虑带有异常值的数据:

>>> from scipy.interpolate import CubicSpline, PchipInterpolator, Akima1DInterpolator
>>> x = np.array([1., 2., 3., 4., 4.5, 5., 6., 7., 8])
>>> y = x**2
>>> y[4] += 101
>>> import matplotlib.pyplot as plt
>>> xx = np.linspace(1, 8, 51)
>>> plt.plot(xx, CubicSpline(x, y)(xx), '--', label='样条')
>>> plt.plot(xx, Akima1DInterpolator(x, y)(xx), '-', label='Akima1D')
>>> plt.plot(xx, PchipInterpolator(x, y)(xx), '-', label='pchip')
>>> plt.plot(x, y, 'o')
>>> plt.legend()
>>> plt.show()

B样条插值#

使用B样条进行插值 B样条形成了一种替代的(形式上等价的)分段多项式表示。这种基底通常比幂基底在计算上更稳定,并且对于包括插值、回归和曲线表示在内的各种应用都很有用。详细信息在 分段多项式部分 中给出,在这里我们通过构造正弦函数的插值来说明它们的用法:

>>> x = np.linspace(0, 3/2, 7)
>>> y = np.sin(np.pi*x)

为了在给定数据数组 xy 的情况下构造插值对象,我们使用 make_interp_spline 函数:

>>> from scipy.interpolate import make_interp_spline
>>> bspl = make_interp_spline(x, y, k=3)

该函数返回一个对象,其接口与 CubicSpline 对象的接口类似。特别是,它可以在数据点处进行求值和求导:

>>> der = bspl.derivative()      # 表示导数的 B样条
>>> import matplotlib.pyplot as plt
>>> xx = np.linspace(0, 3/2, 51)
>>> plt.plot(xx, bspl(xx), '--', label=r'$\sin(\pi x)$ approx')
>>> plt.plot(x, y, 'o', label='data')
>>> plt.plot(xx, der(xx)/np.pi, '--', label=r'$d \sin(\pi x)/dx / \pi$ approx')
>>> plt.legend()
>>> plt.show()
../../_images/1D-4.png

注意,通过在 make_interp_spline 调用中指定 k=3,我们请求了一个三次样条(这是默认值,因此 k=3 可以省略);三次的导数是一个二次的:

>>> bspl.k, der.k
(3, 2)

默认情况下,make_interp_spline(x, y) 的结果等价于 CubicSpline(x, y)。不同之处在于前者允许几种可选功能:它可以构造各种次数的样条(通过可选参数 k)和预定义的节点(通过可选参数 t)。

样条插值的边界条件可以通过 bc_type 参数用于 make_interp_spline 函数和 CubicSpline 构造函数。默认情况下,两者都使用 ‘not-a-knot’ 边界条件。

参数化样条曲线#

到目前为止,我们考虑的是样条*函数*,其中数据 y 被期望显式依赖于独立变量 x——因此插值函数满足 \(f(x_j) = y_j\)。样条*曲线*将 xy 数组视为平面上的点 \(\mathbf{p}_j\) 的坐标,并通过这些点的插值曲线由某个附加参数(通常称为 u)参数化。请注意,这种构造很容易推广到更高维度,其中 \(\mathbf{p}_j\) 是 N 维空间中的点。

样条曲线可以很容易地利用插值函数处理多维数据数组的事实来构建。与数据点对应的参数 u 的值需要由用户单独提供。

参数化的选择取决于问题,不同的参数化可能会产生截然不同的曲线。作为一个例子,我们考虑三个参数化(一个有些困难的)数据集,我们取自 BSpline 文档字符串中列出的参考文献 [1] 的第 6 章:

>>> x = [0, 1, 2, 3, 4, 5, 6]
>>> y = [0, 0, 0, 9, 0, 0, 0]
>>> p = np.stack((x, y))
>>> p
array([[0, 1, 2, 3, 4, 5, 6],
       [0, 0, 0, 9, 0, 0, 0]])

我们将 p 数组的元素视为平面上的七个点的坐标,其中 p[:, j] 给出点 \(\mathbf{p}_j\) 的坐标。

首先,考虑*均匀*参数化,\(u_j = j\)

>>> u_unif = x

其次,我们考虑所谓的*弦长*参数化,这只是连接数据点的直线段的累积长度:

\[u_j = u_{j-1} + |\mathbf{p}_j - \mathbf{p}_{j-1}|\]

对于 \(j=1, 2, \dots\)\(u_0 = 0\)。这里 \(| \cdots |\) 是平面 上连续点 \(p_j\) 之间的距离。

>>> dp = p[:, 1:] - p[:, :-1]      # 点之间的2维向量距离
>>> l = (dp**2).sum(axis=0)        # 点之间2维向量的长度的平方
>>> u_cord = np.sqrt(l).cumsum()   # 2范数的累积和
>>> u_cord = np.r_[0, u_cord]      # 第一个点在零处参数化

最后,我们考虑有时被称为*向心*参数化的情况:\(u_j = u_{j-1} + |\mathbf{p}_j - \mathbf{p}_{j-1}|^{1/2}\)。 由于额外的平方根,连续值 \(u_j - u_{j-1}\) 之间的差异将小于弦长参数化的情况:

>>> u_c = np.r_[0, np.cumsum((dp**2).sum(axis=0)**0.25)]

现在绘制生成的曲线:

>>> from scipy.interpolate import make_interp_spline
>>> import matplotlib.pyplot as plt
>>> fig, ax = plt.subplots(1, 3, figsize=(8, 3))
>>> parametrizations = ['uniform', 'cord length', 'centripetal']
>>>
>>> for j, u in enumerate([u_unif, u_cord, u_c]):
...    spl = make_interp_spline(u, p, axis=1)    # 注意 p 是一个2D数组
...
...    uu = np.linspace(u[0], u[-1], 51)
...    xx, yy = spl(uu)
...
...    ax[j].plot(xx, yy, '--')
...    ax[j].plot(p[0, :], p[1, :], 'o')
...    ax[j].set_title(parametrizations[j])
>>> plt.show()
../../_images/1D-5.png

在给定数据定义的域内任意位置使用线性插值。通过传递构成数据的1-D向量来创建此类的一个实例。此类实例定义了一个``__call__``方法,因此可以像函数一样处理,该函数在已知数据值之间进行插值以获得未知值。在实例化时可以指定边界行为。以下示例演示了其使用,包括线性和三次样条插值:

>>> from scipy.interpolate import interp1d
>>> x = np.linspace(0, 10, num=11, endpoint=True)
>>> y = np.cos(-x**2/9.0)
>>> f = interp1d(x, y)
>>> f2 = interp1d(x, y, kind='cubic')
>>> xnew = np.linspace(0, 10, num=41, endpoint=True)
>>> import matplotlib.pyplot as plt
>>> plt.plot(x, y, 'o', xnew, f(xnew), '-', xnew, f2(xnew), '--')
>>> plt.legend(['data', 'linear', 'cubic'], loc='best')
>>> plt.show()
"此代码生成一个时间序列的X-Y图,Y轴为振幅,X轴为时间。原始时间序列显示为一系列蓝色标记,大致定义某种振荡。一条橙色轨迹显示线性插值,绘制在数据之上,形成原始信号的锯齿状表示。一条虚线绿色三次插值也被绘制,似乎平滑地表示源数据。"

interp1d 的 ‘cubic’ 类型等同于 make_interp_spline,而 ‘linear’ 类型等同于 numpy.interp,同时允许 N 维 y 数组。

interp1d 中的另一组插值是 nearestpreviousnext,它们分别返回沿 x 轴的最近点、前一点或后一点。最近点和后点可以被视为因果插值的特例。 插值滤波器。以下示例展示了它们的使用,使用与前一个示例相同的数据:

>>> from scipy.interpolate import interp1d
>>> x = np.linspace(0, 10, num=11, endpoint=True)
>>> y = np.cos(-x**2/9.0)
>>> f1 = interp1d(x, y, kind='nearest')
>>> f2 = interp1d(x, y, kind='previous')
>>> f3 = interp1d(x, y, kind='next')
>>> xnew = np.linspace(0, 10, num=1001, endpoint=True)
>>> import matplotlib.pyplot as plt
>>> plt.plot(x, y, 'o')
>>> plt.plot(xnew, f1(xnew), '-', xnew, f2(xnew), '--', xnew, f3(xnew), ':')
>>> plt.legend(['data', 'nearest', 'previous', 'next'], loc='best')
>>> plt.show()
"此代码生成一个时间序列的X-Y图,Y轴为振幅,X轴为时间。原始时间序列以一系列蓝色标记显示,大致定义某种振荡。一条橙色轨迹显示最近邻插值,绘制在原始数据之上,呈现出阶梯状外观,原始数据位于每个阶梯的中间。一条绿色轨迹显示前邻插值,看起来与橙色轨迹相似,但原始数据位于每个阶梯的后面。同样,一条红色虚线轨迹显示后邻插值,穿过每个前邻点,但它在每个阶梯的前缘居中。"

一般而言,该库严格遵循IEEE 754语义,其中``NaN``表示*非数值*,即非法数学运算的结果(例如除以零),而非*缺失*。