# Time-domain interpolation using the Fast Fourier Transform

If you read an earlier article about FIR filter impulse response
interpolation using** ****frequency domain zero stuffing**, it is logical to
conclude that if we can interpolate time-domain impulse responses, we should be
able to interpolate time-domain signals using the same frequency-domain zero
stuffing method.

This is true, and the above
interpolation-by-*M *process applied to time signals is sometimes called *exact
interpolation*—because its performance is equivalent to using an ideal,
infinite-stopband attenuation, time-domain interpolation filter—and has made
its way into DSP textbooks, journal articles, and classroom notes of famous DSP
professors.

To establish our notation, let’s say we compute the FFT
of an *N*-point *x*(*n*) time sequence to produce its *X*(*m*)
frequency-domain samples. Next, we stuff (*M*–1)*N *zeros in the
middle of *X*(*m*) to yield the *MN*-length *X*_{int}(*m*)
frequency samples, where *MN *is an integer power of two. Then we perform
an *MN*-point inverse FFT on *X*_{int}(*m*) to obtain the
interpolated-by-*Mx*_{int}(*n*) times samples. Using this
frequency-domain zero stuffing to implement time-domain signal interpolation
involves two important issues upon which we now focus: 1) interpolating real
signals and 2) interpolating analytic signals.

**Computing
Interpolated Real Signals **

The first issue: to ensure the interpolated *x*_{int}(*n*)
time sequence is real only, conjugate symmetry must be maintained in the
zero-stuffed *X*_{int}(*m*) frequency samples. If the *X*(*m*)
sequence has a nonzero sample at *X*_{int}(*N*/2), the *f _{s}*/2
frequency component, we must use the following steps in computing

*X*

_{int}(

*m*) to guarantee conjugate symmetry:

.
• Perform an *N*-point FFT on an *N*-point *x*(*n*)
time sequence, yielding *N *frequency samples, *X*(*m*).

.
• Create an *MN*-point spectral sequence *X*_{int}(*m*)
initially set to all zeros.

.
• Assign *X*_{int}(*m*)= *X*(*m*),
for 0 ≤ *m *≤ (*N*/2)–1.

.
• Assign both *X*_{int}(*N*/2) and *X*_{int}(*MN–N*/2)
equal to *X*(*N*/2)/2. (*This
step, to maintain conjugate symmetry and improve interpolation accuracy, is not
so well known.)*

.
• Assign *X*(*q*), where *MN*–(*N*/2)+1 ≤ *m
*≤ *MN*–1, and

* ^{X}*int

^{(m)= }(

*N*/2)+1 ≤

*q*≤

*N*–1.

• Compute the

*MN*-point inverse FFT of

*X*

_{int}(

*m*) yielding the desired

*MN*-length interpolated

*x*

_{int}(

*n*) sequence.

