# Algorithm for 1d spline interpolation suitable for 8 bit microcontroler

What is a concise, fast, down to earth algorithm for doing (or closely approximating) spline interpolation on a 1d continuous stream of data?

(Edit1: The paragraph below equates to saying "the data is uniform in interval." in an awkward way.)

The data is 1d in that sampling on x is at fixed, regular intervals (likely a power of 2 constant) known well ahead of time. This makes y the only variant, which should allow for quite a bit of simplification and precomputation (LUT?).

Here is a graphical example of approximately what we're trying to do. It's our algo applied to a coarsely sampled sine function.

(Edit2: Note that this is only an example of what our algo should do with coarsely sampled sine data, however, the data we would like to process will be random in the set [0,255], most likely 3 or 4 points)

(source: mathworks.com)

Assume high accuracy is not required, but that we must compute the number of results (red dots) between samples (blue circles) in less than 8ms (~120Hz). Also assume the signal processor available is limited in power and instruction set (PIC/AVR), so has only the following relevant instructions and limitations.

• (signed + carry) Addition/subtraction instructions.
• (unsigned 8x8-to-16) Multiplication instruction. (PIC18, megaAVR only)
• Byte wide Boolean instructions (XOR, OR, NOR, AND, NAND, etc.)
• Single bit left and right logical shifts. (no barrel shifter)
• Can only execute at 2~4 MIPS

• Responses would preferentially be in pseudocode, so they are more generally applicable.
• Cheating is totally OK; it doesn't need to be perfect, just better than liner interpolation.
• Bonus points for alternatives that don't need multiplication.
• More bonus points for alternatives likely to complete in less than 1ms!

This is for an RGB mood lamp x-mas present for my sister and mommy :3, which I would do myself, but the maths for this are apparently beyond me.

Edit 12-21-2016: Better list formatting

• Do you need to calculate the output in order or can the algorithm be such that it puts those original samples sparsely into a large array in memory and filters that array in multiple passes? – Olli Niemitalo Dec 27 '16 at 8:37
• It probably could be the sparse scenario, yes. My intention is to create an endless smooth waveform with no discontinuity's having random peeks and valleys. (see the last picture in my own answer) Right now I'm taking samples from a 4 byte array of random numbers and 4x up sampling (cubic Hermite) to a 32 byte ring buffer. Every 4 samples the 4 byte array gets shifted 1 byte, and a new original sample is loaded on the freshly emptied end. – Charlie Dec 27 '16 at 13:13
• A discrete low-pass filter followed by linear interpolation should also work well. Doesn't have to be a linear-phase filter if your input is noise. – Olli Niemitalo Dec 27 '16 at 17:18

Take a look at the cubic Hermite spline. The interpolated function is continuous at the data points and the first derivative is also continuous. Away from the data points all of the derivatives are continuous.

Let's say that the function $f(x)$ is defined by equally-spaced data points for all $x$ that is an integer. This means you know the values of $f(0), f(1), f(2), ...$

Then separate $x$ into an integer and fractional parts:

$$x \triangleq n+u$$

where

$$n = \lfloor x \rfloor = \operatorname{floor}(x)$$

and

$$u = x - n \quad \text{ , } \quad 0 \le u < 1$$



\begin{align} f(n+u) & = \begin{bmatrix} 1 & u & u^2 & u^3 \\ \end{bmatrix} \cdot \begin{bmatrix} 0 & 1 & 0 & 0 \\ -\tfrac12 & 0 & \tfrac12 & 0 \\ 1 & -\tfrac52 & 2 & -\tfrac12 \\ -\tfrac12 & \tfrac32 & -\tfrac32 & \tfrac12 \\ \end{bmatrix} \cdot \begin{bmatrix} f(n-1) \\ f(n) \\ f(n+1) \\ f(n+2) \end{bmatrix} \\ \\ & = \frac12 \begin{bmatrix} -u^3 +2u^2 - u \\ 3u^3 - 5u^2 + 2 \\ -3u^3 + 4u^2 + u \\ u^3 - u^2 \end{bmatrix}^T \cdot \begin{bmatrix} f(n-1)\\f(n)\\f(n+1)\\f(n+2) \end{bmatrix} \\ \\ & = \frac12 \begin{bmatrix} u ((2-u) u-1) \\ u^2 (3 u-5)+2 \\ u ((4-3 u) u+1) \\ u^2 (u-1) \end{bmatrix}^T \cdot \begin{bmatrix} f(n-1)\\f(n)\\f(n+1)\\f(n+2) \end{bmatrix} \\ \\ & = \tfrac12 \bigg( (u^2(2-u)-u)f(n-1) \ + \ (u^2(3u-5)+2)f(n) \\ & \quad \quad \quad \quad + \ (u^2(4-3u)+u)f(n+1) \ + \ u^2(u-1)f(n+2) \bigg) \\ \end{align}

