R/empi2tf.R

In MatchingPursuit: Processing Time Series Data Using the Matching Pursuit Algorithm

#' Creates a time-frequency map using atoms from the Matching Pursuit algorithm
#'
#' Creates a time-frequency map using atoms from the Matching Pursuit algorithm.
#' The created map can be: 1) displayed on the screen, 2) saved in \code{.png} file,
#' or 3) saved as an \code{.RData} object.
#'
#' @importFrom graphics rasterImage par points text axis mtext layout plot.new plot.window box abline
#' @importFrom grDevices hcl.colors graphics.off pdf dev.off png
#' @importFrom utils tail
#' @importFrom DescTools DrawEllipse
#' @importFrom imager as.cimg resize
#'
#' @param db.file The SQLite file created after executing the \code{empi.execute()} function.
#' In this case, the \code{db.list} parameter must be \code{NULL}.
#'
#' @param db.list The list created after executing the \code{empi.execute()} function.
#' In this case, the \code{db.file} parameter must be \code{NULL}.
#'
#' @param channel Channel from the SQLite file to process.
#'
#' @param mode \code{"sqrt"}, \code{"log"}, or \code{"linear"}. It determines the intensity
#' with which the so-called blobs are displayed on the T-F map.
#'
#' @param freq.divide Specifies how many times the displayed frequency in the T-F map
#' should be decreased. For example, if the sampling frequency is \code{f=256Hz}, the maximum
#' frequency in the T-F map will be \code{f/2/freq.divide} (f/2 is the Nyquist rule).
#'
#' @param increase.factor Factor of increasing the number of pixels in the f-axis, the most
#' sensible are non-negative integers (e.g. 2, 4, 5, 8).
#'
#' @param shortening.factor.x Usually, for better visualization of atoms, a value of 2 will
#' be appropriate.
#'
#' @param shortening.factor.y Usually, for better visualization of atoms, a value of 2 will
#' be appropriate.
#'
#' @param display.crosses Whether small crosses should be displayed in the canters of atoms.
#'
#' @param display.atom.numbers Whether atom numbers should be displayed in the canters of atoms.
#'
#' @param display.grid Whether to draw grid lines.
#'
#' @param crosses.color Colour of small crosses.
#'
#' @param palette Palette from the list returned by \code{hcl.pals()} function or the string
#' \code{"my custom palette"}.
#'
#' @param rev \code{rev} param in \code{hcl.colors()} function.
#'
#' @param out.mode One of the following:
#'   \itemize{
#'      \item \code{"plot"} - draws a T-F map on the screen.
#'      \item \code{"file"} - saves a T-F map to file \code{file.name} (as \code{png} file).
#'      \item \code{"RData"} - saves the T-F map of \code{size} in the \code{file.name}
#'      (as R's matrix), resampling is performed using the function \code{imager::resize()} function.
#'      \item \code{"RData2"} - saves the T-F map of \code{size} in the \code{file.name}
#'      (as R's matrix), resampling is performed using the function using \code{raster::resample()} function.
#'    }
#'
#' @param path Path where \code{png}, \code{RData} or \code{pdf} files will be written. If \code{NULL},
#' files will be written to the cache directory.
#'
#' @param file.name Name of the \code{png} file (if \code{out.mode="file}) or name of the \code{RData}
#' file (if \code{out.mode="RData"} or \code{out.mode="RData2})
#'
#' @param size \code{png} file size in pixels (if \code{out.mode="file"}) or size of the T-F matrix
#' (if \code{out.mode="RData"} of \code{out.mode="RData2"}).
#'
#' @param draw.ellipses Only for testing. User can set it to \code{TRUE} to see the effect.
#' Works properly only if \code{out.mode="plot"}.
#'
#' @param plot.signals Whether the original and reconstructed signals should also be displayed.
#'
#' @param write.atoms If \code{TRUE}, writes all atom plots into \code{Atoms.pdf} file
#' (to the cache directory or user specified one - depending on \code{path} variable)
#'
#' @return Depending on the \code{out.mode} parameter the function returns:
#'    \itemize{
#'      \item Time-Frequency map plotted on the screen
#'      \item Time-Frequency map saved in a .png file
#'      \item Time-Frequency map saved as .RData file
#'   }
#' Regardless of the above, the function returns the following:
#'   \itemize{
#'     \item all the Gabor functions
#'     \item reconstructed signal
#'     \item original signal
#'     \item sampling frequency
#'     \item grid size in t axis
#'     \item grid size in f axis
#'     \item epoch size in samples
#'     \item length of the signal in seconds
#'     \item time-frequency map
#'     \item time-frequency map after resampling
#'     (if \code{out.mode="RData"} or if \code{out.mode="RData2"}, otherwise, \code{NULL} is returned)
#'     \item channel number processed
#'     \item frequency divide
#'   }
#'
#' @export
#'
#' @examples
#' file <- system.file("extdata", "sample1.db", package = "MatchingPursuit")
#'
#' out <- empi2tf(
#'   db.file = file,
#'   channel = 1,
#'   mode = "sqrt",
#'   freq.divide = 4,
#'   increase.factor= 4,
#'   display.crosses = TRUE,
#'   display.atom.numbers = FALSE,
#'   out.mode = "plot",
#' )
#'
empi2tf <- function(
    db.file = NULL,
    db.list = NULL,
    channel,
    mode = "sqrt",
    freq.divide = 1,
    increase.factor = 1,
    shortening.factor.x = 2,
    shortening.factor.y = 2,
    display.crosses = TRUE,
    display.atom.numbers = FALSE,
    display.grid = FALSE,
    crosses.color = "white",
    palette = 'my custom palette',
    rev = TRUE,
    out.mode = "plot",
    path = NULL,
    file.name = NULL,
    size = c(512, 512),
    draw.ellipses = FALSE,
    plot.signals = TRUE,
    write.atoms = FALSE) {

  # Store
  old.par <- par("mfrow", "pty", "mai", "mgp", "las", "xaxs", "yaxs")
  on.exit(par(old.par))

  if (out.mode != "plot" && out.mode != "file" && out.mode != "RData" && out.mode != "RData2")
    stop("Incorrect value for 'out.mode' parameter.'")

  if (is.null(db.file) && is.null(db.list))
    stop("Specify input as SQLite file _OR_ a list returned by the 'empi.execute()' function.")

  if (!is.null(db.file) && !is.null(db.list))
    stop("Specify input as SQLite file _OR_ a list returned by the 'empi.execute()' function.")

  if (is.null(file.name)) {
    if  (out.mode == "file") fn <- "TFmap.png"
    if  (out.mode == "RData" || out.mode == "RData2") fn <- "TFmap.RData"
  } else {
    fn <- file.name
  }

  if (out.mode != "plot") {
    if (is.null(path)) {
      dest.dir <- tools::R_user_dir("MatchingPursuit", "cache")
      dir.create(dest.dir, recursive = TRUE, showWarnings = FALSE)
      file.name <- file.path(dest.dir, fn)
    } else {
      if (!dir.exists(path)) {
        ok <- dir.create(path, recursive = TRUE, showWarnings = FALSE)
        if (!ok && !dir.exists(path)) {
          stop("Cannot create directory '", path, "'.")
        }
      }
      file.name <- file.path(path, fn)
    }
  }

  if (write.atoms) {
    if (is.null(path)) {
      dest.dir <- tools::R_user_dir("MatchingPursuit", "cache")
      dir.create(dest.dir, recursive = TRUE, showWarnings = FALSE)
      atoms.file.name <- file.path(dest.dir, "Atoms.pdf")
    } else {
      if (!dir.exists(path)) {
        ok <- dir.create(path, recursive = TRUE, showWarnings = FALSE)
        if (!ok && !dir.exists(path)) {
          stop("Cannot create directory '", path, "'.")
        }
      }
      atoms.file.name <- file.path(path, "Atoms.pdf")
    }
  }

  if (out.mode != "RData" & out.mode != "RData2")
    tf.map.resampled <- NULL

  if (palette == 'my custom palette') {
    col <-  c(
      "#000f82", "#001385", "#011789", "#011b8d", "#021f91", "#022395", "#032798", "#042b9c",
      "#042fa0", "#0533a4", "#0537a8", "#063bab", "#073faf", "#0743b3", "#0847b7", "#084bbb",
      "#094fbf", "#0a53c2", "#0a57c6", "#0b5bca", "#0b5fce", "#0c63d2", "#0d67d5", "#0d6bd9",
      "#0e6fdd", "#0e73e1", "#0f77e5", "#107be8", "#107fec", "#1183f0", "#1187f4", "#128bf8",
      "#1390fc", "#1693f8", "#1996f4", "#1c9af0", "#1f9dec", "#22a1e8", "#25a4e4", "#28a8e0",
      "#2cabdc", "#2faed8", "#32b2d4", "#35b5d0", "#38b9cc", "#3bbcc8", "#3ec0c4", "#41c3c0",
      "#45c7bc", "#48cab8", "#4bcdb4", "#4ed1b0", "#51d4ac", "#54d8a8", "#57dba4", "#5adfa0",
      "#5ee29c", "#61e598", "#64e994", "#67ec90", "#6af08c", "#6df388", "#70f784", "#73fa80",
      "#77fe7c", "#7bf978", "#7ff575", "#83f171", "#88ed6e", "#8ce96a", "#90e567", "#94e063",
      "#99dc60", "#9dd85c", "#a1d459", "#a5d055", "#aacc52", "#aec74e", "#b2c34b", "#b6bf47",
      "#bbbb44", "#bfb741", "#c3b33d", "#c7af3a", "#ccaa36", "#d0a633", "#d4a22f", "#d89e2c",
      "#dd9a28", "#e19625", "#e59121", "#e98d1e", "#ee891a", "#f28517", "#f68113", "#fa7d10",
      "#ff790d", "#fb750c", "#f8720c", "#f56f0b", "#f26c0b", "#ef690a", "#eb660a", "#e8630a",
      "#e56009", "#e25c09", "#df5908", "#db5608", "#d85308", "#d55007", "#d24d07", "#cf4a06",
      "#cc4706", "#c84306", "#c54005", "#c23d05", "#bf3a04", "#bc3704", "#b83404", "#b53103",
      "#b22e03", "#af2a02", "#ac2702", "#a82402", "#a52101", "#a21e01", "#9f1b00", "#9c1800",
      "#991500")
  } else {
    col <- hcl.colors(128, palette, rev = rev)
  }

  if (!is.null(db.file)) {
    out <- read.empi.db.file(db.file)
  }

  if (!is.null(db.list)) {
    out <- db.list
  }

  if (length(which(out$atoms$channel_id == channel)) == 0)
    stop("There is no channel number ", channel, ".")

  f <- out$f

  # according to the Nyquist criterion
  maxf <- round((f / 2) / freq.divide)

  epochSize <- length(out$t)
  s <- epochSize / f

  # grid size in t
  t <- seq(from = 0, to = s, length.out = epochSize)

  # grid size in f
  y <- seq(from = 0, to = maxf, length.out = maxf * increase.factor)

  # t-f map
  tf.map <- matrix(0, nrow = epochSize, ncol = maxf * increase.factor)

  grid <- expand.grid(x = t, y = y)

  # number of atoms
  rows <- which(out$atoms$channel_id == channel)
  num.atoms <- length(rows)

  # atoms params
  energy <- out$atoms$energy[rows]
  scale <- out$atoms$scale[rows]
  position <- out$atoms$position
  frequency <- out$atoms$frequency[rows]
  original.signal <- out$original.signal[, channel]
  reconstruction <- out$reconstruction[, channel]
  gabors <- out$gabors[[channel]]

  # Signal energy
  o <- round(sum(original.signal^2), 2)
  r <- round(sum(reconstruction^2), 2)

  message("Channel #: ", channel)
  message("Number of atoms: ", length(rows))
  message("Sampling rate: ", f)
  message("Epoch size (in points): ", epochSize)
  message("Signal length (in seconds): ", s)

  message("\nEnergy of the original signal:       ",o)
  message("nEnergy of the reconstructed signal: ",r)
  message("reconstruction / original %:         ", round(r / o * 100, digits = 2), "\n")

  if(write.atoms) {
    graphics.off()
    pdf(atoms.file.name, width = 15, height = 30)
    message("Atom plots saved in '", atoms.file.name, "'")

    nn <- num.atoms
    # mai: c(bottom, left, top, right)
    par(mfrow = c(nn, 1), pty = "m", mai = c(0.05, 4, 0.0, 0.1), mgp = c(0, 0, 0), las = 1)

    for (m in 1:num.atoms) {
      if (m %% 2 == 0) cc = "blue" else cc = "red"
      if (m == 1) {
        plot(gabors[, m], xlab = "", ylab = "", xaxt = "n", yaxt = "n", type = "l", bty = "n", col = cc)
        lab <- seq(from = 0, to = ceiling(tail(t, 1)), length.out = 11)
        axis(
          side = 1, las = 1, cex.axis = 0.9,
          at = seq(from = 0, to = ceiling(tail(t * f, 1)), length.out = 11),
          labels = c(formatC(lab, format = "f", digits = 2))
        )
      } else {
        plot(gabors[, m], xlab = "", ylab = "", xaxt = "n", yaxt = "n", type = "l", bty = "n", col = cc)
      }

      txt <- paste(
        "A", m,
        ", f=",
        round(frequency[m], 2), "Hz",
        ", t=",
        round(position[m], 2), "s",
        ", sd=",
        round(scale[m], 2), "s",
        ", E=", round(energy[m], 3),
        sep = "")
      mtext(txt, side = 2, line = 0, las = 1, cex = 1)
    }
    dev.off()
  } # if(write.atoms)


  # Empty chart on which the ellipses will appear
  if (draw.ellipses && out.mode == "plot") {
    par(mfrow = c(1, 1), pty = "m", mai = c(0.9, 0.9, 0.2, 0.4))
    plot(0, xlim = c(0, tail(t, 1)), ylim = c(0, tail(y, 1)), type = "n", las = 1,
         xlab = "Time [s]", ylab = "Frequency [Hz]", yaxs = "i", xaxs = "i")
  }

  for (n in 1:num.atoms) {

    if (draw.ellipses && out.mode == "plot") {
      ellipse <- DrawEllipse(
        x = position[n],
        y = frequency[n],
        # from Heisenberg rule: delta_t x delta_omega >= 1/2
        radius.x = (scale[n] / 2),
        radius.y = 1 / ((scale[n])),
        col = "lightgray",
        border = "black",
        plot = TRUE,
        nv = 100)

      if (display.crosses) {
        points(position[n], frequency[n] , pch = 3, col = "black", cex = 1)
      }

      if(display.atom.numbers) {
        text(position[n], frequency[n], n, col = "black", cex = 1)
      }
    }

    if (mode == "sqrt")
      A  <- sqrt(energy[n] * f)

    if (mode == "log")
      A  <- log(energy[n] * f)

    if (mode == "linear")
      A  <- energy[n] * f

    # from Heisenberg rule: delta_t x delta_omega >= 1/2
    radius.x <- (scale[n] / 2) / shortening.factor.x
    radius.y <- 1 / ((scale[n])) / shortening.factor.y

    x0 <- position[n]
    y0 <- frequency[n]
    sx <- radius.x
    sy <- radius.y

    # 2D Gaussian
    z <- with(grid,
              A * exp(-((x - x0)^2 / (2 * sx^2) +
                        (y - y0)^2 / (2 * sy^2))))

    # convert to a matrix, because image() requires it
    z.mtx <- matrix(z, nrow = length(t), ncol = length(y))
    tf.map <- tf.map + z.mtx

  } # for (n in 1:num.atoms)

  if (out.mode == "plot") {
    if (plot.signals) {
      grid.matrix <- cbind(c(1, 1, 1, 2, 3))
      layout(grid.matrix, widths = c(1, 1, 1), heights = c(2, 1, 1))
      # mai: c(bottom, left, top, right)
      par(pty = "m", mai = c(0.4, 0.7, 0.2, 0.4), xaxs = "i", yaxs = "i")
    } else {
      par(mfrow = c(1, 1), pty = "m", mai = c(0.9, 0.9, 0.2, 0.4), xaxs = "i", yaxs = "i")
    }

    # Drawing with graphics::image() is very slow, especially for large matrices.
    # Graphics::rasterImage() is much faster.

    z.col <- col[cut(tf.map, breaks = 129)]
    # image() and rasterImage() differ in the orientation of the Y axis.
    # Often you also need to do the following:
    rot90 <- function(m) t(m)[ncol(m):1, ]
    tf.map.rot90 <- rot90(matrix(z.col, nrow(tf.map)))
    plot.new()
    plot.window(range(t), range(y))
    rasterImage(tf.map.rot90, 0, 0, tail(t, 1), tail(y, 1))

    lab <- seq(from = 0, to = ceiling(tail(t, 1)), length.out = 11)
    axis(
      side = 1, las = 1, cex.axis = 0.9,
      at = seq(from = 0, to = ceiling(tail(t, 1)), length.out = 11),
      labels = c(formatC(lab, format = "f", digits = 2))
    )

    lab <- seq(from = 0, to = ceiling(tail(y, 1)), length.out = 11)
    axis(
      side = 2, las = 1,  cex.axis = 0.9,
      at = seq(from = 0, to = ceiling(tail(y, 1)), length.out = 11),
      labels = c(formatC(lab, format = "f", digits = 2))
    )

    box()
    mtext("Time [s]", side = 1, line = 2, cex = 0.8)
    mtext("Frequency [Hz]", side = 2, line = 3.5, cex = 0.8)

    # At the centres of the atoms, the atom numbers
    if (display.atom.numbers) {
      for (n in 1:num.atoms) {
        text(position[n], frequency[n], n, col = "white", cex = 1)
      }
    }

    # We display small crosses in the centres of atoms
    for (n in 1:num.atoms) {
      if (display.crosses) {
        points(position[n], frequency[n], pch = 3, col = crosses.color, cex = 1)
      }
    }

    if (display.grid)
      grid(col = "grey")

    if (plot.signals) {
      xx <- seq(from = 0, to = epochSize / f, length.out =  epochSize)
      plot(x = xx, original.signal, type = "l", xlab = "", ylab = "", xaxs = "i", las = 1, main = "Original signal", panel.first = grid())
      abline(h = 0, col = "blue")

      plot(x = xx, reconstruction, type = "l", xlab = "", ylab = "", xaxs = "i", las = 1, main = "Reconstructed signal", panel.first = grid())
      abline(h = 0, col = "blue")
    }

  } # if (out.mode == "plot")

  if (out.mode == "file") {
    graphics.off()
    png(file.name, width = size[1], height = size[2], pointsize = 18)
    par(pty = "m", mai = c(0, 0, 0, 0))
    graphics::image(x = t, y = y, z = tf.map, col = col)
    dev.off()
    message("PNG file saved in '", file.name, "'")
  }

  if (out.mode == "RData") {
    im <- as.cimg(tf.map)
    im.resampled <- imager::resize(im, size_x = size[1], size_y = size[2], interpolation_type = 3)
    tf.map.resampled <- as.matrix(im.resampled)

    # Rescaling to the range 0-1.
    # Protect against a situation where a zero appears in the denominator.
    if (max(tf.map.resampled) - min(tf.map.resampled) == 0) {
      tf.map.resampled <- matrix(0, size[1], size[2])
    } else {
      tf.map.resampled <- (tf.map.resampled - min(tf.map.resampled)) / (max(tf.map.resampled) - min(tf.map.resampled))
    }

    save(tf.map.resampled, file = file.name)
    message("RData file saved in '", file.name, "'")
  }

  if (out.mode == "RData2") {
    zz  <- tf.map
    rr <- raster::raster(nrow = ncol(zz), ncol = nrow(zz)) # # this is how it should be: nrow = ncol(zz), ncol = nrow(zz)
    rr[] <- t(zz)
    tt <- raster::raster(ncol = size[1], nrow = size[2])
    tt <- raster::resample(rr, tt)
    tf.map.resampled <- matrix(tt@data@values, size[1], size[2])

    # Rescaling to the range 0-1.
    # Protect against a situation where a zero appears in the denominator.
    if (max(tf.map.resampled) - min(tf.map.resampled) == 0) {
      tf.map.resampled <- matrix(0, size[1], size[2])
    } else {
      tf.map.resampled <- (tf.map.resampled - min(tf.map.resampled)) / (max(tf.map.resampled) - min(tf.map.resampled))
    }

    save(tf.map.resampled, file = file.name)
    message("RData file saved in '", file.name, "'")
  } # if (out.mode == "RData2")

  # Restore
  par(old.par)

  list(
    gabor.functions = gabors,
    reconstruction = reconstruction,
    original.signal = original.signal,
    f = f,
    grid.size.t = t,
    grid.size.f = y,
    epochSize = epochSize,
    number.of.secs = s,
    tf.map = tf.map,
    tf.map.resampled = tf.map.resampled,
    channel = channel,
    freq.divide = freq.divide)
}