psd (version 2.1.0)

pspectrum: Adaptive sine multitaper power spectral density estimation

Description

This is the primary function to be used in this package: it returns power spectral density estimates of a timeseries, with an optimal number of tapers at each frequency based on iterative reweighted spectral derivatives. If the object given is a multicolumn object, the cross spectrum (multivariate PSD) will be calculated using the same iterative procedure.

Usage

pspectrum(x, ...)
# S3 method for ts
pspectrum(x, output_column = NULL, ...)
# S3 method for matrix
pspectrum(x, x.frqsamp, ...)
# S3 method for spec
pspectrum(x, ...)
# S3 method for default
pspectrum(
  x,
  x.frqsamp = 1,
  ntap.init = NULL,
  niter = 3,
  output_column = NULL,
  AR = FALSE,
  Nyquist.normalize = TRUE,
  verbose = TRUE,
  no.history = FALSE,
  plot = FALSE,
  ...
)
pspectrum_basic(x, ntap.init = 7, niter = 5, verbose = TRUE, ...)
adapt_message(stage, dvar = NULL)

Arguments

x

vector; series to find PSD estimates for; if this is a multicolumn object, a cross spectrum will be calculated.

...

Optional parameters passed to riedsid2

output_column

scalar integer; If the series contains multiple columns, specify which column contains the output. The default assumes the last column is the output and the others are all inputs.

x.frqsamp

scalar; the sampling rate (e.g. Hz) of the series x; equivalent to frequency.

ntap.init

scalar; the number of sine tapers to use in the pilot spectrum estimation; if NULL then the default in pilot_spec is used.

niter

scalar; the number of adaptive iterations to execute after the pilot spectrum is estimated.

AR

logical; should the effects of an AR model be removed from the pilot spectrum?

Nyquist.normalize

logical; should the units be returned on Hz, rather than Nyquist?

verbose

logical; Should messages be given?

no.history

logical; Should the adaptive history not be saved?

plot

logical; Should the results be plotted?

stage

integer; the current adaptive stage (0 is pilot)

dvar

numeric; the spectral variance; see also vardiff etc

Value

Object with class 'spec', invisibly. It also assigns the object to "final_psd" in the working environment.

Details

See the Adaptive estimation section in the description of the psd-package for details regarding adaptive estimation.

NEW as of version 2.0: use pspectrum to calculate the cross spectrum if x is a multi-column array.

pspectrum_basic is a simplified implementation used mainly for testing.

See Also

psdcore, pilot_spec, riedsid2, prewhiten

Examples

# NOT RUN {
#REX
library(psd)
library(RColorBrewer)

##
## Adaptive multitaper PSD estimation
## (see also the  "psd_overview"  vignette)
##

data(magnet)
Xr <- magnet$raw
Xc <- magnet$clean

# adaptive psd estimation (turn off diagnostic plot)
PSDr <- pspectrum(Xr, plot=FALSE)
PSDc <- pspectrum(Xc, plot=FALSE)

# plot them on the same scale
plot(PSDc, log="dB",
     main="Raw and cleaned Project MAGNET power spectral density estimates",
     lwd=3, ci.col=NA, ylim=c(0,32), yaxs="i")
plot(PSDr, log="dB", add=TRUE, lwd=3, lty=5)
text(c(0.25,0.34), c(11,24), c("Clean","Raw"), cex=1)

## Change sampling, and inspect the diagnostic plot
plot(pspectrum(Xc, niter=1, x.frqsamp=10, plot=TRUE))

## Say we forgot to assign the results: we can recover from the environment with:
PSDc_recovered <- psd_envGet("final_psd")
plot(PSDc_recovered)

# }
# NOT RUN {
#REX
# }