Proper Implementation of Cubic Spline Interpolation

Cubic spline implementation in Matlab

I think you are confused about what you are interpolating. You should interpolate x1 and y1 separately, and afterwards plot them against each other. The following example produces a smooth curve:

x1=[1000 3000 5000 6000 4000 2000 500 0   -1000 -3000 -5000 -6000 -4000 -2000 -500 1   999 2999 4999];
y1=[5000 4999 4990 3500 2500 2499 2498 2497 2496 2495  2494  1000  -1000 -999 -998 -997 -996 -995 -994];

s = [0,cumsum(sqrt(diff(x1).^2+diff(y1).^2))]
N = length(s);

figure();
plot(x1,y1);
hold on

s_fine = interp1(linspace(0,1,N),s,linspace(0,1,5*N));
pcx=interp1(s,x1,s_fine,'spline');
pcy=interp1(s,y1,s_fine,'spline');
plot(pcx,pcy,'r-');

How to perform cubic spline interpolation in python?

Short answer:

from scipy import interpolate

def f(x):
    x_points = [ 0, 1, 2, 3, 4, 5]
    y_points = [12,14,22,39,58,77]

    tck = interpolate.splrep(x_points, y_points)
    return interpolate.splev(x, tck)

print(f(1.25))

Long answer:

scipy separates the steps involved in spline interpolation into two operations, most likely for computational efficiency.

1. The coefficients describing the spline curve are computed,
   using splrep(). splrep returns an array of tuples containing the
   coefficients.

2. These coefficients are passed into splev() to actually
   evaluate the spline at the desired point x (in this example 1.25).
   x can also be an array. Calling f([1.0, 1.25, 1.5]) returns the
   interpolated points at 1, 1.25, and 1,5, respectively.

This approach is admittedly inconvenient for single evaluations, but since the most common use case is to start with a handful of function evaluation points, then to repeatedly use the spline to find interpolated values, it is usually quite useful in practice.

understand implementation of cubic interpolation

The sample codes in the links you posted are all cubic interpolations either from 4 points or from 2 points and 2 derivatives.

The 'mu' is the parameter at which you want to evaluate the y value. The input y0, y1,y2 and y3 are the y coordinates of the 4 input points. You can do the same to evaluate x value at parameter 'mu' by passing x0, x1, x2 and x3 to the function.

You said you want to have a function that will take x value as parameter input and get back y value as output. This is only possible when you represent your curve in the explicit form, i.e., y=f(x). However, explicit form is not capable of representing curves with multiple y values with the same x value. So, they are not good ways of representing curves in general. Parametric form is a better way of representing a curve. For example, a 2D curve C(t) is represented as C(t) = (x(t), y(t)) or a 3D curve can be represented as C(t)=(x(t), y(t), z(t)). Representing a full circle as (x(t), y(t))=(rcos(t),rsin(t)) is a good example of parametric form.

If you want to draw a smooth curve connecting a series of data points, the easiest way is to implement the Catmull Rom spline or Overhauser spline. Both of them are special form of cubic Hermite curve. Read the Wikipedia page about cubic Hermite curve here and go to the Catmull Rom spline section within the same page and you shall be able to know how to create a Catmull Rom spline interpolating the data points and how to evaluate it to get (x, y) value to draw the spline in your canvas.

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

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.")

Cubic Spline Python code producing linear splines

Ok, got this working. The problem was in my implementation. I got it working with a different approach, where the splines are constructed individually instead of continuously. This is fully functioning cubic spline interpolation by method of first constructing the coefficients of the spline polynomials (which is 99% of the work), then implementing them. Obviously this is not the only way to do it. I may work on a different approach and post that if there is interest. One thing that would clarify the code would be an image of the linear system that is solved, but i can't post pics until my rep gets up to 10. If you want to go deeper into the algorithm, see the text book link in the comments above.

import matplotlib.pyplot as plt
from pylab import arange
from math import e
from math import pi
from math import sin
from math import cos
from numpy import poly1d

# need some zero vectors...
def zeroV(m):
    z = [0]*m
    return(z)

