Free Books
   Free Books    Physical Audio Signal Processing  

Windowed Sinc Interpolation

Bandlimited interpolation of discrete-time signals is a basic tool having extensive application in digital signal processing.5.8In general, the problem is to correctly compute signal values at arbitrary continuous times from a set of discrete-time samples of the signal amplitude. In other words, we must be able to interpolate the signal between samples. Since the original signal is always assumed to be bandlimited to half the sampling rate, (otherwise aliasing distortion would occur upon sampling), Shannon's sampling theorem tells us the signal can be exactly and uniquely reconstructed for all time from its samples by bandlimited interpolation.

Considerable research has been devoted to the problem of interpolating discrete points. A comprehensive survey of ``fractional delay filter design'' is provided in [267]. A comparison between classical (e.g., Lagrange) and bandlimited interpolation is given in [407]. The book Multirate Digital Signal Processing [97] provides a comprehensive summary and review of classical signal processing techniques for sampling-rate conversion. In these techniques, the signal is first interpolated by an integer factor $ L$ and then decimated by an integer factor $ M$. This provides sampling-rate conversion by any rational factor $ L/M$. The conversion requires a digital lowpass filter whose cutoff frequency depends on $ \max\{L,M\}$. While sufficiently general, this formulation is less convenient when it is desired to resample the signal at arbitrary times or change the sampling-rate conversion factor smoothly over time.

In this section, a public-domain resampling algorithm is described which will evaluate a signal at any time specifiable by a fixed-point number. In addition, one lowpass filter is used regardless of the sampling-rate conversion factor. The algorithm effectively implements the ``analog interpretation'' of rate conversion, as discussed in [97], in which a certain lowpass-filter impulse response must be available as a continuous function. Continuity of the impulse response is simulated by linearly interpolating between samples of the impulse response stored in a table. Due to the relatively low cost of memory, the method is quite practical for hardware implementation.

In the next section, the basic theory is presented, followed by sections on implementation details and practical considerations.


Theory of Ideal Bandlimited Interpolation

We review briefly the ``analog interpretation'' of sampling rate conversion [97] on which the present method is based. Suppose we have samples $ x(nT)$ of a continuous absolutely integrable signal $ x(t)$, where $ t$ is time in seconds (real), $ n$ ranges over the integers, and $ T$ is the sampling period. We assume $ x(t)$ is bandlimited to $ \pm f_s/2$, where $ f_s=1/T$ is the sampling rate. If $ X(\omega)$ denotes the Fourier transform of $ x(t)$, i.e., $ X(\omega)=\int_{-\infty}^{\infty} x(t) e^{-j\omega t} dt$, then we assume $ X(\omega)=0$ for $ \vert\omega\vert\geq\pi f_s$. Consequently, Shannon's sampling theorem gives us that $ x(t)$ can be uniquely reconstructed from the samples $ x(nT)$ via

$\displaystyle {\hat x}(t) \isdef \sum_{n=-\infty}^{\infty} x(nT) h_s(t-nT) \equiv x(t), \protect$ (5.13)

where

$\displaystyle h_s(t) \isdef$   sinc$\displaystyle (f_st) \isdef \frac{\sin(\pi f_st)}{\pi f_st}. $

To resample $ x(t)$ at a new sampling rate $ f_s^\prime=1/T^\prime$, we need only evaluate Eq.$ \,$(4.13) at integer multiples of $ T^\prime$.

When the new sampling rate $ f_s^\prime$ is less than the original rate $ f_s$, the lowpass cutoff must be placed below half the new lower sampling rate. Thus, in the case of an ideal lowpass, $ h_s(t) = \min\{1,f_s^\prime/f_s\}$   sinc$ (\min\{f_s,f_s^\prime\}t)$, where the scale factor maintains unity gain in the passband.

A plot of the sinc function sinc$ (t) \isdeftext \sin(\pi t)/(\pi t)$ to the left and right of the origin $ t=0$ is shown in Fig.4.21. Note that peak is at amplitude $ 1$, and zero-crossings occur at all nonzero integers. The sinc function can be seen as a hyperbolically weighted sine function with its zero at the origin canceled out. The name sinc function derives from its classical name as the sine cardinal (or cardinal sine) function.

Figure 4.21: The sinc function plotted for seven zero-crossings to the left and right.
\includegraphics[width=3in]{eps/Sinc}

If ``$ \ast $'' denotes the convolution operation for digital signals, then the summation in Eq.$ \,$(4.13) can be written as $ (x\ast h_s)(t)$.

