parabolic_weights: Calculate parabolic weighting factors.
In williamdeleo/psd: Adaptive, Sine-Multitaper Power Spectral Density Estimation

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

The resampled spectrum involves summing weighted tapers; this produces the weighting factors. parabolic_weights_rcpp is the fastest implementation, used by resample_fft_rcpp, but it takes only a single value.

parabolic_weights_rcpp(ntap = 1L)

parabolic_weights(ntap, ...)

## S3 method for class 'tapers'
parabolic_weights(ntap, tap.index = 1L)

## Default S3 method:
parabolic_weights(ntap = 1L)

`ntap`	integer (or `tapers` object); the number of tapers to provide weightings for.
`...`	optional arguments
`tap.index`	integer; if `ntap` is a `tapers` object, the index from which to produce a sequence of weights for

If one has a tapers object, specify the taper.index to produce a sequence of weights up to the value at that index; the user is likely to never need to use this function though.

Weighting factors, w, are calculated as follows:

w = 6 (k^2 - (K-1)^2) / (k (4 * k - 1) (k + 1))

where k is the total number of tapers, and K is the integer sequence [0,K-1].

The sum of tapers should equal 1, within machine precision, when k>0.

A list with the number of tapers, indices of the taper sequence, and the weights w.

A.J. Barbour adapted the original algorithm (R.L. Parker), and authored the optimized versions.

resample_fft_rcpp, psdcore, riedsid2

## Not run: #REX
library(psd)
library(grDevices)
library(RColorBrewer)

##
## Show parabolic weighting factors as a function of maximum tapers
##

# maximum number of tapers
maxx <- 1e3
# sequence in logspace
xseq <- seq(from=1,to=2.8,by=0.2)

# plot palette
pal <- "Spectral"
npal <- switch(pal, RdYlBu=11, Spectral=11, Blues=9)
pal.col <- RColorBrewer::brewer.pal(npal, pal)
cols <- rev(grDevices::colorRampPalette(pal.col)(maxx))

to_df <- function(W){
  # convert parabolic results to data.frame
  with(W, data.frame(taper_seq=as.vector(taper_seq), taper_weights=as.vector(taper_weights)))
}

## a roundabout way of bootstrapping y-axis limits:
#  upper
WgtsU <- parabolic_weights(5)
DfU <- to_df(WgtsU)
#  lower
WgtsL <- parabolic_weights(maxx)
DfL <- to_df(WgtsL)

ylims <- range(pretty(dB(c(DfL$taper_weights, DfU$taper_weights)))) + c(-2,5)

# function for plotting text
TFUN <- function(Df.){
  tx <- max(Df.$taper_seq)
  ty <- mean(Df.$taper_weights)
  text(log10(tx)+0.1, dB(ty), sprintf("%i", tx), col=cols[tx])
}

# function for weighting factors and plotting
WFUN <- function(x){
  message(x)
  Wgts <- parabolic_weights(x)
  Df <- to_df(Wgts)
  lcol <- cols[x]
  lines(dB(taper_weights) ~ log10(taper_seq), Df, type="s", lwd=2, col=lcol)
  TFUN(Df)
}

## Plot parabolic weighting, in dB, colored by maximum num tapers
plot(dB(taper_weights) ~ log10(taper_seq), DfU, type="s", 
     xlim=c(0, log10(maxx)+0.2), 
     ylim=ylims, yaxs="i",
     col=cols[5], lwd=2,  
     main="Multitaper weighting factors by maximum tapers applied",
     xlab="log10 taper sequence", 
     ylab="dB")
TFUN(DfU)
invisible(lapply(round(10**xseq), FUN=WFUN))
WFUN(maxx)

##

## End(Not run)#REX