Nothing
#' Find valid MRS data files recursively from a directory path.
#' @param dir a directory path.
#' @return a vector of valid MRS data files.
#' @export
find_mrs_files <- function(dir) {
# recursively find all files
all_files <- list.files(dir, recursive = TRUE, full.names = TRUE)
# empty vector of found files
mrs_files <- c()
# check each file, catching any errors
for (file in all_files) {
tryCatch(
{
dummy <- read_mrs(file)
mrs_files <- c(mrs_files, file)
},
error = function(cond) {
}
)
}
return(mrs_files)
}
#' Apply a function over specified array axes.
#' @param x an array.
#' @param axes a vector of axes to apply fun over.
#' @param fun function to be applied.
#' @param ... optional arguments to fun.
#' @return array.
#' @examples
#' z <- array(1:1000, dim = c(10, 10, 10))
#' a <- apply_axes(z, 3, fft)
#' a[1,1,] == fft(z[1,1,])
#' a <- apply_axes(z, 3, sum)
#' a[1,1,] == sum(z[1,1,])
#' @export
apply_axes <- function(x, axes, fun, ...) {
ndim <- length(dim(x))
# apply apply function to each axis
for (axis in axes) {
margin <- c(1:ndim)[-axis]
z <- apply(x, margin, fun, ...)
if (length(dim(z)) != ndim) {
# insert the dropped dimension
dim(z) <- append(dim(z), 1, axis - 1)
x <- z
} else {
# permute back to the original order
perm_vec <- 2:ndim
perm_vec <- append(perm_vec, 1, axis - 1)
print(perm_vec)
x <- aperm(z, perm_vec)
}
}
x
}
#' Simulate an ideal pulse excitation profile by smoothing a top-hat function
#' with a Gaussian.
#' @param bw top-hat bandwidth (Hz).
#' @param sigma Gaussian width smoothing parameter (Hz).
#' @param fa intended flip angle of the pulse.
#' @return data frame containing the frequency scale, excitation profile and
#' corresponding flip-angles.
#' @export
sim_th_excit_profile <- function(bw = 1500, sigma = 50, fa = 180) {
y <- c(rep(0, bw / 2), rep(1, bw), rep(0, bw / 2))
if (sigma != 0) y <- mmand::gaussianSmooth(y, sigma)
width <- length(y)
x <- seq(from = -width / 2, to = width / 2, length.out = width)
y[y > 1] <- 1 # acos doesn't like values gt 1
y[y < -1] <- -1 # acos doesn't like values lt -1
angle <- acos(1 - 2 * y) * fa / pi
return(data.frame(freq = x, Mxy = y, fa = angle))
}
#' Repeat an array over a given dimension.
#' @param x array.
#' @param rep_dim dimension to extend.
#' @param n number of times to repeat.
#' @return extended array.
#' @export
rep_array_dim <- function(x, rep_dim, n) {
# add a dimension if needed
if ((rep_dim == 7) && length(dim(x)) == 6) dim(x) <- c(dim(x), 1)
dims <- length(dim(x))
orig_dim <- dim(x)
new_dim <- orig_dim
new_dim[rep_dim] <- new_dim[rep_dim] * n
# make rep_dim the last dimension
perm_vec <- (1:dims)[-rep_dim]
perm_vec <- c(perm_vec, rep_dim)
z <- aperm(x, perm_vec)
# duplicate the data
z <- rep(z, n)
# set the new dimensions
dim(z) <- new_dim[perm_vec]
# reorder
perm_vec <- append((1:(dims - 1)), dims, rep_dim - 1)
aperm(z, perm_vec)
}
#rep_array_dim <- function(x, rep_dim, n) {
# if (dim(x)[rep_dim] != 1) stop("Starting dimension extent does not equal one.")
# z <- replicate(n, x)
# dims <- length(dim(z))
# perm_vec <- 1:dims
# perm_vec[rep_dim] <- dims
# perm_vec[dims] <- rep_dim
# z <- aperm(z, perm_vec)
# # drop the last dimension
# abind::adrop(z, dims)
#}
#' Covert a beta value in the time-domain to an equivalent linewidth in Hz:
#' x * exp(-i * t * t * beta).
#' @param beta beta damping value.
#' @return linewidth value in Hz.
#' @export
beta2lw <- function(beta) {2 * (-beta * log(0.5)) ^ 0.5 / pi}
#' Covert a linewidth in Hz to an equivalent beta value in the time-domain ie:
#' x * exp(-t * t * beta).
#' @param lw linewidth in Hz.
#' @return beta damping value.
#' @export
lw2beta <- function(lw) {(lw * pi / 2) ^ 2 / (-log(0.5))}
alpha2lw <- function(alpha) {alpha / pi}
#' Covert a linewidth in Hz to an equivalent alpha value in the time-domain ie:
#' x * exp(-t * alpha).
#' @param lw linewidth in Hz.
#' @return beta damping value.
#' @export
lw2alpha <- function(lw) {lw * pi}
#' Perform a fft and ffshift on a vector.
#' @param vec_in vector input.
#' @return output vector.
#' @export
ft_shift <- function(vec_in) {pracma::fftshift(stats::fft(vec_in))}
#' Perform an iffshift and ifft on a vector.
#' @param vec_in vector input.
#' @return output vector.
#' @export
ift_shift <- function(vec_in) {pracma::ifft(pracma::ifftshift(vec_in))}
#' Perform a fft and fftshift on a matrix with each column replaced by its
#' shifted fft.
#' @param mat_in matrix input.
#' @return output matrix.
#' @export
ft_shift_mat <- function(mat_in) {
mat_in_ft <- stats::mvfft(mat_in)
mvfftshift(mat_in_ft)
}
#' Perform an ifft and ifftshift on a matrix with each column replaced by its
#' shifted ifft.
#' @param mat_in matrix input.
#' @return output matrix.
#' @export
ift_shift_mat <- function(mat_in) {
mat_in_shift <- mvifftshift(mat_in)
stats::mvfft(mat_in_shift, inverse = TRUE) / NROW(mat_in)
}
#' Perform a fftshift on a matrix, with each column replaced by its shifted
#' result.
#' @param x matrix input.
#' @return output matrix.
#' @export
mvfftshift <- function(x) {
m <- NROW(x)
p <- ceiling(m / 2)
idx <- c((p + 1):m, 1:p)
x[idx,,drop = FALSE]
}
#' Perform an ifftshift on a matrix, with each column replaced by its shifted
#' result.
#' @param x matrix input.
#' @return output matrix.
#' @export
mvifftshift <- function(x) {
m <- NROW(x)
p <- floor(m / 2)
idx <- c((p + 1):m, 1:p)
x[idx,,drop = FALSE]
}
hilbert <- function(x) {
x <- Re(x)
x <- stats::fft(x)
n <- length(x)
h <- vector()
h[1] <- 1
h[floor(n / 2) + 1] <- 1
h[2:round(n / 2 + 0.5)] <- 2
h[floor(n / 2 + 2):n] <- 0
y <- h * x
stats::fft(y, inverse = TRUE) / n
}
re_max <- function(x) {
max(Re(x))
}
re_max_interp <- function(data_pts, interp_f) {
if (is.na(data_pts[1])) return(NA)
data_pts <- stats::spline(Re(data_pts), n = interp_f * length(data_pts))
max(data_pts$y)
}
re_mean <- function(x) {
mean(Re(x))
}
re_sd <- function(x) {
stats::sd(Re(x))
}
is.installed <- function(mypkg) {
is.element(mypkg, utils::installed.packages()[,1])
}
# From : http://mathworld.wolfram.com/LeastSquaresFittingExponential.html
measure_lorentz_amp <- function(y, t, start_pnt = 10, end_pnt = 50) {
# check for masked data
if (is.na(y[1])) return(NA)
# crop to time region and take Mod
y <- Mod(y[start_pnt:end_pnt])
t <- t[start_pnt:end_pnt]
a <- sum((t ^ 2) * y) * sum(y * log(y)) - sum(t * y) * sum(t * y * log(y))
a <- a / (sum(y) * sum(t ^ 2 * y) - sum(t * y) ^ 2)
A <- exp(a)
A
}
measure_td_amp <- function(y, start_pnt = 10, end_pnt = 50) {
# check for masked data
if (is.na(y[1])) return(NA)
# crop to time region and take Mod
y <- Mod(y[start_pnt:end_pnt])
stats::spline(start_pnt:end_pnt, y, xout = 1, method = "natural")$y
}
measure_td_amp_poly <- function(y, start_pnt = 10, end_pnt = 50) {
# check for masked data
if (is.na(y[1])) return(NA)
x <- start_pnt:end_pnt
y <- Mod(y[x])
lm_res <- stats::lm(y ~ poly(x, 3))
new <- data.frame(x = 1)
return(stats::predict(lm_res, new))
}
# return the sum of squares of a complex vector
cplx_sum_sq <- function(x) sum(Re(x) ^ 2) + sum(Im(x) ^ 2)
#' Perform a polynomial fit, subtract and return the standard deviation of the
#' residuals.
#' @param y array.
#' @param degree polynomial degree.
#' @return standard deviation of the fit residuals.
#' @export
calc_sd_poly <- function(y, degree = 1) {
x <- 1:length(y)
model <- stats::lm(y ~ poly(x, degree))
stats::sd(model$residual)
}
#' Compute the vector cross product between vectors x and y. Adapted from
#' http://stackoverflow.com/questions/15162741/what-is-rs-crossproduct-function
#' @param x vector of length 3.
#' @param y vector of length 3.
#' @return vector cross product of x and y.
#' @export
crossprod_3d <- function(x, y) {
i <- 1:3
# Project inputs into 3D, since the cross product only makes sense in 3D.
To3D <- function(x) utils::head(c(x, rep(0, 3)), 3)
x <- To3D(x)
y <- To3D(y)
# Indices should be treated cyclically (i.e., index 4 is "really" index 1, and
# so on). Index3D() lets us do that using R's convention of 1-based (rather
# than 0-based) arrays.
Index3D <- function(i) (i - 1) %% 3 + 1
# The i'th component of the cross product is:
# (x[i + 1] * y[i + 2]) - (x[i + 2] * y[i + 1])
# as long as we treat the indices cyclically.
return(x[Index3D(i + 1)] * y[Index3D(i + 2)] -
x[Index3D(i + 2)] * y[Index3D(i + 1)])
}
# stolen from the interweb
add_alpha <- function(col, alpha = 1) {
if (missing(col))
stop("Please provide a vector of colours.")
apply(sapply(col, grDevices::col2rgb) / 255, 2,
function(x) grDevices::rgb(x[1], x[2], x[3], alpha = alpha))
}
l2_norm_vec <- function(x) {
return(x / norm(x, "2"))
}
crop_range <- function(map, lower, upper) {
# TODO check upper > lower and both are >0<1
map_range <- range(map, na.rm = TRUE)
#upper_lim <- map_range[1]+upper/100*(map_range[2] - map_range[1])
#lower_lim <- map_range[1]+lower/100*(map_range[2] - map_range[1])
upper_lim <- upper
lower_lim <- lower
map <- ifelse(map > upper_lim, upper_lim,map)
ifelse(map < lower_lim, lower_lim, map)
}
# https://stackoverflow.com/questions/13432863/determine-level-of-nesting-in-r
depth <- function(this) ifelse(is.list(this), 1L + max(sapply(this, depth)), 0L)
#' Simulate MRS data with a similar appearance to normal brain (by default).
#' @param acq_paras list of acquisition parameters or an mrs_data object. See
#' \code{\link{def_acq_paras}}.
#' @param type type of spectrum, only "normal" is implemented currently.
#' @param pul_seq pulse sequence function to use.
#' @param xlim range of frequencies to simulate in ppm.
#' @param full_output when FALSE (default) only output the simulated MRS data.
#' When TRUE output a list containing the MRS data, basis set object and
#' corresponding amplitudes.
#' @param amps a vector of basis amplitudes may be specified to modify the
#' output spectrum.
#' @param basis_lb apply additional Gaussian line-broadening to the basis (Hz).
#' @param zero_lip_mm zero the amplitudes of any lipid or macromolecular
#' components based on their name starting with "MM" or "Lip".
#' @param remove_lip_mm remove any lipid or macromolecular basis components
#' based on their name starting with "MM" or "Lip".
#' @param ... extra parameters to pass to the pulse sequence function.
#' @return see full_output option.
#' @export
sim_brain_1h <- function(acq_paras = def_acq_paras(), type = "normal_v2",
pul_seq = seq_slaser_ideal, xlim = c(0.5, 4.2),
full_output = FALSE, amps = NULL,
basis_lb = NULL, zero_lip_mm = FALSE,
remove_lip_mm = FALSE, ...) {
if (type == "normal_v1") {
brain_basis_paras <- get_1h_brain_basis_paras_v1(ft = acq_paras$ft)
if (is.null(amps)) {
amps <- c(0.000000000, 0.009799548, 0.072152490, 0.077845526, 0.045575002,
0.005450371, 0.000000000, 0.028636132, 0.076469056, 0.028382618,
0.069602483, 0.001763720, 0.042031981, 0.013474549, 0.000000000,
0.000000000, 0.015705756, 0.001903173, 0.063409950, 0.051807236,
0.222599530, 0.110431759, 0.023451957, 0.000000000, 0.032474090,
0.007033074, 0.000000000)
}
} else if (type == "normal_v2") {
brain_basis_paras <- get_1h_brain_basis_paras_v3(ft = acq_paras$ft)
if (is.null(amps)) {
# levels are taken from within the range specified by "In Vivo NMR
# Spectroscopy 3rd ed" in normal grey matter
amps <- c("-CrCH2" = 0,
"Ala" = 0.2,
"Asp" = 1.5,
"Cr" = 5,
"GABA" = 1,
"Glc" = 1,
"Gln" = 3,
"Gly" = 1,
"GSH" = 2,
"Glu" = 10,
"GPC" = 1.2,
"Ins" = 5,
"Lac" = 0.5,
"Lip09" = 0,
"Lip13a" = 0,
"Lip13b" = 0,
"Lip20" = 0,
"MM09" = 4,
"MM12" = 4,
"MM14" = 4,
"MM17" = 4,
"MM20" = 4,
"NAA" = 10,
"NAAG" = 2,
"PCh" = 0.5,
"PCr" = 4.5,
"PEth" = 0.8,
"sIns" = 0.4,
"Tau" = 1.5)
}
} else {
stop("invalid spectrum type")
}
basis <- sim_basis(brain_basis_paras, pul_seq, acq_paras, xlim = xlim, ...)
if (zero_lip_mm) {
if (length(grep("^MM", basis$names)) > 0) {
zero_idx <- grep("^MM", basis$names)
amps[zero_idx] <- 0
}
if (length(grep("^Lip", basis$names)) > 0) {
zero_idx <- grep("^Lip", basis$names)
amps[zero_idx] <- 0
}
}
if (remove_lip_mm) {
if (length(grep("^MM", basis$names)) > 0) {
rem_idx <- grep("^MM", basis$names)
amps <- amps[-rem_idx]
basis$names <- basis$names[-rem_idx]
basis$data <- basis$data[, -rem_idx]
}
if (length(grep("^Lip", basis$names)) > 0) {
rem_idx <- grep("^Lip", basis$names)
amps <- amps[-rem_idx]
basis$names <- basis$names[-rem_idx]
basis$data <- basis$data[, -rem_idx]
}
}
names(amps) <- basis$names
if (!is.null(basis_lb)) basis <- lb(basis, basis_lb)
mrs_data <- basis2mrs_data(basis, sum_elements = TRUE, amps = amps)
if (!full_output) {
return(mrs_data)
} else {
return(list(mrs_data = mrs_data, basis = basis, amps = amps))
}
}
#' Get the point spread function (PSF) for a 2D phase encoded MRSI scan.
#' @param FOV field of view in mm.
#' @param mat_size acquisition matrix size (not interpolated).
#' @param sampling can be either "circ" for circular or "rect" for rectangular.
#' @param hamming should Hamming k-space weighting be applied (default FALSE).
#' @param ensure_odd add 1mm to the FOV when required to ensure the output pdf
#' has odd dimensions. Required when using get_mrsi2d_seg.
#' @return A matrix of the PSF with 1mm resolution.
#' @export
get_2d_psf <- function(FOV = 160, mat_size = 16, sampling = "circ",
hamming = FALSE, ensure_odd = TRUE) {
# round FOV to ensure int value
FOV <- round(FOV)
if (ensure_odd) {
if ((FOV %% 2) == 0) FOV <- FOV + 1
}
if (sampling == "circ") {
g <- expand.grid(1:FOV, 1:FOV)
dist <- sqrt((g$Var1 - FOV / 2 - 0.5) ^ 2 + (g$Var2 - FOV / 2 - 0.5) ^ 2)
dist <- matrix(dist, FOV, FOV)
kspace <- 1 * (dist <= mat_size / 2)
} else if (sampling == "rect") {
kspace <- matrix(0, nrow = FOV, ncol = FOV)
start_ind <- FOV / 2 - mat_size / 2 + 1
inds <- seq(from = start_ind, by = 1, length.out = mat_size)
kspace[inds, inds] <- 1
} else {
stop("Incorrect sampling option.")
}
if (hamming) {
ham_mat <- pracma::repmat(signal::hamming(mat_size), mat_size, 1)
ham_mat <- ham_mat * t(ham_mat)
inds <- seq(from = start_ind, by = 1, length.out = mat_size)
kspace[inds, inds] <- kspace[inds, inds] * ham_mat
}
kspace_shift <- apply(kspace, 1, pracma::ifftshift)
kspace_shift <- apply(kspace_shift, 1, pracma::ifftshift)
imspace <- apply(kspace_shift, 1, stats::fft, FALSE)
imspace <- apply(imspace, 1, stats::fft, FALSE)
imspace <- apply(imspace, 1, pracma::fftshift)
imspace <- apply(imspace, 1, pracma::fftshift)
imspace <- imspace / max(Re(imspace))
return(imspace)
}
#' Mask fit result spectra depending on a vector of bool values.
#' @param fit_result fit result object to be masked.
#' @param mask_vec a Boolean vector with the same number of rows as there are
#' rows in the results table.
#' @param amps_only only mask the amplitude and associated error estimate
#' columns.
#' @return a masked fit result object.
#' @export
mask_fit_res <- function(fit_result, mask_vec, amps_only = FALSE) {
if (length(mask_vec) != nrow(fit_result$res_tab)) {
stop("mask_vec and fit_res dimensions do not agree")
}
# mask any NA values too
mask_vec[is.na(mask_vec)] <- TRUE
if (amps_only) {
fit_result$res_tab[mask_vec, 6:(fit_result$amp_cols * 2 + 5)] <- NA
} else {
fit_result$res_tab[mask_vec, 6:ncol(fit_result$res_tab)] <- NA
}
return(fit_result)
}
#' @importFrom RNifti readNifti
#' @export
RNifti::readNifti
#' @importFrom RNifti writeNifti
#' @export
RNifti::writeNifti
# prob not as fast as functions in matrixStats but better than
# apply
colSdColMeans <- function(x, na.rm) {
if (na.rm) {
n <- colSums(!is.na(x)) # thanks @flodel
} else {
n <- nrow(x)
}
colVar <- colMeans(x * x, na.rm = na.rm) - (colMeans(x, na.rm = na.rm)) ^ 2
return(sqrt(colVar * n / (n - 1)))
}
# Rotation matrix from axis and angle
# https://en.wikipedia.org/wiki/Rotation_matrix#cite_note-5
rotate_vec <- function(vec_in, ax, theta) {
ct <- cos(theta)
st <- sin(theta)
rotate_mat <- matrix(nrow = 3, ncol = 3)
rotate_mat[1, 1] = ct + ax[1] * ax[1] * (1 - ct)
rotate_mat[1, 2] = ax[1] * ax[2] * (1 - ct) - ax[3] * st
rotate_mat[1, 3] = ax[1] * ax[3] * (1 - ct) + ax[2] * st
rotate_mat[2, 1] = ax[2] * ax[1] * (1 - ct) + ax[3] * st
rotate_mat[2, 2] = ct + ax[1] * ax[1] * (1 - ct)
rotate_mat[2, 3] = ax[2] * ax[3] * (1 - ct) - ax[1] * st
rotate_mat[3, 1] = ax[3] * ax[1] * (1 - ct) - ax[2] * st
rotate_mat[3, 2] = ax[3] * ax[2] * (1 - ct) + ax[1] * st
rotate_mat[3, 3] = ct + ax[3] * ax[3] * (1 - ct)
vec_out <- rotate_mat %*% vec_in
as.numeric(vec_out / sum((vec_out) ^ 2) ^ 0.5)
}
cross <- function(a, b) {
vec_out <- c(NA, NA, NA)
vec_out[1] <- a[2] * b[3] - a[3] * b[2]
vec_out[2] <- a[3] * b[1] - a[1] * b[3]
vec_out[3] <- a[1] * b[2] - a[2] * b[1]
vec_out <- vec_out / (sum(vec_out ^ 2) ^ 0.5)
vec_out
}
#' Create a logical circular mask spanning the full extent of an n x n matrix.
#' @param d diameter of the mask.
#' @param n number of matrix rows and columns.
#' @param offset offset the mask centre in matrix dimension units.
#' @return logical n x n mask matrix.
#' @export
circ_mask <- function(d, n, offset = 1) {
g <- expand.grid(1:n, 1:n)
dist <- sqrt((g$Var1 - n / 2 - offset) ^ 2 + (g$Var2 - n / 2 - offset) ^ 2)
dist <- matrix(dist, n, n)
return(dist <= d / 2)
}
# equation for following function taken from:
# https://math.stackexchange.com/questions/2645689/what-is-the-parametric-equation-of-a-rotated-ellipse-given-the-angle-of-rotatio
#' Create an elliptical mask stored as a matrix of logical values.
#' @param xN number of pixels in the x dimension.
#' @param yN number of pixels in the y dimension.
#' @param x0 centre of ellipse in the x direction in units of pixels.
#' @param y0 centre of ellipse in the y direction in units of pixels.
#' @param xr radius in the x direction in units of pixels.
#' @param yr radius in the y direction in units of pixels.
#' @param angle angle of rotation in degrees.
#' @return logical mask matrix with dimensions fov_yN x fov_xN.
#' @export
elliptical_mask <- function(xN, yN, x0, y0, xr, yr, angle) {
g <- pracma::meshgrid(0:(xN - 1) - xN / 2, 0:(yN - 1) - yN / 2)
phi <- angle / 180 * pi
mask <- (((g$X - x0) * cos(phi) + (g$Y - y0) * sin(phi)) ^ 2 / xr ^ 2 +
(((g$X - x0) * sin(phi) - (g$Y - y0) * cos(phi))) ^ 2 / yr ^ 2) <= 1
return(mask)
}
#' Create a rectangular mask stored as a matrix of logical values.
#' @param xN number of pixels in the x dimension.
#' @param yN number of pixels in the y dimension.
#' @param x0 centre of rectangle in the x direction in units of pixels.
#' @param y0 centre of rectangle in the y direction in units of pixels.
#' @param xw width in the x direction in units of pixels.
#' @param yw width in the y direction in units of pixels.
#' @param angle angle of rotation in degrees.
#' @return logical mask matrix with dimensions fov_yN x fov_xN.
#' @export
rectangular_mask <- function(xN, yN, x0, y0, xw, yw, angle) {
g <- pracma::meshgrid(0:(xN - 1) - xN / 2, 0:(yN - 1) - yN / 2)
phi <- angle / 180 * pi
x_mask <- (((g$X - x0) * cos(phi) + (g$Y - y0) * sin(phi)) ^ 2 /
(xw / 2) ^ 2 <= 1)
y_mask <- (((g$X - x0) * sin(phi) - (g$Y - y0) * cos(phi)) ^ 2 /
(yw / 2) ^ 2 <= 1)
mask <- x_mask & y_mask
return(mask)
}
#' Create a two dimensional Gaussian window function stored as a matrix.
#' @param xN number of pixels in the x dimension.
#' @param yN number of pixels in the y dimension.
#' @param x0 centre of window function in the x direction in units of pixels.
#' Note, only integer values are applied.
#' @param y0 centre of window function in the y direction in units of pixels.
#' Note, only integer values are applied.
#' @param xw the reciprocal of the standard deviation of the Gaussian window in
#' x direction.
#' @param yw the reciprocal of the standard deviation of the Gaussian window in
#' y direction.
#' @return matrix with dimensions fov_yN x fov_xN.
#' @export
gausswin_2d <- function(xN, yN, x0, y0, xw, yw) {
x_gaus_vec <- signal::gausswin(xN, xw)
x_gaus_vec <- pracma::circshift(x_gaus_vec, x0)
x_gaus_mat <- pracma::repmat(x_gaus_vec, yN, 1)
y_gaus_vec <- signal::gausswin(yN, yw)
y_gaus_vec <- pracma::circshift(y_gaus_vec, y0)
y_gaus_mat <- pracma::repmat(y_gaus_vec, xN, 1)
y_gaus_mat <- t(y_gaus_mat)
return(x_gaus_mat * y_gaus_mat)
}
#' Check if an object is defined, which is the same as being not NULL.
#' @param x object to test for being NULL.
#' @return logical value.
is.def <- function(x) !is.null(x)
convolve_td <- function(x, y) {
x <- stats::fft(x * Conj(y) * length(x))
return(x)
}
# Directly from the matrixcalc package
# I don't like copying code like this, but CRAN don't like packages
# that have lots of dependencies so...
hankel.matrix <- function(n, x) {
if (!is.numeric(n))
stop("argument n is not numeric")
if (n != trunc(n))
stop("argument n ix not an integer")
if (n < 2)
stop("argument n is less than 2")
if (!is.vector(x))
stop("argument x is not a vector")
m <- length(x)
if (m < n)
stop("length of argument x is less than n")
y <- x
H <- matrix(0, nrow = n, ncol = n)
H[1, ] <- y[1:n]
for (i in 2:n) {
y <- c(y[2:m], y[1])
H[i, ] <- y[1:n]
}
return(H)
}
# Directly from the MASS package
ginv <- function(X, tol = sqrt(.Machine$double.eps)) {
if (length(dim(X)) > 2L || !(is.numeric(X) || is.complex(X)))
stop("'X' must be a numeric or complex matrix")
if (!is.matrix(X))
X <- as.matrix(X)
Xsvd <- svd(X)
if (is.complex(X))
Xsvd$u <- Conj(Xsvd$u)
Positive <- Xsvd$d > max(tol * Xsvd$d[1L], 0)
if (all(Positive))
Xsvd$v %*% (1/Xsvd$d * t(Xsvd$u))
else if (!any(Positive))
array(0, dim(X)[2L:1L])
else Xsvd$v[, Positive, drop = FALSE] %*% ((1/Xsvd$d[Positive]) *
t(Xsvd$u[, Positive, drop = FALSE]))
}
#' Matrix exponential function taken from complexplus package to reduce the
#' number of spant dependencies.
#' @param x a square complex matrix.
#' @return the matrix exponential of x.
#' @export
matexp <- function(x) {
d <- dim(x)
n <- d[1]
p <- d[2]
Ar <- Re(x)
Ai <- Im(x)
E <- rbind(cbind(Ar, -Ai), cbind(Ai, Ar))
eE <- expm::expm(E)
eA <- eE[1:n, 1:n] + (0 + 1i) * eE[(1:n) + n, 1:n]
eA <- matrix(Imzap(eA), ncol = n)
return(eA)
}
#' Complex rounding function taken from complexplus package to reduce the number
#' of spant dependencies.
#' @param x a scalar or vector, real or complex.
#' @param tol a tolerance, 10^-6 by default. Prevents possible numerical
#' problems. Can be set to 0 if desired.
#' @export
Imzap <- function(x, tol = 1e-06) {
if (all(abs(Im(z <- zapsmall(x))) <= tol))
suppressWarnings(as.double(x))
else x
}
# a trick to allow nice interpolation to work with data containing NAs
interpolate_nas <- function(map, factor, sigma) {
map_box_interp <- mmand::rescale(map, factor, mmand::boxKernel())
map_gaus <- mmand::gaussianSmooth(map, sigma)
map_gaus[is.na(map_gaus)] <- 0 # set any remaing NA values to 0
map_replace <- map
map_replace[is.na(map_replace)] <- map_gaus[is.na(map_replace)]
map_interp <- mmand::rescale(map_replace, factor, mmand::mnKernel())
map_interp[is.na(map_box_interp)] <- NA
return(map_interp)
}
# is this dataset SVS?
is_svs <- function(mrs_data) identical(dim(mrs_data$data)[2:4], c(1L, 1L, 1L))
#' Calculate the mean of adjacent blocks in a vector.
#' @param x input vector.
#' @param block_size number of adjacent elements to average over.
#' @return vector data averaged in blocks.
#' @export
mean_vec_blocks <- function(x, block_size) {
if ((length(x) %% block_size) != 0) {
warning("Block size does not fit into vector length without truncation.")
}
new_dyns <- floor(length(x) / block_size)
x_out <- x[seq(1, new_dyns * block_size, block_size)]
for (n in 2:block_size) {
x_out <- x_out + x[seq(n, new_dyns * block_size, block_size)]
}
return(x_out / block_size)
}
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.