Digital Signal Processing Logic

Fast Fourier Transform (FFT)

["Numerical Methods and Analysis", Buchanan and Turner]. (09 Nov 1994)

An algorithm for computing the Fourier transform of a set of discrete data values. Given a finite set of data points, for example a periodic sampling taken from a real-world signal, the FFT expresses the data in terms of its component frequencies. It also solves the essentially identical inverse problem of reconstructing a signal from the frequency data.

The FFT is a mainstay of numerical analysis. Gilbert Strang described it as "the most important algorithm of our generation". The FFT also provides the asymptotically fastest known algorithm for multiplying two polynomials. Versions of the algorithm (in C and Fortran) can be found on-line from the GAMS server

Steve Thackery Says:

The highest frequency you can extract is half the sampling frequency. In other words, if you sample at 10kHz then the highest frequency you can pull out is 5kHz (this is known as the Nyquist frequency).

The lowest frequency you can extract is 1 / duration of sound bite. So, if your soundbite is 100ms long, then the lowest frequency you can extract is 10Hz. With a 2 second long soundbite your lowest extractable frequency is 0.5Hz, and so on.

The rest of the steps in the FFT output are multiples (harmonics) of that lowest frequency. So, with our 100ms soundbite the first item in the FFT output represents 10Hz, the next represents 20Hz, and so on up to the Nyquist frequency. (This is just a better way of explaining what I said in my previous attempt).

Most FFT engines require you to present them with a soundbite which is exactly a power of 2 samples long. In other words, the soundbite must consist of 4, 8, 16, ........ 256, 512, 1024, 2048, etc samples. If your soundbite isn't the right length, you can either truncate it down to the previous power of 2, or you can pad it with zeroes up to the next power of two. For example, if your soundbite is (say) 700 samples long you must either chop it back to 512, or pad it to 1024. Generally padding is the better bet if you can spare the processing time. Both can cause complications if you are using a "window" (which is too heavy to go into now).

The industry standard sampling rates we tend to use (e.g. 11025Hz, 22050Hz, etc), combined with the need to adjust the length of the soundbite to be a power of 2, means that the frequencies in the FFT output are usually fractional numbers rather than integers. For example, a soundbite of 1024 samples, which was recorded at 11025Hz, gives us a bottom frequency of 10.766Hz, the next frequency is 21.533Hz, and so on up to 5512.5Hz.

Finally, to repeat: filter out everything above the Nyquist frequency *before* you do the sampling. In other words, the filtering needs to be done in the analogue domain. For example, with a 11025Hz sampling frequency, you would filter out everything above 5.5kHz. The reason for this is that frequencies above the Nyquist frequency get "mirrored" back into the 0 to 5kHz range and give you misleading results.

Robert A. LaBudde, PhD, PAS, Dpl. ACAFS says:

If you only need a feew points in the spectrum (< 10), you can get by with a direct DFT without the fancy FFT logic.

You will need:

1. Table of 8-bit sines and cosines incremental frequencies (can store in EEROM).
2. 8 x 8 multiply routine (available from MicroChip app note).
3. A 16-bit add routine (ditto).
4. Directly compute the coefficients via

           a(k) = sum(0,n-1) {x(k) * sin(2 pi k/n) }
           b(k) = sum(0,n-1) {x(k) * cos(2 pi k/n) }
   

5. You compute the power at each frequency by

           p(k) = a(k)*a(k) + b(k)*b(k)
   

6. If you want the rms power, you also need a sqrt routine /techref/method/math/sqrt.htm to calculate

           s(k) = sqrt( p(k) )
   

   .

Another approach would be to switch a more microprocessor friendly analysis, such as use of z-transforms or the Goertzel algorithm instead of Fourier transforms.

[one may also need to worry] about scaling.

See also:

• PIC spectrum analyzer source code /microchip/fft/picspect.htm (Microchip Application Note AN542 doesn't work) http://www.edtn.com/embapps/emba035.htm
• http://www.fftw.org/ GPL'd fast FFT in C

Comments:

  .