Equation Eq.$ \,$(4.13) can be interpreted as a superpositon of shifted and scaled sinc functions $ h_s$. A sinc function instance is translated to each signal sample and scaled by that sample, and the instances are all added together. Note that zero-crossings of sinc$ (z)$ occur at all integers except $ z=0$. That means at time $ t=nT$, (i.e., on a sample instant), the only contribution to the sum is the single sample $ x(nT)$. All other samples contribute sinc functions which have a zero-crossing at time $ t=nT$. Thus, the interpolation goes precisely through the existing samples, as it should.

A plot indicating how sinc functions sum together to reconstruct bandlimited signals is shown in Fig.4.22. The figure shows a superposition of five sinc functions, each at unit amplitude, and displaced by one-sample intervals. These sinc functions would be used to reconstruct the bandlimited interpolation of the discrete-time signal $ x = [\ldots ,0,1,1,1,1,1,0,\ldots ]$. Note that at each sampling instant $ t=nT$, the solid line passes exactly through the tip of the sinc function for that sample; this is just a restatement of the fact that the interpolation passes through the existing samples. Since the nonzero samples of the digital signal are all $ 1$, we might expect the interpolated signal to be very close to $ 1$ over the nonzero interval; however, this is far from being the case. The deviation from unity between samples can be thought of as ``overshoot'' or ``ringing'' of the lowpass filter which cuts off at half the sampling rate, or it can be considered a ``Gibbs phenomenon'' associated with bandlimiting.

Figure 4.22: Bandlimited reconstruction of the signal $ x(t)$ from its samples $ x = [\ldots ,0,1,1,1,1,1,0,\ldots ]$. The dots show the signal samples, the dashed lines show the component sinc functions, and the solid line shows the unique bandlimited reconstruction from the samples obtained by summing the component sinc functions.
\includegraphics{eps/SincSum}

A second interpretation of Eq.$ \,$(4.13) is as follows: to obtain the interpolation at time $ t$, shift the signal samples under one sinc function so that time $ t$ in the signal is translated under the peak of the sinc function, then create the output as a linear combination of signal samples where the coefficient of each signal sample is given by the value of the sinc function at the location of each sample. That this interpretation is equivalent to the first can be seen as a result of the fact that convolution is commutative; in the first interpretation, all signal samples are used to form a linear combination of shifted sinc functions, while in the second interpretation, samples from one sinc function are used to form a linear combination of samples of the shifted input signal. The practical bandlimited interpolation algorithm presented below is based on the second interpretation.

From Theory to Practice

The summation in Eq.$ \,$(4.13) cannot be implemented in practice because the ``ideal lowpass filter'' impulse response $ h_s(t)$ actually extends from minus infinity to infinity. It is necessary in practice to window the ideal impulse response so as to make it finite. This is the basis of the window method for digital filter design [115,362]. While many other filter design techniques exist, the window method is simple and robust, especially for very long impulse responses. In the case of the algorithm presented below, the filter impulse response is very long because it is heavily oversampled. Another approach is to design optimal decimated ``sub-phases'' of the filter impulse response, which are then interpolated to provide the ``continuous'' impulse response needed for the algorithm [358].

Figure 4.23 shows the frequency response of the ideal lowpass filter. This is just the Fourier transform of $ h_s(t)$.

Figure 4.23: Frequency response of the ideal lowpass filter.
\includegraphics{eps/IdealLowpass}

If we truncate $ h_s(t)$ at the fifth zero-crossing to the left and the right of the origin, we obtain the frequency response shown in Fig.4.24. Note that the stopband exhibits only slightly more than 20 dB rejection.

Figure 4.24: Frequency response of the ideal lowpass filter after rectangularly windowing the ideal (sinc) impulse response at the fifth zero crossing to the left and right of the time origin. The vertical axis is in units of decibels (dB), and the horizontal axis is labeled in units of spectral samples between plus and minus half the sampling rate.
\includegraphics{eps/SincTruncatedF}

If we instead use the Kaiser window [221,438] to taper $ h_s(t)$ to zero by the fifth zero-crossing to the left and the right of the origin, we obtain the frequency response shown in Fig.4.25. Note that now the stopband starts out close to $ -80$ dB. The Kaiser window has a single parameter which can be used to modify the stop-band attenuation, trading it against the transition width from pass-band to stop-band.

Figure 4.25: Frequency response of the ideal lowpass filter Kaiser windowed at the fifth zero crossing to the left and right.
\includegraphics{eps/SincKaiseredF}

Implementation