Is this enough information for how to do this in your PIC? You need to be able to separate into integer and fractional parts and you need to be able to multiply.

in my opinion, Olli's method [now moved into its own answer] is not the best way of looking at it for the OP's case of implementing this simply in a PIC. [his formulation] separates the four data points and computes four coefficients that are attached to powers of $u$. that's the way to do it if your fractional ordinate is any arbitrary value that is $0 \le u < 1$. but the OP has only a few values like $u=0, \tfrac14, \tfrac12, \tfrac34$. or maybe 8 multiples of $\tfrac18$.

so my recommendation is to compute the vaules of these four polynomials:

$$c_{-1} = \tfrac12 (-u^3 +2u^2 - u) \\ c_0 = \tfrac12 (3u^3 - 5u^2 + 2) \\ c_1 = \tfrac12 (-3u^3 + 4u^2 + u) \\ c_2 = \tfrac12 (u^3 - u^2)$$

and do that for every fractional value of $u$ (such as $u=0, \tfrac14, \tfrac12, \tfrac34$) that you will use many, many times.

then the code in the PIC only needs to implement is a dot product between the 4 data points and the selected set of coefficients:

$$f(x) = f(n+u) = c_{-1} f(n-1) + c_0 f(n) + c_1 f(n+1) + c_2 f(n+2)$$

since $c_{-1}$ and $c_2$ can be shown to always be negative for $0 < u < 1$, then put into the table their absolute values and subtract their terms:

$$f(x) = f(n+u) = c_0 f(n) + c_1 f(n+1) - (-c_2) f(n+2) - (-c_{-1}) f(n-1)$$

the coefficients stored will be 256 times bigger than their actual value (and stored as 8-bit unsigned integers) then, after multiplying and accumulating your answer (that is 256 times too big), you add 128 (for rounding) and shift right 8 bits (which is the same as taking the answer out of the higher-order byte of the product).

• but when you have negative numbers and no 2's comp multiply, you have to bias your numbers to be multiplied unless you're certain they are non-negative. is all of your plotting data positive? if not, add a bias (that you know of) and bump it up. since the coefficients are calculated in advance, you can use their abs value and add or subtract the result. otherwise, you better learn about and understand offset binary arithmetic (sometimes called "excess-K"). – robert bristow-johnson Dec 22 '16 at 2:46
• okay, so you need not offset them. then the coefficients you get from evaluating those 4 polynomials of $u$, with $u=$ 0, 0.25, 0.50, 0.75, when the coefficient is negative, then use the abs value and subtract after the multiplication. also, those coefs have fractional values, so you must represent them as some $2^n$ times their actual value. then after adding up the terms, you must shift right by $n$ bits to undo the $2^n$ factor. that is essentially how fixed-point arithmetic is done. – robert bristow-johnson Dec 22 '16 at 3:21
• @Olli Niemitalo -Just FYI, the code you've provided (which I assume is C) for the Hermite spline algo attempts to index an array with a negative number (y[-1]). This is almost always going to NOT point to valid address/data, because no valid element exists before the first element in an array. I would advise incrementing all those indexes by one, so y[1], y[2], y[3], y[4]. Or at least making note of this problem. (Q.E.D., my compiler throws an error because of this.) – Charlie Dec 23 '16 at 13:20
• ANSI C does not bounds check indices of array references. And you can define a pointer into the middle of an array, and have negative indices w.r.t. that pointer in the middle and the reference WILL point to valid data. – robert bristow-johnson Dec 23 '16 at 17:51
• how many points between your given points do you want to interpolate? if it's a small number, like 4 or 8, you should precompute those coefficients (that are the cubic polynomials of $u$ above) in advance. – robert bristow-johnson Dec 24 '16 at 5:27

