Extracting exact frequencies from FFT output

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:

Quadratic Interpolation of Spectral Peaks

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.