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?

Question

bingqing zou 2018-6-1

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/403512-why-are-the-results-of-two-spectral-density-estimation-methods-periodogram-and-pwelch-inconsistent

评论： MIn SUN 2021-9-30

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!

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Honglei Chen 2018-6-1

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/403512-why-are-the-results-of-two-spectral-density-estimation-methods-periodogram-and-pwelch-inconsistent#answer_322846

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

5 个评论
显示 3更早的评论隐藏 3更早的评论

CARMEN OLIVAS 2019-4-5

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 2020-4-18

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

请先登录，再进行评论。

Answer 2

Pit Hoffmann 2021-4-26

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/403512-why-are-the-results-of-two-spectral-density-estimation-methods-periodogram-and-pwelch-inconsistent#answer_684960

在 MATLAB Online 中打开

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)');