This is a different way to do cubic Hermite interpolation than the one explained in Robert's answer. In his notation, we can also write:

\begin{align}f(n+u) =\, &u^3\left(-\frac{1}{2}f(n-1) + \frac{3}{2}f(n) - \frac{3}{2}f(n+1) + \frac{1}{2}f(n+2)\right)\\ +\, &u^2\left(f(n-1) - \frac{5}{2}f(n) + 2f(n+1) - \frac{1}{2}f(n+2)\right)\\ +\, &u\left(\frac{1}{2}f(n+1) - \frac{1}{2}f(n-1)\right)\\ +\, &f(n)\end{align}

My code has different variable names but does the calculation in essentially the same order. When you put the Hermite code to real use, it will sometimes address one sample (y[-1]) before the first sample in your data and one sample (y[2]) after the last sample in your data. I normally make those extra "safety" samples available in the memory just outside the array. Another warning is that in worst case cubic Hermite interpolation overshoots the original input range, say from maximum values [-128, 127] to maximum values [-159.875, 158.875] for worst-case inputs [127, -128, -128, 127] and [-128, 127, 127, -128]. This is floating point code but can be converted to fixed point.

// x = 0..1 is the fractional position.
// Interpolating between y[0] and y[1], using also y[-1] and y[2].
float c0 = y[0];
float c1 = 1/2.0*(y[1]-y[-1]);
float c2 = y[-1] - 5/2.0*y[0] + 2*y[1] - 1/2.0*y[2];
float c3 = 1/2.0*(y[2]-y[-1]) + 3/2.0*(y[0]-y[1]);
return ((c3*x+c2)*x+c1)*x+c0;


Try to implement linear interpolation first if you have trouble:

// x = 0..1 is the fractional position.
// Interpolating between y[0] and y[1].
return (y[1]-y[0])*x+y[0];


Here's vintage 1998, Pentium optimized fixed point assembly cubic Hermite interpolation code for 32-bit x86 architecture:

;8192-times oversampling Hermite interpolation of signed 8-bit integer data.
;ESI.ECX = position in memory, 32.32-bit unsigned fixed point, lowest 19 bits ignored.
;EAX = output, 24.8-bit signed fixed point.

data:
ipminus1        dd      0
ip1             dd      0
ip2             dd      0

code:
movsx   EBP, byte [ESI-1]
movsx   EDX, byte [ESI+1]
movsx   EBX, byte [ESI+2]
movsx   EAX, byte [ESI]
sal     EBX, 8
sal     EDX, 8
mov     dword [ip2], EBX
sal     EAX, 8
mov     dword [ip1], EDX
mov     EBX, EAX
sub     EAX, EDX
sal     EBP, 8

mov     [ipminus1], EBP
lea     EAX, [EAX*4+EDX]
mov     EDX, ECX
sub     EAX, EBX
shr     EDX, 19
sub     EAX, EBP
lea     EBP, [EBX*4+EBX]

imul    EAX, EDX

sar     EAX, 32-19+1
sar     EBP, 1
sub     EAX, EBP
mov     EBP, [ip1]
sub     EBP, [ipminus1]

imul    EAX, EDX

sar     EBP, 1
sar     EAX, 32-19

imul    EAX, EDX

sar     EAX, 32-19


The above methods are useful if you need to interpolate at "random" positions. If you need to evaluate the interpolation polynomial at equidistant points, there's the forward difference method. There's an article about it in Dr Dobb's. You can do it without any multiplications in the inner loop, and also the rest of the multiplications are constant multiplications that in fixed point arithmetic can be done by shifts, additions and subtractions. Here's C/C++ demonstration code using floating point numbers:

