用Python做一元插值

interp1d

最常用的插值函数就是interp1d，按照字面意思理解就是插值一个一维函数。其必不可少的输入参数，就是将要被插值的函数的自变量和因变量，输出为被插值后的函数，示例如下

import numpy as np
import matplotlib.pyplot as plt
import scipy.interpolate as si
x = np.arange(0, 10, 0.1)
y = np.sin(x)
f = si.interp1d(x, y)
xnew = np.arange(0, 99)/10
ynew = f(xnew)  #调用经由interp1d返回的函数
plt.plot(x, y, 'o', xnew, ynew, '-')
plt.show()

结果为

而所谓插值，要求只能在特定的两个值之间插入，而对于超出定义域范围的值，是无法插入的

>>> f(15)
# 报错

在无声明的情况下，插值方法默认是线性插值linear，如有其他需求，可变更kind参数来实现，可选插值方法如下：

• 样条插值：其0、1、2、3阶插值参数分别为zero、slinear、quadratic、cubic
• 返回单点：next和previous用于返回上一个或下一个值
• 最邻近插值：nearest采取向下取整；nearest-up采用向上取整。

为了对比这些不同插值算法的差异，可以生成一组数据来看一下

x = np.arange(10)
y = np.sin(x)
plt.scatter(x[1:-1],y[1:-1])

xNew = np.arange(1,9,0.1)

ks = ['zero', 'slinear', 'quadratic', 'cubic']
cs = ['r', 'g', 'b', 'gray']

for i in range(4):
    f = si.interp1d(x,y,kind=ks[i])
    plt.plot(xNew, f(xNew), c=cs[i])

plt.show()

下图中，红、绿、蓝、灰分别代表0到3次插值，可见，尽管只有10个点，但分段的二次函数已经描绘出了三角函数的形状，其插值效果还是不错的。

splrep

interp1d返回的是一个插值对象，而非可读数组。相比之下，splrep可以直接返回样条插值参数，更加直观。

通过splrep得到的插值参数，可以通过splev来生成插值结果。

tck = si.splrep(x,y)    # 默认三次插值
'''
tck返回一个元组，包含三个元素，分别是插值节点、插值系数和
(array([0., 0., 0., 0., 2., 3., 4., 5., 6., 7., 9., 9., 9., 9.]),
 array([-5.3791e-19,  7.6306e-01,  1.3851e+00,  1.6502e-01,
       -8.9345e-01, -1.1320e+00, -3.3192e-01,  1.1549e+00,
        1.1340e+00,  4.1212e-01,  0.0000e+00,  0.0000e+00,
        0.0000e+00,  0.0000e+00]),
 3)
'''
yNew = si.splev(xNew, tck)

在splrep中，可通过参数k来调节插值次数，取值范围是1到5，下面是测试代码，由于得到的图和interp1d的图极其相似，所以就不贴上来了。

plt.scatter(x[1:-1],y[1:-1])
for i in range(4):
    tck = si.splrep(x,y,k=i+1)
    yNew = si.splev(xNew, tck)
    plt.plot(xNew, yNew, c=cs[i])

plt.show()