Effacer les filtres
Effacer les filtres

Power spectral density of gaussian white noise

74 vues (au cours des 30 derniers jours)
TheBeginner
TheBeginner le 8 Nov 2013
Réponse apportée : Jeremy le 18 Juin 2015
Hi, I just wanted to check that the matlab function "pwelch" gives a correct estimates of the PSD of a gaussian white noise. I plot the estimate of the PSD and also the variance, which is supposed to be equal to the mean of PSD. I always have a bias that I don't understand.
Here's my code :
sigma = 1000;
L = 200000;
%PSD
noise= randn(1,L) * sigma;
[PSD_noise,vect_freq_noise] = pwelch(noise,75000,[],4000);
%Display
plot(vect_freq_noise,2*PSD_noise,'-g'), %PSD
%Theoritical value of the PSD
hold on, plot([vect_freq_noise(1) vect_freq_noise(end)],[sigma^2 sigma^2],'r');
%Value of the mean
plot([vect_freq_noise(1) vect_freq_noise(end)],[mean(2*PSD_noise) mean(2*PSD_noise)],'b');
Any idea where's the mistake?
Thank you
  1 commentaire
TheBeginner
TheBeginner le 8 Nov 2013
I have also tried to calculate myself the PSD :
sigma = 1000;
L = 200000;
noise= randn(1,L) * sigma;
PSD = abs(fft(noise)).^2/length(noise);
vect_freq = linspace(0,pi,length(PSD));
plot(vect_freq,PSD,'-.g');
hold on, plot([vect_freq(1) vect_freq(end)],[sigma^2 sigma^2],'r');
hold on, plot([vect_freq(1) vect_freq(end)],[mean(PSD) mean(PSD)],'-.b');
And it works fine. I think there's a normalization factor in pwelch that I can't quite figure out.

Connectez-vous pour commenter.

Réponses (1)

Jeremy
Jeremy le 18 Juin 2015
The "density" in PSD means that the power is normalized to something, usually 1 Hz, but in this case it is the Nyquist frequewncy since there was sampling rate input into pwelch. The energy of white noise will be spread over all frequencies so you need to look at the integral of the signal:
sum(PSD_noise*.0016) % this should equal omega^2

Catégories

En savoir plus sur Parametric Spectral Estimation dans Help Center et File Exchange

Tags

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by