What is the difference between cubic interpolation and cubic "Spline" interpolation?. How to use it for upsampling purpose?

Question

After considering a couple of advices and suggestions for upsampling techniques here, I finally converged to use the cubic interpolation technique to estimate the voltage values corresponding to intermediate samples present between the original or previous samples. I know that spline interpolation is basically used for getting smoother curves, but what makes it different from the normal cubic interpolation technique as both of them use a 3rd degree polynomial to estimate intermediate values?.

Another implementation issue is that, for example, If I have some voltages corresponding to some samples say

    V = 3.674, 6.791, 8.888, 9.667.....
    Sample = 2, 3, 4, 5.....

Now, If we have to find the voltage information corresponding to the intermediate sample 3.5, then using cubic polynomial method, we arrive at 4 equations and 4 unknowns obtained by using the information provided by the neighboring samples closest to sample 3.5.

    V(2) = a + 2b + 4c + 8d
    V(3) = a + 3b + 9c + 16d
    V(4) = a + 4b + 16c + 64d
    V(5) = a + 5b + 125c + 125d

So, solving these equations I arrived at

   a = -4.428
   b = 4.3756
   c = -0.06299
   d = -0.04966

Using these values , we can calculate the voltage value at the sample 3.5 as

  V(3.5) = a + 3.5b + 12.25c + 42.875d
  V(3.5) = 7.186 Volts

Now my question - Is this method of interpolation suitable for large sampling rates?. How can I use this technique for upsampling a signal say from 10 sps to 100 sps for N = 1024 sample points?. I know that I have to develop a function that performs this cubic interpolation task between original samples (in C++), but I am just wondering how to implement it for upsampling for a continuous series of samples . Any suggestions, ideas or advices regarding the topic and its implementation would be appreciated. Thanks!

I understand that cubic interpolation can operate on 4 data points and the more sophisticated technique I can think of is cubic spline. In case I am using the normal cubic interpolation, how about I loop through the "N" sample points i.e. 1024, for a condition below the "input sampling rate" i.e. 10 sps considering 4 data points each and then performing the interpolation function based on the up sampling factor between each of those 4 consecutive data points (Meaning - Interpolating/Estimating 10 values between each of those 4 data points) and then the function considers the next 4 data points to perform the same operation and it goes on until a 100 samples has bee acquired i.e. the output sampling rate!

Answer 1

The difference between cubic interpolation as described in your question and cubic spline interpolation is that in cubic interpolation you use 4 data points to compute the polynomial. There are no constraints on the derivatives. Cubic spline interpolation computes a third order polynomial only from two data points with the additional constraint that the first and second derivative at the interpolation points are continuous. So if you have 4 points, then you compute 3 different polynomials (between points 1-2, 2-3, and 3-4), and these polynomials are smoothly connected in the sense that their first and second derivatives are equal at the given data points. This is not the case with standard polynomial interpolation.

But if you want an efficient implementation I wouldn't necessarily use polynomial or spline interpolation but I would go for a standard polyphase implementation of a low-pass interpolation filter. If $L$ is the upsampling factor you get $L$ parallel polyphase filters operating at the lower sampling rate, followed by $L$ upsamplers, the results of which are added. The low-pass filter can be designed using the window method (e.g. a Kaiser window). There is a lot of material available on the polyphase implementation of interpolators. This Master's thesis is very accessible. Just have a look at single-stage polyphase implementations using linear-phase FIR filters and ignore the rest.

EDIT:

If you want or need to use cubic interpolation, here's what you need to do. First of all, some coefficients in your equations are wrong, that's why also the result is not correct. These would be the correct equations:

    V(2) = a + 2b + 4c + 8d
    V(3) = a + 3b + 9c + 27d
    V(4) = a + 4b + 16c + 64d
    V(5) = a + 5b + 25c + 125d

However, you don't want to do it this way because if you always use the sample index as the $x$-coordinate you'll end up with a very large range for the coefficients of the system of linear equations (because you need to take the third power), and your coefficient matrix will sooner or later become ill-conditioned (asking for numerical problems). Second, and even more importantly, you would get a different system of equations for every new segment of the input signal, which would be very inefficient. You can do the interpolation by solving a $4\times 4$ system of linear equations with the same coefficient matrix each time. Just the left-hand side vector (containing the values of the given data points, V in your example) will be different for each segment.

