The FFT and Digital Sampling

March 29, 2018

The Fast Fourier transform (FFT) is a computer algorithm developed by James Cooley and John Tukey [1]. The algorithm computes the coefficients of a Fourier Series representation of a sequence.

(1)   \begin{equation*}  x(t) = \bar{x} + \sum_{i=1}^{\infty}{[C_{i}\cos(2 \pi f_{i} t)]} + S_{i}\sin(2 \pi f_{i} t)]} \end{equation*}

Equation 1

The number of samples, N, used in the FFT must be an integer power of 2. Therefore, 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; however,  in practice, it is limited by the sample rate (SR). The usable frequencies in the FFT are limited to those less than one-half of the SR. Half of the sampling rate of a system is denoted by the Nyquist frequency,  fq.

Digital Sampling and Potential Discontinuity

In practice, the number of samples is often denoted by N = 1k, 2k, 4k, and so on, where k = 1024.  It is customary to number the samples starting with 0, so the samples are labeled xn,  n = 0, … , 7.

For example, consider the point sample N = 8 in Figure 2.15.  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 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 Figure 2.15). As this is not the case for an arbitrary random signal, a potential error is introduced. The discontinuity between the end of the digitized signal and the beginning of the presumed periodic repetition is covered in the Windowing lesson.

The usable frequencies in this example are shown in Figure 2.16, 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.

Digitally sampled time signal

Figure 2.15. A digitally sampled time signal, where the digital sample is denoted by red dots and x(t) by a blue dashed line.

detailsbehindpsd-figure2-125pct

Figure 2.16. A set of cosine and sine waves used to represent the signal in Figure 2.15. 

Determining the Fourier Series Coefficients

Each coefficient is determined by multiplying equation 1 by the associated cosine or sine wave. Then, the average value of these products over the N samples is computed. The result is equation 2, where  f j = j / T,  j = 1, …, N/2 – 1.

(2)   \begin{equation*}  \frac{1}{N} \sum_{n=0}^{N-1}{[C_{i}\cos(2 \pi f_{j} t_{n}) * x_{n}]} = \frac{1}{2}C_{j} \end{equation*}

(3)   \begin{equation*}  \frac{1}{N} \sum_{n=0}^{N-1}{[C_{i}\sin(2 \pi f_{j} t_{n}) * x_{n}]} = \frac{1}{2}C_{j} \end{equation*}

Equation 2

The right side of equation 2 is simplified. 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 one-half. This assumes the averages are taken over an integer number of cycles, which is the case in the FFT.

The left side of Equation 2 requires many multiplications to complete. However, in the FFT algorithm, N is an integer power of 2 and 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 [2].

Computing Fourier Series Coefficients with Long-Form Equation

Continuing with the computation of the Fourier series coefficients using the long hand version of Equation 2 gives:

(4)   \begin{equation*} x(t)=x0, x1, x2, x3, x4, x5, x6, x7 = 1, 2, 0, -1, -3, -2, 1, -1 \end{equation*}

(5)   \begin{equation*} \bar{x}=(x0 + x1 + x2 + x3 + x4 + x5 + x6 + x7)/N = -0.375 \end{equation*}

(computed as C0 at f0 = 0 in the FFT)

(6)   \begin{equation*} C_{1}=\frac{2}{N}(x0+\sqrt{2}x1-\sqrt{2}x3-x4-\sqrt{2}x5+\sqrt{2}x7)=1.707; f_{1}=\frac{1}{T}=250Hz \end{equation*}

(7)   \begin{equation*} S_{1}=\frac{2}{N}(\sqrt{2}x1+x2+\sqrt{2}x3-\sqrt{2}x5-x6-\sqrt{2}x7)=-0.457 \end{equation*}

(8)   \begin{equation*} C_{2}=\frac{2}{N}(x0-x2+x4-x6)=-0.75; f_{2}=\frac{2}{T}=500Hz \end{equation*}

(9)   \begin{equation*} S_{2}=\frac{2}{N}(x1-x3+x5-x7)=-0.5 \end{equation*}

(10)   \begin{equation*} C_{3}=\frac{2}{N}(x0-\sqrt{2}x1+\sqrt{2}x3-x4+\sqrt{2}x5-\sqrt{2}x7)=0.293; f_{3}=\frac{3}{T}=750Hz \end{equation*}

(11)   \begin{equation*} S_{3}=\frac{2}{N}(\sqrt{2}x1-x2+\sqrt{2}x3-\sqrt{2}x5+x6-\sqrt{2}x7)=0.957 \end{equation*}

These computations illustrate 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 equation 2.

Computing the PSD of the signal

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).

(12)   \begin{equation*} PSD_{x}(0) = \frac {C_{0}^2}{\Delta f} \end{equation*}

(13)   \begin{equation*} PSD_{x}(f_{j}) = \frac {0.5(C_{j}^2 + S_{j}^2)}{\Delta f}; j=1, \dots, N/2-1 \end{equation*}

Equation 3

The PSD for the example in Figure 2.15 is shown in Figure 2.17, 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  Δ / 2).

detailsbehindpsd-figure3-125pct

Figure 2.17. The PSD of the signal in Figure 2.15. 

The FFT computes values for N frequencies,  f 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 idea is discussed further in the Aliasing lesson.

Recovering the Mean Square, Standard Deviation, and RMS

The mean square, standard deviation, and RMS values of the time signal can be (approximately) recovered from the PSD. The PSD is the mean square value of the Fourier series components. Therefore, the mean square of the signal is found by summing the PSD values and multiplying by the frequency bandwidth, Δf  = 1 / T , (Eq. 4).

(14)   \begin{equation*} \bar{x^2} = \Delta f * \sum_{j=0}^{N/2-1} PSD_{x}(f_{j}) \end{equation*}

Equation 4

The RMS value is given by the root mean square of x.

(15)   \begin{equation*} x_{rms} = \sqrt{\frac{1}{n}(x_{1}^2 + x_{2}^2,\dots , + x_{n}^2)} \end{equation*}

The standard deviation is given by:

(16)   \begin{equation*} \sigma = \sqrt{\bar{x^2} - (\bar{x})^2} \end{equation*}

Table 1 compares the values computed from the signal in Figure 2.15 and the PSD in Figure 2.17.  There are small differences because 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

Table 1. Comparison of the values computed from the signal and the PSD.

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.

References

Cooley, J.W. & Tukey, J.W. (1965). An algorithm for the machine calculation of complex Fourier series. Mathematics of Computation, 19, 297-301.

Smith, S. W. (1997). The Scientist and Engineer’s Guide to Digital Signal Processing. (1st ed.). California Technical Pub.