高效的 numpy 零阶保持

Question

有没有一种有效的方法使用零阶保持对 numpy 数组进行重新采样？理想情况下，签名类似于 numpy.interp ？

我知道 scipy .interpolate.interp1d，但我确信矢量化替代方案可用于处理此类情况。

Answer 1

numpy < 1.12

由于您不会插入任何新值，因此最有效的方法是保持原始数组不变并用浮点数对其进行索引。这实际上是零阶保持。

>>> import numpy as np
>>> A = np.array(range(10))
>>> [A[i] for i in np.linspace(0, 9, num=25)]
[0, 0, 0, 1, 1, 1, 2, 2, 3, 3, 3, 4, 4, 4, 5, 5, 6, 6, 6, 7, 7, 7, 8, 8, 9]

Numpy >= 1.12

numpy v1.11 和 中已弃用浮点数索引在 v1.12 中删除（2017 年 1 月）。上面显示的代码在当前 numpy 版本中引发 IndexError 异常。

您可以通过使用包装器访问数组来重现旧 numpy 版本的浮点索引行为，将浮点索引动态转换为整数。当考虑内存效率时，这将避免使用 numpy.interp 或 scipy.interpolate.interp1d。

Answer 2

Numpy 数组不再允许非整数索引，因此接受的解决方案不再有效。 scipy.interpolate.interp1d 很强大，但在许多情况下不是最快的。这是一个强大的代码解决方案，可以时使用 np.searchsorted，不能时则恢复到 scipy 版本。

import numpy as np
from scipy.interpolate import interp1d

def zero_order_hold(x, xp, yp, left=np.nan, assume_sorted=False):
    r"""
    Interpolates a function by holding at the most recent value.

    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.
    yp: 1-D sequence of float or complex
        The y-coordinates of the data points, same length as xp.
    left: int or float, optional, default is np.nan
        Value to use for any value less that all points in xp
    assume_sorted : bool, optional, default is False
        Whether you can assume the data is sorted and do simpler (i.e. faster) calculations

    Returns
    -------
    y : float or complex (corresponding to fp) or ndarray
        The interpolated values, same shape as x.

    Notes
    -----
    #.  Written by DStauffman in July 2020.

    Examples
    --------
    >>> import numpy as np
    >>> xp = np.array([0., 111., 2000., 5000.])
    >>> yp = np.array([0, 1, -2, 3])
    >>> x = np.arange(0, 6001, dtype=float)
    >>> y = zero_order_hold(x, xp, yp)

    """
    # force arrays
    x  = np.asanyarray(x)
    xp = np.asanyarray(xp)
    yp = np.asanyarray(yp)
    # find the minimum value, as anything left of this is considered extrapolated
    xmin = xp[0] if assume_sorted else np.min(xp)
    # check that xp data is sorted, if not, use slower scipy version
    if assume_sorted or np.all(xp[:-1] <= xp[1:]):
        ix = np.searchsorted(xp, x, side='right') - 1
        return np.where(np.asanyarray(x) < xmin, left, yp[ix])
    func = interp1d(xp, yp, kind='zero', fill_value='extrapolate', assume_sorted=False)
    return np.where(np.asanyarray(x) < xmin, left, func(x))

这通过了一系列测试用例：

import unittest

import numpy as np
from scipy.interpolate import interp1d