#include <stdio.h>
#include <math.h>

// Forward difference cubic Hermite interpolation

const float x[4] = {-1, 2, -3, 4}; // Input data

int main() {
const float *y = &x[1]; // Interpolate between the middle two values
const int m = 4; // Parameter: Interpolate 2^m values for each input value.
// Cubic Hermite specific:
float c0 = y[0];
float c1 = 1/2.0*(y[1]-y[-1]);
float c2 = y[-1] - 5/2.0*y[0] + 2*y[1] - 1/2.0*y[2];
float c3 = 1/2.0*(y[2]-y[-1]) + 3/2.0*(y[0]-y[1]);
// The rest works for any cubic polynomial:
float diff0 = 3*pow(2, 1 - 3*m)*c3;
float diff1 = pow(2, 1 - 2*m)*c2 + 3*pow(2, 1 - 3*m)*c3;
float diff2 = pow(2, -m)*c1 + pow(2, -2*m)*c2 + pow(2, -3*m)*c3;
float poly = c0;
for (int k = 0; k < (1<<m)+1; k++) {
printf("%d, %f\n", k, poly);
poly += diff2;
diff2 += diff1;
diff1 += diff0;
}
}


Compared to Robert's method, this is less work in total, especially if hardware multiplication is slow or unavailable. A possible advantage of Robert's method is the balanced workload per output sample. Here there is also serial dependency. For PIC it is not a problem, but with processor architectures that have more parallel execution pipelines it becomes a bottleneck. That potential problem can be alleviated by parallelizing the calculation to groups of say four output samples with independent update of their [diff1, diff2, poly] state vectors, like in this (C/C++ code):

#include <stdio.h>
#include <math.h>

// Parallelized forward difference cubic Hermite interpolation

const float x[4] = {-1, 2, -3, 4}; // Input data

struct state {
float diff1;
float diff2;
float poly;
};

int main() {
const float *y = &x[1]; // Interpolate between the middle two values
const int m = 4; // Parameter: Interpolate 2^m values for each input value.
const int n = 2; // Parameter: 2^n parallel state vectors.
// Cubic Hermite specific:
float c0 = y[0];
float c1 = 1/2.0*(y[1]-y[-1]);
float c2 = y[-1] - 5/2.0*y[0] + 2*y[1] - 1/2.0*y[2];
float c3 = 1/2.0*(y[2]-y[-1]) + 3/2.0*(y[0]-y[1]);
// The rest works for any cubic polynomial:
state states[1<<n];
float diff0 = 3*pow(2, 1 - 3*m)*c3;
float diff1 = pow(2, 1 - 2*m)*c2 + 3*pow(2, 1 - 3*m)*c3;
float diff2 = pow(2, -m)*c1 + pow(2, -2*m)*c2 + pow(2, -3*m)*c3;
states[0].poly = c0;
printf("%d, %f\n", 0, states[0].poly);
for (int k = 1; k < (1<<n); k++) {
states[k].poly = states[k-1].poly + diff2;
printf("%d, %f\n", k, states[k].poly);
diff2 += diff1;
diff1 += diff0;
}
diff0 = 3*pow(2, 3*(n-m) + 1)*c3;
for (int k = 0; k < (1<<n); k++) {
// These are polynomials in k so could also be evaluated by forward difference, avoiding multiplicaton
states[k].diff1 = pow(2, 2*(n-m) + 1)*c2 + pow(2, 1 - 3*m)*(3*(1<<3*n)*c3 + 3*(1<<2*n)*c3*k);
states[k].diff2 = pow(2, n - m)*c1 + pow(2, - 2*m)*((1<<2*n)*c2 + (1<<n+1)*c2*k) + pow(2, - 3*m)*((1<<3*n)*c3 + 3*(1<<2*n)*c3*k + 3*(1<<n)*c3*k*k);
}
for (int i = 1; i < 1<<(m-n); i++) {
for (int k = 0; k < (1<<n); k++) {
states[k].poly += states[k].diff2;
states[k].diff2 += states[k].diff1;
states[k].diff1 += diff0;
printf("%d, %f\n", (i<<n)+k, states[k].poly);
}
}
printf("%d, %f\n", 1<<m, states[0].poly + states[0].diff2);
}

