Discrepancy in (number of) frequencies between psd::pspectrum vs stats::spectrum (and others)

[gmgeorg]

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)