# Time-domain interpolation using the Fast Fourier Transform

If you read an earlier article about FIR filter impulse responseinterpolation using**frequency domain zero stuffing**, it is logical toconclude that if we can interpolate time-domain impulse responses, we should beable to interpolate time-domain sigÂnals using the same frequency-domain zerostuffing method .

This is true, and the aboveinterpolation-by-*M * process applied to time signals is someÂtimes called *exactinterpolation * â€”because its performance is equivalent to using an ideal,infinite-stopband attenuation, time-domain interpolation filterâ€”and has madeits way into DSP textbooks, journal articles, and classroom notes of famous DSPprofessors.

To establish our notation, letâ€™s say we compute the FFTof an *N* -point *x* (*n* ) time sequence to produce its *X* (*m* )frequency-domain samples. Next, we stuff (*M* â€“1)*N * zeros in themiddle of *X* (*m* ) to yield the *MN* -length *X* _{int } (*m* )frequency samples, where *MN * is an integer power of two. Then we performan *MN* -point inverse FFT on *X* _{int } (*m* ) to obtain theinterpolated-by-*Mx* _{int } (*n* ) times samples. Using thisfrequency-domain zero stuffing to implement time-domain signal interpolationinvolves two important issues upon which we now focus: 1) interpolating realsignals and 2) interpolating analytic signals.

**ComputingInterpolated Real Signals **

The first issue: to ensure the interpolated *x* _{int } (*n* )time sequence is real only, conjugate symmetry must be maintained in thezero-stuffed *X* _{int } (*m* ) frequency samples. If the *X* (*m* )sequence has a nonzero sample at *X* _{int } (*N* /2), the *f _{s } * /2freÂquency 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. (*Thisstep, to maintain conjugate symmetry and improve interpolation accuracy, is notso well known.)*

. â€¢Assign *X* (*q* ), where *MN* â€“(*N* /2)+1 â‰¤ *m* â‰¤ *MN* â€“1, and

* ^{X } * int

^{(m )= }(

*N*/2)+1 â‰¤

*q*â‰¤

*N*â€“1.

â€¢Computethe

*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 * amÂplitude loss induced byinterpolation.

Whew! Ourmathematical 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 interpoÂlate 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 thenice properties of the above algorithm is that every *M* th *x* _{int } (*n* )sample coincides with the original *x* (*n* ) samples. In practice, dueto 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 inÂverse FFT of *X* _{int } (*m* ).

Hereâ€™s the second issue regarding time-domain real signalinterpolation, and itâ€™s a big deal indeed. This *exact * interpolationalgorithm provides correct results only if the original *x* (*n* )sequence is periodic within its full time interÂval.

If *X* (*m* ) exhibits any spectral leakage, like the signalswith which we usuÂally work, the interpolated *x* _{int } (*n* )sequence can contain noticeable amplitude errors in its beginning and endingsamples, 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 correctinterpolation; (b) interpolation error.

In that figure, squares (both white and black) represent the 48-pointinterpolated *x* _{int } (*n* ) seÂquence, 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 becausetheyâ€™re hidden behind the squares.) The interpolation error, the differencebetween the correct samples and *x* _{int } (*n* ), is given in Figure13â€“74(b).

These interpolation errors exist because *X* _{int } (*m* )is not identical to the spectrum obtained had we sampled *x* (*n* ) ata rate of *Mf* and performed an *MN* -point forward FFT. Thereâ€™s noclosed-form mathematical expression enÂabling us to predict these errors. Theerror depends on the spectral compoÂnent magnitudes and phases within *x* (*n* ),as well as *N * and *M* . Reference [72] reports a reduction ininterpolation errors for two-dimensional image interÂpolation when, in aneffort to reduce spectral leakage in *X* (*m* ), simple winÂdowing isapplied to the first few and last few samples of *x* (*n* ).

With theadvent of fast hardware DSP chips and pipelined FFT techÂniques, the abovetime-domain interpolation algorithm may be viable for a number of applications,such as computing selectable sample rate time seÂquences of a test signal thathas a fixed spectral envelope shape; providing inÂterpolation, by selectablefactors, of signals that were filtered in the frequency domain using the fastconvolution method; or digital image re-sampling.

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

**ComputingInterpolated Analytic Signals **

We can use the frequency-domain zero stuffing scheme to generate aninterpolated-by-*M * analytic time signal based upon the real *N* -pointtime seÂquence *x* (*n* ), if *N * is even. The process is asfollows:

. â€¢Â Â Â Â Â Â Â 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 * amÂplitude loss induced byinterpolation.

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

*Used with thepermission of the publisher, Prentice Hall, this on-going series of articles onEmbedded.com is based on copyrighted material from “UnderstandingDigital Signal Processing, Second Edition” by Richard G. Lyons. The book canbe purchased on line. *

*RichardLyons * * is a consultingsystems engineer and lecturer with Besser Associates. As a lecturer with Besserand an instructor for the University of California Santa Cruz Extension, Lyonshas delivered digitasl signal processing seminars and training course attechnical conferences as well at companies such as Motorola, Freescale, LockheedMartin, Texas Instruments, Conexant, Northrop Grumman, Lucent, Nokia, Qualcomm,Honeywell, National Semiconductor, General Dynamics andInfinion. *