#INPUT: n; x0, x1, ... ,xn; a0 = f(x0), a1 =f(x1), ... , an = f(xn).
def cubic_spline(n, xn, a):
"""function cubic_spline(n,xn, a, xd) interpolates between the knots
   specified by lists xn and a. The function computes the coefficients
   and outputs the ranges of the piecewise cubic splines."""        

    h = zeroV(n-1)

    # alpha will be values in a system of eq's that will allow us to solve for c
    # and then from there we can find b, d through substitution.
    alpha = zeroV(n-1)

    # l, u, z are used in the method for solving the linear system
    l = zeroV(n+1)
    u = zeroV(n)
    z = zeroV(n+1)

    # b, c, d will be the coefficients along with a.
    b = zeroV(n)     
    c = zeroV(n+1)
    d = zeroV(n)    

    for i in range(n-1):
        # h[i] is used to satisfy the condition that 
        # Si+1(xi+l) = Si(xi+l) for each i = 0,..,n-1
        # i.e., the values at the knots are "doubled up"
        h[i] = xn[i+1]-xn[i]  

    for i in range(1, n-1):
        # Sets up the linear system and allows us to find c.  Once we have 
        # c then b and d follow in terms of it.
        alpha[i] = (3./h[i])*(a[i+1]-a[i])-(3./h[i-1])*(a[i] - a[i-1])

    # I, II, (part of) III Sets up and solves tridiagonal linear system...
    # I   
    l[0] = 1      
    u[0] = 0      
    z[0] = 0

    # II
    for i in range(1, n-1):
        l[i] = 2*(xn[i+1] - xn[i-1]) - h[i-1]*u[i-1]
        u[i] = h[i]/l[i]
        z[i] = (alpha[i] - h[i-1]*z[i-1])/l[i]

    l[n] = 1
    z[n] = 0
    c[n] = 0

    # III... also find b, d in terms of c.
    for j in range(n-2, -1, -1):      
        c[j] = z[j] - u[j]*c[j+1]
        b[j] = (a[j+1] - a[j])/h[j] - h[j]*(c[j+1] + 2*c[j])/3.
        d[j] = (c[j+1] - c[j])/(3*h[j]) 

    # Now that we have the coefficients it's just a matter of constructing
    # the appropriate polynomials and graphing.
    for j in range(n-1):
        cub_graph(a[j],b[j],c[j],d[j],xn[j],xn[j+1])

    plt.show()

def cub_graph(a,b,c,d, x_i, x_i_1):
    """cub_graph takes the i'th coefficient set along with the x[i] and x[i+1]'th
       data pts, and constructs the polynomial spline between the two data pts using
       the poly1d python object (which simply returns a polynomial with a given root."""

    # notice here that we are just building the cubic polynomial piece by piece
    root = poly1d(x_i,True)
    poly = 0
    poly = d*(root)**3
    poly = poly + c*(root)**2
    poly = poly + b*root
    poly = poly + a

    # Set up our domain between data points, and plot the function
    pts = arange(x_i,x_i_1, 0.001)
    plt.plot(pts, poly(pts), '-')
    return

If you want to test it, here's some data you can use to get started, which comes from the
function 1.6e^(-2x)sin(3*pi*x) between 0 and 1:

# These are our data points
x_vals = [0, 1./6, 1./3, 1./2, 7./12, 2./3, 3./4, 5./6, 11./12, 1]

# Set up the domain
x_domain = arange(0,2, 1e-2)

fx = zeroV(10)

# Defines the function so we can get our fx values
def sine_func(x):
    return(1.6*e**(-2*x)*sin(3*pi*x))

for i in range(len(x_vals)):
    fx[i] = sine_func(x_vals[i])

# Run cubic_spline interpolant.
cubic_spline(10,x_vals,fx)

Dividing cubic spline computation

You can start the interpolation at any time, using the available points (presumably in the correct order !). Cubic spline interpolation is a very stable process and when you add more points, most of the curve remains unchanged.

If your concern is that you want to avoid redoing the whole computation several times, I guess that it is enough to work on several sections with some overlap (say 20 points) and discard the results of the 10 extreme points of all sections.

---

Related Topics