Python - 线性插值实现

AIHGF

October 9, 2021

402views

No comments

5501 words

Python

插值问题描述：

已知坐标 (x0, y0) 和 (x1, y1)，求区间 [x0, x1] 内某一点位置 x 所对应的 y 值.

计算公式如：

$$ y = y_0 + (x-x_0) \frac{y_1-y_0}{x_1 - x_0} = y_0 + \frac{(x-x_0)y_1 - (x-x_0)y_0}{x_1 - x_0} $$

Python 中 numpy 库和scipy 库提供了线性插值的功能.

1. numpy.interp

https://numpy.org/doc/stable/reference/generated/numpy.interp.html

使用示例：

import numpy as np 
import matplotlib.pyplot as plt

xp = [1, 2, 3]
fp = [3, 2, 0]
y = np.interp(2.5, xp, fp)
#1.0

y = np.interp([0, 1, 1.5, 2.72, 3.14], xp, fp)
#array([3.  , 3.  , 2.5 , 0.56, 0.  ])

UNDEF = -99.0
y = np.interp(3.14, xp, fp, right=UNDEF)
#-99.0

#sine 函数插值
x = np.linspace(0, 2*np.pi, 10)
y = np.sin(x)
xvals = np.linspace(0, 2*np.pi, 50)
yinterp = np.interp(xvals, x, y)

plt.plot(x, y, 'o')
plt.plot(xvals, yinterp, '-x')
plt.show()

#周期 x 坐标的插值
x = [-180, -170, -185, 185, -10, -5, 0, 365]
xp = [190, -190, 350, -350]
fp = [5, 10, 3, 4]
y = np.interp(x, xp, fp, period=360)
#array([7.5 , 5.  , 8.75, 6.25, 3.  , 3.25, 3.5 , 3.75])

#复数插值Complex interpolation:
x = [1.5, 4.0]
xp = [2,3,5]
fp = [1.0j, 0, 2+3j]
y = np.interp(x, xp, fp)
#array([0.+1.j , 1.+1.5j])

函数定义如：

https://github.com/numpy/numpy/blob/v1.23.0/numpy/lib/function_base.py#L1456-L1594

@array_function_dispatch(_interp_dispatcher)
def interp(x, xp, fp, left=None, right=None, period=None):
    """
    One-dimensional linear interpolation for monotonically increasing sample points.
    Returns the one-dimensional piecewise linear interpolant to a function
    with given discrete data points (`xp`, `fp`), evaluated at `x`.
    Parameters
    ----------
    x : array_like
        The x-coordinates at which to evaluate the interpolated values.
    xp : 1-D sequence of floats
        The x-coordinates of the data points, must be increasing if argument
        `period` is not specified. Otherwise, `xp` is internally sorted after
        normalizing the periodic boundaries with ``xp = xp % period``.
    fp : 1-D sequence of float or complex
        The y-coordinates of the data points, same length as `xp`.
    left : optional float or complex corresponding to fp
        Value to return for `x < xp[0]`, default is `fp[0]`.
    right : optional float or complex corresponding to fp
        Value to return for `x > xp[-1]`, default is `fp[-1]`.
    period : None or float, optional
        A period for the x-coordinates. This parameter allows the proper
        interpolation of angular x-coordinates. Parameters `left` and `right`
        are ignored if `period` is specified.
        
    Returns
    -------
    y : float or complex (corresponding to fp) or ndarray
        The interpolated values, same shape as `x`.
    Raises
    ------
    ValueError
        If `xp` and `fp` have different length
        If `xp` or `fp` are not 1-D sequences
        If `period == 0`
    See Also
    --------
    scipy.interpolate
    Warnings
    --------
    The x-coordinate sequence is expected to be increasing, but this is not
    explicitly enforced.  However, if the sequence `xp` is non-increasing,
    interpolation results are meaningless.
    Note that, since NaN is unsortable, `xp` also cannot contain NaNs.
    A simple check for `xp` being strictly increasing is::
        np.all(np.diff(xp) > 0)
    """

    fp = np.asarray(fp)

    if np.iscomplexobj(fp):
        interp_func = compiled_interp_complex
        input_dtype = np.complex128
    else:
        interp_func = compiled_interp
        input_dtype = np.float64

    if period is not None:
        if period == 0:
            raise ValueError("period must be a non-zero value")
        period = abs(period)
        left = None
        right = None

        x = np.asarray(x, dtype=np.float64)
        xp = np.asarray(xp, dtype=np.float64)
        fp = np.asarray(fp, dtype=input_dtype)

        if xp.ndim != 1 or fp.ndim != 1:
            raise ValueError("Data points must be 1-D sequences")
        if xp.shape[0] != fp.shape[0]:
            raise ValueError("fp and xp are not of the same length")
        # normalizing periodic boundaries
        x = x % period
        xp = xp % period
        asort_xp = np.argsort(xp)
        xp = xp[asort_xp]
        fp = fp[asort_xp]
        xp = np.concatenate((xp[-1:]-period, xp, xp[0:1]+period))
        fp = np.concatenate((fp[-1:], fp, fp[0:1]))

    return interp_func(x, xp, fp, left, right)