The implementation below provides signal evaluation at an arbitrary time, where time is specified as an unsigned binary fixed-point number in units of the input sampling period (assumed constant).

Figure 4.26: Time register format.
\includegraphics[scale=0.8]{eps/TimeRegisterFormat}

Figure 4.27: Waveforms and parameters in the interpolator.
\includegraphics[scale=0.8]{eps/Waveforms}

Figure 4.26 shows the time register $ t$, and Figure 4.27 shows an example configuration of the input signal and lowpass filter at a given time. The time register is divided into three fields: The leftmost field gives the number $ n$ of samples into the input signal buffer, the middle field is an initial index $ l$ into the filter coefficient table $ h(l)$, and the rightmost field is interpreted as a number $ \epsilon $ between 0 and $ 1$ for doing linear interpolation between samples $ l$ and $ l+1$ (initially) of the filter table. The concatenation of $ l$ and $ \epsilon $ are called $ P\in[0,1)$ which is interpreted as the position of the current time between samples $ n$ and $ n+1$ of the input signal.

Let the three fields have $ {n_n}$, $ n_l$, and $ n_\eta$ bits, respectively. Then the input signal buffer contains $ N=2^{n_n}$ samples, and the filter table contains $ L=2^n_l$ ``samples per zero-crossing.'' (The term ``zero-crossing'' is precise only for the case of the ideal lowpass; to cover practical cases we generalize ``zero-crossing'' to mean a multiple of time $ t_c=0.5/f_c$, where $ f_c$ is the lowpass cutoff frequency.) For example, to use the ideal lowpass filter, the table would contain $ h(l)=$sinc$ (l/L)$.

Our implementation stores only the ``right wing'' of a symmetric finite-impulse-response (FIR) filter (designed by the window method based on a Kaiser window [362]). Specifically, if $ h(l)$, $ l\in[-LN_z,LN_z]$, denotes a length $ 2LN_z+1$ symmetric finite impulse response, then the right wing of $ h$ is defined as the set of samples $ h(l)$ for $ l\in[0,LN_z]$. By symmetry, the left wing can be reconstructed as $ h(-l)=h(l)$, $ l=1,2,\ldots,LN_z$.

Our implementation also stores a table of differences $ \hbar(l) = h(l+1) - h(l)$ between successive FIR sample values in order to speed up the linear interpolation. The length of each table is $ {\hat N}=LN_z+1$, including the endpoint definition $ \hbar({\hat N})=0$.

Consider a sampling-rate conversion by the factor $ \rho = f_s^\prime/f_s$. For each output sample, the basic interpolation Eq.$ \,$(4.13) is performed. The filter table is traversed twice--first to apply the left wing of the FIR filter, and second to apply the right wing. After each output sample is computed, the time register is incremented by $ 2^{n_l+n_\eta}/\rho$ (i.e., time is incremented by $ 1/\rho$ in fixed-point format). Suppose the time register $ t$ has just been updated, and an interpolated output $ y(t)$ is desired. For $ \rho\geq1$, output is computed via

\begin{eqnarray*} v & \gets & \sum_{i=0}^{\mbox{$h$\ end}} x(n-i) \left[h(l+iL) ... ...$\ end}} x(n+1+i) \left[h(l+iL) + \epsilon \hbar(l+iL)\right], \end{eqnarray*}

where $ x(n)$ is the current input sample, and $ \epsilon \in[0,1)$ is the interpolation factor. When $ \rho<1$, the initial $ P$ is replaced by $ P^\prime=\rho P$, $ 1-P$ becomes $ \rho-P^\prime = \rho(1-P)$, and the step-size through the filter table is reduced to $ \rho L$ instead of $ L$; this lowers the filter cutoff to avoid aliasing. Note that $ \epsilon $ is fixed throughout the computation of an output sample when $ \rho\geq1$ but changes when $ \rho<1$.

When $ \rho<1$, more input samples are required to reach the end of the filter table, thus preserving the filtering quality. The number of multiply-adds per second is approximately $ (2N_z+1)\max\{f_s,f_s^\prime\}$. Thus the higher sampling rate determines the work rate. Note that for $ \rho<1$ there must be $ \lceil{N_zf_s/f_s^\prime}\rceil $ extra input samples available before the initial conversion time and after the final conversion time in the input buffer. As $ \rho\to 0$, the required extra input data becomes infinite, and some limit must be chosen, thus setting a minimum supported $ \rho$. For $ \rho\geq1$, only $ N_z$ extra input samples are required on the left and right of the data to be resampled, and the upper bound for $ \rho$ is determined only by the fixed-point number format, viz., $ \rho_{\mbox{max}}= 2^{n_l+n_\eta}$.

