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

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!

• You could cross-check with this answer. – Matt L. Sep 6 '20 at 16:16
• also, I'm not an audio expert, but you very rarely interpolate any time-series signals that are harmonic in nature using quadratic interpolation; that's something I've so far mostly found in image processing. Doing quadratic resampling on audio sounds like a pretty bad idea, considering it's not a linear way of interpolation, and thus you'll get new frequency components that weren't even there in the original audio. Why did you choose that? – Marcus Müller Sep 6 '20 at 16:39
• I didn't choose anything. Sorry, I didn't mean to indicate that I thought this was a good choice. Just giving context for how I happened upon this. I just see other projects using it for upsampling and I got curious about this particular implementation that I see popping up from time to time. I already understand the implications of using it for audio. – accumulator Sep 6 '20 at 16:58
• Marcus, text-book interpolation will do some combination of zero-filling and sample-dropping, also causing new frequency components to arise, to be suppressed by a linear filter. So how is cubic interpolation fundamentally different (asides from having few taps and a restricted space to choose tap weights from)? – Knut Inge Sep 8 '20 at 15:17

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.

• Thank you so much for the analysis, Olli! That’s definitely enough for me to move on from this. I learned a lot from this exploration. Cheers! – accumulator Sep 8 '20 at 17:11
• so would it make a difference if a 3rd-order Lagrange interpolation was used vs. a 3rd-order Hermite interpolation when the input is a straight-line ramp function of non-zero slope? – robert bristow-johnson Sep 9 '20 at 4:41
• @robertbristow-johnson Order $N$ Lagrange will, by definition, exactly fit any polynomial up to degree $N$, given its samples. 3rd order Hermite will exactly fit any polynomial up to degree $2$, given its samples. So both 3rd order Lagrange and 3rd oder Hermite exactly fit any straight line or any quadratic polynomial. 3rd order Lagrange further fits exactly any cubic polynomial. (Uniform sampling is assumed.) – Olli Niemitalo Sep 9 '20 at 8:27
• this all looks a lot like your pink elephant paper (or Duane's and my old interpolation paper). seems like a goofy asymmetric polynomial interpolator. looks like they're trying to get a wider bandwidth; less low-pass filtering in the baseband. and it looks like you repeated a lotta work that you (and Duane and me) did back in the previous century. good work, Olli! – robert bristow-johnson Sep 9 '20 at 10:21

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.