R/makeFunctions.R

In jrwishart/mwaved: Multichannel Wavelet Deconvolution with Additive Long Memory Noise

# Auxiliary file that contains the functions used in the simulations. 

#' @title Multichannel Gamma density blur
#' 
#' @param n An integer specifying the number of observations in each channel.
#' @param shape A numeric vector of length m where each element specifies the shape parameter for the Gamma density used to blur each channel.
#' @param scale A numeric vector of length m where each element specifies the scale parameter for the Gamma density used to blur each channel.
#' 
#' @description Create blur of the regular smooth type generated from Gamma densities
#' 
#' @details Function creates a matrix of dimension n by m which contains a normalised (unit energy in each column) set of blur functions of regular smooth type, generated by a Gamma density. These Gamma densities are generated using the \code{\link{dgamma}} base R function and then normalised to have unit energy.
#' 
#' @return A numeric n by m matrix of normalised Gamma blur.
#' 
#' @seealso \code{\link{dgamma}} for details of the Gamma density
#' @examples
#' n <- 1024
#' m <- 3
#' shape <- seq(from = 0.5, to = 1, length = m)
#' scale <- rep(0.25,m)
#' 
#' # Plot the smooth (gamma) blur
#' x <- (1:n)/n
#' blurMat <- gammaBlur(n, shape, scale)
#' matplot(x, blurMat, type = 'l', main = paste('Set of Gamma', m,'Gamma blur densities.'))
#' 
#' # Plot a LIDAR signal and its multichannel smooth blurred version
#' signal <- makeLIDAR(n)
#' matplot(x, signal, type = 'l', main = 'LIDAR test signal')
#' blurredSignal <- blurSignal(signal, blurMat)
#' matplot(x, blurredSignal, type = 'l', main = 'Smooth blurred LIDAR test signals')
#' @seealso \code{\link{boxcarBlur}}, \code{\link{blurSignal}}
#' @export
gammaBlur <- function(n, shape, scale) {
  if (length(shape) != length(scale)){
    warning("Dimension mismatch: Length of vectors for scale and shape parameters do not match")
  } 
  m <- length(shape)
  G <- sapply((1:n)/n, dgamma, shape = shape, scale = scale)
  
  if( is.null(dim(G)) ){
    G <- G/max(G)
    G <- G/sum(G)
    G <- as.matrix(G)
  } else {
    G <- G/(apply(G, 1, max))
    G <- G/(apply(G, 1, sum))
    G <- t(G)
  }
	return(G) 
}

#' @title Multichannel box car blur
#' 
#' @param n An integer specifying the number of observations in each channel
#' @param width A numeric vector of length m where each element specifies the box car widths to blur in each channel.
#' 
#' @description Create blur of the box car type.
#' 
#' @details Creates a matrix of dimension n by m which contains a normalised (unit energy in each column) set of box car blur functions. The generated box-car functions are governed by an input vector \code{BA} which specifies the box car width. For the method to adhere to the theory described in Wishart (2014), the box car widths should be a Badly Approximable (BA) number.
#' 
#' @return A n by m matrix of normalised box car blur.
#' 
#' @examples
#' n <- 1024
#' width <- 1/sqrt(c(89,353))
#' 
#' # Plot the box car blur
#' blurMat <- boxcarBlur(n, width)
#' x <- (1:n)/n
#' matplot(x, blurMat, type = 'l', main = paste('Set of box car blur functions'))
#' 
#' # Plot a LIDAR signal and its multichannel box car blurred version
#' signal <- makeLIDAR(n)
#' matplot(x, signal, type = 'l', main = 'LIDAR test signal')
#' blurredSignal <- blurSignal(signal, blurMat)
#' matplot(x, blurredSignal, type = 'l', main = 'Box car blurred LIDAR test signals')
#' @seealso \code{\link{gammaBlur}}, \code{\link{blurSignal}}
#' @export
boxcarBlur <- function(n, width) {
  m <- length(width)
  
  box.car <- function(n, a) {
    left <- seq(from = 0, to = 1 - 1/n, length = n)
    left <- (left < a/2) * 1/(a*n)
    aux <- c(0,rev(left[-1]))
    y <- left + aux
  }
  
  G <- sapply(width, function(i) box.car(n, i))
   
  return(G) 
}


#' @name makeSignals
#' @rdname makeSignals
#' @title Generate test signals for simulation
#' @description Generates some of the test signals used the standard nonparametric deconvolution literature.
#' @param n An integer specifying the length of the desired signal.
#' @return A numeric vector of length n giving the desired test signal.
#' @examples
#' n <- 1024
#' x <- (1:n)/n
NULL

