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.

Quadratic 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