class Test_zero_order_hold(unittest.TestCase):
    r"""
    Tests the zero_order_hold function with the following cases:
        Subsample high rate
        Supersample low rate
        xp Not sorted
        x not sorted
        Left extrapolation
        Lists instead of arrays

    Notes
    -----
    #.  Uses scipy.interpolate.interp1d as the gold standard (but it's slower)
    """
    def test_subsample(self):
        xp = np.linspace(0., 100*np.pi, 500000)
        yp = np.sin(2 * np.pi * xp)
        x  = np.arange(0., 350., 0.1)
        func = interp1d(xp, yp, kind='zero', fill_value='extrapolate', assume_sorted=True)
        y_exp = func(x)
        y = zero_order_hold(x, xp, yp)
        np.testing.assert_array_equal(y, y_exp)
        y = zero_order_hold(x, xp, yp, assume_sorted=True)
        np.testing.assert_array_equal(y, y_exp)

    def test_supersample(self):
        xp = np.array([0., 5000., 10000., 86400.])
        yp = np.array([0, 1, -2, 0])
        x  = np.arange(0., 86400.,)
        func = interp1d(xp, yp, kind='zero', fill_value='extrapolate', assume_sorted=True)
        y_exp = func(x)
        y = zero_order_hold(x, xp, yp)
        np.testing.assert_array_equal(y, y_exp)
        y = zero_order_hold(x, xp, yp, assume_sorted=True)
        np.testing.assert_array_equal(y, y_exp)

    def test_xp_not_sorted(self):
        xp    = np.array([0, 10, 5, 15])
        yp    = np.array([0, 1, -2, 3])
        x     = np.array([10, 2, 14,  6,  8, 10, 4, 14, 0, 16])
        y_exp = np.array([ 1, 0,  1, -2, -2,  1, 0,  1, 0,  3])
        y     = zero_order_hold(x, xp, yp)
        np.testing.assert_array_equal(y, y_exp)

    def test_x_not_sorted(self):
        xp    = np.array([0, 5, 10, 15])
        yp    = np.array([0, -2, 1, 3])
        x     = np.array([10, 2, 14,  6,  8, 10, 4, 14, 0, 16])
        y_exp = np.array([ 1, 0,  1, -2, -2,  1, 0,  1, 0,  3])
        y     = zero_order_hold(x, xp, yp)
        np.testing.assert_array_equal(y, y_exp)

    def test_left_end(self):
        xp    = np.array([0, 5, 10, 15, 4])
        yp    = np.array([0, 1, -2, 3, 0])
        x     = np.array([-4, -2, 0, 2, 4, 6])
        y_exp = np.array([-5, -5, 0, 0, 0, 1])
        y     = zero_order_hold(x, xp, yp, left=-5)
        np.testing.assert_array_equal(y, y_exp)

    def test_lists(self):
        xp    = [0, 5, 10, 15]
        yp    = [0, 1, 2, 3]
        x     = [-4, -2, 0, 2, 4, 6, 20]
        y_exp = [-1, -1, 0, 0, 0, 1, 3]
        y     = zero_order_hold(x, xp, yp, left=-1)
        np.testing.assert_array_equal(y, y_exp)

前两个测试的速度测试表明，numpy 解决方案对于高速率数据的二次采样大约快 100 倍，对于低速率数据的超采样大约快 3 倍。

Answer 3

聚会有点晚了，但这是我的想法：

from numpy import zeros, array, sign

def signal_zoh(x,y,epsilon = 0.001):
    """
        Fills in the data from a Zero-Order Hold (stair-step) signal
    """
    deltaX = array(x[1:],dtype='float') - x[:-1]
    fudge = min(deltaX) *epsilon
    retX = zeros((len(x)*2-1,))
    retY = zeros((len(y)*2-1,))
    retX[0::2] = x
    retX[1::2] = x[1:]+fudge*sign(deltaX)
    retY[0::2] = y
    retY[1::2] = y[:-1]
    return retX,retY

Answer 4

这是 numpy 免费版本，具有相同的签名。数据需要按递增顺序排列 - b/c 当您通过“巧妙”使用列表作为嵌套函数默认值时，它们会被转储（加速 100 倍）：

def interp0(x, xp, yp):
    """Zeroth order hold interpolation w/ same
    (base)   signature  as numpy.interp."""
    from collections import deque

    def func(x0, xP=deque(xp), yP=deque(yp)):
      if x0 <= xP[0]:
        return yP[0]
      if x0 >= xP[-1]:
        return yP[-1]    
      while x0 > xP[0]:
        xP.popleft()     # get rid of default
        y = yP.popleft() # data as we go.
      return y

return map(func, x)