#' @rdname makeSignals
#' @examples
#' signal <- makeLIDAR(n)
#' plot(x, signal, main = 'LIDAR test signal', type = 'l')
#' @export
makeLIDAR <- function(n) {
  x <- (1:n)/n
  y <- 0.7 * (1 * (x > 0.15) * (x < 0.65) + 1 * (x > 0.28) * (x < 0.48) + (133.33 * 
    x - 106.66) * (x > 0.8) * (x < 0.815) + (-133.33 * x + 110.6639) * (x > 
    0.815) * (x < 0.83) + (133.33 * x - 119.997) * (x > 0.9) * (x < 0.915) + 
    (-133.33 * x + 123.9969) * (x > 0.915) * (x < 0.93))
  y
}

#' @rdname makeSignals
#' @examples
#' signal <- makeBumps(n)
#' plot(x, signal, main = 'Bumps test signal', type = 'l')
#' @export
makeBumps <- function(n) {
    x <- 1:n/n
    pos <- c(0.1, 0.13, 0.15, 0.23, 0.25, 0.4, 0.44, 0.65, 0.76, 0.78, 0.81)
    hgt <- c(4, 5, 3, 4, 5, 4.2, 2.1, 4.3, 3.1, 5.1, 4.2)
    wth <- c(0.005, 0.005, 0.006, 0.01, 0.01, 0.03, 0.01, 0.01, 0.005, 0.008, 
        0.005)
    y <- rep(0, n)
    for (j in 1:length(pos)) {
        y <- y + hgt[j]/(1 + abs((x - pos[j])/wth[j]))^4
    }
    y <- y * 1.8
  y
}

#' @rdname makeSignals
#' @examples
#' signal <- makeDoppler(n)
#' plot(x, signal, main = 'Doppler test signal', type = 'l')
#' @export
makeDoppler <- function(n) {
    x <- (1:n)/n
    y <- (sqrt(x * (1 - x))) * sin((2 * pi * 1.05)/(x + 0.05))
    y <- y * 2.4
    y
}

#' @rdname makeSignals
#' @examples
#' signal <- makeCusp(n)
#' plot(x, signal, main = 'Cusp test signal', type = 'l')
#' @export
makeCusp <- function(n) {
  x <- (1:n)/n
  y <- 2.4 * sqrt(abs(x - 0.37))
  y
}

#' @rdname makeSignals
#' @examples
#' signal <- makeBlocks(n)
#' plot(x, signal, main = 'Blocks test signal', type = 'l')
#' @export
makeBlocks <- function(n) {
    x <- (1:n)/n
    pos <- c(0.1, 0.13, 0.15, 0.23, 0.25, 0.4, 0.44, 0.65, 0.76, 0.78, 0.81)
    hgt <- c(4, (-5), 3, (-4), 5, (-4.2), 2.1, 4.3, (-3.1), 2.1, (-4.2))
    y <- rep(0, n)
    for (j in 1:length(pos)) {
        y <- y + (1 + sign(x - pos[j])) * (hgt[j]/2)
    }
    y
}

#' @rdname makeSignals
#' @examples
#' signal <- makeHeaviSine(n)
#' plot(x, signal, main = 'HeaviSine test signal', type = 'l')
#' @export
makeHeaviSine <- function(n) {
  x <- (1:n)/n
  y <- 4 * sin(4 * pi * x) - sign(x - 0.3) - sign(0.72 - x)
  y
}

#' @title Blur an input signal 
#' 
#' @description An input signal is blurred by a set of functions to obtain a blurred multichannel signal.
#' 
#' @param signal A numeric vector of length n of the signal of interest
#' @param G An n by m matrix of blur functions to be applied to the signal of interest.
#' 
#' @details Applies the convolution operator to the signal of interest with each column of G to create a multichannel blurred signal. This operation is done in the Fourier domain using the base R fft transforms.
#' @examples
#' n <- 1024
#' m <- 3
#' blur <- gammaBlur(n, shape = seq(from = 0.5, to = 1, length = m), scale = rep(0.25, m))
#' x <- (1:n)/n
#' signal <- makeLIDAR(n)
#' par(mfrow = c(2,1))
#' plot(x, signal, type = 'l', main = 'Direct LIDAR signal')
#' indirectSignal <- blurSignal(signal, blur)
#' matplot(x, indirectSignal, type = 'l', main = 'Set of blurred LIDAR signals')
#' 
#' @seealso \code{\link{gammaBlur}}, \code{\link{boxcarBlur}}
#' @export
blurSignal <- function(signal, G) {
  n <- length(signal)
  G <- as.matrix(G)
  if( dim(G)[1] != n){
    stop("Dimension mismatch: Number of rows in G needs to match length of signal")
  }
  m <- dim(G)[2]
  signal <- matrix(rep(signal, m), ncol = m, nrow = n)
  y <- Re(mvfft(mvfft(signal) * mvfft(G), inverse = TRUE))/n
  y
}

