How to Interpolate the Peak Location of a DFT or FFT if the Frequency of Interest is Between Bins

How to interpolate the peak location of a DFT or FFT if the frequency of interest is between bins

by Matt Donadio


Problem

If the actual frequency of a signal does not fall on the center frequency of a DFT (FFT) bin, several bins near the actual frequency will appear to have a signal component. In that case, we can use the magnitudes of the nearby bins to determine the actual signal frequency.

This is a summary of some methods found at Some Frequency Estimation Algorithms. (Other methods of sinusoidal parameter estimation which do not rely on the DFT are found on that page; those methods also work well, and are sometimes faster than these.)

The following notation is used below:

  • k =  index of the max (possibly local) magnitude of an DFT.
  • X[i]  =  bin "i" of an DFT |X[i]| =  magnitude of DFT at bin "i".
  • k'  =  the interpolated bin location.

Solutions

The first two methods appeared in comp.dsp posts.

Quadradic Method ( see post) :

  1. y1 =  |X[k - 1]|
  2. y2  =  |X[k]|
  3. y3 =  |X[k + 1]|
  4. d = (y3 - y1) / (2 * (2 * y2 - y1 - y3))
  5. k'  =  k + d

Barycentric method (see post) :

  1. y1  = |X[k - 1]|
  2. y2 =  |X[k]|
  3. y3 =  |X[k + 1]|
  4. d = (y3 - y1) / (y1 + y2 + y3)
  5. k'  =  k + d

Quinn's First Estimator:

  1. ap = (X[k + 1].r * X[k].r + X[k + 1].i * X[k].i)  /  (X[k].r * X[k].r + X[k].i * X[k].i)
  2. dp = -ap  / (1.0 - ap)
  3. am = (X[k - 1].r * X[k].r + X[k - 1].i * X[k].i)  /  (X[k].r * X[k].r + X[k].i * X[k].i)
  4. dm = am / (1.0 - am)
  5. if (dp > 0) AND (dm > 0) then
       d = dp
    else
       d = dm
    end
  6. k' = k + d

Quinn's Second Estimator:

  1. tau(x) = 1/4 * log(3x^2 + 6x + 1) - sqrt(6)/24 * log((x + 1 - sqrt(2/3))  /  (x + 1 + sqrt(2/3)))
  2. ap = (X[k + 1].r * X[k].r + X[k+1].i * X[k].i)  /  (X[k].r * X[k].r + X[k].i * X[k].i)
  3. dp = -ap / (1 - ap)
  4. am = (X[k - 1].r * X[k].r + X[k - 1].i * X[k].i)  /  (X[k].r * X[k].r + X[k].i * X[k].i)
  5. dm = am / (1 - am)
  6. d = (dp + dm) / 2 + tau(dp * dp) - tau(dm * dm)
  7. k' = k + d

Jain's Method :

  1. y1 = |X[k-1]|
  2. y2 = |X[k]|
  3. y3 = |X[k+1]|
  4. if y1 > y3 then
       a = y2  /  y1
       d = a  /  (1 + a)
       k' = k - 1 + d
    else
       a = y3  /  y2
       d = a  /  (1 + a)
       k' = k + d
    end

Of the above methods, Quinn's second estimator has the least RMS error.

References

Posted 1999/04/19, Released 1999/05/03, Revised 1999/05/05

Phase Estimator?

I got good results with Quinn's Second Estimator, so now have the frequency. The phase of each bin is atan(X[k].i,X[k].r). Is there any estimator that calculates the phase corresponding to the estimated frequency?

Thanks.