Search code examples
pythonnumpyscipy

How to interpolate over a piecewise-constant set of points with scipy or numpy?


I am having troubles in interpolating a set of piecewise-constant points. To better explain, consider the following code:

import numpy as np  
from scipy.interpolate import UnivariateSpline
import matplotlib.pyplot as plt


N = 30
Ts = 0.2

times = np.arange(0.0, N * Ts, Ts)
signal = np.concatenate((1.5 * np.ones(10), -0.5 * np.ones(10), np.ones(10)))
interp = UnivariateSpline(times[::3], signal[::3])

signal_interp = interp(times)

plt.plot(times, signal)
plt.plot(times, signal_interp)
plt.show()

I expected the plot of signal_interp to overlap the plot of signal but that is not happening.

I have found something here but the iterp1d is considered legacy.

I also tried to look here but I could not find any guidance for my specific case.


Solution

  • interp is a function of times, not a function of signal.

    You meant
    signal_interp = interp(times)
    not
    signal_interp = interp(signal)

    Which gives

    enter image description here

    It still doesn't "overlap" the initial signal. But splines cannot take the shape of piecewise-constant function. They are smooth by nature. It is what they do. They are used precisely for not being piecewise-constant. The reason why they appeared to be able to produce a piece-wise constant behavior with your initial code, is because what you were plotting were just 3 points of the spline, repeated 10 times each (since the abscissa used was signal, which contains only 3 values, repeated 10 times).

    Nevertheless, you can still more, by reducing smooth factor s

    interp = UnivariateSpline(times[::3], signal[::3], s=0.1)
    plt.plot(times, signal)
    plt.plot(times, interp(times))
    

    enter image description here

    np interp

    For a simple linear interpolation between each pair of consecutive point (so the less smooth possible),

    signal_interp = np.interp(times, times[::3], signal[::3])
    plt.plot(times, signal)
    plt.plot(times, signal_interp)
    

    enter image description here