Merge pull request #5 from jkennel/master

Multivariate psd methods thanks to @jkennel

.travis.yml

@@ -19,6 +19,9 @@ r_packages:
  - fftw
  - covr
 
+r_github_packages:
+   - gagolews/stringi
+
 notifications:
   email:
     on_success: change

DESCRIPTION

@@ -34,7 +34,7 @@ Suggests:
     RSEIS,
     rbenchmark,
     reshape2,
-    testthat
+    testthat (>= 2.1.0)
 Encoding: UTF-8
 VignetteBuilder: knitr
 LinkingTo: Rcpp, RcppArmadillo

NAMESPACE

@@ -21,6 +21,7 @@ S3method(parabolic_weights,default)
 S3method(parabolic_weights,tapers)
 S3method(pgram_compare,amt)
 S3method(pilot_spec,default)
+S3method(pilot_spec,matrix)
 S3method(pilot_spec,ts)
 S3method(plot,tapers)
 S3method(points,tapers)
@@ -29,8 +30,10 @@ S3method(prewhiten,ts)
 S3method(print,summary.tapers)
 S3method(print,tapers)
 S3method(psdcore,default)
+S3method(psdcore,matrix)
 S3method(psdcore,ts)
 S3method(pspectrum,default)
+S3method(pspectrum,matrix)
 S3method(pspectrum,spec)
 S3method(pspectrum,ts)
 S3method(riedsid,default)
@@ -52,13 +55,15 @@ S3method(varddiff,spec)
 export(.spec_confint)
 export(adapt_message)
 export(as.tapers)
+export(coherence)
 export(colvec)
 export(constrain_tapers)
 export(create_poly)
 export(ctap_loess)
 export(ctap_simple)
 export(dB)
 export(data.frame.tapers)
+export(det_vector)
 export(get_adapt_history)
 export(get_psd_env_name)
 export(get_psd_env_pointer)
@@ -74,9 +79,11 @@ export(new_adapt_history)
 export(normalize)
 export(ones)
 export(parabolic_weights)
+export(parabolic_weights_field)
 export(parabolic_weights_rcpp)
 export(parabolic_weights_rcpp2)
 export(pgram_compare)
+export(phase)
 export(pilot_spec)
 export(prewhiten)
 export(psd_envAssign)
@@ -92,8 +99,10 @@ export(pspectrum_basic)
 export(rcpp_ctap_simple)
 export(resample_fft_rcpp)
 export(resample_fft_rcpp2)
+export(resample_mvfft)
 export(riedsid)
 export(riedsid2)
+export(riedsid_rcpp)
 export(rowvec)
 export(spec_confint)
 export(spec_details)

R/RcppExports.R

