Blogs Rick Lyons FFT Interpolation Based on FFT Samples: A Detective Story With a Surprise

   Blogs    Rick Lyons  

FFT Interpolation Based on FFT Samples: A Detective Story With a Surprise Ending

https://www.dsprelated.com/showarticle/1156.php

This blog presents several interesting things I recently learned regarding the estimation of a spectral value located at a frequency lying between previously computed FFT spectral samples. My curiosity about this FFT interpolation process was triggered by reading a spectrum analysis paper written by three astronomers [1].

My fixation on one equation in that paper led to the creation of this blog.

Background

 

This article is available in PDF format for easy printing

The notion of FFT interpolation is straightforward to describe. That is, for example, given an N = 16 sample x(n) time-domain sequence shown in Figure 1(a), performing an N = 16 point FFT on x(n) produces the |X(m)| magnitude of samples shown by the red dots in Figure 1(b). The blue dashed curve in Figure 1(b) is the magnitude of the discrete-time Fourier transform (DTFT) of x(n), what I like to call the "true spectrum" of x(n).

 

Figure 1: (a) An x(n) time sequence; (b) the spectral magnitude of x(n).

The process of FFT interpolation is estimating an X(k) spectral value at, for example, a non‑integer frequency of k = 4.4688 lying between the m = 4 and m = 5 FFT samples. The magnitude of X(k) is marked by the bold 'x' in Figure 1(b).

An X(k)Estimation Equation That Does Not Work

Reference [1] presents the following equation for computing a single X(k) spectral value based on N FFT samples:

where:

• N = the number of x(n) input samples and the number of FFT samples
• X(m) = the original FFT samples, the FFT of x(n)
• m = the frequncy-domain integer index of X(m) where 0 ≤ m ≤ N-1
• k = the non-integer "FFT bin-like" frequency value of the desired spectral value
• j = sqrt(-1)

To avoid overwhelming the reader with too many algebra equations, I've included the derivation of Eq. (1) in Appendix A of the downloadable PDF file associated with this blog.

Not having seen Eq. (1) before I decided to model it using MATLAB software. Plugging my Figure 1(b) complex-valued X(m) FFT samples and k = 4.4688 into Eq. (1) I computed an XRef.[1](4.4688) value of:

  XRef.[1](4.4688) = -1.1023 -j1.1053.

Next, padding the x(n) time sequence with lots of zero-valued samples I computed a finer-granularity spectrum shown by the blue dashed curve in Figure 1(b). From that finer-resolution spectrum I learned the true value of XTrue(4.4688) is:

  XTrue(4.4688) = 0.3537 -j1.0489

which is not even close to the above XRef.[1](4.4688) spectral value computed using Eq. (1)! So what's going on here(!)?

Knowing that my home-grown MATLAB code is about as reliable as the paper towel dispenser in the Men's Room, I asked our DSP colleague Neil Robertson to implement Eq. (1) using my above x(n) sequence and a value of k = 4.4688. Neil graciously agreed and his MATLAB code produced the same incorrect XRef.[1](4.4688) = ‑1.1023 ‑j1.1053 that I obtained.

So now it seems likely that Eq. (1) is not valid. But if that is true, how could the three authors of Reference [1] have missed that error? That was the mystery I set out to solve.

We'll revisit Eq. (1) later in this blog and learn something interesting in our quest to understand that equation.

An X(k) Estimation Equation That Does Work

Knowing that Eq. (1) is questionable I started looking through my DSP books to find a valid FFT interpolation formula. In an old out-of-print book written by two DSP pioneers [2], I found the following expression for computing the spectral value at non-integer frequency k based on N X(m) FFT samples:

Equation (2) produces correct FFT interpolation results and its derivation is given in Appendix B of the downloadable PDF file. Using [2]'s z-transform derivation approach I performed a modified derivation to obtain the following computationally simpler FFT interpolation equation that also works just fine:

 

My derivation of Eq. (3) is presented in Appendix C of the downloadable PDF file.

Quiz of the Day

Equation (3) is a correct way to compute an interpolated X(k) value, but look at it carefully. Before continuing, do you see a way to reduce the number of arithmetic operations needed to compute a single X(k) value? Did you notice that neither the numerator nor the 1/N factor inside the summation in Eq. (3) are functions of m? That fact allows us to move them both outside the summation giving us:

 

Equation (4) has a significantly reduced computational workload relative to that of Eq. (3) as we shall see. Don't feel bad if you didn't foresee those beneficial simplifications leading to Eq. (4). Neither did I—while writing this blog I found those slick algebraic moves in [3].

So now, Eq. (4) is my preferred method for computing a single interpolated X(k) spectral value based on N FFT samples. Now that we have a workable Eq. (4), let's briefly revisit the suspicious Eq. (1).

Here's What's Wrong With Eq. (1)

If you download the PDF file associated with this blog, in its Appendix A you'll see the published derivation of the troublesome Eq. (1). One intermediate step in that derivation, my Eq. (A‑5) repeated here as Eq. (5), produces the following expression for our desired X(k) as:

 

