Algorithm for 1d spline interpolation suitable for 8 bit microcontroler

Question

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

---

Additional Notes:

• 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

Answer 1

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) $$

$ n = \lfloor x \rfloor$ is the sole integer such that

$$ x-1 < \lfloor x \rfloor \le x < \lfloor x \rfloor + 1 $$

and the fractional part

$$ 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}^\mathrm{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}^\mathrm{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) \\ \\ & = \tfrac12 \bigg(\Big(\big(-f(n-1) + 3f(n) - 3f(n+1) + f(n+2)\big) u \\ & \quad \quad \quad \quad + \big(2f(n-1) - 5f(n) + 4f(n+1) - f(n+2)\big)\Big)u \\ & \quad \quad \quad \quad + \big(-f(n-1) + f(n+1)\big)\bigg)u \ + \ f(n)\\ \end{align}$$

where $[\cdot]^\mathrm{T}$ means the matrix transpose and $\lfloor\cdot\rfloor=\operatorname{floor}(\cdot)$ means the floor() function, the largest integer that does not exceed the argument. The last expression is using Horner's Method, which might be useful (perhaps more useful with floating-point).

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

Answer 2

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.

Here is a link to the spreadsheet. (Opens in browser)

---

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.

Answer 3

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(-\tfrac{1}{2}f(n-1) + \tfrac{3}{2}f(n) - \tfrac{3}{2}f(n+1) + \tfrac{1}{2}f(n+2)\right)\\ +\, &u^2\left(f(n-1) - \tfrac{5}{2}f(n) + 2f(n+1) - \tfrac{1}{2}f(n+2)\right)\\ +\, &u\left(\tfrac{1}{2}f(n+1) - \tfrac{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

add     EAX, [ip2]

lea     EBP, [EBX*4+EBX]
imul    EAX, EDX
sar     EAX, 32-19+1

add     EBP, [ip2]

sar     EBP, 1

add     EAX, [ip1]

add     EAX, [ip1]

add     EDI, 8

sub     EAX, EBP

mov     EBP, [ip1]

add     EAX, [ipminus1] 
sub     EBP, [ipminus1]
imul    EAX, EDX
sar     EBP, 1

sar     EAX, 32-19

add     EAX, EBP
imul    EAX, EDX
sar     EAX, 32-19

add     EAX, EBX

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.0y[0] + 2y[1] - 1/2.0y[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 = 3pow(2, 1 - 3m)c3;
  float diff1 = pow(2, 1 - 2m)c2 + 3pow(2, 1 - 3m)c3;
  float diff2 = pow(2, -m)c1 + pow(2, -2m)c2 + pow(2, -3m)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.0y[0] + 2y[1] - 1/2.0y[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 = 3pow(2, 1 - 3m)c3;
  float diff1 = pow(2, 1 - 2m)c2 + 3pow(2, 1 - 3m)c3;
  float diff2 = pow(2, -m)c1 + pow(2, -2m)c2 + pow(2, -3m)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 = 3pow(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 - 3m)(3(1<<3n)c3 + 3(1<<2n)c3k);
    states[k].diff2 = pow(2, n - m)c1 + pow(2, - 2m)((1<<2n)c2 + (1<<n+1)c2k) + pow(2, - 3m)((1<<3n)c3 + 3(1<<2n)c3k + 3(1<<n)c3k*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);
}

Answer 4

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.