3.7.1 一维插值

一维数据的插值运算可以通过interpld( )完成。其调用形式如下，它实际上不是函数而是一个类：

    interp1d(x, y, kind='linear', ...)

其中，x和y参数是一系列已知的数据点，kind参数是插值类型，可以是字符串或整数，它给出插值的B样条曲线的阶数，可以有如下候选值：

●'zero'、'nearest'：阶梯插值，相当于0阶B样条曲线。

●'slinear'、'linear'：线性插值，用一条直线连接所有的取样点，相当于一阶B样条曲线，'slinear'使用扩展库中的相关函数计算，而'linear'则直接使用Python编写的函数进行运算，其结果一样。

●'quadratic'、'cubic'：二阶和三阶B样条曲线，更高阶的曲线可以直接使用整数值指定。

interpld对象可以计算x的取值范围之内任意点的函数值。它可以像函数一样直接调用，和NumPy的ufunc函数一样能对数组中的每个元素进行计算，并返回一个新的数组。

下面的程序演示了kind参数以及与其对应的插值曲线，结果如图3-38所示。程序中我们使用循环对相同的数据进行4种不同阶数的插值运算。❶首先使用数据点创建一个interpld对象f，通过kind参数指定其阶数。❷调用f( )计算出一系列的插值结果。本例中，决定插值曲线的数据点一共有11个，插值之后的曲线数据点有101个。

图3-38 interpld的各阶插值

高次interp1d( )插值的运算量很大，因此对于点数较多的数据，建议使用后面介绍的UnivariateSpline( )。

    from scipy import interpolate
    
    x = np.linspace(0, 10, 11)
    y = np.sin(x)
    
    xnew = np.linspace(0, 10, 101)
    pl.plot(x,y,'ro')
    for kind in ['nearest', 'zero', 'slinear', 'quadratic']:
        f = interpolate.interp1d(x,y,kind=kind) ❶
        ynew = f(xnew) ❷
        pl.plot(xnew, ynew, label=str(kind))
    
    pl.legend(loc='lower right')

1．外推和Spline拟合

上节所介绍的interpld类要求其参数x是一个递增的序列，并且只能在x的取值范围之内进行内插计算，不能用它进行外推运算，即计算x的取值范围之外的数据点。UnivariateSpline类的插值运算比interpld更高级，它支持外推和拟合运算，其调用形式如下：

    UnivariateSpline(x, y, w=None, bbox=[None, None], k=3, s=None)

●x、y是保存数据点的X-Y坐标的数组，其中x必须是递增序列。

●w是为每个数据点指定的权重值。

●k为样条曲线的阶数。

●s是平滑系数，它使得最终生成的样条曲线满足条件：。即当s>0时，样条曲线并不一定通过各个数据点。为了让曲线通过所有数据点，必须将s参数设置为0。此外，还可以使用InterpolatedUnivariateSpline类，它与UnivariateSpline的唯一区别就是它通过所有的数据点，相当于将s设置为0。

下面的程序演示了使用UnivariateSpline对数据进行插值、外推以及样条曲线拟合：

❶如图3-39（上）所示，UnivariateSpline能够进行外推运算，虽然输入数据中没有X轴大于10的点，但是它能计算出X轴在0到12的插值结果。在X轴大于10的部分，样条曲线仍然呈现出和正弦波类似的形状，越远离输入数据范围，误差会越大，因此外推的范围是有限的。由于s参数为0，因此插值曲线经过所有的数据点。

图3-39 使用UnivariateSpline进行插值：外推（上）和数据拟合（下）

❷图3-39（下）则显示了s参数不为零时的结果，对于带噪声的输入数据，选择合适的s参数能够使得样条曲线接近无噪声时的波形，可以把它看作使用样条曲线对数据进行拟合运算。

    x1 = np.linspace(0, 10, 20)
    y1 = np.sin(x1)
    sx1 = np.linspace(0, 12, 100)
    sy1 = interpolate.UnivariateSpline(x1, y1, s=0)(sx1) ❶
    
    x2 = np.linspace(0, 20, 200)
    y2 = np.sin(x2) + np.random.standard_normal(len(x2))*
0.2
    sx2 = np.linspace(0, 20, 2000)
    spline2 = interpolate.UnivariateSpline(x2, y2, s=8) ❷
    sy2 = spline2(sx2) 
    
    pl.figure(figsize=(8, 5))
    pl.subplot(211)
    pl.plot(x1, y1, ".", label=u"数据点")
    pl.plot(sx1, sy1, label=u"spline曲线")
    pl.legend( )
    
    pl.subplot(212)
    pl.plot(x2, y2, ".", label=u"数据点")
    pl.plot(sx2, sy2, linewidth=2, label=u"spline曲线")
    pl.plot(x2, np.sin(x2), label=u"无噪声曲线")
    pl.legend( )

