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

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.

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.