parabolic_weights: parabolic_weights_field

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

parabolic_weights_field

Description

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.

Usage

parabolic_weights_rcpp(ntap = 1L)

parabolic_weights_field(ntap)

parabolic_weights(ntap, ...)

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

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

Arguments

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

Details

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 \equiv \frac{6 (n^2 - K^2)}{n (4 * n - 1) (n + 1)}

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

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

Value

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

Author(s)

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

See Also

resample_fft_rcpp, psdcore, riedsid2

Examples

## 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

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