Borrar filtros
Borrar filtros

Power spectral density of gaussian white noise

74 visualizaciones (últimos 30 días)
TheBeginner
TheBeginner el 8 de Nov. de 2013
Respondida: Jeremy el 18 de Jun. de 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 comentario
TheBeginner
TheBeginner el 8 de Nov. de 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.

Iniciar sesión para comentar.

Respuestas (1)

Jeremy
Jeremy el 18 de Jun. de 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

Categorías

Más información sobre Parametric Spectral Estimation en Help Center y File Exchange.

Etiquetas

Community Treasure Hunt

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

Start Hunting!

Translated by