As shown below, if $ n_c$ denotes the word-length of the stored impulse-response samples, then one may choose $ n_l=n_c/2$, and $ n_\eta=n_c/2$ to obtain $ n_c-1$ effective bits of precision in the interpolated impulse response.

Note that rational conversion factors of the form $ \rho=L/m$, where $ L=2^n_l$ and $ m$ is an arbitrary positive integer, do not use the linear interpolation feature (because $ \epsilon \equiv0$). In this case our method reduces to the normal type of bandlimited interpolator [97]. With the availability of interpolated lookup, however, the range of conversion factors is boosted to the order of $ 2^{n_l+n_\eta}/m$. E.g., for $ \rho\approx1$, $ n_l=9,n_\eta=8$, this is about $ 5.1$ decimal digits of accuracy in the conversion factor $ \rho$. Without interpolation, the number of significant figures in $ \rho$ is only about $ 2.7$.

The number $ N_z$ of zero-crossings stored in the table is an independent design parameter. For a given quality specification in terms of aliasing rejection, a trade-off exists between $ N_z$ and sacrificed bandwidth. The lost bandwidth is due to the so-called ``transition band'' of the lowpass filter [362]. In general, for a given stop-band specification (such as ``80 dB attenuation''), lowpass filters need approximately twice as many multiply-adds per sample for each halving of the transition band width.

As a practical design example, we use $ N_z=13$ in a system designed for high audio quality at $ 20$% oversampling. Thus, the effective FIR filter is $ 27$ zero crossings long. The sampling rate in this case would be $ 50$ kHz.5.9In the most straightforward filter design, the lowpass filter pass-band would stop and the transition-band would begin at $ 20$ kHz, and the stop-band would begin (and end) at $ 25$ kHz. As a further refinement, which reduces the filter design requirements, the transition band is really designed to extend from $ 20$ kHz to $ 30$ kHz, so that the half of it between $ 25$ and $ 30$ kHz aliases on top of the half between $ 20$ and $ 25$ kHz, thereby approximately halving the filter length required. Since the entire transition band lies above the range of human hearing, aliasing within it is not audible.

Using $ 512$ samples per zero-crossing in the filter table for the above example (which is what we use at CCRMA, and which is somewhat over designed) implies desiging a length $ 27\times 512 = 13824$ FIR filter having a cut-off frequency near $ \pi/512$. It turns out that optimal Chebyshev design procedures such as the Remez multiple exchange algorithm used in the Parks-McLellan software [362] can only handle filter lengths up to a couple hundred or so. It is therefore necessary to use an FIR filter design method which works well at such very high orders, and the window method employed here is one such method.

It is worth noting that a given percentage increase in the original sampling rate (``oversampling'') gives a larger percentage savings in filter computation time, for a given quality specification, because the added bandwidth is a larger percentage of the filter transition bandwidth than it is of the original sampling rate. For example, given a cut-off frequency of $ 20$ kHz, (ideal for audio work), the transition band available with a sampling rate of $ 44$ kHz is about $ 2$ kHz, while a $ 48$ kHz sampling rate provides a $ 4$ kHz transition band. Thus, a $ 10$% increase in sampling rate halves the work per sample in the digital lowpass filter.

Choice of Table Size and Word Lengths

It is desirable that the stored filter impulse response be sampled sufficiently densely so that interpolating linearly between samples does not introduce error greater than the quantization error. It is shown in [462] that this condition is satisfied when the filter impulse-response table contains at least $ L=2^{1+n_c/2}$ entries per ``zero-crossing'', where $ n_c$ is the number of bits allocated to each table entry. (A later, sharper, error bound gives that $ L=2^{n_c/2}$ is sufficient.) It is additionally shown in [462] that the number of bits in the interpolation between impulse-response samples should be near $ n_c/2$ or more. With these choices, the linear interpolation error and the error due to quantized interpolation factors are each about equal to the coefficient quantization error. A signal resampler designed according to these rules will typically be limited primarily by the lowpass filter design, rather than by quantization effects.

Summary of Windowed Sinc Interpolation

The digital resampling method described in this section is convenient for bandlimited interpolation of discrete-time signals at arbitrary times and/or for arbitrary changes in sampling rate. The method is well suited for software or hardware implementation, and widely used free software is available.

Next Section:
Delay-Line Interpolation Summary
Previous Section:
Thiran Allpass Interpolators