@@ -125,7 +125,74 @@ resample_fft_rcpp2 <- function(fftz, tapers, verbose = TRUE, dbl = TRUE, tapcap
 #' @param ntaper scalar or vector; number of tapers to apply optimization
 #' 
 #' @return kopt vector
+#' @export
+#' 
 riedsid_rcpp <- function(PSD, ntaper) {
     .Call('_psd_riedsid_rcpp', PACKAGE = 'psd', PSD, ntaper)
 }
 
+#' @title parabolic_weights_field
+#' @rdname parabolic_weights
+#' 
+#' @param ntap the maximum number of tapers
+#' 
+#' @export
+#' 
+parabolic_weights_field <- function(ntap) {
+    .Call('_psd_parabolic_weights_field', PACKAGE = 'psd', ntap)
+}
+
+#' @title Resample an fft using varying numbers of sine tapers
+#' 
+#' @description
+#' Produce an un-normalized psd based on an fft and a vector of optimal sine 
+#' tapers.
+#' 
+#' @details
+#' To produce a psd estimate with our adaptive spectrum estimation method,
+#' we need only make one fft calculation initially and then apply the weighting
+#' factors given by \code{\link{parabolic_weights}}, which this function 
+#' does.
+#' 
+#' @param fftz complex; a matrix representing the dual-length \code{\link{fft}}; see also the \code{dbl} argument
+#' @param tapers integer; a vector of tapers
+#' @param verbose logical; should messages be given?
+#' @param dbl logical; should the code assume \code{fftz} is dual-length or single-length?
+#' @param tapcap integer; the maximum number of tapers which can be applied; note that the length is
+#' automatically limited by the length of the series.
+#' 
+#' @return list that includes the auto and cross-spectral density, and the
+#' number of tapers
+#' 
+#' @seealso \code{\link{riedsid}}
+#' 
+#' @examples
+#' fftz <- complex(real=1:8, imaginary = 1:8)
+#' taps <- 1:4
+#' try(resample_mvfft(fftz, taps))
+#' 
+#' @export
+resample_mvfft <- function(fftz, tapers, verbose = TRUE, dbl = TRUE, tapcap = 10000L) {
+    .Call('_psd_resample_mvfft', PACKAGE = 'psd', fftz, tapers, verbose, dbl, tapcap)
+}
+
+#' @title
+#' det_vector
+#'
+#' @description
+#' Determinant for an array
+#'
+#' @param x \code{numeric array} values to evaluate
+#'
+#' @return vector of determinants
+#'
+#'
+#' @export
+det_vector <- function(x) {
+    .Call('_psd_det_vector', PACKAGE = 'psd', x)
+}
+
+solve_tf <- function(x) {
+    .Call('_psd_solve_tf', PACKAGE = 'psd', x)
+}
+

R/func_coherence_phase.R

@@ -0,0 +1,64 @@
+#' coherence
+#'
+#' calculate coherency from the spectra and cross-spectra
+#'
+#' @param pgram \code{numeric array} must be multivariate
+#'
+#'
+#' @return list of coherency
+#'
+#' @export
+#'
+#'
+coherence <- function(pgram) {
+
+    dims = dim(pgram)
+    n_ser <- dims[2]
+    n_r <- dims[1]
+
+    if (n_ser == 1) {
+      warning('Cannot calculate coherency for a single periodogram')
+      return(NA)
+    }
+
+    coh <- matrix(NA_real_, nrow = n_r, ncol = (n_ser - 1))
+
+    for (i in 2L:(n_ser)) {
+      coh[, i-1]   <- Re(Mod(pgram[, 1, i])^2 / (pgram[, 1, 1] * pgram[, i, i]))
+    }
+
+    return(coh)
+  }
+
+#' phase
+#'
+#' calculate phase from the spectra and cross-spectra
+#'
+#' @param pgram \code{numeric array} must be multivariate
+#'
+#'
+#' @return list of phase
+#'
+#' @export
+#'
+phase <- function(pgram) {
+
+  dims = dim(pgram)
+  n_ser <- dims[2]
+  n_r <- dims[1]
+
+  if (n_ser == 1) {
+    warning('Cannot calculate phase for a single periodogram')
+    return(NA)
+  }
+
+  phase <- matrix(NA_real_, nrow = n_r, ncol = (n_ser - 1))
+
+  for (i in 2L:(n_ser)) {
+
+    phase[, i-1] <- Arg(pgram[, 1, i])
+
+  }
+
+  return(phase)
+}

R/func_pilot.R

@@ -62,10 +62,19 @@ pilot_spec <- function(x, ...) UseMethod("pilot_spec")
 #' @export
 pilot_spec.ts <- function(x, ...){
   stopifnot(is.ts(x))
+  pilot_spec.default(x, ...)  
+}
+
+
+#' @rdname pilot_spec
+#' @aliases pilot_spec.matrix
+#' @export
+pilot_spec.matrix <- function(x, ...){
   frq <- stats::frequency(x)
-  pilot_spec.default(as.vector(x), x.frequency=frq, ...)  
+  pilot_spec(stats::ts(x, frequency=frq), ...)
 }
 
+
 #' @rdname pilot_spec
 #' @export
 pilot_spec.default <- function(x, x.frequency=NULL, ntap=NULL, remove.AR=NULL, plot=FALSE, verbose=FALSE, fast = FALSE, ...){
@@ -84,7 +93,6 @@ pilot_spec.default <- function(x, x.frequency=NULL, ntap=NULL, remove.AR=NULL, p
   if (REMAR) remove.AR <- max(1, min(100, abs(remove.AR)))
 
   xprew <- prewhiten(x, x.fsamp=x.frequency, AR.max=remove.AR, detrend=TRUE, impute=TRUE, plot=FALSE, verbose=verbose)
-
   ## Remove and AR model
   if (REMAR){
     # AR fit