abarbour / psd Goto Github PK

The preprocessing flags are not reflective of the actual preprocessing in the final product, e.g., in the title of the plot.spec method.

I run into a namespace issue when using psd::pspectrum:

data(magnet, package = "psd")
x <- magnet$clean
spec <- psd::pspectrum(x, plot = TRUE)

gives error on

Error in psdcore.default(ntaper = 7, AR = FALSE, preproc = FALSE, refresh = TRUE,  : 
  !is.null(ops) is not TRUE

Same issue is described here at end of document:
https://towardsdatascience.com/attach-packages-mindfully-in-r-b9c8151b3fb4

For a 1D map can we have some functions to return a boolean value (TRUE or FALSE)
the function may be named

is.chaotic()

or

is.quasiPeriodical()

Judging by width of 'peaks' might work. A better idea may be to calculate square of errors when fitting a smoothed (e.g. linear) model to the P[dB] ~ log(frequency) plot.
Is this achievable?

Thanks a lot for releasing the multivariate support for psd. Incorporating this into ForeCA was straightforward as psd returns the 3D array power spectra. However, I am running into a subtle inconsistency in the calculated frequencies of the periodogram. psd::pspectrum frequencies differ from various other packages including the baseline stats::spectrum. I am aware that different implementations handle frequency 0 and the (a)symmetry of frequencies around 0 differently (returning only the non-negative side of the spectrum). However, it seems that psd frequencies are actually computed with a different denominator ([n/2]-1) compared to other implementations ([n+1 / 2] - 1?).

The psd spectra still are great estimates for ForeCA, it's just that if possible I'd like to avoid the discrepancy in the first place as the psd::pspectrum return values have a lower dimensional shape compared to various other spectrum estimates (sapa::SDF, astsa::mvspec and stats::spectrum all agree on the (number of) frequencies).

Is there any way to pass a flag to psectrum so it calculates according to same frequencies as seen in other implementations? The only hacky solution I can think of is to add a row of 0s to the data and compute pspectrum on that augmented dataset.

See below for a reproducible example

set.seed(10)
library(astsa)
library(psd)

nn <- 20  # also try 18 for odd number of observations; guess related to the [(n-1)/2] - 1 vs [n/2] - 1 discrepancy
X <- matrix(rnorm(nn), ncol = 2)

num_non_zero_freqs <- function(freqs) {
  print(c(length(freqs), length(freqs[freqs > 0])))
  return(freqs)
}

num_non_zero_freqs(psd::pspectrum(X)$freq)
# hacky fix to make pspectrum output consistent
num_non_zero_freqs(psd::pspectrum(rbind(0, X, 0), ntap.init=3, Nyquist.normalize=TRUE)$freq)
num_non_zero_freqs(stats::spectrum(X)$freq)
num_non_zero_freqs(astsa::mvspec(X)$freq)

The first and last values of the spectra are strangely very low. The solution thus far is to extrapolate values, which is precarious; however, those points are very nearly useless for most applications of power spectral density estimates.

See ?psdcore

Technically speaking, is it possible to calculate the optimal number of tapers at each frequency for a cross-spectrum of two univariate timeseries?