当曲线为三阶曲线时，UnivariateSpline.roots( )可以用于计算曲线与y=0的交点横坐标。下面显示了图3-39（下）中的曲线与X轴的6个交点的横坐标：

    print np.array_str( spline2.roots( ), precision=3 )
    [  3.288   6.329   9.296  12.578  15.75   18.805]

如果需要计算曲线与任意横线的交点，可以事先将曲线的拟合数据在Y轴方向上进行平移。但若要计算与多条y=c横线的交点，则需要对同样的数据进行平移和拟合多次。

❶下面的root_at( )通过直接修改拟合曲线的参数，实现曲线的平移，从而可以计算与任意横线的交点。❷将roots_at( )动态添加为UnivariateSpline类的方法。❸对多条横线求交点，并进行绘图，其结果如图3-40所示。

图3-40 计算Spline与水平线的交点

    def roots_at(self, v): ❶
        coeff = self.get_coeffs( )
        coeff -= v
        try:
            root = self.roots( )
            return root
        finally:
            coeff += v
    
    interpolate.UnivariateSpline.roots_at = roots_at ❷
    
    pl.plot(sx2, sy2, linewidth=2, label=u"spline曲线")
    
    ax = pl.gca( )
    for level in [0.5, 0.75, -0.5, -0.75]:
        ax.axhline(level, ls=":", color="k")
        xr = spline2.roots_at(level) ❸
        pl.plot(xr, spline2(xr), "ro")

2．参数插值

前面介绍的插值函数都需要X轴的数据是按照递增顺序排列的，就像一般的y=f(x)函数曲线一样。数学上还有一种参数曲线，它使用参数t和两个函数x=f(t)、y=g(t)来定义二维平面上的一条曲线，例如圆形、心形等曲线都是参数曲线。参数曲线的插值可以通过splprep( )和splev( )实现，这组函数支持高维空间的曲线的插值，这里以二维曲线为例介绍其用法。

❶首先调用splprep( )，其第一个参数为一组一维数组，每个数组是各点在对应轴上的坐标。s参数为平滑系数，与UnivariateSpline的含义相同。splprep( )返回两个对象，其中tck是一个元组，它包含了插值曲线的所有信息。t是自动计算出的参数曲线的参数数组。

❷调用splev( )进行插值运算，其第一个参数为一个新的参数数组，这里将t的取值范围等分200份，第二个参数为splprep( )返回的第一个对象。实际上，参数数组t是正规化之后的各个线段长度的累计，因此t的范围为0到1。

其结果如图3-41所示，图中比较了平滑系数为0和1e-4时的插值曲线。当平滑系数为0时，插值曲线通过所有的数据点。

图3-41 使用参数插值连接二维平面上的点

    x = [ 4.913,  4.913,  4.918,  4.938,  4.955,  4.949,  4.911,
          4.848,  4.864,  4.893,  4.935,  4.981,  5.01 ,  5.021]
    
    y = [ 5.2785,  5.2875,  5.291 ,  5.289 ,  5.28  ,  5.26  ,  5.245 ,
          5.245 ,  5.2615,  5.278 ,  5.2775,  5.261 ,  5.245 ,  5.241]
    
    pl.plot(x, y, "o")
    
    for s in (0, 1e-4):
        tck, t = interpolate.splprep([x, y], s=s) ❶
        xi, yi = interpolate.splev(np.linspace(t[0], t[-1], 200), tck) ❷
        pl.plot(xi, yi, lw=2, label=u"s=%g" % s)
    
    pl.legend( )

3．单调插值

前面介绍的几种插值方法不能保证数据点的单调性，即曲线的最值可能出现在数据点之外的地方。PchipInterpolator类（别名pchip）使用单调三次插值，能够保证曲线的所有最值都出现在数据点之上。下面的程序用pchip( )对数据点进行插值，并绘制其一阶导数曲线，由图3-42的导数曲线可知，所有最值点处的导数都为0。

图3-42 单调插值能保证两个点之间的曲线为单调递增或递减

    x = [0, 1, 2, 3, 4, 5]
    y = [1, 2, 1.5, 2.5, 3, 2.5]
    xs = np.linspace(x[0], x[-1], 100)
    curve = interpolate.pchip(x, y)
    ys = curve(xs)
    dys = curve.derivative(xs)
    pl.plot(xs, ys, label=u"pchip")
    pl.plot(xs, dys, label=u"一阶导数")
    pl.plot(x, y, "o")
    pl.legend(loc="best")
    pl.grid( )
    pl.margins(0.1, 0.1)
上一章目录下一章