2. scipy.interpolate.interp1d

使用示例如：

from scipy.interpolate import interp1d

x = [1, 2, 3]
y = [3, 2, 0]
f = interp1d(x,y,fill_value=(3,0),bounds_error=False) # 线性内插
out = f([0, 1, 1.5, 2.72, 3.14])
#array([3. , 3. , 2.5 , 0.56, 0. ])

fe = interp1d(x,y, fill_value='extrapolate') # 线性内插+外插
out = fe([0, 1, 1.5, 2.72, 3.14])
#array([ 4. , 3. , 2.5 , 0.56, -0.28])

函数定义如：

https://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.interp1d.html

scipy.interpolate.interp1d(x, y, kind='linear', axis=-1, copy=True, bounds_error=None, fill_value=nan, assume_sorted=False)
'''
一维函数插值

x和y是用于逼近函数y=f(x)的数组

返回一个函数，函数调用后通过插值方式找到新数据点的返回值。

参数：
x: 1-D数值数组，一般是升序排列的x数据点

y: N-D数值数组，与x数据点对应的y坐标，插值维的长度必须与x长度相同

kind: 字符串或整数，给出插值的样条曲线的阶数，线性插值用'linear'

bounds_error: 布尔值，越界是否报错，除非fill_value='extrapolate'，否则默认越界时报错

fill_value: 数组或'extrapolate'，指定不在x范围内时的填充值或填充方法. 当为'extrapolate'时，返回的函数会对落在x范围外的值进行外插.
'''

Python - 线性插值实现

AIHGF • 2021 年 10 月 09 日

插值问题描述：

已知坐标 (x0, y0) 和 (x1, y1)，求区间 [x0, x1] 内某一点位置 x 所对应的 y 值.

计算公式如：

$$ y = y_0 + (x-x_0) \frac{y_1-y_0}{x_1 - x_0} = y_0 + \frac{(x-x_0)y_1 - (x-x_0)y_0}{x_1 - x_0} $$

Python 中 numpy 库和scipy 库提供了线性插值的功能.

1. numpy.interp

https://numpy.org/doc/stable/reference/generated/numpy.interp.html

使用示例：

import numpy as np 
import matplotlib.pyplot as plt

xp = [1, 2, 3]
fp = [3, 2, 0]
y = np.interp(2.5, xp, fp)
#1.0

y = np.interp([0, 1, 1.5, 2.72, 3.14], xp, fp)
#array([3.  , 3.  , 2.5 , 0.56, 0.  ])

UNDEF = -99.0
y = np.interp(3.14, xp, fp, right=UNDEF)
#-99.0

#sine 函数插值
x = np.linspace(0, 2*np.pi, 10)
y = np.sin(x)
xvals = np.linspace(0, 2*np.pi, 50)
yinterp = np.interp(xvals, x, y)

plt.plot(x, y, 'o')
plt.plot(xvals, yinterp, '-x')
plt.show()

#周期 x 坐标的插值
x = [-180, -170, -185, 185, -10, -5, 0, 365]
xp = [190, -190, 350, -350]
fp = [5, 10, 3, 4]
y = np.interp(x, xp, fp, period=360)
#array([7.5 , 5.  , 8.75, 6.25, 3.  , 3.25, 3.5 , 3.75])

#复数插值Complex interpolation:
x = [1.5, 4.0]
xp = [2,3,5]
fp = [1.0j, 0, 2+3j]
y = np.interp(x, xp, fp)
#array([0.+1.j , 1.+1.5j])

函数定义如：

https://github.com/numpy/numpy/blob/v1.23.0/numpy/lib/function_base.py#L1456-L1594

