Nothing
# 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
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.