.
• Finally, if desired, multiply *x*_{int}(*n*)
by *M *to compensate for the 1/*M *amplitude loss induced by
interpolation.

Whew! Our
mathematical notation makes this signal interpolation scheme look complicated,
but it’s really not so bad. **Table 13–8 below**
shows the frequency-domain *X*_{int}(*m*) sample assignments,
where 0 ≤ *m *≤ 15, to interpolate an *N *= 8-point *x*(*n*)
sequence by a factor of *M *=2.

**Table 13–8 ***X_{INT}(m)
Assignments for Interpolation by Two *

One of the
nice properties of the above algorithm is that every *M*th *x*_{int}(*n*)
sample coincides with the original *x*(*n*) samples. In practice, due
to our finite-precision computing, the imaginary parts of our final *x*_{int}(*n*)
may have small non-zero values. As such, we take *x*_{int}(*n*)
to the be real part of the inverse FFT of *X*_{int}(*m*).

Here’s the second issue regarding time-domain real signal
interpolation, and it’s a big deal indeed. This *exact *interpolation
algorithm provides correct results only if the original *x*(*n*)
sequence is periodic within its full time interval.

If *X*(*m*) exhibits any spectral leakage, like the signals
with which we usually work, the interpolated *x*_{int}(*n*)
sequence can contain noticeable amplitude errors in its beginning and ending
samples, as shown in **Figure 13–74(a) below**
where an *N *= 24 sample *x*(*n*) sequence is interpolated by *M
*= 2.

**
**

**Figure 13–74 **Interpolation results, N = 24, M = 2: (a)
interpolated x_{int}(n), original x(n), and correct
interpolation; (b) interpolation error.

In that figure, squares (both white and black) represent the 48-point
interpolated *x*_{int}(*n*) sequence, white squares are the original *x*(*n*)
samples, and circles represent the exactly correct interpolated sample values.
(In the center portion of the figure, the circles are difficult to see because
they’re hidden behind the squares.) The interpolation error, the difference
between the correct samples and *x*_{int}(*n*), is given in Figure
13–74(b).

These interpolation errors exist because *X*_{int}(*m*)
is not identical to the spectrum obtained had we sampled *x*(*n*) at
a rate of *Mf*and performed an *MN*-point forward FFT. There’s no
closed-form mathematical expression enabling us to predict these errors. The
error depends on the spectral component magnitudes and phases within *x*(*n*),
as well as *N *and *M*. Reference [72] reports a reduction in
interpolation errors for two-dimensional image interpolation when, in an
effort to reduce spectral leakage in *X*(*m*), simple windowing is
applied to the first few and last few samples of *x*(*n*).

With the advent of fast hardware DSP chips and pipelined FFT techniques, the above time-domain interpolation algorithm may be viable for a number of applications, such as computing selectable sample rate time sequences of a test signal that has a fixed spectral envelope shape; providing interpolation, by selectable factors, of signals that were filtered in the frequency domain using the fast convolution method; or digital image re-sampling.

One scenario
to consider is using the efficient 2*N*-Point Real FFT technique to
compute the forward FFT of the real-valued *x*(*n*). Of course, the
prudent engineer would conduct a literature search to see what algorithms are
available for efficiently performing inverse FFTs when many of the
frequency-domain samples are zeros.

**Computing
Interpolated Analytic Signals **

We can use the frequency-domain zero stuffing scheme to generate an
interpolated-by-*M *analytic time signal based upon the real *N*-point
time sequence *x*(*n*), if *N *is even. The process is as
follows:

.
• Perform an *N*-point FFT on an *N*-point real *x*_{r}(*n*)
time sequence, yielding *N *frequency samples, *X*_{r}(*m*).

.
• Create an *MN*-point spectral sequence *X*_{int}(*m*)
initially set to all zeros, where *MN *is an integer power of two.

.
• Assign *X*_{int}(0) = *X*_{r}(0),
and *X*_{int}(*N*/2)
= *X*_{r}(*N*/2).

.
• Assign *X*_{int}(*m*) = 2*X*_{r}(*m*),
for 1 ≤ *m *≤ =(*N*/2)–1.

.
• Compute the *MN*-point inverse FFT of *X*_{int}(*m*)
yielding the desired *MN*-length interpolated analytic (complex) *x*_{c,int}(*n*)
sequence.

.
• Finally, if desired, multiply *x*_{c,int}(*n*)
by *M *to compensate for the 1/*M *amplitude loss induced by
interpolation.

The complex *x*_{c,int}(*n*)
sequence will also exhibit amplitude errors in its beginning and ending
samples.

*Used with the
permission of the publisher, Prentice Hall, this on-going series of articles on
Embedded.com is based on copyrighted material from "Understanding
Digital Signal Processing, Second Edition" by Richard G. Lyons. The book can
be purchased on line. *

*Richard
Lyons** is a consulting
systems engineer and lecturer with Besser Associates. As a lecturer with Besser
and an instructor for the University of California Santa Cruz Extension, Lyons
has delivered digitasl signal processing seminars and training course at
technical conferences as well at companies such as Motorola, Freescale, Lockheed
Martin, Texas Instruments, Conexant, Northrop Grumman, Lucent, Nokia, Qualcomm,
Honeywell, National Semiconductor, General Dynamics and
Infinion.*

## Loading comments... Write a comment