@array_function_dispatch(_interp_dispatcher)
def interp(x, xp, fp, left=None, right=None, period=None):
    """
    One-dimensional linear interpolation for monotonically increasing sample points.
    Returns the one-dimensional piecewise linear interpolant to a function
    with given discrete data points (`xp`, `fp`), evaluated at `x`.
    Parameters
    ----------
    x : array_like
        The x-coordinates at which to evaluate the interpolated values.
    xp : 1-D sequence of floats
        The x-coordinates of the data points, must be increasing if argument
        `period` is not specified. Otherwise, `xp` is internally sorted after
        normalizing the periodic boundaries with ``xp = xp % period``.
    fp : 1-D sequence of float or complex
        The y-coordinates of the data points, same length as `xp`.
    left : optional float or complex corresponding to fp
        Value to return for `x < xp[0]`, default is `fp[0]`.
    right : optional float or complex corresponding to fp
        Value to return for `x > xp[-1]`, default is `fp[-1]`.
    period : None or float, optional
        A period for the x-coordinates. This parameter allows the proper
        interpolation of angular x-coordinates. Parameters `left` and `right`
        are ignored if `period` is specified.
        
    Returns
    -------
    y : float or complex (corresponding to fp) or ndarray
        The interpolated values, same shape as `x`.
    Raises
    ------
    ValueError
        If `xp` and `fp` have different length
        If `xp` or `fp` are not 1-D sequences
        If `period == 0`
    See Also
    --------
    scipy.interpolate
    Warnings
    --------
    The x-coordinate sequence is expected to be increasing, but this is not
    explicitly enforced.  However, if the sequence `xp` is non-increasing,
    interpolation results are meaningless.
    Note that, since NaN is unsortable, `xp` also cannot contain NaNs.
    A simple check for `xp` being strictly increasing is::
        np.all(np.diff(xp) > 0)
    """

    fp = np.asarray(fp)

    if np.iscomplexobj(fp):
        interp_func = compiled_interp_complex
        input_dtype = np.complex128
    else:
        interp_func = compiled_interp
        input_dtype = np.float64

    if period is not None:
        if period == 0:
            raise ValueError("period must be a non-zero value")
        period = abs(period)
        left = None
        right = None

        x = np.asarray(x, dtype=np.float64)
        xp = np.asarray(xp, dtype=np.float64)
        fp = np.asarray(fp, dtype=input_dtype)

        if xp.ndim != 1 or fp.ndim != 1:
            raise ValueError("Data points must be 1-D sequences")
        if xp.shape[0] != fp.shape[0]:
            raise ValueError("fp and xp are not of the same length")
        # normalizing periodic boundaries
        x = x % period
        xp = xp % period
        asort_xp = np.argsort(xp)
        xp = xp[asort_xp]
        fp = fp[asort_xp]
        xp = np.concatenate((xp[-1:]-period, xp, xp[0:1]+period))
        fp = np.concatenate((fp[-1:], fp, fp[0:1]))

    return interp_func(x, xp, fp, left, right)

2. scipy.interpolate.interp1d

使用示例如：

from scipy.interpolate import interp1d

x = [1, 2, 3]
y = [3, 2, 0]
f = interp1d(x,y,fill_value=(3,0),bounds_error=False) # 线性内插
out = f([0, 1, 1.5, 2.72, 3.14])
#array([3. , 3. , 2.5 , 0.56, 0. ])

fe = interp1d(x,y, fill_value='extrapolate') # 线性内插+外插
out = fe([0, 1, 1.5, 2.72, 3.14])
#array([ 4. , 3. , 2.5 , 0.56, -0.28])

函数定义如：

https://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.interp1d.html

scipy.interpolate.interp1d(x, y, kind='linear', axis=-1, copy=True, bounds_error=None, fill_value=nan, assume_sorted=False)
'''
一维函数插值

x和y是用于逼近函数y=f(x)的数组

返回一个函数，函数调用后通过插值方式找到新数据点的返回值。

参数：
x: 1-D数值数组，一般是升序排列的x数据点

y: N-D数值数组，与x数据点对应的y坐标，插值维的长度必须与x长度相同

kind: 字符串或整数，给出插值的样条曲线的阶数，线性插值用'linear'

bounds_error: 布尔值，越界是否报错，除非fill_value='extrapolate'，否则默认越界时报错

fill_value: 数组或'extrapolate'，指定不在x范围内时的填充值或填充方法. 当为'extrapolate'时，返回的函数会对落在x范围外的值进行外插.
'''

Python - 线性插值实现

1. numpy.interp

2. scipy.interpolate.interp1d

※相关文章推荐※

※最新文章推荐※

Leave a Comment Cancel reply

Python - 线性插值实现

1. numpy.interp

2. scipy.interpolate.interp1d