Interpolation

As mentioned before, interpolation schemes are often used to increase the accuracy of table lookup. Here we will give a somewhat simplified account of the effects of table sizes and interpolation schemes on the result of table lookup.

To speak of error in table lookup, we must view the wavetable as a sampled version of an underlying function--and when we ask for a value of the underlying function which lies between the points of the wavetable, the error is the difference between the result of the wavetable lookup and the ``ideal" value of the function at that point. The most revealing study of wavetable lookup error assumes that the underlying function is the SINUSOID of chapter 1. We can then understand what happens to other wavetables by considering them as superpositions (sums) of sinusoids.

The accuracy of lookup from a wavetable containing a real sinusoid depends on two factors: the quality of the interpolation scheme, and the period of the sinusoid. In general, the longer the period of the sinusoid, the more accurate the result.

In the case of a synthetic wavetable, we might know its sinusoidal components from having specified them--in which case the issue becomes one of choosing a wavetable size appropriately, when calculating the wavetable, to match the interpolation algorithm and meet the desired standard of accuracy. In the case of recorded sounds, the accuracy analysis might lead us to adjust the sample rate of the recording, either at the outset or else by resampling later.

Interpolation error for a sinusoidal wavetable can have two components: first, the continuous signal (the theoretical result of reading the wavetable continuously in time, as if the output sample rate were infinite) might not be a pure sinusoid; and second, the amplitude might be wrong. (It is possible to get phase errors as well, but only through carelessness; and we won't worry about that here.)

In this treatment we'll only consider polynomial interpolation schemes such as rounding, linear interpolation, and cubic interpolation. These schemes amount to evaluating polynomials (of degree zero, one, and three, respectively) in the interstices between points of the wavetable. The idea is that, for any index $x$, we choose a nearby ``good" point ${x_0}$, and let the output be calculated by some polynomial:

$${a_0} + {a_1} (x - {x_0}) + {a_2} {{({x - x_0})}^ 2 } + \cdots + {a_n} {{({x - x_0})}^ n } .$$

Usually we choose the polynomial which passes through the $n+1$ nearest points of the wavetable. For 1-point interpolation (a zero-degree polynomial) this means letting $a_0$ equal the nearest point of the wavetable. For two-point interpolation, we draw a line segment between the two points of the wavetable on either side of the desired point $x$. We can let $x_0$ be the point at the left of $x$ (which we write as $\lfloor x \rfloor$) and then the formula for linear interpolation is:

$$y[{x_0}] + (y[{x_0} + 1]- y[{x_0}]) \cdot (x - {x_0}).$$

or in other words,

$${a_0} = y[{x_0}] ,$$

$${a_1} = (y[{x_0} + 1]- y[{x_0}]).$$

In general, you can fit exactly one polynomial of degree $n-1$ through any $n$ points as long as their $x$ values are all different.

Figure 2.11 shows the effect of using linear (two-point) interpolation to fill in a sinusoid of period 6. At the top are three traces: the original sinusoid, the linearly-interpolated result of using 6 points per period to represent the sinusoid, and finally, another sinusoid, of slightly smaller amplitude, which better matches the six-segment waveform. The error introduced by replacing the original sinusoid by the linearly interpolated version has two components: first, a (barely perceptible) change in amplitude, and second, a (very perceptible) distortion of the wave shape.

Figure 2.11: Linear interpolation of a sinusoid: (a) the original sinusoid, the interpolated sinusoid, and the best sinusoidal fit back to the interpolated version; (b) the error.

The bottom graph in the figure shows the difference between the segmented waveform and the best-fitting sinusoid. This is a residual signal all of whose energy lies in overtones of the original sinusoid. As the number of points increases, the error gets smaller in magnitude. Since the error is the difference between a sinusoid and a sequence of approximating line segments, the magnitude of the error is roughly proportional to the square of the phase difference between each pair of points, or in other words, inversely proportional to the square of the number of points in the wavetable. Put another way, wavetable error decreases by 12 dB each time the table doubles in size. (This approximation is only good for tables with 4 or more points.)

Four-point (cubic) interpolation works similarly. The formula is:

$$-f (f-1)(f-2)/6 \cdot y[{x_0}-1] + (f+1)(f-1)(f-2)/2 \cdot y[{x_0}]$$

$$- (f+1) f (f-2) / 2 \cdot y[{x_0}+1] + (f+1) f (f-1) / 6 \cdot y[{x_0}+2] ,$$

where $f = x - {x_0}$ is the fractional part of the index. For tables with 4 or more points, doubling the number of points on the table tends to improve the RMS error by 24 dB. Table 2.5 shows the calculated RMS error for sinusoids at various periods for 1, 2, and 4 point interpolation. (A slightly different quantity is measured in [Moo90, p.164]. There, the errors in amplitude and phase are also added in, yielding slightly more pessimistic results. See also [Har87].)

Table 2.1: RMS error for table lookup using 1, 2, and 4 point interpolation at various table sizes.

period	interpolation points
	1	2	4
2	-1.2	-17.1	-20.2
3	-2.0	-11.9	-15.5
4	-4.2	-17.1	-24.8
8	-10.0	-29.6	-48.4
16	-15.9	-41.8	-72.5
32	-21.9	-53.8	-96.5
64	-27.9	-65.9	-120.6
128	-34.0	-77.9	-144.7

The allowable input domain for table lookup depends on the number of points of interpolation. In general, if using $k$-point interpolation into a table with $N$ points, the inputs may range over an interval of $N + 1 - k$ points. If $k=1$ (i.e., no interpolation at all), the domain is from $0$ to $N$ (including the endpoint at $0$ but excluding the one at $n$) if input values are truncated (and we will always use truncation in this book when doing non-interpolated table lookup). The domain is from -$1/2$ to $N-1/2$ if, instead, we round the input to the nearest integer instead of interpolating. In either case, the domain stretches over a length of $N$ points.

For two-point interpolation, the inputs must lie between the first and last points, that is, between $0$ and $N-1$. So the $N$ points suffice to define the function over a domain of length $N-1$. For four-point interpolation, we cannot get values for inputs between 0 and 1 (not having the required two points to the left of the input) and neither can we for the space between the last two points ($N-2$ and $N-1$. So in this case the domain reaches from $1$ to $N-2$ and has length $N-3$.

Periodic waveforms stored in wavetables require special treatment at the ends of the table. For example, suppose we wish to store a pure sinusoid of length $N$. For a noninterpolating table, it suffices to set, for example,

$$x[n] = cos(2 \pi n / N) , n = 0, \cdots, N-1 .$$

For two-point interpolation, we need $N+1$ points:

$$x[n] = cos(2 \pi n / N) , n = 0, \cdots, N ;$$

in other words, we must repeat the first ($n=0$) point at the end, so that the last segment from $N-1$ to $N$ reaches back to the beginning value.

For four-point interpolation, the cycle must be adjusted to start at the point $n=1$, since we can't get properly interpolated values out for inputs less than one. If, then, one cycle of the wavetable is arranged from $1$ to $N$, we must supply extra points for $0$ (copied from $N$), and also $N+1$ and $N+2$, copied from $1$ and $2$, to make a table of length $N+3$. For the same sinusoid as above, the table should contain:

$$x[n] = cos(2 \pi (n-1) / N) , n = 0, \cdots, N+2 .$$