The FFT and Digital Sampling
March 29, 2018
The FFT is a computer algorithm developed by Cooley and Tukey  to efficiently compute the coefficients of the Fourier Series representation of a sequence of digital samples of a signal. The number of samples, N, used in the FFT must be an integer power of 2. That is, N = 2p, where p is a positive integer. This minimizes the number of multiplications (and therefore the computation time) needed to compute the coefficients of the Fourier Series.
The Fourier Series represents a time signal x(t) by the sum of cosine and sine waves:
where x̅ is the mean value. In theory this sum has an infinite number of terms but in practice it is limited by the Sample Rate (SR). The usable frequencies in the FFT are limited to those less than ½ SR ≡ fq, called the Nyquist frequency.
As an example consider the N = 8 point sample of a signal shown in Fig. 8. (In practice the number of samples is often N = 1k, 2k, 4k, or more, where k = 1024.) It is customary to number the samples starting with 0, so the samples are labeled xn , n = 0, …, 7. In this example the samples are spaced in time by Δt = 0.5 msec, so SR = 1 / Δt = 2000 samples/sec.
The total time duration of this digital sample is T = N *Δt = 4.0 msec. However, there is no sample at T = 4.0 msec. This is because the N samples are assumed to start at t = 0. The Fourier Series assumes the signal is periodic in time T, so the sample at T would, by definition, be the same value as the first sample (illustrated by the dashed point in Fig. 1). Since this is not the case for an arbitrary random signal, a potential error is introduced. Special attention must be paid to this issue, which is covered in the Windowing lesson.
The usable frequencies in this example are shown in Fig. 2 where the amplitude coefficients Ci and Si are to be determined by the FFT. Summing all of the cosine and sine waves with their determined coefficients will give a result which exactly duplicates the signal at the sampled points.
Each coefficient is determined by multiplying Eq. 1 by the associated cosine or sine wave and then computing the average value of these products over the N samples. The result is
where f j = j / T, j = 1, …, N/2 – 1
The right hand side is greatly simplified because the average value of the product of two different cosine and sine waves is zero, while the average value of the product of a cosine or sine wave with itself is 1/2. (This assumes the averages are taken over an integer number of cycles, which is the case in the FFT.) The left hand side requires a large number of multiplications to complete. However, in the FFT algorithm the requirement that N is an integer power of 2 means that a large number of these computations are redundant and can be avoided. More information on the details of the FFT computation algorithm can be found in reference .
Continuing with the computation of the Fourier Series coefficients using the long hand version of Eq. 2 gives
x(t) = x0, x1, x2, x3, x4, x5, x6, x7 = 1, 2, 0, -1, -3, -2, 1, -1
x̅ = (x0 + x1 + x2 + x3 + x4 + x5 + x6 + x7) / N = – 0.375 (computed as C0 at f0 = 0 in the FFT)
C1 = (2/N) (x0 + √ 2 x1 – √ 2 x3 – x4 – √ 2 x5 + √ 2 x7) = 1.707 f1 = 1 / T = 250 Hz
S1 = (2/N) (√ 2 x1 + x2 + √ 2 x3 – √ 2 x5 – x6 – √ 2 x7) = – 0.457
C2 = (2/N) (x0 – x2 + x4 – x6) = – 0.75 f2 = 2 / T = 500 Hz
S2 = (2/N) (x1 – x3 + x5 – x7) = – 0.5
C3 = (2/N) (x0 – √ 2 x1 + √ 2 x3 – x4 + √ 2 x5 – √ 2 x7) = 0.293 f3 = 3 / T = 750 Hz
S3 = (2/N) (√ 2 x1 – x2 + √ 2 x3 – √ 2 x5 + x6 – √ 2 x7) = 0.957
This illustrates the redundancy in the calculations. It should be noted that some FFT algorithms (such as those in Excel) do not multiply by (2/N) in their computations but only return the summation values in Eq. 2.
The values of the PSD are the mean square amplitudes of the FFT at each frequency divided by the frequency spacing, Δf = 1 / T (Eq. 3). The value at f0 = 0 is computed differently because it represents the mean value given by C0 (since cos(0) = 1).
The PSD for the example in Fig. 8 is shown in Fig. 310 (where A is used for the units of x). Each PSD value is associated with a frequency band centered at the associated FFT frequency. The frequency bandwidths are equal to Δf (except for f = 0, where the bandwidth is Δf / 2).
The FFT actually computes values for N frequencies, f j = j /T, j = 0, …, N–1, but only the first N/2 values are needed for the PSD. In addition, the signal x(t) must be low-pass filtered before the digital sampling in order to assure there is no frequency content above the Nyquist frequency. Otherwise, errors will occur in the FFT computation. This is discussed further in the Aliasing lesson.
The mean square, standard deviation, and rms values of the time signal can be (approximately) recovered from the PSD. Since the PSD is the mean square value of the Fourier Series components, the mean square of the signal is found by summing the PSD values and multiplying by the frequency band width, Δf = 1 / T , (Eq. 4). The rms value is given by by xrms = .
The standard deviation is given by . Table 1 compares these values computed from the signal in Fig. 1 and the PSD in Fig. 10. The small differences are due to the fact that the FFT uses the presumed signal value at T, whereas the calculation from the signal samples does not. This is discussed further in the Windowing lesson.
|Using xn||Using PSD|
|Mean Square Value, x̅2̅||2.625||2.609|
|Mean Square Value, x̅||– 0.375||– 0.375|
|Standard Deviation, σ||1.576||1.571|
[Note: In many random vibration signals the mean value is assumed to be zero so xrms = σ, but that is not the case in this example.]
- Cooley, J. W., Tukey, J. W., “An algorithm for the machine calculation of complex Fourier series,” Math. Comput. 19, 297–301, 1965.
- Smith, S. W., The Scientist and Engineer’s Guide toDigital Signal Processing, California Technical Pub; 1st ed., 1997.