[R] noise poser spectral density

[Yudi Pawitan]

At 10:12 AM 5/11/00 +0100, you wrote:
>I know there are some time series experts on this list, and I am hoping
>one of you might help me with this problem. I am trying to understand the
>relationship between No (noise power spectral density) and Gaussian white
>noise sigma^2.
>
>I want to check numerically the theoretical relationship (from Selin,1965)
>No = sigma^2/(2*W) = sigma^2 * delta t, where
>No is power spectral density
>delta t is sampling interval in sec = 1/(2*W)
>W is the bandwidth: the power spectral density (psd) is flat (equal to No)
>for freqs between -W and W Hz
>
>I generate a vector x of n samples from N(0, sigma^2).
>n<-2048
>sigma<-10
>x<- rnorm(n=n,sd=sigma)
>
>I checked R's fft against the examples in Bracewell 1978 p. 364, e.g.
>{1,0,0,0} gives fft 1/4*{1,1,1,1}
>To make R's fft give this result I do fft(x)/n. Then
>
>psd<-Mod(fft(x)/n)^2

the periodogram should be

   Mod(fft(x))^2/n

so the average periodogram (equal to No in the
white noise case) is sigma^2 (set delta t to 1)

>
>The psd varies wildly across freqs, as expected. I do an average over
>freqs
>mean(psd)
>It turns out that this is roughly equal to sigma^2/n. But
>n = duration/(delta t), so I find:
>			mean(psd) = [sigma^2 * (delta t)]/duration
>However theory says:	mean(psd)=sigma^2*(delta t)
>The numerical result is too small by a factor of duration.
>
>Is the theory wrong?
>Should I just multiply mean(psd) by duration?
>
>I am confused. Thanks for any help.
>
>Bill Simpson
>
>-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
.-.-
>r-help mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
>Send "info", "help", or "[un]subscribe"
>(in the "body", not the subject !)  To: r-help-request at stat.math.ethz.ch
>_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
._._
>
Yudi Pawitan     yudi at stat.ucc.ie
Department of Statistics UCC
Cork, Ireland
Ph   353-21-490 2906
Fax 353-21-427 1040 
-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-help mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !)  To: r-help-request at stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._