Search code examples
numpyscipyinterpolationarray-broadcasting

Linear interpolation of two 2D arrays

In a previous question (fastest way to use numpy.interp on a 2-D array) someone asked for the fastest way to implement the following:

np.array([np.interp(X[i], x, Y[i]) for i in range(len(X))])

assume X and Y are matrices with many rows so the for loop is costly. There is a nice solution in this case that avoids the for loop (see linked answer above).

I am faced with a very similar problem, but I am unclear on whether the for loop can be avoided in this case:

np.array([np.interp(x, X[i], Y[i]) for i in range(len(X))])

In other words, I want to use linear interpolation to upsample a large number of signals stored in the rows of two matrices X and Y. I was hoping to find a function in numpy or scipy (scipy.interpolate.interp1d) that supported this operation via broadcasting semantics but I so far can't seem to find one.

Other points:

  • If it helps, the rows X[i] and x are pre-sorted in my application. Also, in my case len(x) is quite a bit larger than len(X[i]).

  • The function scipy.signal.resample almost does what I want, but it doesn't use linear interpolation...

Solution

  • This is a vectorized approach that directly implements linear interpolation. First, for each x value and each i, j compute the weight w expressing how much of the interval (X[i, j], X[i, j+1]) is to the left of x.

    • If the entire interval is to the left of x, the weight of that interval is 1.
    • If none of the subinterval is to the left, the weight is 0
    • Otherwise, the weight is a number between 0 and 1, expressing the proportion of that interval to the left of x.

    Then the value of PL interpolant is computed as Y[i, 0] + sum of differences dY[i, j] multiplied by the corresponding weight. The logic is to follow by how much the interpolant changes from interval to interval. The differences dY = np.diff(Y, axis=1) show how much it changes over the entire interval. Multiplication by the weight prorates that change accordingly.

    Setup, with some small data arrays

    import numpy as np
X = np.array([[0, 2, 5, 6, 9], [1, 3, 4, 7, 8]])
Y = np.array([[3, 5, 2, 4, 1], [8, 6, 9, 5, 4]])
x = np.linspace(1, 8, 20)

    The computation

    dX = np.diff(X, axis=1)
dY = np.diff(Y, axis=1)
w = np.clip((x - X[:, :-1, None])/dX[:, :, None], 0, 1)
y = Y[:, [0]] + np.sum(w*dY[:, :, None], axis=1)

    Demonstration

    This is only to show that the interpolation is correct. Blue points: original data, red ones are computed.

    import matplotlib.pyplot as plt
plt.plot(x, y[0], 'ro')
plt.plot(X[0], Y[0], 'bo')
plt.plot(x, y[1], 'rd')
plt.plot(X[1], Y[1], 'bd')
plt.show()

    PL interpolants