How the PSD is Computed
March 29, 2018
When computing the PSD of a signal, it is essential to determine the frequency spectrum. This is done using the Fourier series, named after the mathematician Joseph Fourier who invented the method in 1807 . It is based on the principle that any signal x of length T can be constructed by the sum of a series of sine and cosine waves, each having an integer number of cycles in T. Mathematically this is written as:
where f1 = 1 / T, fi = i f1, and Ci and Si are the Fourier coefficients to be determined. (x is the mean value of x.) This is illustrated graphically in Figure 2.6.
In practice, the signal x is digitized into a sequence of N numbers, xn, n = 1 to N, with a time interval of Δt between samples, as illustrated in Figure 2.7. The value of each Fourier coefficient is determined by multiplying both sides of Eq. (1) by the corresponding sine or cosine wave and computing the mean value over N = T/Δt. For example, the value of C37 is determined by Equation 2.
The following summation is used to compute the mean value and tn = n Δt:
This procedure works because the different sine and cosine functions are independent of each other. That is, the mean value of the product of any two different ones is equal to zero, as illustrated in Figure 2.3. The terms in Equation 2 have a mean value of zero except for the C37 term. Then:
In Equation 3, the mean-square value of the cosine wave is 1/2 has been used (as illustrated in Figure 2.2). In general, the coefficients are determined by the following equation (Equation 4):
Historically, the computation of the summations in Equation 4 was very time consuming until 1965, when Cooley and Tukey published an efficient algorithm for this known as the FFT (fast Fourier transform) . Their algorithm uses a sequence of N digital samples where N is a power of 2 (e.g. 2, 4, 8, 16, 32, etc.). The highest frequency that can be resolved in the FFT is fq = 1 / (2 Δt ) Hz, called the Nyquist frequency, so there are N / 2 frequency values. The signal x must be conditioned by a low pass filter before it is digitized to make sure that frequencies above fq have been removed.
The FFT computes the values of the Fourier coefficients Ci and Si. The mean square amplitude is given by ( Ci2 + Si2 ) / 2. Since the frequency resolution is Δf = 1 / T = 1 / ( N Δt ), the PSD is normalized to 1Hz bandwidth by dividing the mean square amplitude by Δf.
A number of FFT computational details must be addressed as part of the PSD creation process. These FFT details include digital time sampling, aliasing, windowing, frequency resolution, and averaging.
Cooley, J.W., & Tukey, J.W. (1965). An algorithm for the machine calculation of complex Fourier series. Mathematics of Computation, 19, 297-301.
Fourier, J. (1807). Mémoire sur la propagation de la chaleur dans les corps solides. Nouveau bulletin des sciences par la société philomathique de Paris, 1, 112–116. https://www.biodiversitylibrary.org/bibliography/45848#/.