How the PSD is Computed
March 29, 2018
The essential procedure in computing the PSD of a signal is determining the frequency spectrum. This is done using a 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.) Graphically this is illustrated in Fig. 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 Fig. 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
where the summation, , computes 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 Fig. 3. All of the terms on the right hand side of Eq. (2) have a mean value of zero except for the C37 term. Then
where the fact that the mean-square value of the cosine wave is 1/2 has been used (as illustrated in Fig. 2). In general, the coefficients are determined by
Historically, the computation of the summations in Eq. (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 1 Hz 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.
- Fourier, J., “Mémoire sur la propagation de la chaleur dans les corps solides,” Présenté le 21 décembre 1807 à l’Institut national – Nouveau Bulletin des sciences par la Société philomatique de Paris I (6), 112–116, 1808.
- Cooley, J. W., Tukey, J. W., “An algorithm for the machine calculation of complex Fourier series,” Math. Comput. 19, 297–301, 1965.