Digital Analysis of Geophysical Signals and Waves

3.6 Calculating the Discrete Fourier Transform with the Fast Fourier Transform

The famous fast Fourier transform (FFT) algorithm was published in 1965 by two Americans, James W. Cooley and John W. Tukey (Cooley and Tukey, 1965). Cooley was an applied mathematician with IBM Corporation and Tukey was a statistician who is credited with introducing the terms "bit" and "software."  The FFT revolutionized digital analysis by greatly reducing the computing time to  calculate a discrete Fourier transform.  It is based on an algorithm for the rapid calculation of the summations

The computation time is especially reduced when the digital sequence consists of N= 2ⁿ samples; i.e., if the number of points is a power of 2.  In this case the number of numerical computations is proportional to nN instead of N² for the direct evaluation of the transform.  This means, for example, that the computation savings are over 99% if N = 1024 = 2¹⁰.  There are many clear descriptions of the algorithm (e.g., Brigham, 1974; Mesko, 1984) which we will not discuss.  Since the FFT is a discrete Fourier transform, its properties are those of the DFT.

We discovered in Figure 3.11f that the discrete frequency values in the frequency domain can be obtained by the multiplication of the Dirac comb with "teeth" spacing of Δf.  This means that the corresponding convolution in the time (or space) domain produces a discrete periodic sequence as in Figure 3.11g.  Thus, the discrete Fourier transform, in essence, treats the finite, discrete signal as one period of an infinite periodic sequence.

Assuming that there are N samples, we can compute the discrete Fourier transform using equation 3.6.1, i.e.,

Practical application requires first choosing the sampling interval Δt and the length of the data, T = NDt.  Choice of Dt must be made with full knowledge of the highest frequency content in the signal and the consequences of aliasing (Appendix B).  The length of the data (i.e., the number of points, N) requires that the necessary frequency resolution, Δf = 1/NΔt be achieved (Section 2.4). 

The discrete Fourier transform of N points will provide values at N distinct frequency values. As we have said earlier, it is clear from the symmetry properties (Section 3.3) that only N/2 of these are independent values.  The frequency values are the zero frequency (d.c.) component, and + 1/NΔt, + 2/NΔt, ..., +(N/2-1)/NΔt, +(N/2)/NΔt.  The latter frequency is, of course, the Nyquist frequency since (N/2)/NΔt = 1/2Δt. The Nyquist frequency is the point of symmetry for the frequency values that are produced by the DFT as illustrated in Figure 3.9b.  Explicitly, when the digital signal, s_k is real the DFT frequency domain results have the following symmetry property

Here, the asterisk when used as a superscript indicates a complex conjugate.  This means that there is a simple change of sign in the imaginary part and the real part retains its sign.  This follows the symmetry discussed in Section 3.3 for the Fourier transform of a real function.  That is, the imaginary part (and phase) is an odd function of frequency and the real part (and amplitude) is an even function.

Figure 3.9 illustrates the symmetry resulting from the discrete Fourier transform including the FFT. As dictated by equation 3.6.3, the values of the real part of the discrete Fourier transform for n>N/2 are actually the even, negative frequency values.  The odd symmetry of the imaginary part of the discrete Fourier transform values results in the negative frequency results graphing as shown in Figure 3.9b.  Often the DFT results are shifted so the results are displayed from –f to +f like the continuous Fourier transform results are graphed (Section 3.4).  However, the format in Figure 3.9a and b is commonly followed since the results are plotted sequentially as a function of n = 0, 1, 2, ..., N-1.  This is the order that the numbers are listed in FFT output.  Always remember that the values for n>N/2 are negative frequencies in reversed order beginning with -(N/2-1) / Lx. The symmetry is described above and plotted in Figure 3.9.  Also, it is important to know that the frequency domain data must be "loaded" into an FFT program in the sequence shown in Figure 3.9b in order to yield the desired time domain sequence shown in Figure 3.9a after discrete inverse Fourier transform.