Interpolation

When we use phase truncation on a phasor whose frequency is low relative the tablesize and sampling rate, consecutive samples may not change. To counter this, we can use different methods of interpolation to more accurately represent the waveforms when the periods of the phasor are very small compared to the size of the table.

It is essentially the ratio between the period of the phasor and size of the table…

Linear Interpolation

Linear interpolation involves using the fractional part of the calculated phase to calculate the “distance” between the two points on which we lie. If, for instance, our calculated index is 32.76, we would need to look at both indexes 32 and 33. We would then take their respective values, `t[32]` and `t[33]`, calculate the difference between the latter index and the earlier, multiply by the fractional part, then add it to the index as we would in phase truncation:

 // calculate the fractional part
frac = index - int(index);

// get the difference and wrap the index (can also use a modulo)
while(index > tablesize)
    index = index - tablesize;

if (index == tablesize-1) 
    diff = t[0] - t[index]; // wrap
else
    diff = t[index+1] - t[index]; // no need to wrap

// get the interpolated output
out = t[index] + (diff * frac);

where `t` is the table filled with the sine wave and `index` is the calculated index with a phasor multiplied by the size of the table minus 1.

Hermite (Cubic) Interpolation

A more accurate interpolation technique is the Hermite or cubic method. What this method essentially does is to find the coefficients a, b, c, and d for the third degree polynomial f(x) and its derivative, f'(x), at x = 0 and x = 1.

\begin{aligned}& f(x) = ax^3 + bx^2 + cx + d \\ & f'(x) = 3ax^2 + 2bx + c \end{aligned}

Calculating the Coefficients

Given four points, (x0, y0), (x1, y1), (x2, y2), (x3, y3), we can calculate the coefficients with the following:

\begin{aligned} & a = -\frac{y_{0}}{2} + \frac{3y_{1}}{2} - \frac{3y_{2}}{2} + \frac{y_{3}}{2} \\ & b = y_{0} - \frac{5y_{1}}{2} + 2y_{2} - \frac{y_{3}}{2} \\ & c = -\frac{y_{0}}{2} + \frac{y_{2}}{2} \\ & d = y_{1} \end{aligned}
This is called the Catmull-Rom spline.

Calculating these coefficients and plugging them into the third degree polynomial (the first equation) gives us the cubic interpolative values between x = 0 and x = 1.

Putting it Together

Given an index i, we need to get the indices i-1, i, i+1, and i+2 to use them in our cubic interpolator.

float cubicInterpolator(float idx, float *table)
{
    long trunc = (long) idx; // truncate the index but don't overwrite
    float frac = idx - trunc; // get the fractional part

    // get the indices
    int x0 = wrap(trunc-1, tableSize-1);
    int x1 = wrap(trunc, tableSize-1);
    int x2 = wrap(trunc+1, tableSize-1);
    int x3 = wrap(trunc+2, tableSize-1);

    // calculate the coefficients
    float a = -0.5*table[x0] + 1.5*table[x1] - 1.5*table[x2] + 0.5*table[x3];
    float b = table[x0] - 2.5*table[x1] + 2.0*table[x2] - 0.5*table[x3];
    float c = -0.5*table[x0] + 0.5*table[x2];
    float d = table[x1];

    return (a*(frac*frac*frac) + b*(frac*frac) + c*frac + d);
}

Phase Truncation, Linear Interpolation, and Cubic Interpolation Compared

If we push this to an extreme to show the differences, this plot shows samples 32 to 64 of the same 100Hz sine wave lookup sampled at 8kHz using either phase truncation (purple), linear interpolation (green), or cubic interpolation (blue) from a tablesize of 16:

Listen

Listen to each of the different interpolation methods below. Be particularly aware of the differences at extremely high and low frequencies.