2
$\begingroup$

I would to calculate the PSD of a signal using FFT however the result do not match with periodogram command. What did was as follow :

len=length(signal);
NFFT=2^10;
signal_fft=fftshift(fft(signal,NFFT));
px=signal_fft.*conj(signal_fft)/(NFFT*len);
f=-[NFFT/2:NFFT/2-1]/NFFT*Fs;
plot(f,px);

However I am wondering :

Should I use Fs instead of NFFT in the denominator for calculation signal_fft? Since most websites mentioned it should be Fs, but according to this link the number of FFT points are used. However NFFT makes more sense since that is number of points I have.

The reason is that if I take 1024 point FFT the resolution is different from $2^{20}$ points FFT so NFFT should be used somewhere to adjust for this change.

$\endgroup$
  • $\begingroup$ Duplicate of what ? I read the forum, no question is exactly like this. I am asking for more clarification about calculating PSD. $\endgroup$ – justin Jul 19 '16 at 14:53
  • $\begingroup$ Duplicate of post linked at the top of this question. $\endgroup$ – jojek Jul 19 '16 at 14:59
  • $\begingroup$ I understand I am the reading the post right now. But in my question I am referring to something else I added the link that I was referring to and what was used in that website to calculate psd. I need to see if somebody can explain about that. I am asking if I have N points FFT should I normalize the psd by NFFT or by Fs. $\endgroup$ – justin Jul 19 '16 at 15:05
  • $\begingroup$ You never normalize by sampling frequency. If you don't do any windowing, then it is very straightforward. Otherwise please see answer 1 and answer 2. $\endgroup$ – jojek Jul 19 '16 at 15:17
  • $\begingroup$ Thanks no I am not using any windowing. However in your duplicate link you posted as well as Mathwork website it mentioned that psd should be scaled by Fs X length(signal) however in this website gaussianwaves.com/2014/07/… it says it should be scaled by number of FFT points X length of signal. $\endgroup$ – justin Jul 19 '16 at 15:22
8
$\begingroup$

There is only one correct way of scaling DFT when calculating PSD with RMS values. Given input signal $x$ and its DFT $X$, the exact formula is:

$$\mathrm{PSD}=\frac{2\cdot \hat{X}}{f_s\cdot S} $$

where:

$\hat{X}=|X|^2=X\cdot X^*$ - squared spectrum magnitude

$S=\sum_{i=1}^{N} w_{i}^{2}$ - scaling factor defined as sum of squared samples of window function

$f_s$ - sampling frequency


Scaling coefficient $2$ takes into account removal of energy at negative frequencies (we drop this side of PSD).

Sampling frequency $f_s$ happens to be in the equation since it is a part of the Effective Noise Bandwidth, and PSD is simply PS (Power Spectrum) divided by ENBW. So in the end, the unit of PSD is $x^2/\mathrm{Hz}$, where $x$ is the quantity that you want to measure, i.e. Volts.

Scaling factor $S$ is compensating for any loss of energy when applying a window function. Obviously for rectangular window (or no window in fact) this doesn't change anything and one, might simply use $N$ instead (length of input sequence).


You did also mention that you were planning to use high number of DFT points ($2^{20}$) since it "makes more sense". Well, it doesn't. If your input sequence has length $N$ (assume 1024), then you can't get better frequency resolution than:

$$\Delta f = \frac{f_s}{1024}$$

Even if you increase number of DFT points to $2^{20}$ by zero-padding, those spectral samples in between will be only interpolated window spectrum - useless garbage. For more detailed explanation please see my other answer.

Another thing to remember. Zero-padding doesn't change the scaling factor! Samples are added after the input sequence had been multiplied by the window, hence the total energy will not change.


Since you wanted to match the calculation with MATLAB's periodogram function, here comes the code. It is a modified version of this one. I tried to comment the code.

