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