• +1 because it is, indeed helpful. It is a lot more readable as it's own answer too. Thanks for your help. – Charlie Dec 26 '16 at 9:00

OK, I'm now using (abusing?) this answer as a checkpoint for the progress I'm making. Eventually, this will fill out and become a "true" answer and this header can be removed... please bear with me.

Precomputed Constants for $u$ at 1/4 fractions.

This is related to the accepted answer; the Cubic Hermite spline case. It's here, because it needs to be formatted correctly to be remotely legible.

0      ::::    0            1           0           0
0.25   ::::   -0.0703125    0.8671875   0.2265625   -0.0234375
0.5    ::::   -0.0625       0.5625      0.5625      -0.0625
0.75   ::::   -0.0234375    0.2265625   0.8671875   -0.0703125

x256   ::::    0            256         0           0
x256   ::::   -18           222         58          -6
x256   ::::   -16           144         144         -16
x256   ::::   -6            58          222         -18


Edit: Thank you Robert. You were correct, of course, there was an error. The error was in the first columns polynomial. I was cubing $u$ in the second term when I should have squared it. The table is now correct, the spreadsheet will follow.

I have a *.ods spreadsheet that I used to make this that I will relinquish on request.

So, after bashing my head on the (wonderful) answers provided so far for the last week, I digressed to a tangential algorithm, the Centripetal Catmull–Rom spline. The Wiki page has Python code which isn't very hard to get working. The code provided does almost exactly what I was asking, only with a TON of extra baggage that is not needed. I spent the better part of the night cutting and simplifying the algo, and it's getting close to perfect now.

The only thing it needs now is...

• The vectoring needs to be unfolded so elements can be processed one-by-one.
• It needs the rest of the constants precomputed down.
• A linear interpolation stage will need to be wrapped around this.

Edit: After a day or two of messing with it, I have fully ported and partly simplified the Catmull-Rom spline algo from the Python code to working PIC XC8 code. What's more, it is reasonably fast, even though it calculates the constants in real time. On a PIC18 chip (w/hardware multiplier) operating at 1 MIPS, it takes ~2.5ms to output one new point. This is about 3x faster than the absolute minimum required for 120Hz operation. 1 MIPS is mostly a worst case as that is a paltry pace for most PICs, especially PIC18's. This is perfectly functional for my needs, and more or less solves my problem/question to my satisfaction.

Here is the relevant code.

    unsigned char j;
unsigned char l = 0;
for(j = 0; j < 16; j++)
{
// (global) unsigned char y[4] = {0};
y[0] = y[1];
y[1] = y[2];
y[2] = y[3];
y[3] = randchar(); // Wrapper, limits standard rand to [0,255]
// Debug for overshoot worst case. (y[] should alternate 2xMAX, 2xMIN)
//y[3] = y[0];

//further limit our starting points to prevent overshoot
if (y[3] > (255-16)){y[3]=(255-16);}
if (y[3] < 12){y[3]=12;}

unsigned char k;
const static unsigned char c0 = 64; // amount of fixed point shift.
for(k = c0; k < c0*2; k = k+(c0/16)) {
signed int A1 = (((c0 - k) * y[0] + k * y[1]) / c0);
signed int A2 = ((((c0*2) - k) * y[1] + (k - c0) * y[2]) / c0);
signed int A3 = ((((c0*3) - k) * y[2] + (k - (c0*2)) * y[3]) / c0);

signed int B1 = ((((c0*2) - k) / 2 * A1 + k / 2 * A2) / c0);
signed int B2 = ((((c0*3) - k) / 2 * A2 + (k - c0) / 2 * A3) / c0);

// (global) unsigned char buff[256] = {0};
buff[l] = ((((c0*2) - k) * B1 + (k - c0) * B2) + (c0*16))/ c0;
l++;
}
}


Notes:

