pspectrum.R

set.seed(123)

library(oce)

## Loading required package: mapproj
## Loading required package: maps
## Loading required package: proj4

library(psd)

## Loading required package: fftw
## Loaded psd (0.4.1) -- Adaptive multitaper spectrum estimation.

Fs <- 1000
t <- seq(0, 0.296, 1/Fs)
x <- cos(2 * pi * t * 200) + rnorm(n=length(t))
xts <- ts(x, frequency=Fs)
par(mar=c(3, 3, 1, 1), mgp=c(2, 0.7, 0))
s <- spectrum(xts, spans=c(3,2), main="", log='no')
grid()
w <- pwelch(xts, plot=FALSE)
lines(w$freq, w$spec, col="red")

S <- pspectrum(x, x.frqsamp=Fs)

## Stage  0 est. (pilot) 
##  environment  ** .psdEnv **  refreshed
## detrending (and demeaning)
## Stage  1 est. (Ave. S.V.R. -13.9 dB) 
## Stage  2 est. (Ave. S.V.R. -32.5 dB) 
## Stage  3 est. (Ave. S.V.R. -36.4 dB) 
## Normalized  single-sided PSD  (PSD)  to single-sided PSD for sampling-freq.  1000

lines(S$freq, S$spec, col='blue')

S <- pspectrum(x, x.frqsamp=Fs, ntap.init=1, niter=2)

## Stage  0 est. (pilot) 
##  environment  ** .psdEnv **  refreshed
## detrending (and demeaning)
## Stage  1 est. (Ave. S.V.R. -14.3 dB) 
## Stage  2 est. (Ave. S.V.R. -30.8 dB) 
## Normalized  single-sided PSD  (PSD)  to single-sided PSD for sampling-freq.  1000

lines(S$freq, S$spec, col='darkgreen')

legend("topright", lwd=par("lwd"), title="Random + 200Hz",
       col=c("black", "red", "blue", "darkgreen"),
       legend=c("spectrum", "pwelech", "pspectrum default", "pspectrum modified"))