riedsid: Constrained, optimal tapers using the Riedel &...

In abarbour/psd: Adaptive, Sine-Multitaper Power Spectral Density and Cross Spectrum Estimation

Constrained, optimal tapers using the Riedel & Sidorenko–Parker method

Description

Estimates the optimal number of tapers at each frequency of given PSD, using a modified Riedel-Sidorenko MSE recipe (RS-RLP).

Usage

riedsid(PSD, ...)

## S3 method for class 'spec'
riedsid(PSD, ...)

## Default S3 method:
riedsid(
  PSD,
  ntaper = 1L,
  tapseq = NULL,
  Deriv.method = c("local_qls", "spg"),
  constrained = TRUE,
  c.method = NULL,
  verbose = TRUE,
  ...
)

riedsid2(PSD, ...)

## S3 method for class 'spec'
riedsid2(PSD, ...)

## Default S3 method:
riedsid2(
  PSD,
  ntaper = 1L,
  constrained = TRUE,
  verbose = TRUE,
  fast = FALSE,
  riedsid_column = 0L,
  ...
)

Arguments

PSD	vector or class 'amt' or 'spec'; the spectral values used to optimize taper numbers
...	optional arguments passed to constrain_tapers
ntaper	scalar or vector; number of tapers to apply optimization
tapseq	vector; representing positions or frequencies (same length as PSD)
Deriv.method	character; choice of gradient estimation method
constrained	logical; apply constraints with constrain_tapers; FALSE turns off constraints
c.method	string; constraint method to use with constrain_tapers, only if constrained=TRUE
verbose	logical; should messages be printed?
fast	logical; use faster method?
riedsid_column	scalar integer; which column to use in multivariate optimization. If the value is 0 the maximum number of tapers for all columns is chosen. If the value is < 0 the minimum number of tapers for all columns is chosen. If the value is 1, 2, 3, etc. the number of tapers is based on the column selected.

Details

The optimization is as follows. First, weighted derivatives of the input PSD are computed. Using those derivatives the optimal number of tapers is found through the RS-RLP formulation. Constraints are then placed on the practicable number of tapers.

riedsid2 is a new (faster) implementation which does not allow for multiple constraint methods; this is the preferred function to use.

Taper constraints

The parameter c.method provides an option to change the method of taper constraints. A description of each may be found in the documentation for constrain_tapers.

Once can use constrained=FALSE to turn off all taper constraints; this could lead to strange behavior though.

Spectral derivatives

The parameter Deriv.method determines which method is used to estimate derivatives.

• "local_qls" (default) uses quadratic weighting and local least-squares estimation; this can be slower than "spg".

• "spg" uses splineGrad; then, additional arguments may be passed to control the smoothness of the derivatives (e.g spar in smooth.spline).

Value

Object with class 'tapers'

Warning

The "spg" can become numerically unstable, and it's not clear when it will be the preferred over the "local_qls" method, other than for efficiency's sake.

Author(s)

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

See Also

constrain_tapers, resample_fft_rcpp, psdcore, pspectrum

Examples

## Not run: #REX
library(psd)

##
## Riedel-Sidorenko-Parker taper optimization
##

set.seed(1234)
# some params
nd <- 512 # num data
ntap <- 10 # num tapers
nrm <- 40 # sharpness of the peaks rel 2*variance
#
# create a pseudo spectrum
# with broad peaks
x <- 0:(nd-1)
riex <- rnorm(nd) + nrm*abs(cos(pi*x/180) + 1.2)
riex <- riex + 8*nrm*dcauchy(x, nd/3)
riex <- riex + 5*nrm*dnorm(x, nd/2)
# some flat regions
riex[riex<25] <- 25
ried <- dB(riex, invert=TRUE)

# optimize tapers
rtap <- riedsid(riex, ntaper=ntap) # deprecation warning
rtap2 <- riedsid2(riex, ntaper=ntap)
rtap3 <- riedsid2(riex, ntaper=ntap, fast=TRUE)

# plot
op <- par(no.readonly = TRUE)
par(mfrow=c(2,1), mar=rep(1.3,4), mai=rep(0.6,4))
# ... the mock spectrum
plot(riex, type="h", xaxs="i", ylim=c(0,200), 
     main='Pseudo-spectrum') 

# ... tapers
plot(rtap2, col=NA, xaxs="i",
     main='Original and Optimized tapers', 
     ylim=c(0,max(c(ntap, rtap,rtap2,rtap3)))) 
# original tapers:
abline(h=ntap, lty=2)
# optimized tapers
lines(rtap, col="red")
# 2 and 2-fast
lines(rtap2, lwd=3, col="blue")
lines(rtap3, col="cyan")
par(op)


## End(Not run)#REX

abarbour/psd documentation built on Aug. 15, 2023, 8:56 a.m.