A cubic interpolation function: folkloric copypasta or clever trade-off?

Question

I've been reading on interpolation methods recently and I have come across an implementation of cubic interpolation that is leaving my head scratching. Every other variant and example of cubic interpolation I've come across has included lots of fractions and exponents. This implementation I found on Paul Bourke's website uses only additions and subtractions to figure the coefficients. I've seen this implementation used in several codebases but unfortunately without any documentation.

double Cubic(double x, double a, double b, double c, double d) {
    double A = d - c - a + b;
    double B = a - b - A;
    double C = c - a;
    double D = b;
    return A * (x * x * x) + 
           B * (x * x) +
           C * x +
           D;
}

If I may, let me first walk through how this implementation might have come to be just to make sure I at least understand the algebra going on here. But I'm curious as to why it can be made so simple.

As I understand it, cubic interpolation takes four sample points, say $a$, $b$, $c$, and $d$, and interpolates between b and c using the function:

$$ f(x) = Ax^3 + Bx^2 + Cx + D \tag{1}\label{1} $$

Where $0 \le x \le 1$.

A major condition I often see (one that's obviously important for me) is that $f(x)$ must run through the sample points. So, in this case, samples $b$ and $c$ correspond to $x = 0$ and $x = 1$ respectively.

$$ f(0) = b \tag{2}\label{2} $$

$$ f(1) = c \tag{3}\label{3} $$

There is another condition that says that the derivative at samples $b$ and $c$ are the same as the slope between their surrounding points. I believe this is called a Catmull-Rom spline?

$$ f'(x) = 3Ax^2 + 2Bx + C \tag{4}\label{4} $$

$$ f'(0) = \frac{c - a}{2} \tag{5}\label{5} $$

$$ f'(1) = \frac{d - b}{2} \tag{6}\label{6} $$

Now, if I follow that path, I get Paul Breeuwsma's solution here and that all makes sense to me.

However, in order to get Bourke's implementation, I have to multiply the slopes at $b$ and $c$ by 2! That is:

$$ \begin{align} f'(0) = c - a \tag{5a}\label{5b}\\ f'(1) = d - b \tag{6a}\label{6b} \end{align} $$

If I follow from that, I get the below for the coefficients $A$, $B$, $C$, and $D$ in order to match Bourke's.

First, $D$:

$$ \require{cancel} $$

$$ \begin{align} f(0) &= \cancel{A(0)^3} + \cancel{B(0)^2} + \cancel{C(0)} + D \\ f(0) &= D \\ D = b \tag{7}\label{7} \end{align} $$

Now, $C$:

$$ \begin{align} f'(0) &= \cancel{3A(0)^2} + \cancel{2B(0)} + C \\ f'(0) &= C \\ C = c - a \tag{8}\label{8} \end{align} $$

Now, $B$:

$$ \begin{align} f(1) &= A(1)^3 + B(1)^2 + C(1) + D \\ f(1) &= A + B + C + D \\ c &= A + B + c - a + b \\ B = a - b - A \tag{9}\label{9} \end{align} $$

Finally, $A$:

$$ \begin{align} f'(1) &= 3A(1)^2 + 2B(1) + C \\ d - b &= 3A + 2(a - b - A) + c - a \\ d - b &= 3A + 2a - 2b - 2A + c - a \\ d - b &= A + a - 2b + c \\ A = d - c + -a + b \tag{10}\label{10} \end{align} $$

Ok, so... why the heck? This looks like nothing I've seen while reading up on this.

By doubling the slope at $b$ and $c$, besides having a fast and definitely appealing implementation, is this a trade-off or are there additional conditions that were able to be made that I missed? It seems like this could cause the interpolant to shoot around much more but I haven't yet coded up a comparison to look at and listen to. I've been scouring articles, blog posts, and papers and I can't seem to map anything to this. If anything, other implementations end up quite hairy and difficult for me to understand. Is there a name for this kind of interpolation? Where did it come from?

I know this was a long question. Thanks for reading!

Answer 1

A reasonable, analytically derived cubic interpolation method will give a straight line if the inputs are on a straight line. With $a = 0$, $b = 1$, $c = 2$, $d = 3$ you get from cubic_bourke $A = d - c - a + b = 3 - 2 - 0 + 1 = 2$ as the coefficient for $x^3$ whereas for a straight line you would have $A = 0$ (and $B = 0$ as the coefficient for $x^2$). So the method doesn't give a straight line for straight-line input. Instead it will, as you say, wobble. In the frequency response that will manifest itself as a low multiplicity of the roots at multiples of the sampling frequency.

The impulse response is not symmetrical. Normally splines have a symmetrical impulse response. The asymmetry will result in phase non-linearity:

Figure 1. Impulse response of cubic_bourke.

The magnitude frequency response shows a rather nice treble response for audio purposes, but gives a lot of aliasing, worse than linear interpolation:

Figure 2. Magnitude frequency rseponse of cubic_bourke (red), linear interpolation (blue) and cubic Hermite spline (green). Half the input sampling frequency is at $\omega = \pi$.

Cubic Hermite (Catmull-Rom) will give better treble response up to 15 kHz for a 44.1 kHz input sampling frequency, and lower aliasing, than cubic_bourke.

So I say scrap it.

Answer 2

I'm going to chalk this up as an optimization with a trade-off. I wrote a program to plot the results of the two interpolation methods, the one I inlined above (Bourke) and the other one that I originally derived (Breeuwsma).

The interpolating functions used:

# https://www.paulinternet.nl/?page=bicubic

def cubic_breeuwsma(x, a, b, c, d): 
  return b + 0.5 * x*(c - a + x*(2.0*a - 5.0*b + 4.0*c - d + x*(3.0*(b - c) + d - a)))

# http://paulbourke.net/miscellaneous/interpolation/

def cubic_bourke(x, a, b, c, d): 
  A = d - c - a + b 
  B = a - b - A 
  C = c - a 
  D = b 
  return A * (x**3) + B * x**2 + C * x + D 

My observations from playing with various wav files:

• Bourke does indeed overshoot minima and maxima and "wobble" through steep slopes (especially visible in the 8x figure)
• Still, they are actually more similar than I thought

I guess I can see why a lot of projects "approximate" with this method but I hope to still find out how someone came up with the idea of simplifying/overestimating the slope of the surrounding points. Probably just a curiosity to most folks but this was pretty fascinating to compare. My kingdom for a documentation comment...

Attached are some plots of 2x, 4x, and 8x upsampling comparing the two methods using the beginning of a kick drum wave file, 16-bit @ 44.1kHz.