Extracting exact frequencies from FFT output

Say I pass 512 samples into my FFT

My microphone spits out data at 10KHz, so this represents 1/20s.

(So the lowest frequency FFT would pick up would be 40Hz).

The FFT will return an array of 512 frequency bins - bin 0: [0 - 40Hz) - bin 1: [40 - 80Hz) etc

So if my original sound contained energy at say 115Hz, how can I accurately retrieve this frequency?

That is going to lie in bin #2, but very close to bin #3. so I would expect both bins to contain something nonzero.

Question: how about the bins either side of this? Would they be guaranteed to be zero if there are no other frequencies close in the original signal?

Main question: is there some algorithm for deciphering the original frequency given the relative bin strengths?

-
You are aware of being able to decompose a complex number into its magnitude and phase, yes? – J. M. Nov 8 '10 at 11:37
Yes, but I can't see how this helps – P-i- Nov 8 '10 at 12:26
This paper, by Barry Kulp, gives a technique for using two adjacent bin values for interpolating the exact frequency: tedknowlton.com/resume/FFT_Bin_Interp.html. Hope it's of use. – mtrw Dec 22 '10 at 22:15

Remember that the FFT is circular. Inputs which contain an integer number of cycles will come out clean as a single point, in the corresponding bin. Those which do not, act as if they are multiplied by a rectangular pulse in the time domain, which creates convolution by a sinc function in the frequency domain. Since sinc has unlimited support, your supposition that all bins except the closest two would be zero is incorrect.

Finding a closed-form analytic solution may be impossible, in which case your best bet would be to start with the center frequency for the two strongest bins and use binary search to find the frequency in-between that most closely corresponds to your actual spectrum.

-

Generally you need to apply a windowing function to your signal before performing the FFT, to get clean spikes for the frequency components (without "skirts"). (The only exception is if your frequency components fit exactly into the length of the FFT, so that they each fill exactly one bin.)

If you use a Gaussian windowing function, and then fit a parabola to the highest three points in your FFT, you can get a theoretically exact result for the frequency:

-

I just noticed this old question and thought I'd expand on J.M.'s comment (that is, what I think he/she was hinting at).

First of all, a small remark on the "bins" you talk about: The frequencies associated with the coefficients should be thought of as lying in the center of their "bin," and you're also off by a factor of 2 unless I'm missing something. So the first bin would be 0-10Hz, the second 10-29Hz, the third 29-49Hz, and so on. However, as others have pointed out, anything but a sine/cosine wave with an exactly matching frequency will appear in more than one "bin" anyway, so it's better to drop the notion of a "bin" and just think of frequencies.

This is because one way to interpret the FFT is that it decomposes the signal into a sum of cosine waves with certain amplitudes (given by the absolute values of the coefficients) and time offsets (given by the phase). For example, if your input is a 39Hz sine wave with an amplitude of 1, the third coefficient will have an absolute value of 1 (or 512, depending on the algorithm) and a phase of $-\pi/2$ (or $\pi/2$). This translates to a cosine wave that is shifted to the right by $\pi$, giving a sine wave, and then stretched appropriately in both directions.

If you slowly shift your FFT window to the right, the phase will decrease, wrap around (because a phase of $2\pi$ is the same as a phase of $0$), and return to $-\pi/2$ after 1/39s or 256 samples.

The same also works for frequencies that are not exact multiples of 19.5Hz. So you can determine how many samples it takes for the phase to wrap around, and this will give you a frequency that is as accurate as you want: To increase accuracy, just let it wrap around once more. By the way, you don't actually need to compute the phase; you can just check when the real or imaginary part (whichever you like) crosses zero.

-

You can try a zero-padded FFT ( add zeros to the end of your signal and then take FFT).

-

In the end I derived a formula that calculates the exact frequency from the rate at which each bin rotates.

To understand bin rotation, just look at a single bin. So let's consider theta = 2pi/12, so exp(i.n.theta) for n in {0, ..., 11} will give us our twelve 12th-roots of unity; if you look on the complex plane it will make a 12 pointed star.

Now if the signal matches that frequency, Ie 12 samples makes one wave, then repeatedly projecting that signal onto that star ( proj[k] = sig[k]*star[k%12] ) will stretch out one spine (the one corresponding to the maximum positive amplitude of sig[k], the peak of the sine wave if you like) and will squash or negatively stretch (depending on whether the sine wave is vertically centred, i.e. DC offset = 0) the diametrically opposite spine.

If you keep a running vector sum of the last 12 projections, the result will come out in this direction, by symmetry. And the angle of the result gives the bin phase, and length gives the bin magnitude.

And the phase remains constant. So no bin rotation.

But now if you change the input frequency to {13 samples = one wave}, you can see that this vector sum is now going to rotate 1/12 of a cycle If a peak corresponds to an angle of 0, the next peak will no longer line up, it will correspond to 1/12 of a revolution

And this is the bin rotation.

While I wouldn't expect the reader to derive the answer from this basis, it will hopefully give some intuitive understanding of what is going on.

-