[Rd] Proposed speedup of spec.pgram from spectrum.R

[Jonathan Kennel]

Hello,

I propose two small changes to spec.pgram to get modest speedup when
dealing with input (x) having multiple columns.  With plot = FALSE, I
commonly see ~10-20% speedup, for a two column input matrix and the speedup
increases for more columns with a maximum close to 45%.  In the function as
it currently exists, only the upper right triangle of pgram is necessary
and pgram is not returned by the function, therefore, calculating the lower
left portion (ie when j < i) is not required.  We change two nested loops
to index from i:ncol(x) instead of 1L:ncol(x) :

newspec.pgram <-
    function (x, spans = NULL, kernel = NULL, taper = 0.1,
              pad = 0, fast = TRUE,
              demean = FALSE, detrend = TRUE,
              plot = TRUE, na.action = na.fail, ...)
{
    ## Estimate spectral density from (smoothed) periodogram.
    series <- deparse(substitute(x))
    x <- na.action(as.ts(x))
    xfreq <- frequency(x)
    x <- as.matrix(x)
    N <- N0 <- nrow(x)
    nser <- ncol(x)
    if(!is.null(spans)) # allow user to mistake order of args
        kernel <- {
            if(is.tskernel(spans)) spans else
            kernel("modified.daniell", spans %/% 2)
        }
    if(!is.null(kernel) && !is.tskernel(kernel))
        stop("must specify 'spans' or a valid kernel")
    if (detrend) {
        t <- 1L:N - (N + 1)/2
        sumt2 <- N * (N^2 - 1)/12
        for (i in 1L:ncol(x))
            x[, i] <- x[, i] - mean(x[, i]) - sum(x[, i] * t) * t/sumt2
    }
    else if (demean) {
x <- sweep(x, 2, colMeans(x), check.margin=FALSE)
    }
    ## apply taper:
    x <- spec.taper(x, taper)
    ## to correct for tapering: Bloomfield (1976, p. 194)
    ## Total taper is taper*2
    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)))
    N <- nrow(x)
    Nspec <- floor(N/2)
    freq <- seq.int(from = xfreq/N, by = xfreq/N, length.out = Nspec)
    xfft <- mvfft(x)
    pgram <- array(NA, dim = c(N, ncol(x), ncol(x)))
    for (i in 1L:ncol(x)) {
        for (j in i:ncol(x)) { # N0 = #{non-0-padded}
            pgram[, i, j] <- xfft[, i] * Conj(xfft[, j])/(N0*xfreq)
            ## value at zero is invalid as mean has been removed, so
interpolate:
            pgram[1, i, j] <- 0.5*(pgram[2, i, j] + pgram[N, i, j])
        }
    }
    if(!is.null(kernel)) {
for (i in 1L:ncol(x)) for (j in i:ncol(x))
    pgram[, i, j] <- kernapply(pgram[, i, j], kernel, circular = TRUE)
df <- df.kernel(kernel)
bandwidth <- bandwidth.kernel(kernel)
    } else { # raw periodogram
df <- 2
bandwidth <- sqrt(1/12)
    }
    df <- df/(u4/u2^2)
    df <- df  * (N0 / N) ## << since R 1.9.0
    bandwidth <- bandwidth * xfreq/N
    pgram <- pgram[2:(Nspec+1),,, drop=FALSE]
    spec <- matrix(NA, nrow = Nspec, ncol = nser)
    for (i in 1L:nser) spec[, i] <- Re(pgram[1L:Nspec, i, i])
    if (nser == 1) {
        coh <- phase <- NULL
    } else {
        coh <- phase <- matrix(NA, nrow = Nspec, ncol = nser * (nser - 1)/2)
        for (i in 1L:(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])
            }
        }
    }
    ## correct for tapering
    for (i in 1L: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 = N, orig.n = N0,# "n.orig" =
"n..."
             series = series, snames = colnames(x),
             method = ifelse(!is.null(kernel), "Smoothed Periodogram",
                             "Raw Periodogram"),
             taper = taper, pad = pad, detrend = detrend, demean = demean)
    class(spg.out) <- "spec"
    if(plot) {
plot(spg.out, ...)
        return(invisible(spg.out))
    } else return(spg.out)
}



If for some reason the entire pgram array is desired (for future changes to
spec.pgram) we can minimize calls to kernapply by calculating the lower
left portion as the complex conjugate:

for (i in 1L:ncol(x)) {
    for (j in 1L:ncol(x)) {
        if(i <= j) {
            pgram[, i, j] <- kernapply(pgram[, i, j], kernel, circular =
TRUE)
        } else {
            pgram[, i, j] <- Conj(pgram[, j, i])
        }
    }
}


Cheers,
Jonathan

	[[alternative HTML version deleted]]