1# File src/library/stats/R/spectrum.R 2# Part of the R package, https://www.R-project.org 3# 4# Copyright (C) 1999-2020 The R Core Team 5# Copyright (C) 1994-9 W. N. Venables and B. D. Ripley 6# 7# This program is free software; you can redistribute it and/or modify 8# it under the terms of the GNU General Public License as published by 9# the Free Software Foundation; either version 2 of the License, or 10# (at your option) any later version. 11# 12# This program is distributed in the hope that it will be useful, 13# but WITHOUT ANY WARRANTY; without even the implied warranty of 14# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 15# GNU General Public License for more details. 16# 17# A copy of the GNU General Public License is available at 18# https://www.R-project.org/Licenses/ 19 20## based on code by Martyn Plummer, plus kernel code by Adrian Trapletti 21spectrum <- function (x, ..., method = c("pgram", "ar")) 22{ 23 switch(match.arg(method), 24 pgram = spec.pgram(x, ...), 25 ar = spec.ar(x, ...) 26 ) 27} 28 29## spec.taper based on code by Kurt Hornik 30spec.taper <- function (x, p = 0.1) 31{ 32 if (any(p < 0) || any(p > 0.5)) 33 stop("'p' must be between 0 and 0.5") 34 a <- attributes(x) 35 x <- as.matrix(x) 36 nc <- ncol(x) 37 if (length(p) == 1L) 38 p <- rep(p, nc) 39 else if (length(p) != nc) 40 stop("length of 'p' must be 1 or equal the number of columns of 'x'") 41 nr <- nrow(x) 42 for (i in 1L:nc) { 43 m <- floor(nr * p[i]) 44 if(m == 0) next 45 w <- 0.5 * (1 - cos(pi * seq.int(1, 2 * m - 1, by = 2)/(2 * m))) 46 x[, i] <- c(w, rep_len(1, nr - 2 * m), rev(w)) * x[, i] 47 } 48 attributes(x) <- a 49 x 50} 51 52spec.ar <- function(x, n.freq, order = NULL, plot = TRUE, 53 na.action = na.fail, method = "yule-walker", ...) 54{ 55 ## can be called with a ts or a result of an AR fit. 56 if(!is.list(x)) { 57 series <- deparse1(substitute(x)) 58 x <- na.action(as.ts(x)) 59 xfreq <- frequency(x) 60 nser <- NCOL(x) 61 x <- ar(x, is.null(order), order, na.action=na.action, method=method) 62 } else { ## result of ar() 63 cn <- match(c("ar", "var.pred", "order"), names(x)) 64 if(anyNA(cn)) 65 stop("'x' must be a time series or an ar() fit") 66 series <- x$series 67 xfreq <- x$frequency 68 if(is.array(x$ar)) nser <- dim(x$ar)[2L] else nser <- 1 69 } 70 order <- x$order 71 if(missing(n.freq)) n.freq <- 500 72 freq <- seq.int(0, 0.5, length.out = n.freq) 73 if (nser == 1) { 74 coh <- phase <- NULL 75 var.p <- as.vector(x$var.pred) 76 spec <- 77 if(order >= 1) { 78 cs <- outer(freq, 1L:order, function(x, y) cos(2*pi*x*y)) %*% x$ar 79 sn <- outer(freq, 1L:order, function(x, y) sin(2*pi*x*y)) %*% x$ar 80 var.p/(xfreq*((1 - cs)^2 + sn^2)) 81 } else rep.int(var.p/xfreq, length(freq)) 82 } else .NotYetImplemented() 83 spg.out <- list(freq = freq*xfreq, spec = spec, coh = coh, phase = phase, 84 n.used = nrow(x), series = series, 85 method = paste0("AR (", order, ") spectrum ") 86 ) 87 class(spg.out) <- "spec" 88 if(plot) { 89 plot(spg.out, ci = 0, ...) 90 invisible(spg.out) 91 } else spg.out 92} 93 94spec.pgram <- 95 function (x, spans = NULL, kernel = NULL, taper = 0.1, 96 pad = 0, fast = TRUE, 97 demean = FALSE, detrend = TRUE, 98 plot = TRUE, na.action = na.fail, ...) 99{ 100 ## Estimate spectral density from (smoothed) periodogram. 101 series <- deparse1(substitute(x)) 102 x <- na.action(as.ts(x)) 103 xfreq <- frequency(x) 104 x <- as.matrix(x) 105 N <- N0 <- nrow(x) 106 nser <- ncol(x) 107 if(!is.null(spans)) # allow user to mistake order of args 108 kernel <- { 109 if(is.tskernel(spans)) spans else 110 kernel("modified.daniell", spans %/% 2) 111 } 112 if(!is.null(kernel) && !is.tskernel(kernel)) 113 stop("must specify 'spans' or a valid kernel") 114 if (detrend) { 115 t <- 1L:N - (N + 1)/2 116 sumt2 <- N * (N^2 - 1)/12 117 for (i in 1L:ncol(x)) 118 x[, i] <- x[, i] - mean(x[, i]) - sum(x[, i] * t) * t/sumt2 119 } 120 else if (demean) { 121 x <- sweep(x, 2, colMeans(x), check.margin=FALSE) 122 } 123 ## apply taper: 124 x <- spec.taper(x, taper) 125 ## to correct for tapering: Bloomfield (1976, p. 194) 126 ## Total taper is taper*2 127 u2 <- (1 - (5/8)*taper*2) 128 u4 <- (1 - (93/128)*taper*2) 129 if (pad > 0) { 130 x <- rbind(x, matrix(0, nrow = N * pad, ncol = ncol(x))) 131 N <- nrow(x) 132 } 133 NewN <- if(fast) nextn(N) else N 134 x <- rbind(x, matrix(0, nrow = (NewN - N), ncol = ncol(x))) 135 N <- nrow(x) 136 Nspec <- floor(N/2) 137 freq <- seq.int(from = xfreq/N, by = xfreq/N, length.out = Nspec) 138 xfft <- mvfft(x) 139 pgram <- array(NA, dim = c(N, ncol(x), ncol(x))) 140 for (i in 1L:ncol(x)) { 141 for (j in 1L:ncol(x)) { # N0 = #{non-0-padded} 142 pgram[, i, j] <- xfft[, i] * Conj(xfft[, j])/(N0*xfreq) 143 ## value at zero is invalid as mean has been removed, so interpolate: 144 pgram[1, i, j] <- 0.5*(pgram[2, i, j] + pgram[N, i, j]) 145 } 146 } 147 if(!is.null(kernel)) { 148 for (i in 1L:ncol(x)) for (j in 1L:ncol(x)) 149 pgram[, i, j] <- kernapply(pgram[, i, j], kernel, circular = TRUE) 150 df <- df.kernel(kernel) 151 bandwidth <- bandwidth.kernel(kernel) 152 } else { # raw periodogram 153 df <- 2 154 bandwidth <- sqrt(1/12) 155 } 156 df <- df/(u4/u2^2) 157 df <- df * (N0 / N) ## << since R 1.9.0 158 bandwidth <- bandwidth * xfreq/N 159 pgram <- pgram[2:(Nspec+1),,, drop=FALSE] 160 spec <- matrix(NA, nrow = Nspec, ncol = nser) 161 for (i in 1L:nser) spec[, i] <- Re(pgram[1L:Nspec, i, i]) 162 if (nser == 1) { 163 coh <- phase <- NULL 164 } else { 165 coh <- phase <- matrix(NA, nrow = Nspec, ncol = nser * (nser - 1)/2) 166 for (i in 1L:(nser - 1)) { 167 for (j in (i + 1):nser) { 168 coh[, i + (j - 1) * (j - 2)/2] <- 169 Mod(pgram[, i, j])^2/(spec[, i] * spec[, j]) 170 phase[, i + (j - 1) * (j - 2)/2] <- Arg(pgram[, i, j]) 171 } 172 } 173 } 174 ## correct for tapering 175 for (i in 1L:nser) spec[, i] <- spec[, i]/u2 176 spec <- drop(spec) 177 spg.out <- 178 list(freq = freq, spec = spec, coh = coh, phase = phase, 179 kernel = kernel, df = df, 180 bandwidth = bandwidth, n.used = N, orig.n = N0,# "n.orig" = "n..." 181 series = series, snames = colnames(x), 182 method = if(is.null(kernel)) "Raw Periodogram" else "Smoothed Periodogram", 183 taper = taper, pad = pad, detrend = detrend, demean = demean) 184 class(spg.out) <- "spec" 185 if(plot) { 186 plot(spg.out, ...) 187 invisible(spg.out) 188 } else 189 spg.out 190} 191 192plot.spec <- 193 function (x, add = FALSE, ci = 0.95, log = c("yes", "dB", "no"), 194 xlab = "frequency", ylab = NULL, 195 type = "l", ci.col = "blue", ci.lty = 3, 196 main = NULL, sub = NULL, 197 plot.type = c("marginal", "coherency", "phase"), ...) 198{ 199 spec.ci <- function (spec.obj, coverage = 0.95) 200 { 201 ## A utility function for plot.spec which calculates the confidence 202 ## interval (centred around zero). We use a conditional argument to 203 ## ensure that the ci always contains zero. 204 205 if (coverage < 0 || coverage >= 1) 206 stop("coverage probability out of range [0,1)") 207 tail <- (1 - coverage) 208 df <- spec.obj$df 209 upper.quantile <- 1 - tail * pchisq(df, df, lower.tail = FALSE) 210 lower.quantile <- tail * pchisq(df, df) 211 1/(qchisq(c(upper.quantile, lower.quantile), df)/df) 212 } 213 214 plot.type <- match.arg(plot.type) 215 log <- match.arg(log) 216 m <- match.call() 217 if(plot.type == "coherency") { 218 ## need stats:: for non-standard evaluation 219 m[[1L]] <- quote(stats::plot.spec.coherency) 220 m$plot.type <- m$log <- m$add <- NULL 221 return(eval(m, parent.frame())) 222 } 223 if(plot.type == "phase") { 224 ## need stats:: for non-standard evaluation 225 m[[1L]] <- quote(stats::plot.spec.phase) 226 m$plot.type <- m$log <- m$add <- NULL 227 return(eval(m, parent.frame())) 228 } 229 if(is.null(ylab)) 230 ylab <- if(log == "dB") "spectrum (dB)" else "spectrum" 231 if(is.logical(log)) 232 log <- if(log) "yes" else "no" 233 if(missing(log) && getOption("ts.S.compat")) log <- "dB" 234 log <- match.arg(log) 235 ylog <- "" 236 if(log=="dB") x$spec <- 10 * log10(x$spec) 237 if(log=="yes") ylog <- "y" 238 dev.hold(); on.exit(dev.flush()) 239 if(add) { 240 matplot(x$freq, x$spec, type = type, add=TRUE, ...) 241 } else { 242 matplot(x$freq, x$spec, xlab = xlab, ylab = ylab, type = type, 243 log = ylog, ...) 244 if (ci <= 0 || !is.numeric(x$df) || log == "no") { 245 ## No confidence limits 246 ci.text <- "" 247 } else { 248 ## The position of the error bar has no meaning: only the width 249 ## and height. It is positioned in the top right hand corner. 250 ## 251 conf.lim <- spec.ci(x, coverage = ci) 252 if(log=="dB") { 253 conf.lim <- 10*log10(conf.lim) 254 conf.y <- max(x$spec) - conf.lim[2L] 255 conf.x <- max(x$freq) - x$bandwidth 256 lines(rep(conf.x, 2), conf.y + conf.lim, col=ci.col) 257 lines(conf.x + c(-0.5, 0.5) * x$bandwidth, rep(conf.y, 2), 258 col=ci.col) 259 ci.text <- paste0(", ", round(100*ci, 2), "% C.I. is (", 260 paste(format(conf.lim, digits = 3), 261 collapse = ","), 262 ")dB") 263 } else { 264 ci.text <- "" 265 conf.y <- max(x$spec) / conf.lim[2L] 266 conf.x <- max(x$freq) - x$bandwidth 267 lines(rep(conf.x, 2), conf.y * conf.lim, col=ci.col) 268 lines(conf.x + c(-0.5, 0.5) * x$bandwidth, rep(conf.y, 2), 269 col=ci.col) 270 } 271 } 272 if (is.null(main)) 273 main <- paste(if(!is.null(x$series)) paste("Series:", x$series) 274 else "from specified model", 275 x$method, sep = "\n") 276 if (is.null(sub) && is.numeric(x$bandwidth)) 277 sub <- paste0("bandwidth = ", format(x$bandwidth, digits = 3), 278 ci.text) 279 title(main = main, sub = sub) 280 } 281 invisible(x) 282} 283 284## based on code in Venables & Ripley 285plot.spec.coherency <- 286 function(x, ci = 0.95, 287 xlab = "frequency", ylab = "squared coherency", ylim=c(0,1), 288 type = "l", main = NULL, ci.col="blue", ci.lty = 3, ...) 289{ 290 nser <- NCOL(x$spec) 291 ## Formulae from Bloomfield (1976, p.225) 292 gg <- 2/x$df 293 se <- sqrt(gg/2) 294 z <- -qnorm((1-ci)/2) 295 if (is.null(main)) 296 main <- paste(paste("Series:", x$series), 297 "Squared Coherency", sep = " -- ") 298 if(nser == 2) { 299 plot(x$freq, x$coh, type=type, xlab=xlab, ylab=ylab, ylim=ylim, ...) 300 coh <- pmin(0.99999, sqrt(x$coh)) 301 lines(x$freq, (tanh(atanh(coh) + z*se))^2, lty=ci.lty, col=ci.col) 302 lines(x$freq, (pmax(0, tanh(atanh(coh) - z*se)))^2, 303 lty=ci.lty, col=ci.col) 304 title(main) 305 } else { 306 dev.hold(); on.exit(dev.flush()) 307 opar <- par(mfrow = c(nser-1, nser-1), mar = c(1.5, 1.5, 0.5, 0.5), 308 oma = c(4, 4, 6, 4)) 309 on.exit(par(opar), add = TRUE) 310 plot.new() 311 for (j in 2:nser) for (i in 1L:(j-1)) { 312 par(mfg=c(j-1,i, nser-1, nser-1)) 313 ind <- i + (j - 1) * (j - 2)/2 314 plot(x$freq, x$coh[, ind], type=type, ylim=ylim, axes=FALSE, 315 xlab="", ylab="", ...) 316 coh <- pmin(0.99999, sqrt(x$coh[, ind])) 317 lines(x$freq, (tanh(atanh(coh) + z*se))^2, lty=ci.lty, col=ci.col) 318 lines(x$freq, (pmax(0, tanh(atanh(coh) - z*se)))^2, 319 lty=ci.lty, col=ci.col) 320 box() 321 if (i == 1) { 322 axis(2, xpd = NA) 323 title(ylab=x$snames[j], xpd = NA) 324 } 325 if (j == nser) { 326 axis(1, xpd = NA) 327 title(xlab=x$snames[i], xpd = NA) 328 } 329 mtext(main, 3, 3, TRUE, 0.5, 330 cex = par("cex.main"), font = par("font.main")) 331 } 332 } 333 invisible() 334} 335 336plot.spec.phase <- 337 function(x, ci = 0.95, 338 xlab = "frequency", ylab = "phase", ylim=c(-pi, pi), 339 type = "l", main = NULL, ci.col = "blue", ci.lty = 3, ...) 340{ 341 nser <- NCOL(x$spec) 342 ## Formulae from Bloomfield (1976, p.225) 343 gg <- 2/x$df 344 if (is.null(main)) 345 main <- paste(paste("Series:", x$series), 346 "Phase spectrum", sep = " -- ") 347 if(nser == 2) { 348 plot(x$freq, x$phase, type=type, xlab=xlab, ylab=ylab, ylim=ylim, ...) 349 coh <- sqrt(x$coh) 350 cl <- asin( pmin( 0.9999, qt(ci, 2/gg-2)* 351 sqrt(gg*(coh^{-2} - 1)/(2*(1-gg)) ) ) ) 352 lines(x$freq, x$phase + cl, lty=ci.lty, col=ci.col) 353 lines(x$freq, x$phase - cl, lty=ci.lty, col=ci.col) 354 title(main) 355 } else { 356 dev.hold(); on.exit(dev.flush()) 357 opar <- par(mfrow = c(nser-1, nser-1), mar = c(1.5, 1.5, 0.5, 0.5), 358 oma = c(4, 4, 6, 4)) 359 on.exit(par(opar), add = TRUE) 360 plot.new() 361 for (j in 2:nser) for (i in 1L:(j-1)) { 362 par(mfg=c(j-1,i, nser-1, nser-1)) 363 ind <- i + (j - 1) * (j - 2)/2 364 plot(x$freq, x$phase[, ind], type=type, ylim=ylim, axes=FALSE, 365 xlab="", ylab="", ...) 366 coh <- sqrt(x$coh[, ind]) 367 cl <- asin( pmin( 0.9999, qt(ci, 2/gg-2)* 368 sqrt(gg*(coh^{-2} - 1)/(2*(1-gg)) ) ) ) 369 lines(x$freq, x$phase[, ind] + cl, lty=ci.lty, col=ci.col) 370 lines(x$freq, x$phase[, ind] - cl, lty=ci.lty, col=ci.col) 371 box() 372 if (i == 1) { 373 axis(2, xpd = NA) 374 title(ylab=x$snames[j], xpd = NA) 375 } 376 if (j == nser) { 377 axis(1, xpd = NA) 378 title(xlab=x$snames[i], xpd = NA) 379 } 380 mtext(main, 3, 3, TRUE, 0.5, 381 cex = par("cex.main"), font = par("font.main")) 382 } 383 } 384 invisible() 385} 386