spectral_properties: Calculate properties of multitaper power spectral density...
In williamdeleo/psd: Adaptive, Sine-Multitaper Power Spectral Density Estimation

Description Usage Arguments Details Value Author(s) See Also Examples

Various spectral properties may be computed from the vector of tapers, and if necessary the sampling frequency.

spectral_properties(x, ...)

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

## S3 method for class 'tapers'
spectral_properties(x, ...)

## Default S3 method:
spectral_properties(x, f.samp = 1, n.freq = NULL,
  p = 0.95, db.ci = FALSE, ...)

`x`	object to calculate spectral properties for; or a vector of number of tapers
`...`	additional arguments
`f.samp`	numeric; the sampling frequency (e.g. Hz) of the series the tapers are for
`n.freq`	integer; the number of frequencies of the original spectrum (if `NULL` the length of the tapers object is assumed to be the number)
`p`	numeric; the coverage probability, bound within [0,1)
`db.ci`	logical; should the uncertainty confidence intervals be returned as decibels?

Parameter Details:

Uncertainty

See spec_confint for details.

Resolution

The frequency resolution depends on the number of tapers (K), and is found from

\frac{K \cdot f_N}{N_f}

where f_N is the Nyquist frequency and N_f is the number of frequencies estimated.

Degrees of Freedom

There are two degrees of freedom for each taper K:

ν = 2 K

Bandwidth

The bandwidth of a multitaper estimate depends on the number of tapers. Following Walden et al (1995) the effective bandwidth is \approx 2W where

W = \frac{K + 1}{2 N}

and N is the number of terms in the series, which makes N \cdot W the approximate time-bandwidth product.

A list with the following properties (and names):

taper: the number of tapers
stderr.chi .upper, .lower, .median: results returned from spec_confint
resolution: effective spectral resolution
dof: degrees of freedom; will be slightly inaccurate for single-taper periodograms
bw: effective bandwidth of the spectrum

A.J. Barbour

spec_confint, psd-package

## Not run: #REX
library(psd)

##
## Spectral properties from the number of tapers used
## (portions extracted from overview vignette)
##

#
# Theoretical uncertainties from Chi^2 distribution
#
sp <- spectral_properties(as.tapers(1:50), p=0.95, db.ci=TRUE)
par(las=1)
plot(stderr.chi.upper ~ taper, sp, type="s",
       ylim=c(-10,20), yaxs="i", xaxs="i",
       xlab=expression("number of tapers ("* nu/2 *")"), ylab="dB",
       main="Spectral uncertainties")
lines(stderr.chi.lower ~ taper, sp, type="s")
lines(stderr.chi.median ~ taper, sp, type="s", lwd=2)
lines(stderr.chi.approx ~ taper, sp, type="s", col="red",lwd=2)

#
# An example using the Project MAGNET dataset
#
data(magnet)
tapinit <- 15 # tapers
dt <- 1 # 1/km

# remove mean/trend (not really necessary but good practice; also, done internally)
ats <- prewhiten(ts(magnet$clean, deltat=dt), plot=FALSE)$prew_lm

# normal and adaptive multitaper spectra
Pspec <- psdcore(ats, dt, tapinit)
Aspec <- pspectrum(ats, dt, tapinit, niter=3, plot=FALSE)

# calculate spectral properties
spp <- spectral_properties(Pspec$taper, db.ci=TRUE)
spa <- spectral_properties(Aspec$taper, db.ci=TRUE)

# function to create polygon data, and create them
pspp <- create_poly(Pspec$freq, dB(Pspec$spec), spp$stderr.chi.approx)
psppu <- create_poly(Pspec$freq, dB(Pspec$spec), spp$stderr.chi.upper)
pspa <- create_poly(Aspec$freq, dB(Aspec$spec), spa$stderr.chi.approx)
pspau <- create_poly(Aspec$freq, dB(Aspec$spec), spa$stderr.chi.upper)

##
## Project MAGNET uncertainties
##
plot(c(0,0.5),c(-8,35),col="white",
       main="Project MAGNET Spectral Uncertainty (p > 0.95)",
       ylab="", xlab="spatial frequency, 1/km", yaxt="n", frame.plot=FALSE)
lines(c(2,1,1,2)*0.01,c(5,5,8.01,8.01)-8)
text(.05, -1.4, "3.01 dB")
polygon(psppu$xx, (psppu$yy), col="light grey", border="black", lwd=0.5)
polygon(pspp$xx, (pspp$yy), col="dark grey", border=NA)
text(0.15, 6, "With adaptive\ntaper refinement", cex=1.2)
polygon(pspau$xx, (pspau$yy)-10, col="light grey", border="black", lwd=0.5)
polygon(pspa$xx, (pspa$yy)-10, col="dark grey", border=NA)
text(0.35, 22, "Uniform tapering", cex=1.2)

##
## Project MAGNET resolution
##
frq <- Aspec$freq
relp <- dB(1/spa$resolution)
par(las=1)
plot(frq, relp,
     col="light grey",
     ylim=dB(c(1,5)),
     type="h", xaxs="i", yaxs="i",
     ylab="dB", xlab="frequency, 1/km",
     main="Project MAGNET Spectral Resolution and Uncertainty")
lines(frq, relp)
lines(frq, spp$stderr.chi.upper+relp, lwd=1.5, lty=3)
lines(frq, spa$stderr.chi.upper+relp, lwd=3, lty=2)
abline(h=dB(sqrt(vardiff(Aspec$spec))), lwd=1.5, lty=2, col="red")

##

## End(Not run)#REX