#' @title Generate multichannel noise
#' @description Generate a matrix of multichannel (possibly long memory) noise variables
#' @param n An integer specifying the number of observations per channel.
#' @param alpha A vector specifying the dependence level in each channel.  
#' @param sigma A vector giving the noise levels (standard deviation) for each channel in the multichannel model (see details).
#' @param ... Additional arguments to pass to the \code{fracdiff} package to tightly control the long memory noise.
#' @details Generates a n by m matrix of noise variables. Long memory variables can be generated by the use of the optional \code{fracdiff} package (if installed). The dependence is specified using the \code{alpha} parameter where \code{alpha} = 2 - 2H where H = Hurst parameter. Long memory is ensured when alpha is between 0 and 1 (H between 1/2 and 1). If \code{alpha} is a single element and \code{sigma} has more than one element (multichannel), then the same dependence level of \code{alpha} is used amongst all of the channels. Otherwise the size of \code{alpha} and \code{sigma} should be the same size. 
#' @seealso \code{\link{sigmaSNR}}
#' @examples
#' n <- 1024
#' m <- 3
#' signal <- makeLIDAR(n)
#' blur <- gammaBlur(n, c(0.5, 0.75, 1), rep(1, m))
#' X <- blurSignal(signal, blur)
#' SNR <- 10*1:3
#' sigma <- sigmaSNR(X, SNR)
#' E <- multiNoise(n, sigma, alpha = c(0.5, 0.75, 1))
#' matplot(X + E, type = 'l')
#' @export
multiNoise <- function(n, sigma = 1, alpha = length(sigma), ...) {
  # Simulate independent Gaussian noise by default
  if ( length(alpha) == 1 && alpha %% 1 == 0 ){
    m <- length(sigma)
    noise <- matrix(rnorm(n * m, sd = sigma), ncol = m, byrow = TRUE)
  } else {
    if ( requireNamespace("fracdiff", quietly = TRUE) ){
      m <- length(sigma)
      # Repeat dependence on all channels if alpha is a single element.
      if (length(alpha) == 1 & m > 1){
        alpha <- rep(alpha, m)
      }
      if ( length(sigma) != length(alpha) ){
        stop("Dimension mismatch: length of sigma should equal length of alpha")
      }
      d <- (1 - alpha)/2
      noise <- sapply(1:m, function(x,y) y[x] * fracdiff::fracdiff.sim(n, d = d[x], ...)$series, y = sigma)
    } else {
      # ignore the alpha parameter
      m <- length(sigma)
      noise <- matrix(rnorm(n * m, sd = sigma), ncol = m, byrow = TRUE)
    }
  }
  noise
}

#' @title Determine noise scale levels from specified \bold{S}ignal to \bold{N}oise \bold{R}atios
#' @description Compute the noise scale levels for each channel using the \bold{S}ignal to \bold{N}oise \bold{R}atios
#' @param signal Noisefree multichannel input signal
#' @param SNR A numeric vector specifying the desired \bold{S}ignal to \bold{N}oise \bold{R}atio for each channel.
#' @details The output noise scale levels (theoretical standard deviation for the process noise process in each channel) is governed by the blurred \bold{S}ignal-to-\bold{N}oise \bold{R}atio (SNR) measured in decibels (dB) where,
#' \deqn{SNR  = 10 log_{10} (\frac{||k*f||^2}{\sigma^2)}}
#' and k*f is the blurred signal, \eqn{||\cdot||} is the norm operator and \eqn{\sigma} is the standard deviation of the noise. Roughly speaking, noise levels are considered high, medium and low for the cases 10 dB, 20 dB and 30 dB respectively.
#' @return A numeric vector with m elements giving the scales (standard deviation of the noise in each channel) to achieve the desired SNR.
#' @seealso \code{\link{multiNoise}} \code{\link{multiSigma}}
#' @examples
#' n <- 1024
#' m <- 3
#' signal <- makeLIDAR(n)
#' blur <- gammaBlur(n, c(0.5, 0.75, 1), rep(1, m))
#' X <- blurSignal(signal, blur)
#' SNR <- 10*1:3
#' sigma <- sigmaSNR(X, SNR)
#' E <- multiNoise(n, sigma)
#' sigmaEst <- multiSigma(E)
#' @export
sigmaSNR <- function(signal, SNR) {
  Y <- as.matrix(signal)
  n <- dim(Y)[1]
  m <- dim(Y)[2]
  if( m > 1 && length(SNR) == 1 ){
    # Repeat SNR across all channels
    sigma <- sqrt(apply(Y^2, 2, mean)) * 10^(-SNR/20)
  } else {
    if( length(SNR) >= m ) {
      if( length(SNR) == m ){
        sigma <- sqrt(apply(Y^2, 2, mean)) * 10^(-SNR/20)
      } else {
        warning('Dimension mismatch: Length of SNR too long and must match number of columns of input signal. Only first m elements used.')
        sigma <- sqrt(apply(Y^2, 2, mean)) * 10^(-SNR[1:m]/20)
      }
    } else {
      # Repeat first element of SNR on all channels
      warning("Dimension mismatch: Length of SNR must match number of columns of input signal. Only first element of SNR used")
      sigma <- sqrt(apply(Y^2, 2, mean)) * 10^(-SNR[1]/20)
    }
  }
  sigma
}