% Parameters and signal
Fs = 1000; 
dt = 1/Fs; 
t = 0:dt:1-dt;
x = cos(2*pi*100*t)';
N = length(x);
% Add a bit of noise
x = x + 0.001*randn(N, 1);
% Use window function of your choice. Rectangular window means no window.
win = rectwin(N); 
% The actual number of points used when calculating the DFT might be higher
% this can help in speeding up the computation (next power of 2), but won't
% increase the frequency resolution (only frequency spacing is affected).
% See: https://dsp.stackexchange.com/a/32161/8202
NFFT = N; % N=2^20
xw = x.*win; % Apply windowing here (periodogram does it on it's own)
xdft = fft(xw, NFFT);
psdx = xdft.*conj(xdft); % psdx = abs(xdft).^2
psdx = psdx(1:NFFT/2+1); % Take only first half of the PSD
% Calculate the scaling factor. 
% Note that in case of rectangular window this is simply N.
S = sum(win.^2);
psdx = (2/(Fs*S)) * psdx;
freq = 0:Fs/NFFT:Fs/2;
% Calculate the corresponding PSD using MATLAB's in-built
[pxx, ~] = periodogram(x, rectwin(N), NFFT, Fs);

% Plot the PSD's
plot(freq, 10*log10(psdx), 'r')
grid on
hold on
plot(freq, 10*log10(pxx), 'b--')
title('PSD calculation comparison')
xlabel('Frequency [Hz]')
ylabel('Power/Frequency [dB/Hz]')
legend({'DFT', 'periodogram'})

% Plot error between PSD's
figure
diff = pxx - psdx;
plot(freq, diff)
grid on
xlabel('Frequency [Hz]')
ylabel('Difference')
title('Difference between PSDs')

Please note the perfect match between both PSD's:

enter image description here

And here is the difference:

enter image description here

Probably you can notice that there is a bigger error for 1'st and last bin (Nyquist in fact). That is because we multiplied all frequency bins by scaling factor $2$ to compensate for negative frequencies. This shouldn't be done for DC component (1'st bin) and for Nyquist bin (if it's present) - both of them are unique.

But in a real life, Nyquist bin will be removed by anti-aliasing filter and people will remove mean from the signal anyway. Also DC scaling for RMS values is pointless, especially when doing windowing.

$\endgroup$
  • $\begingroup$ if you increase N (i.e. take more number of samples) the power of peak doesn't remain same. How do you explain that? $\endgroup$ – reluctant_programmer May 14 '18 at 16:00
  • $\begingroup$ @reluctant_programmer: Where exactly are you changing the N? $\endgroup$ – jojek May 14 '18 at 16:45
  • $\begingroup$ If you change the line NFFT = N to NFFT = 2^20 then you will get the exact same peak at -3dB. You must be doing something wrong. $\endgroup$ – jojek May 14 '18 at 17:24
2
$\begingroup$

I'll try to explain how one arrives at the correct scaling factor. In order to keep things simple I will assume a rectangular window. Other windows can be taken into account by an additional factor that is determined by the window weights (as shown in jojek's answer).

The periodogram for estimating the power spectral density (PSD) of a signal $x(t)$ is given by

$$P_x(f)=\frac{1}{T}\left|\int_0^Tx(t)e^{-j2\pi ft}dt\right|^2\tag{1}$$

Using samples $x[n]:= x(nT_s)$, $n=0,1,\ldots,N-1$, and approximating the integral in $(1)$ by a sum, we get

$$P_x(f)\approx \frac{1}{NT_s}\left|\sum_{n=0}^{N-1}x[n]e^{-j2\pi fnT_s}\cdot T_s\right|^2=\frac{T_s}{N}\left|\sum_{n=0}^{N-1}x[n]e^{-j2\pi fnT_s}\right|^2\tag{2}$$

Evaluating $(2)$ at equidistant frequencies $f_k=kf_s/K$, $k=0,1,\ldots K-1$, with $f_s=1/T_s$, and $K\ge N$ gives

$$P_x\left(\frac{kf_s}{K}\right)\approx \frac{T_s}{N}\left|\sum_{n=0}^{N-1}x[n]e^{-j2\pi kn/K}\right|^2=\frac{T_s}{N}\left|X[k]\right|^2\tag{3}$$

where $X[k]$ is the length $K$ DFT of the length $N$ sequence $x[n]$, where $x[n]$ is zero-padded to length $K$ if $K>N$.

Eq. $(3)$ shows that the correct scaling factor is $T_s/N$ (for a rectangular window), and for the two-sided PSD. For the one-sided PSD you just need to add a factor of $2$, except for DC and Nyquist. Furthermore, the scaling factor is independent of the DFT size $K$, but it only depends on the number of time domain samples $N$.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.