ENH: Allow interp to process multiple signals by passing multidimensional fp

[tim994018]

Proposed new feature or change:

Consider multiple 1D signals with identical x-coordinate sample points. It is sometimes desired to resample all of these signals to the same interpolated x-coordinates via linear interpolation. Currently, this requires multiple calls to numpy.interp, one for each signal. This is wasteful, because it requires re-calculation of interpolation weights.

For example, consider original time sample coordinates of t, and updated time sample coordinates of t_new. Also, consider three signals sampled in this way: s0, s1, and s2. To resample all three with current methods:

s0_resamp = np.interp(t_new, t, s0)
s1_resamp = np.interp(t_new, t, s1)
s2_resamp = np.interp(t_new, t, s2)

It would be helpful to do this in a single call, like this:

s = np.row_stack((s0, s1, s2))
s_new = np.interp(t_new, t, s)

In this case, the number of rows in s_new would be 3 (the number of 1D signals). The number of columns in s_new would be len(t_new).

In general, fp could be any number of dimensions. Interpolation could be performed along the final dimension by default, with an additional optional input to specify a different dimension for interpolation.

Note: Matlab's interp1 function offers this functionality, as noted in the function description: "If you have multiple sets of data that are sampled at the same point coordinates, then you can pass v as an array."

[WarrenWeckesser]

FYI, I implemented the function linear_interp1d(x, xp, fp) as a gufunc in ufunclab.

In [1]: import numpy as np

In [2]: from ufunclab import linear_interp1d

In [3]: xp = np.array([0, 1, 2, 3, 5, 8])

In [4]: fp = np.array([[1.0, 1.0, 2.0, 2.0, 3.0, 3.0],
   ...:                [9.0, 8.0, 7.0, 6.0, 5.0, 4.0],
   ...:                [0.0, 2.0, 4.0, 5.0, 4.5, 4.0]])

Interpolate at a single point:

In [5]: linear_interp1d(1.5, xp, fp)
Out[5]: array([1.5, 7.5, 3. ])

Interpolate at a collection of x values. The x input must be reshaped to shape (4, 1) to broadcast correctly with fp:

In [6]: x = np.array([0.5, 1.0, 1.5, 7])

In [7]: linear_interp1d(x[:,np.newaxis], xp, fp)
Out[7]: 
array([[1.        , 8.5       , 1.        ],
       [1.        , 8.        , 2.        ],
       [1.5       , 7.5       , 3.        ],
       [3.        , 4.33333333, 4.16666667]])

[ngoldbaum]

If I'm reading the SciPy docs correctly, scipy.interpolate.interp1d accepts an array of 1D signals.

[tim994018]

Thanks ngoldbaum ... "scipy.interpolate.interp1d" is marked as legacy so I'd prefer to not rely on it. Its documentation points users to np.interp.

[ev-br]

Drive-by comment: a possible non-legacy scipy incantation would be make_interp_spline(xp, fp, k=1) (xnew).