Python scipy.interpolate插值

我们采集到的数据都是以离散的点的形式存在的，只有在采样点上才有具体的值，在其他区域都没有值数据。此时就需要插值分析，将采样点的数值根据一定的算法，推算出其他未采样区域的数值。

在讲scipy.interpolate类方法插值函数之前我们先讲两种常见的插值方法：待定系数法和拉格朗日法插值。

待定系数法插值：待定系数法插值在我们拥有n个插值节点时构造一个n次多项式，

 然后可以构造非齐次线性方程组，

 在高数或线性代数里，我们学过范德蒙德行列式，我们可以根据上述非齐次线性方程组构造出它的系数矩阵，

再根据解线性方程组的克拉默（克莱姆） 法则，线性方程组的解确定且唯一，由此我们便可以得到我们的插值函数。

由python代码编写如下：

import numpy as np
import matplotlib.pyplot as plt
from scipy import interpolate
x0=np.linspace(-10,10,21)
y0=np.sin(x0)
plt.plot(x0,y0,'o')
A=np.vander(x0)
p=np.linalg.inv(A)@y0
x1=np.linspace(-10,10,101)
plt.plot(x1,np.polyval(p,x1))

 拉格朗日插值法：相对于待定系数法更为简便，首先构造一组基函数，

并令其满足以下条件，

得到n次拉格朗日插值多项式，由方程组

解的唯一性，n+1个节点的n次拉格朗日插值多项式存在且唯一。

python实现代码如下:

import numpy as np
from scipy import interpolate
import pylab as pl
x=np.linspace(-15,15,11)
y=np.sin(x)
pl.plot(x,y,'o')
x1=np.linspace(-15,15,101)
p=interpolate.lagrange(x,y)
pl.plot(x1,np.polyval(p,x1))

 

此外还要注意当插值次数太高产生的龙格现象。龙格现象_百度百科在计算方法中，有利用多项式对某一函数的近似逼近，计算相应的函数值。一般情况下，多项式的次数越多，需要的数据就越多，而预测也就越准确。插值次数越高，插值结果越偏离原函数的现象称为龙格现象。https://baike.baidu.com/item/%E9%BE%99%E6%A0%BC%E7%8E%B0%E8%B1%A1/5473475

import numpy as np
import matplotlib.pyplot as plt
from scipy import interpolate
x=np.linspace(-10,10,11)
y=np.sin(x)
x1=np.linspace(-10,10,101)
y1=np.sin(x1)
plt.plot(x1,y1,label='sinx')
f1=interpolate.interp1d(x,y,kind='quadratic')
f2=interpolate.interp1d(x,y,kind='cubic')
f3=interpolate.interp1d(x,y,kind=7)
f4=interpolate.interp1d(x,y,kind=9)
plt.plot(x1,f1(x1),label='2')
plt.plot(x1,f2(x1),label='3')
plt.plot(x1,f3(x1),label='5')
plt.plot(x1,f4(x1),label='9')
plt.legend()
plt.show()

 现在我们开始讲scipy.interpolate类函数插值方法

一维插值：使用interpolate.interp1d函数插值

官方文档：

scipy.interpolate.interp1d — SciPy v1.8.1 Manual

使用时主要定义kind的类型来进行不同类别的插值

一次插值kind='inear'

二次插值kind='quadratic'

三次插值kind='cubic'

此外还有0阶插值nearest和zero方法

此两种方法使用如下：

import numpy as np
from scipy import interpolate
import matplotlib.pyplot as plt
x=np.linspace(0,10,11)
y=np.sin(x)
x1=np.linspace(0,10,101)
y1=np.sin(x1)
plt.plot(x1,y1)
f1=interpolate.interp1d(x,y,kind='nearest')
plt.plot(x1,f1(x1))
f2=interpolate.interp1d(x,y,kind='zero')
plt.plot(x1,f2(x1))
plt.plot(x,y)

可构成阶梯状曲线

 二维插值：使用interpolate.interp2d函数插值

对模糊图像处理：

示例函数如下：

import numpy as np
import matplotlib as mpl
import pylab as pl
from scipy import interpolate

def z(x,y):
    return (x+y)*np.exp(-5*(x**2+y**2))
x,y=np.mgrid[-1:1:15j,-1:1:15j]
newfunc=interpolate.interp2d(x,y,z(x,y),kind='cubic')
xnew=np.linspace(-1,1,200)
ynew=np.linspace(-1,1,200)
fnew=newfunc(xnew,ynew)
pl.subplot(121)
im1=pl.imshow(z(x,y),origin='lower',extent=[-1,1,-1,1],cmap=mpl.cm.hot)
pl.colorbar(im1)
pl.subplot(122)
im2=pl.imshow(fnew,origin='lower',extent=[-1,1,-1,1],cmap=mpl.cm.hot)
pl.colorbar(im2)
pl.show()

 二维插值的三维图片展示

函数仍未上示例函数

import numpy as np
import matplotlib as mpl
from scipy import interpolate
import matplotlib.pyplot as plt

def func(x, y):
    return (x + y) * np.exp(-5.0 * (x ** 2 + y ** 2))

x,y=np.mgrid[-1:1:20j,-1:1:20j]  
fvals = func(x, y)  
fig = plt.figure(figsize=(18,10))  
ax = plt.subplot(121, projection='3d')
surf = ax.plot_surface(x, y, fvals, rstride=1, cstride=1, cmap=mpl.cm.coolwarm, linewidth=0.5,
                       antialiased=True)
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('z')  
plt.colorbar(surf, shrink=0.5, aspect=5,pad=0.1)  
newfunc = interpolate.interp2d(x, y, fvals, kind='cubic')  
xnew = np.linspace(-1, 1, 200) 
ynew = np.linspace(-1, 1, 200) 
fnew = newfunc(xnew, ynew)  
xnew, ynew = np.meshgrid(xnew, ynew)
ax2 = plt.subplot(122, projection='3d')
surf2 = ax2.plot_surface(xnew, ynew, fnew, rstride=1, cstride=1, cmap=mpl.cm.coolwarm, linewidth=0.5, antialiased=True)
ax2.set_xlabel('xnew')
ax2.set_ylabel('ynew')
ax2.set_zlabel('znew')
plt.colorbar(surf2, shrink=0.5, aspect=5,pad=0.1)  
plt.show()