# The following functions are contained in acvs.R
# These functions are loosley based on the lisp code provided with SAPA
# http://lib.stat.cmu.edu/sapaclisp/
# http://lib.stat.cmu.edu/sapaclisp/acvs.lisp
# acvs <- function(timeSeries, centreDataP=T, acsP=F)
# acvs_ws <- function(timeSeries, centreDataP=T, acsP=F)
# biasedACVStoUnbiasedACVS <- function(biasedACVS)
################################################################################
#libraries required
##library(waveslim);
################################################################################
# Requires:
##source("~/PWLisp/basicStatistics.R");
##source("~/PWLisp/utilities.R");
##source("~/PWLisp/tapers.R");
#source("nonparametric.R"); caused probs circular
################################################################################
#coded to match the lisp acvs function in the SAPA code
#Note: see Equations (191a) and (190b) of the SAPA book
acvs <- function(timeSeries, centreDataP=TRUE, acsP=FALSE) {
N <- length(timeSeries);
ntot <- nextPowerOf2(2*N);
N_Xntot <- N*ntot;
if(centreDataP) {
timeSeries <- timeSeries - mean(timeSeries);
}
scratch <- c(timeSeries, array(0, ntot-N));
scratch <- fft(scratch);
scratch <- abs(scratch)^2;
scratch <- fft(scratch);
sampVar <- Re(scratch[1])/N_Xntot;
scaleFactor <- if(acsP) {N_Xntot*sampVar;} else {N_Xntot;}
result <- Re(scratch[1:N])/ scaleFactor;
return(list(acvs=result, var=sampVar));
}
#uses my.acf provided in waveslim
acvs_ws <- function(timeSeries, centreDataP=TRUE, acsP=FALSE) {
if(centreDataP) {
timeSeries <- timeSeries - mean(timeSeries);
}
#requires waveslim
res <- my.acf(timeSeries);
if(acsP) {
res <- res / res[1];
}
return(list(acvs=res, var=res[1]));
}
biasedACVStoUnbiasedACVS <- function(biasedACVS) {
result <- biasedACVS*length(biasedACVS) / (length(biasedACVS):1);
return(result);
}
## i'm not sure what this is for but acvs2 allows
## for matrix valued fft used in dpss.R
acvsB <- function(timeSeries, centreDataP=TRUE, acsP=FALSE) {
N <- length(timeSeries);
ntot <- nextPowerOf2(2*N);
N_Xntot <- N*ntot;
if(centreDataP) {
timeSeries <- timeSeries - mean(timeSeries);
}
scratch <- c(timeSeries, array(0, ntot-N));
scratch <- fft(scratch);
scratch <- abs(scratch)^2;
scratch <- fft(scratch);
sampVar <- Re(scratch[1])/N_Xntot;
scaleFactor <- if(acsP) {N_Xntot*sampVar;} else {N_Xntot;}
result <- Re(scratch[1:N])/ scaleFactor;
return(list(result=result, sampleVariance=sampVar));
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.