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?
73 views (last 30 days)
Show older comments
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/190958/image.png)
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 Comments
Answers (2)
Honglei Chen
on 1 Jun 2018
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 Comments
CARMEN OLIVAS
on 5 Apr 2019
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
on 18 Apr 2020
Because the default window of pwelch is Hamming window which is different with Rectangular window.
Pit Hoffmann
on 26 Apr 2021
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)');
1 Comment
MIn SUN
on 30 Sep 2021
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
See Also
Categories
Find more on Parametric Spectral Estimation in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!