You compute a cubic polynomial for each segment (i.e. for each range between two given data points) by considering the data points defining the segment and the two adjacent data points, just as in your example. If you define the cubic polynomial as

$$P(x)=a_0+a_1x+a_2x^2+a_3x^3$$

and if you define the $x$-axis values of the 4 current input data points as $-1$, $0$, $1$, and $2$, you get the following system of equations:

$$\begin{align}y_0&=a_0-a_1+a_2-a_3\\ y_1&=a_0\\ y_2&=a_0+a_1+a_2+a_3\\ y_3&=a_0+2a_1+4a_2+8a_3\end{align}$$

where $y_i$, $k=0,1,\ldots, 3$ are the current data points, and the current segment is between the points $y_1$ and $y_2$. The system above is readily solved as

$$\begin{align} a_0&=y_1\\ a_1&=-\frac13 y_0 - \frac12 y_1 + y_2 - \frac16 y_3\\ a_2&=\frac12 (y_0 + y_2) - y_1\\ a_3&=\frac12 \left[(y_1 - y_2) + \frac13 (y_3 - y_0) \right]\\ \end{align}$$

So you just need to evaluate the above 4 expressions for the polynomial coefficients for each set of 4 data values. As soon as you have the coefficients, you can compute the interpolated output values in the current segment by evaluating the polynomial at $x_i=i/L$, $i=1,2,\ldots,L-1$, where $L$ is the upsampling factor ($L=10$ in your case). Then you forget the oldest input value, read the new input value, and you have a new set of 4 values from which you compute a new polynomial for the current segment, etc.

Here is a little Matlab/Octave script implementing this method:

x = randn(20,1);        % some input data
L = 10;                 % upsampling factor

% append points at both ends for first and last segment
N = length(x);
x2 = [x(1);x(:);x(N)];
Nseg = N-1;         % # segments

% matrix to compute polynomial coefficients from the current 4 data points: a = A*y
A = [0 1 0 0;-1/3 -.5 1 -1/6;.5 -1 .5 0;-1/6 .5 -.5 1/6];

% x-axis values of interpolated values in the current segment
mu = (0:L-1)/L;

% allocate output vector
Ny = Nseg*L + 1;
y = zeros(Ny,1);

% loop over all segments
for i = 1:Nseg,
    xt = x2(i:i+3);
    a = A*xt;
    yt = a(1) + a(2)*mu + a(3)*mu.^2 + a(4)*mu.^3;
    Iy = (i-1)*L + (1:L);
    y(Iy) = yt;
end

y(Ny) = x(N);

% plot result
tx = 0:N-1;
ty = (0:Ny-1)*(N-1)/(Ny-1);
plot(ty,y,'r-'), hold on
stem(ty,y,'r'), hold on
stem(tx,x), hold off

Answer 2

For those seeking a more DSP/FPGA friendly solution, just substitute 'a' in first equation and rearrange to get: $$ \begin{array}{ll} C_0(x) = &-\frac{1}{6}x^3+\frac{1}{2}x^2-\frac{1}{3}x \\ C_1(x) = &\frac{1}{2}x^3-x^2-\frac{1}{2}x+1 \\ C_2(x) = &-\frac{1}{2}x^3 + \frac{1}{2}x^2+x \\ C_3(x) = &\frac{1}{6}x^3 - \frac{1}{6}x \end{array} $$

And interpolated value: $$ P(x) = \sum_{i=0}^{3} C_i y_i \textrm{ for } x=[0\ldots1) $$

Pre compute required coefficients for each phase and multiply selected phase with data points.

Example coefficients for upsampling factor 10:

-0.0000  1.0000 -0.0000  0.0000
-0.0285  0.9405  0.1045 -0.0165
-0.0480  0.8640  0.2160 -0.0320
-0.0595  0.7735  0.3315 -0.0455
-0.0640  0.6720  0.4480 -0.0560
-0.0625  0.5625  0.5625 -0.0625
-0.0560  0.4480  0.6720 -0.0640
-0.0455  0.3315  0.7735 -0.0595
-0.0320  0.2160  0.8640 -0.0480
-0.0165  0.1045  0.9405 -0.0285

---

Matt L. got error in his equations (but not in a code). a2 should be like this $$a_2 = \frac{1}{2}(y_0+y_2)-y1$$