有没有一种有效的方法来使用零阶保持重新采样一个 numpy 数组？理想情况下具有像numpy.interp这样的签名？
我知道scipy.interpolate.interp1d，但我确信矢量化替代方案可用于处理此类情况。
问问题
3388 次
4 回答
4
麻木 < 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]
麻木> = 1.12
在 numpy v1.11 中不推荐使用浮点索引，并在 v1.12（2017 年 1 月）中删除。上面显示的代码在当前 numpy 版本中引发IndexError异常。
您可以通过使用包装器访问数组，即时将浮点索引转换为整数，从而重现旧 numpy 版本的浮点索引行为。numpy.interp当内存效率是一个问题时，这将避免使用or预先插入中间值的需要scipy.interpolate.interp1d。
于 2012-01-17T03:02:15.153 回答
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 倍。
于 2020-07-15T20:51:07.427 回答
0
聚会有点晚了，但这是我想出的：
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
于 2014-03-31T18:47:40.263 回答
0
这是一个 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)
于 2014-12-13T16:22:54.790 回答