Why are the results of two spectral density estimation methods（periodogram and pwelch） inconsistent? If the amplitude of the power spectrum density of random signals is concerned, which spectral estimation is accurate?

I want to calculate the power spectral density of the noise voltage. The magnitude of the power spectral density is very important to me. I want to use the two calculation methods （periodogram and pwelch） in the Matlab example. It is found that the results obtained from these two methods are inconsistent with the same signal.

periodogram method：

Fs = 1000;
t = 0:1/Fs:1-1/Fs;
x = cos(2*pi*100*t)+ randn(size(t));
N = length(x);
xdft = fft(x);
xdft = xdft(1:N/2+1);
psdx = (1/(Fs*N)) * abs(xdft).^2;
psdx(2:end-1) = 2*psdx(2:end-1);
freq = 0:Fs/length(x):Fs/2;
subplot(2,1,1);
plot(freq,10*log10(psdx))
grid on
title('Periodogram Using FFT')
xlabel('Frequency (Hz)')
ylabel('Power/Frequency (dB/Hz)')

subplot(2,1,2);
periodogram(x,rectwin(length(x)),length(x),Fs);
grid on
title('Periodogram PSD Estimate');
xlabel('Frequency (Hz)');
ylabel('Power/Frequency (dB/Hz)');

pwelch method：

fs = 1000;
t = 0:1/fs:1-1/fs;
x = cos(2*pi*100*t) + randn(size(t));
[pxx,f] = pwelch(x,500,300,500,fs);

plot(f,10*log10(pxx))

xlabel('Frequency (Hz)')
ylabel('PSD (dB/Hz)')

The programs are both the examples in MATLAB，I found them on this website。 But the result is different. The result of periodogram is -3.214dB, and the result of pwelch is -6.616dB. The difference between them is 3dB. Why? How does the amplitude of the power spectral density be calculated? Is there an exact answer? I would be grateful if someone could point me in the right direction. Thank you in advance!

Answer 1

[Honglei Chen]

Your bin widths are different. In Welch's method, you are doing in only 500 points FFT. However, the power is given by multiplying PSD with the frequency bin width. If you do that, the power is preserved.

HTH

[bingqing zou]

Thank you for your answer. I don't understand what do you mean bin width. Can you tell me more about it in more detail? Is bin width means the frequency resolution or the width of the window? Thank you!

[Honglei Chen]

Yes. In periodogram, you sample the [0 Fs] with 1000 points, so the bin width, i.e., the frequency resolution for each sampling point, is 1 Hz. In pwelch, since you are using only 500 points, each bin becomes 2 Hz.

HTH

[bingqing zou]

Thank you for your clear answer! But if I change the frequency resolution of the two methods to 1Hz, the result is still different, however， the gap is narrowed. The difference is about 1.2dB. The result is obtained by running the program many times.Why are the results of two spectral density estimation methods（periodogram and pwelch）still inconsistent?

fs = 1000; N=1000; t = 0:1/fs:1-1/fs; x = cos(2*pi*100*t)+ randn(size(t));

subplot(2,1,1); periodogram(x,rectwin(length(x)),length(x),fs); grid on title('Periodogram PSD Estimate'); xlabel('Frequency (Hz)'); ylabel('Power/Frequency (dB/Hz)');

[pxx,f] = pwelch(x,1000,600,1000,fs); subplot(2,1,2); plot(f,10*log10(pxx)) grid on xlabel('Frequency (Hz)') ylabel('PSD (dB/Hz)')

[CARMEN OLIVAS]

In pwelch, the window must be a vector. Try with:

[pxx,f] = pwelch(x,rectwin(length(x)),600,1000,fs); subplot(2,1,2); plot(f,10*log10(pxx)) grid on xlabel('Frequency (Hz)') ylabel('PSD (dB/Hz)')

(:

[tom a]

Because the default window of pwelch is Hamming window which is different with Rectangular window.

Answer 2

[Pit Hoffmann]

For those who are still interested in this question: The code below gives the exact same solution for using FFT, Periodogram and Pwelch. Note that I changed 'noverlap' to 0. As the 'window' already has the length of the signal there is no need/opportunity for an overlap. The window length would have to be less to use an overlap. Apart from that is it recomennded to have 2^n coordinates of a signal while using FFT to avoid truncation or trailing zeros [Link]. This would be done be changing 'Fs' (e.g. Fs = 1024 (2^10)).

Fs = 1000;
t = 0:1/Fs:1-1/Fs;
x = cos(2*pi*100*t)+ randn(size(t));
figure;
%% FFT
N = length(x);
xdft = fft(x);
xdft = xdft(1:N/2+1);
psdx = (1/(Fs*N)) * abs(xdft).^2;
psdx(2:end-1) = 2*psdx(2:end-1);
freq = 0:Fs/length(x):Fs/2;
subplot(3,1,1);
plot(freq,10*log10(psdx));
grid on; 
title('Periodogram Using FFT');
xlabel('Frequency (Hz)');
ylabel('Power/Frequency (dB/Hz)');
%% Periodogram
subplot(3,1,2);
periodogram(x,rectwin(length(x)),length(x),Fs);
grid on; 
title('Periodogram PSD Estimate');
xlabel('Frequency (Hz)');
ylabel('Power/Frequency (dB/Hz)');
%% Pwelch
subplot(3,1,3);
[pxx,f] = pwelch(x,rectwin(length(x)),0,length(x),Fs);
plot(f,10*log10(pxx));
grid on; 
title('Pwelch PSD Estimate');
xlabel('Frequency (Hz)');
ylabel('Power/Frequency (dB/Hz)');

[MIn SUN]

I disaggre with it.....

pwelch will make up the loss from adding windows, I validate on my code

so

if you add hanning window ,you shall make 1.63 time for amplitude before fft

then you will see what got from fft method is same as the one from pwelch