Search code examples
pythoninterpolationsamplingsplinecubic

Cubic spline interpolation factors required to pass through original datapoints (scipy, python)

I am using scipy.interpolate.interp1d for cubic spline interpolation of a signal. While I believe interpolated signal should pass through all original data points, this is not the case when interpolating with certain factors.

e.g. if there are N samples, with N-1 spaces between samples and an interpolation factor of f, we can insert x points between the samples N*f == (N-1)*x + N. If x is not a whole number, the interpolated signal cannot pass through the original data points. As expected this is the case, code using scipy below with N = 4 and interpolation factor f of 3 or 4.

My question is A) is this correct or am I doing something wrong? and B) Is the formula above where x is a whole number a sufficient check that the original data samples will appear in the interpolated signal (or maybe there are edge cases).

Many Thanks

import scipy.interpolate
import numpy as np

# produce random data and interp
x = np.linspace(0, 2, 4)
np.random.seed(123)
y = np.random.random(4)

interp_f = scipy.interpolate.interp1d(x, y, kind='cubic')

# upsample factor 4
x_f4 = np.linspace(0, 2, 16)
y_f4 = interp_f(x_f4)

# upsample factor 3
x_f3 = np.linspace(0, 2, 12)
y_f3 = interp_f(x_f3)

print("Sample 2 in raw data: {0:.10f}, Sample 6 in interp f4: {1:.10f}, Sample 4 in interp f3: {2:.10f}".format(y[1], y_f4[5], y_f3[4]))
# Sample 2 in raw data: 0.2861393350, Sample 6 in interp f4: 0.2861393350, Sample 5 in interp f3: 0.2657521625
Solution

  • First, it is true, as you wrote, that a cubic interpolation pass through its original points. You can verify this by:

    all(y == interp_f(x))  # gives True

    Your formula for up-sampling seems a bit mixed up. It is easy to see, if we walk through an example:

    • Assume we have the interval [0, w], where w = 2.
    • Having n = 5 samples gives (n-1) intervals of width d = w/(n-1) = .5.
    • To up-sample by a factor of f=4, we except the new interval width to be d_4 = d / f = 0.125.
    • Thus d_4 = w / (n_4 - 1) needs to hold as well, where n_4 are the number samples of the up-sampled signal.
    • Exploiting 1 / d_4 = f * (n-1) / w should be equal to (n_4 - 1) / w results in n_4 = f * (n-1) + 1 = 17

    As long as f is a positive integer, the original samples will be included in the up-sampled signal (due to d_4 = d / f). We can verify our formula by:

    n, f = 5, 4
n_4 = f * (n-1) + 1
x = np.linspace(0, 2, n)
x_4 = np.linspace(0, 2, n_4)
if all(x_4[::f] == x):  # Every fourth sample is equal to the original
        print("Up-sampling works.")