# How the PSD is Computed

March 29, 2018

Back to: Random Testing

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 [1]. 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:

(1)

Equation 1

where *f _{1}* = 1 /

*T*,

*f*=

_{i}*i*

*f*, and

_{1}*C*and

_{i}*S*are the Fourier coefficients to be determined. (x is the mean value of

_{i}*x*.) This is illustrated graphically in Figure 2.6.

In practice, the signal *x* is digitized into a sequence of *N* numbers, *x _{n}*,

*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

*C*

_{37}is determined by Equation 2.

(2)

(3)

Equation 2

The following summation is used to compute the mean value and *t _{n}* =

*n*Δ

*t:*

(4)

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 *C*_{37} term. Then:

(5)

(6)

(7)

Equation 3

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

(8)

(9)

(10)

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) [2]. 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 *f _{q}* = 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

*f*have been removed.

_{q}The FFT computes the values of the Fourier coefficients *C _{i}* and

*S*. The mean square amplitude is given by (

_{i}*C*

_{i}^{2}+

*S*

_{i}^{2}) / 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.*

(11)

Equation 5

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.

**References**

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#/.