R/psdLS.R

In krehfeld/nest: Non-equally spaced time series analysis

#' Compute spectrum via the Lomb-Scargle periodogram
#' 
#' Compute a Lomb-Scargle periodogram
#' 
#' 
#' @usage psdLS(x, y = NULL, spans = NULL, kernel = NULL, taper = 0.1, pad = 0,
#' fast = TRUE, type = "lomb",demean = FALSE, detrend = TRUE, na.action =
#' na.fail, n.outfreq = floor(nrow(x)/2), ...)
#' @param x %% ~~Describe \code{x} here~~
#' @param y %% ~~Describe \code{y} here~~
#' @param spans %% ~~Describe \code{spans} here~~
#' @param kernel %% ~~Describe \code{kernel} here~~
#' @param taper %% ~~Describe \code{taper} here~~
#' @param pad %% ~~Describe \code{pad} here~~
#' @param fast %% ~~Describe \code{fast} here~~
#' @param type %% ~~Describe \code{type} here~~
#' @param demean %% ~~Describe \code{demean} here~~
#' @param detrend %% ~~Describe \code{detrend} here~~
#' @param na.action %% ~~Describe \code{na.action} here~~
#' @param n.outfreq %% ~~Describe \code{n.outfreq} here~~
#' @param \dots %% ~~Describe \code{\dots} here~~
psdLS <-
function (x, y = NULL, spans = NULL, kernel = NULL, taper = 0.1, 
    pad = 0, fast = TRUE, type = "lomb", demean = FALSE, detrend = TRUE, na.action = na.fail,n.outfreq = floor(nrow(x)/2), ...) 
{# Lomb-Scargle-Fourier-Transform adapted from the cts package
    if (NCOL(x) == 2) {
        series <- deparse(substitute(x))
        ti <- x[, 1]
        x <- x[, 2]
    }
    else {
        series <- deparse(substitute(y))
        ti <- x
        x <- y
    }
    x <- na.action(as.ts(x))
   # xfreq <- frequency(x)
	xfreq<-1/mean(diff(ti)) # this didn't work before

    x <- as.matrix(x)
    N <- nrow(x) #<-----------------------------
    nser <- ncol(x)
    if (!is.null(spans)) 
        if (is.tskernel(spans)) 
            kernel <- spans
        else kernel <- kernel("modified.daniell", spans%/%2)
    if (!is.null(kernel) && !is.tskernel(kernel)) 
        stop("must specify spans or a valid kernel")
    if (detrend) {
        t <- ti - mean(ti)
        sumt2 <- sum(t^2)
        for (i in 1:ncol(x)) x[, i] <- x[, i] - mean(x[, i]) - 
            sum(x[, i] * t) * t/sumt2
    }
    else if (demean) {
        x <- sweep(x, 2, colMeans(x))
    }
    x <- spec.taper(x, taper)
    u2 <- (1 - (5/8) * taper * 2)
    u4 <- (1 - (93/128) * taper * 2)
    if (pad > 0) {
        x <- rbind(x, matrix(0, nrow = N * pad, ncol = ncol(x)))
        N <- nrow(x)
    }
    NewN <- if (fast) 
        nextn(N)
    else N
    x <- rbind(x, matrix(0, nrow = (NewN - N), ncol = ncol(x))) # pad with zeroes
    N <- nrow(x)
	#Nspec <- floor(N/2) #<--------------- number of output frequencies were previously predetermined
    	Nspec <- n.outfreq # use new param
	#freq <- seq(from = xfreq/N, by = xfreq/N, length = Nspec) #<----------freq vector old
	freq <- xfreq*seq(from = 0, to = 0.5, length = Nspec+1)[-1] #<----------freq new

    #pgram <- array(NA, dim = c(N, ncol(x), ncol(x)))
    pgram <- array(NA, dim = c(Nspec, ncol(x), ncol(x))) #make sure both f and p have the same length
   # freq.temp <- seq(from = xfreq/N, by = xfreq/N, length = N)
	freq.temp<-freq

#######################################
    if (type == "lomb") {
        for (i in 1:ncol(x)) {
            for (j in 1:ncol(x)) {
                for (k in 1:length(freq.temp)) {
		#print(k)
                  tao <- atan(sum(sin(2 * pi * freq.temp[k] * 
                    ti))/sum(cos(2 * pi * freq.temp[k] * ti)))/(2 * 
                    freq.temp[k])

                  pgram[k, i, j] <- 0.5 * ((sum(x[1:length(ti)] * 
                    cos(2 * pi * freq.temp[k] * (ti - tao))))^2/sum((cos(2 * 
                    pi * freq.temp[k] * (ti - tao)))^2) + (sum(x[1:length(ti)] * 
                    sin(2 * pi * freq.temp[k] * (ti - tao))))^2/sum((sin(2 * 
                    pi * freq.temp[k] * (ti - tao)))^2))
                  #pgram[1, i, j] <- 0.5 * (pgram[2, i, j] + pgram[N, 
                  #  i, j]) # take this out of the loop
                }
		pgram[1, i, j] <- 0.5 * (pgram[2, i, j] + pgram[Nspec,i, j])
            }
        }
    }
    if (type == "ft") {
        for (i in 1:ncol(x)) {
            for (j in 1:ncol(x)) {
                for (k in 1:length(freq.temp)) pgram[k, i, j] <- ((sum(x * 
                  cos(2 * pi * freq.temp[k] * ti)))^2 + (sum(x * 
                  sin(2 * pi * freq.temp[k] * ti)))^2)/N
            }
        }
    }
    if (!is.null(kernel)) {
        for (i in 1:ncol(x)) for (j in 1:ncol(x)) pgram[, i, 
            j] <- kernapply(pgram[, i, j], kernel, circular = TRUE)
        df <- df.kernel(kernel)/(u4/u2^2)
        bandwidth <- bandwidth.kernel(kernel) * xfreq/N
    }
    else {
        df <- 2/(u4/u2^2)
        bandwidth <- sqrt(1/12) * xfreq/N
    }

    pgram <- pgram[1:Nspec, , , drop = FALSE] 
    
    spec <- matrix(NA, nrow = Nspec, ncol = nser)
    for (i in 1:nser) spec[, i] <- Re(pgram[1:Nspec, i, i])
    if (nser == 1) {
        coh <- phase <- NULL
    }
    else {
        coh <- phase <- matrix(NA, nrow = Nspec, ncol = nser * 
            (nser - 1)/2)
        for (i in 1:(nser - 1)) {
            for (j in (i + 1):nser) {
                coh[, i + (j - 1) * (j - 2)/2] <- Mod(pgram[, 
                  i, j])^2/(spec[, i] * spec[, j])
                phase[, i + (j - 1) * (j - 2)/2] <- Arg(pgram[, 
                  i, j])
            }
        }
    }
    for (i in 1:nser) spec[, i] <- spec[, i]/u2
    spec <- drop(spec)
    spg.out <- list(freq = freq, spec = spec, coh = coh, phase = phase, 
        kernel = kernel, df = df, bandwidth = bandwidth, n.used = nrow(x), 
        series = series, snames = colnames(x), method = ifelse(!is.null(kernel), 
            "Smoothed Lomb-Scargle Periodogram", "Raw Lomb-Scargle Periodogram"), 
        taper = taper, pad = pad, detrend = detrend, demean = demean)
    class(spg.out) <- "spec"

return(spg.out)
}