5
$\begingroup$

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!

$\endgroup$
1
  • $\begingroup$ it looks like your interpolation polynomial is the Lagrange interpolation polynomial. consider Hermite (or "osculating") polynomials, instead. they will preserve continuity of as many derivatives at the splice points as possible. $\endgroup$ Feb 6, 2015 at 18:42

2 Answers 2

12
$\begingroup$

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
$\endgroup$
12
  • 1
    $\begingroup$ Note that the polyphase method is an optimization, and not necessary for low-pass filter kernel or windowed-Sinc interpolation. The optimization may even be wrong for systems with fast math but slow tiny data caches. Thus, the usual rule is to avoid premature optimization. $\endgroup$
    – hotpaw2
    Feb 5, 2015 at 17:35
  • 1
    $\begingroup$ Yes, it's just one possible implementation which is very practical when cost is determined by the number of multiplications. $\endgroup$
    – Matt L.
    Feb 5, 2015 at 17:57
  • $\begingroup$ On one of my iOS devices, I was surprised to find that a large polyphase table lookup was slower than an earlier vectorized multiple transcendental function computation (lots of multiplies). I surmised that the cache miss penalty was vast compared to a stream of pipelined FPMACs. $\endgroup$
    – hotpaw2
    Feb 5, 2015 at 18:32
  • $\begingroup$ I am restricted to the use of cubic interpolation or cubic spline interpolation technique. Once this is done, then I can pass it through a filter. Therefore, I am not looking for alternatives and looking forward to suggestions related to the cubic interpolation method. Thanks! $\endgroup$
    – PsychedGuy
    Feb 6, 2015 at 6:20
  • 1
    $\begingroup$ @DigitalGeeK: You just shift the point of reference (x-axis value $0$) to the left end of the current segment. The right end has value $x=1$, and the two points to the left and to the right of the current segment have $x$-axis values $-1$ and $2$. This results in simple algebra, but this choice is completely arbitrary. $\endgroup$
    – Matt L.
    Feb 7, 2015 at 13:59
3
$\begingroup$

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$$
$\endgroup$
1
  • $\begingroup$ Error corrected, thanks. $\endgroup$
    – Matt L.
    Jan 26, 2020 at 11:07

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.