Normalization while computing Power Spectral Density

Question

Method to compute power spectral density:-

F = fft (s);

PSD = (1/N) * F * conj(F);

Where "s" is the input signal which is given to me in the form of an array.

I also know the sampling rate (Fs).

I want to know what should be the value of the Normalizing Factor "N".

Answer 1

There are many different definitions for a power spectral density function, and correspondingly different possibilities for the scaling factor. Section 13.4 of Numerical recipes in C lists several common definitions such as:

• defined for discrete positive, zero, and negative frequencies, and its sum over these is the function mean squared amplitude
• defined for zero and discrete positive frequencies only, and its sum over these is the function mean square amplitude
• defined in the Nyquist interval from -f_c to f_c and its integral over this range is the function mean squared amplitude
• defined from 0 to f_c, and its integral over this range is the function mean square amplitude

The right definition and scaling factor would thus be specific to your application. As an illustration of the impact these different definitions can have on the scaling factor, I've listed below some specific implementations which use different definitions.

Since I mentioned the Numerical recipes book, we can start be looking at the definition chosen for the purpose of showing a sample implementation of the PSD (not suggesting that it is the correct definition). In this case the second definition listed above (ie. "defined for zero and discrete positive frequencies only, and its sum over these is the function mean square amplitude") has been used, which leads to the normalization:

len = length(F);
N   = 0.5*len^2;
PSD = (1/N) * F(1:len/2) * conj(F(1:len/2));

Octave's pwelch on the other hand uses a different definition of the power spectral density (namely the last one listed above), which leads to a different normalization approximated by:

len = length(F);
N   = 0.5*len*Fs; % where Fs is the sampling rate
PSD = (1/N) * F(1:len/2) * conj(F(1:len/2));

Answer 2

N is just the number of points in the FFT. So if your FFT has, say, 2048 points, then you need to scale the magnitude of the FFT output bins by a factor of 1 / 2048.