4

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

methode
  • 5,348
  • 2
  • 31
  • 42
VikramBishnoi
  • 98
  • 1
  • 8

2 Answers2

4

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 -fc to fc and its integral over this range is the function mean squared amplitude
  • defined from 0 to fc, 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));
SleuthEye
  • 14,379
  • 2
  • 32
  • 61
1

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.

Paul R
  • 208,748
  • 37
  • 389
  • 560
  • My input array consists of only real numbers. Then, will the number of points in the FFT be (length_input_array)/2 ? – VikramBishnoi Jul 01 '15 at 06:30
  • You need to decide what N should be, based on required resolution, any real-time constraints, and the nature of your input data. For real input data use N / 2 output points and scale the magnitude by 2 / N. – Paul R Jul 01 '15 at 06:31
  • Don't I need to consider sampling frequency? – VikramBishnoi Jul 01 '15 at 10:56
  • No - you're just normalising magnitude - sampling frequency is only relevant when converting FFT bin indices to frequencies, See: http://stackoverflow.com/questions/4364823/how-do-i-obtain-the-frequencies-of-each-value-in-a-fft/4371627#4371627 – Paul R Jul 01 '15 at 11:24
  • http://in.mathworks.com/help/signal/ug/psd-estimate-using-fft.html Here they have written psdx = (1/(Fs*N)) * abs(xdft).^2; Can you please explain me the reason behind it? – VikramBishnoi Jul 01 '15 at 12:40
  • It just depends what units you want on your PSD Y axis and whether you care about absolute magnitude measurement or not. But yes, if you want absolute `W/Hz` then scale by 1/Fs as well as 1/N. These are just linear scaling factors though, and there will be other such factors that you will need to take into account (e.g. window gain, calibration constants, etc), so I wouldn't worry about it too much - you can just combine all such factors into a single calibration constant. – Paul R Jul 01 '15 at 13:11