tests/testthat/test-pwelch.R

In ravetools: Signal and Image Processing Toolbox for Analyzing Intracranial Electroencephalography Data

test_that("pwelch", {
  rave_pwelch <- function(x,
                          fs,
                          window = 64,
                          noverlap = 8,
                          nfft = 256,
                          col = "black",
                          xlim = NULL,
                          ylim = NULL,
                          main = "Welch periodogram",
                          plot = TRUE,
                          log = "xy",
                          spec_func = stats::spectrum,
                          cex = 1,
                          ...)
  {
    x <- as.vector(x)
    x_len <- length(x)
    nfft <- max(min(nfft, length(x)), window)
    window <- hanning(window)
    window_len <- length(window)
    step <- max(floor(window_len - noverlap + 0.99), 1)
    # print(c(1, x_len - window_len + 1, by = step))
    offset <- seq(1, x_len - window_len + 1, by = step)
    N <- length(offset)
    re <- sapply(seq_len(N), function(i) {
      a <- detrend_naive(x[offset[i] - 1 + seq_len(window_len)])
      a <- fftwtools::fftw_r2c(postpad(a$Y * window, nfft))
      Mod(a)^2
    })
    NN <- ceiling((nfft + 1) / 2)
    freq <- seq(0, fs / 2, length.out = NN)
    spec <- rowMeans(re[seq_len(NN), , drop = FALSE]) / window_len #(window_len/2)^2
    res <- list(freq = freq,
                spec = spec,
                method = "Welch")
    return(invisible(res))
  }


  x <- rnorm(10000)
  a <- pwelch(
    x,
    fs = 100,
    window = 64,
    noverlap = 6,
    nfft = 256,
    plot = FALSE,
    window_family = hanning
  )
  b <- rave_pwelch(
    x,
    fs = 100,
    window = 64,
    noverlap = 6,
    nfft = 256,
    log = ""
  )

  expect_equal(a$freq, b$freq)
  expect_equal(a$spec, b$spec)

  # if(FALSE) {
  #   x <- rnorm(600000)
  #   microbenchmark::microbenchmark(
  #     rave = {
  #       rave_pwelch(x, 2000, log = '')
  #     },
  #     ravetools = {
  #       pwelch(x, 2000, plot = FALSE)
  #     }, times = 20
  #   )
  # }

})