Equation (5) works just fine for computing our desired X(k)! So now it's interesting to see what went wrong in Eq. (1). The authors of [1] applied two algebraic approximations to Eq. (5) to produce the faulty Eq. (1):

 

Examples of the errors of those two approximations are shown in Figure 2.

 

Figure 2: Approximations' errors, versus m, when k = 2.4688 and N = 64: (a) Approximation# 1 error; (b) Approximation# 2 error.

Analyzing the errors induced by those two approximations solves the mystery of why Eq. (1) is flawed. By the way, I've modeled Eq. (1) for values of N = 16 out to N = 16384 with no improvement in its accuracy for larger N. Next we discuss the computational costs of this blog's FFT interpolation equations.

Computational Workload Comparisons

It's enlightening to compare the necessary arithmetic operations needed to implement the correct FFT interpolation equations given in this blog.

Assuming we have previously computed and stored all the constant factors in our equations that are not functions of k or m, Table 1 shows the reduction in the necessary computations needed by Eqs. (4) and (5) relative to Eqs. (2) and (3). (In Table 1, a "Trigonometric Computation" means computing the sine or cosine of some given angle.)

 

We remind the reader that Trigonometric Computations and Real Divides are far more computationally costly than Real Adds and Real Multiplies.

Here's the "Surprise Ending" To This Blog

While struggling to create Table 1 it finally hit me! The fastest way to compute a single interpolated FFT spectral sample at frequency k is to merely compute a brute force N‑point DFT of the original x(n) time samples using:

where, again, 0 ≤ k ≤ N‑1, and k is not an integer. Assuming the x(n) input in Eq. (6) is real-only, the computational workload comparison of Eq. (6) to my favorite Eq. (4) is given in Table 2.

 

From Table 2 we see that Eq. (6) requires significantly fewer arithmetic operations—notably eliminating 2N of those costly 'Real Divides'—than Eq. (4).

Conclusions

I presented an X(k) FFT interpolation equation published in the literature, Eq. (1), for computing a single X(k) spectral value based on N FFT samples, where frequency k is not an integer. Equation (1) is not valid and I showed why that is so. (I remain troubled by this and perhaps someone can prove me wrong.) Next I presented four different expressions, Eq. (2) through Eq. (5), for computing an X(m) sample and compared their computational costs. Finally I realized the most efficient method for computing a single X(m) sample is to simply perform an N‑point DFT using Eq. (6).

I apologize for this blog being so lengthy. I didn't have the time to make it shorter.

 

Postscript [April 26, 2018]: My original typographical error pointed out in DerekR's April 24th 2018 post has been corrected in the above blog and corrected in the updated "Version 2" PDF file associated with this blog.

References

[1] S. Ransom, S. Eikenberry, and J. Middleditch, "Fourier Techniques For Very Long Astrophysical Time-Series Analysis", The Astronomical Journal, Sept. 2002, Section 4.1, pp. 1796. [http://iopscience.iop.org/article/10.1086/342285/pdf]

[2] L. Rabiner and B. Gold, Theory and Application of Digital Signal Processing, Prentice Hall, Upper Saddle River, New Jersey, 1975, p. 54.

[3] J. Proakis and D. Manolakis, Digital Signal Processing-Principles, Algorithms, and Applications, 3rd Ed., Prentice Hall, Upper Saddle River, New Jersey, 1996, pp. 408.

 

---

Previous post by Rick Lyons:
   An Efficient Linear Interpolation Scheme
Next post by Rick Lyons:
   Two Easy Ways To Test Multistage CIC Decimation Filters

---

You might also like... (promoted content)

 

Tutorial: Running a RISC-V Processor on the Arty A7

 

Expert C Programming Techniques for Embedded Developers

 

How to reduce the bill of material costs with Digital Signal Processing

Comments

• Comments
• Write a Comment
• Select to add a comment

[ - ]

Comment by ScottRansom●May 11, 2018

 

Hi Rick,

So I'm the primary author of the paper that you are discussing and I can guarantee you that Fourier Interpolation works extremely well!  It has been used to help discover many hundreds of new pulsars.

One of the big issues here, I think, is that your Eqn 1 (my Eqn 29) was only derived to apply when N  is very large (N larger than hundreds or thousands).  But Eqn 29 isn't the interesting one.  Eqn 30 in my paper is:

 

If N is large (again, that is a requirement), then you can interpolate your FFT at a non-integer frequency "r" to whatever precision you like, by summing ~m local Fourier amplitudes (integer freqs "k" surrounding freq "r").  For single-precision work, I usually use m of between 30-40.

The result of using this type of Fourier interpolation is exactly equivalent to zero-padding the original time series (if you do them equivalently).  But if you only need small portions of the spectrum interpolated, Fourier interpolation via the above equation is hugely faster as it only needs a small number of operations.

Note that Eqn 30 is effectively a complex FIR filter, which means that we can also do the interpolation via correlations using the convolution theorem in the Fourier domain (using FFTs of small portions of the large-N FFT).  And, in fact, that is exactly how my pulsar search software uses Eqn 30 during large-scale Fourier searches.

Let me know if you have any questions,

Scott Ransom

http://www.cv.nrao.edu/~sransom