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

and then decimated by an integer factor

. This provides sampling-rate conversion by any rational factor

. The conversion requires a digital lowpass
filter whose cutoff
frequency depends on

. 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.
We review briefly the ``analog interpretation'' of
sampling rate conversion
[
97] on which the present method is based. Suppose we have
samples

of a continuous absolutely integrable
signal 
,
where

is time in seconds (real),

ranges over the integers, and

is the
sampling period. We assume

is bandlimited to

, where

is the
sampling rate. If

denotes the
Fourier transform of

,
i.e.,

, then we assume

for

. Consequently,
Shannon's sampling
theorem gives us that

can be uniquely reconstructed from the
samples

via
 |
(5.13) |
where
To
resample 
at a new
sampling rate

, we need only
evaluate Eq.

(
4.13) at integer multiples of

.
When the new sampling rate

is less than the original rate

,
the lowpass cutoff must be placed below half the new lower sampling rate.
Thus, in the case of an ideal lowpass,

sinc

, where the scale factor maintains unity gain
in the
passband.
A plot of the
sinc function
sinc

to the left and right of the origin

is shown in Fig.
4.21.
Note that peak is at amplitude

, 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}](http://www.dsprelated.com/josimages_new/pasp/img1149.png) |
If ``

'' denotes the
convolution operation for digital signals, then
the summation in Eq.

(
4.13) can be written as

.
Equation Eq.

(
4.13) can be interpreted as a superpositon of
shifted and scaled sinc functions

. 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

occur at all integers except

. That means at time

, (
i.e., on a sample instant), the only contribution to the
sum is the single sample

. All other samples contribute sinc
functions which have a zero-crossing at time

. 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 ]$](http://www.dsprelated.com/josimages_new/pasp/img58.png)
. Note that at each sampling
instant

, 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

, we might expect the
interpolated signal to be very close to

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
from its samples
. 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.
 |
A second interpretation of Eq.

(
4.13) is as follows: to obtain the
interpolation at time

, shift the signal samples under
one sinc
function so that time

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.
The summation in Eq.

(
4.13) cannot be implemented in practice because
the ``ideal lowpass
filter''
impulse response 
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

.
Figure 4.23:
Frequency response of the ideal lowpass filter.
 |
If we truncate

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.
 |
If we instead use the
Kaiser window [
221,
438] to
taper

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
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.
 |
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}](http://www.dsprelated.com/josimages_new/pasp/img1160.png) |
Figure 4.27:
Waveforms and parameters in the interpolator.
![\includegraphics[scale=0.8]{eps/Waveforms}](http://www.dsprelated.com/josimages_new/pasp/img1161.png) |
Figure
4.26 shows the time register

, 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

of
samples into the input signal buffer, the middle field is an initial
index

into the filter coefficient table

, and the rightmost
field is interpreted as a number

between 0 and

for doing
linear interpolation between samples

and

(initially) of the
filter table. The concatenation of

and

are called

which is interpreted as the position of the current time
between samples

and

of the input signal.
Let the three fields have

,

, and

bits,
respectively. Then the input signal buffer contains

samples, and the filter table contains

``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

, where

is the lowpass cutoff frequency.) For example, to use the ideal
lowpass
filter, the table would contain
sinc
.
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

,
![$ l\in[-LN_z,LN_z]$](http://www.dsprelated.com/josimages_new/pasp/img1174.png)
, denotes a length

symmetric
finite
impulse response, then the
right wing
of

is defined
as the set of samples

for
![$ l\in[0,LN_z]$](http://www.dsprelated.com/josimages_new/pasp/img1176.png)
. By symmetry, the
left wing can be reconstructed as

,

.
Our implementation also stores a table of
differences

between successive FIR sample values in order to
speed up the linear interpolation. The length of each table is

, including the endpoint definition

.
Consider a
sampling-rate conversion by the factor

.
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

(
i.e., time is incremented by

in
fixed-point format). Suppose the time register

has just been
updated, and an interpolated output

is desired. For

, output is computed via
where

is the current input sample, and

is the
interpolation factor. When

, the initial

is replaced by

,

becomes

, and the
step-size through the filter table is reduced to

instead of

; this lowers the filter cutoff to avoid
aliasing. Note that

is fixed throughout the computation of an output sample when

but changes when

.
When

, 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

.
Thus the higher
sampling rate determines the work rate. Note that for

there must be

extra input samples
available before the initial conversion time and after the final conversion
time in the input buffer. As

, the required extra input
data becomes infinite, and some limit must be chosen, thus setting a
minimum supported

. For

, only

extra input samples are required on
the left and right of the data to be
resampled, and the upper bound for

is determined only by the fixed-point number format,
viz.,

.
As shown below, if

denotes the word-length of the stored
impulse-response samples, then one may choose

, and

to obtain

effective bits of precision in the
interpolated impulse response.
Note that rational conversion factors of the form

, where

and

is an arbitrary positive integer, do not use the linear
interpolation feature (because

). 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

.
E.g., for

,

, this is about

decimal digits of
accuracy in the conversion factor

. Without interpolation, the
number of significant figures in

is only about

.
The number

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

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

in a system designed for
high audio quality at

%
oversampling. Thus, the effective FIR
filter is

zero crossings long. The
sampling rate in this case would
be

kHz.
5.9In the most straightforward
filter design, the lowpass filter pass-band
would stop and the transition-band would begin at

kHz, and the
stop-band would begin (and end) at

kHz. As a further refinement,
which reduces the filter design requirements, the transition band is really
designed to extend from

kHz to

kHz, so that the half of it
between

and

kHz aliases on top of the half between

and

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

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

FIR filter
having a cut-off frequency near

. 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

kHz, (ideal for audio work), the transition band
available with a sampling rate of

kHz is about

kHz, while a

kHz sampling rate provides a

kHz transition band. Thus, a

% increase in sampling rate
halves the work per sample in
the digital lowpass filter.
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

entries per ``zero-crossing'', where

is the
number of bits allocated to each table entry. (A later, sharper,
error bound gives that

is sufficient.) It is
additionally shown in [
462] that the number of bits in the interpolation
between impulse-response samples should be near

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.
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 SummaryPrevious Section: Thiran Allpass Interpolators