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.

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

• Some Frequency Estimation Algorithms
• Eric Jacobsen’s Frequency Estimation Page
• Post: Inverse interpolation?
• Post: How do I interpolate between bins of an FFT?
• Post: Quinn on Quinn’s interpolation method

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