How to Interpolate Frequency Peaks

dspGuru

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.

Quadradic Method (see post) :

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

Barycentric method (see post) :

y1 = |X[k - 1]|
y2 = |X[k]|
y3 = |X[k + 1]|
d = (y3 - y1) / (y1 + y2 + y3)
k' = k + d

Quinn's First Estimator:

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)
dp = -ap / (1.0 - ap)
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)
dm = am / (1.0 - am)
if (dp > 0) AND (dm > 0) then
d = dp
else
d = dm
end
k' = k + d

Quinn's Second Estimator:

tau(x) = 1/4 * log(3x^2 + 6x + 1) - sqrt(6)/24 * log((x + 1 - sqrt(2/3)) / (x + 1 + sqrt(2/3)))
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)
dp = -ap / (1 - ap)
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)
dm = am / (1 - am)
d = (dp + dm) / 2 + tau(dp * dp) - tau(dm * dm)
k' = k + d

Jain's Method :

y1 = |X[k-1]|
y2 = |X[k]|
y3 = |X[k+1]|
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