• The arrays y[] and buff[] will have to be defined somewhere.
• The arrays don't necessarily have to be global. Especially y[].
• j times k need to equal the length of buff[].
• All the math is integers only. (well... fixed point)
• The only core operators are addition/subtraction, multiplication, and division by powers of two. This should make it pretty damn fast, and simple.
• Lastly, there is still some room for simplification.

Here is a plot resulting from running the above Python code.

And here is a plot for the new C code, run on the actual PIC, for RGB LED PWM output. Note that it looks jagged because it does not (yet) have a linear interpolation stage applied to it.

• the other thing you have to do to correctly simulate fixed-point is use the floor() function to force each scaled coefficient to an integer value. – robert bristow-johnson Dec 22 '16 at 17:34
• also, instead of full screen shots, try to use the snipping tool to select just the content. or you can also save to a .png file and upload that to SE. – robert bristow-johnson Dec 22 '16 at 17:36
• @robertbristow-johnson You are, of course, totally right about using floor() on interpolation output. I was working on a revised version that did exactly that. However, it became beside the point when I realized that the whole damn thing is actually wrong( :/) It is only an odd coincidence that it outputs exactly what you would expect to see. Of course :sigh: octave would have it's own internal liner interpolation for plotting, which is obviously covering up any flaws in my own code. Thanks to this, I've more or less given up on trying to use Octave to work this out, and switched to just C. – Charlie Dec 23 '16 at 13:39
• you won't need a linear interpolation stage, if you just have enough different values of $u$ between the integer points. – robert bristow-johnson Dec 26 '16 at 17:46
• Agreed. I think the trade off between the speed vs. accuracy when comparing linear interpolation vs. spline interpolation is worth investigating a little. Excellent optimization could be had with the right balance of the two. Even so, the asm cubic Hermite I have going now is probably plenty fast enough for just about any use case. 1 iteration (generation one fractional point.) only takes ~104 machine cycles worst case, 32 every 4th iteration. So... ~85 cycles on average. – Charlie Dec 27 '16 at 9:02

Depends

Splines are good but I'm pretty sure you need division for that, which will be awkward on your PIC.

If both the original data and the interpolated are sampled at uniform intervals, than this turns simply into an up sampling problem. The way your picture looks, you just need to up-sample by a factor of 4. This can be done easily with a polyphase FIR filter that only requires multiplies and adds. On the downside, there is latency, i.e. your interpolated data will be delayed with respect to your original data. I don't know if that's okay or not.

If your output data is really just a sine wave and you simply don't know the frequency and phase (or its time variant), you could wrap a phase locked loop around it.

• BTW, i think implementing Hermite cubic spline in assembly code on a DSP is in the same ballpark of difficulty as implementing it in a PIC. are you writing C code to run in the PIC? and you don't need division at all to implement cubic Hermite splines for equally spaced (along the x-axis) data points (which is what uniform sampling is). it's a perfectly linear operation and is an approximation to $\operatorname{sinc}(\cdot)$-based interpolation. – robert bristow-johnson Dec 21 '16 at 19:54
• @Hilmar 1: Division on the pic is indeed awkward, unless it's by a power of two integer, then it resolves down to a simple right shift by the equivalent number of bits. 2: Yes, starting data comes in uniform intervals. I'll update the question to say as much. 3: Latency is not an issue, the starting data will be buffered, then new points computed ahead of time during an anticipated idle period. 8ms is still the longest it should take. 4: It's, unfortunately, NOT just a sine wave. The data will be a three-tuple stream of random data in the set [0, 255]. i.e. randomly chosen colors in RGB form. – Charlie Dec 22 '16 at 1:26
• @robert. 1: Authoring the algo in ASM probably is the same difficulty bracket for either platform, agreed. Note that a DSP would execute orders of magnitude faster, however. So simplification/optimization is essential. 2: I could go either way on the language, C or ASM. I can also inline ASM in the C. sudocode would probably be better, so everyone can benefit. 3: the data is uniform on x, so It's good that we don't need division for that. 4: sinc() interpolation sounds interesting, C has sin() built in as part of <math.h>, elaboration on the details of this point would be awesome. – Charlie Dec 22 '16 at 1:35