Is there an error in the "Power Spectral Density Estimates Using FFT" article?

In the article "Power Spectral Density Estimates Using FFT" found at https://www.mathworks.com/help/signal/ug/power-spectral-density-estimates-using-fft.html, should the line that says: 'psdx = (1/(Fs*N)) * abs(xdft).^2;' instead say: 'psdx = 1/Fs*N * abs(xdft).^2;' ? This would make it consistent with the "Power Spectum Estimate" block found inside the "doc_phasenoise" Simulink model that is opened by typing "doc_phasenoise" at the MATLAB command line. (See http://www.mathworks.com/help/releases/R2016a/comm/ug/rf-impairments.html )

Answer 1

[Kushagr Gupta]

I understand that the question is regarding whether the variable 'N' should be in the numerator or the denominator of the power spectral density calculation.

The equation (4) in this link explains the relation between PSD and FFT, which points out that N is in the denominator.

Also, the model 'doc_phasenoise' contains N in the denominator inside the "spectral averaging" subsystem.

Hence, I believe it is right to say:

psdx = (1/(Fs*N)) * abs(xdft).^2;

[David]

Thank you for your response Kushagr. I agree that it is right to say psdx = (1/(Fs*N)) * abs(xdft).^2; , where N is the FFT length . However, I believe the N in the denominator inside the “spectral averaging” subsystem in the model ‘doc_phasenoise’ is not the FFT length, but instead the number of frames that are averaged together in the spectral-estimation process, so I do not think this is the explanation. Scaling by the FFT length is effectively done prior to actually performing the FFT, in the “Window Scaling” subsystem. This can be realized by pretending that the window is a constant unit gain, which results in reducing the signal level by a factor of the frame length (which is equivalent to the FFT length). With the actual window, the frame-length scaling still occurs, and there is some additional scaling to account for the shape of the actual window.

I believe the full explanation lies in the fact that two scalings must actually be performed in order to obtain the PSD:

First, the result of the FFT needs to be divided by the FFT length in order to obtain the correct magnitude in the frequency domain. As mentioned above, in the model ‘doc_phasenoise’ this is effectively done by reducing the magnitude of the time-domain signal by a factor of the frame length in the “Window Scaling” subsystem. [Note that if this scaling were to be done after squaring the magnitude of the FFT output, the scaling factor would have to be 1/(FFT_length)^2 instead of 1/(FFT_length).]

The second scaling is to convert the power spectrum into the power spectral density, resulting in units of Watts/Hz. This is done by dividing the power spectrum by Delta_freq, the size of the frequency bins in the frequency domain. Because the FFT assumes Nyquist sampling, Delta_freq = 1/(T*FFT_length), where T is the sample interval. Thus, dividing by Delta_freq is equivalent to multiplying by (T*FFT_length), so this scaling is taken care of by the “Product2” block in the ‘doc_phasenoise’ model.

In the article "Power Spectral Density Estimates Using FFT,” both of the above two scalings are combined into a single scaling factor. Since the FFT-length scaling is done after squaring the magnitude of the FFT output, FFT_length^2 is needed in the denominator for the first scaling. For the second scaling, Delta_freq is needed in the denominator. Combining these two results in a single scaling factor of (1/ FFT_length^2) * (1/ Delta_freq). Since Delta_freq = 1/(T*FFT_length), this becomes (1/ FFT_length^2) *(T*FFT_length), which reduces to T/ FFT_length. But T = 1/Fs, resulting and a single scaling factor of 1/(Fs*FFT_length) that performs both of the above-described scalings, which is what is used in the article.

Thanks again for your answer.

[Shuo Hao]

great expression, many thanks.