Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Improvement: calculate the cross-spectrum of two timeseries #6

Closed
sbcalaff opened this issue Feb 21, 2020 · 6 comments · Fixed by #5
Closed

Improvement: calculate the cross-spectrum of two timeseries #6

sbcalaff opened this issue Feb 21, 2020 · 6 comments · Fixed by #5

Comments

@sbcalaff
Copy link

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

@abarbour
Copy link
Owner

abarbour commented Feb 26, 2020

Yes it is. Recently (9ec9d84) @jkennel expanded the code to do just this while maintaining computational efficiently. I haven't yet had the time to issue a formal release to CRAN, but will do so soon(ish). In the meantime I encourage you to test out what has been added with the development version:

library(remotes)
install_github('abarbour/psd')

Here's a silly example:

library(psd)
set.seed(1113)
n <- 1000
inp <- rnorm(n) # represents the 'input' to the system
resp <- inp + runif(n) # represents the 'response' of the system
io <- cbind(inp, resp)

p <- pspectrum(io)

# Because io is a multi-dimensional array, 'p' should contain coherence, admittance and phase 

layout(matrix(1:3,1))
with(p,{
  plot(freq, coh, ylim=c(0,1), type='l', main='Coherence')
  plot(freq, Mod(transfer), ylim=c(0,1), type='l', main='Admittance (Gain)')
  plot(freq, phase, ylim=pi*c(-1,1), type='l', main='Phase')
})

@abarbour abarbour linked a pull request Feb 26, 2020 that will close this issue
@gmgeorg
Copy link

gmgeorg commented Jun 15, 2020

Are you planning on releasing this version to CRAN anytime soon?

background: I am the author of the ForeCA package, which relies heavily on multivariate (cross) spectral density estimation. I have relied previously on the sapa R package, but unfortunately ifultools package -- and thus sapa -- have been removed from CRAN recently. While looking for a good alternative for the excellent sapa::SDF() I came across psd and this issue in particular. I tried the test code above with the dev version here (1.2.1) and it does exactly what I need for ForeCA (return a 3D complex-valued array with the multivariate cross-spectrum).

@abarbour
Copy link
Owner

abarbour commented Jun 15, 2020 via email

@gmgeorg
Copy link

gmgeorg commented Jun 18, 2020

Apparently my grace period ended earlier, so it has been removed already last week (despite ifultools having fixed the issues; its just that the maintainer handover on their side was not handled in time I think, hence CRAN removed it).

the current ForeCA version on github works w/o SDF by just using astsa package; but it seems that psd would be a closer/better replacement of sapa::SDF so I was checking in if I should release the ForeCA update now, or just wait if psd is getting on CRAN anytime soon.

@abarbour
Copy link
Owner

abarbour commented Jun 18, 2020 via email

@abarbour
Copy link
Owner

Version 2.0 is on CRAN now (https://cran.r-project.org/package=psd). Note that the capabilities are all there, but in the interest of expedience (e.g., to meet the 'kitagawa' package deadline), we released it with with a few known issues that need to be taken care of. These are all minor ones related to methods like plot, as.data.frame, etc, as well as some documentation. So look out for minor version bumps soon. Regardless, comments and issues welcome!

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Projects
None yet
Development

Successfully merging a pull request may close this issue.

3 participants