psd: pspectrum – R documentation

Become an expert in R — Interactive courses, Cheat Sheets, certificates and more!

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 class 'ts'
pspectrum(x, output_column = NULL, ...)

## S3 method for class 'matrix'
pspectrum(x, x.frqsamp, ...)

## S3 method for class 'spec'
pspectrum(x, ...)

## Default S3 method:
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

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.

Value

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

Author(s)

A.J. Barbour adapted original by R.L. Parker

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)


## End(Not run)#REX

psd

Adaptive, Sine-Multitaper Power Spectral Density and Cross Spectrum Estimation

v2.1.0

GPL (>= 2)

Authors

Andrew J. Barbour [aut, cre] (<https://orcid.org/0000-0002-6890-2452>), Jonathan Kennel [aut] (<https://orcid.org/0000-0003-4474-6886>), Robert L. Parker [aut]

Initial release

2020-06-28