R/grf.R

In obreschkow/cooltools: Practical Tools for Scientific Computations and Visualizations

#' Gaussian Random Field generator
#'
#' @importFrom graphics rasterImage curve
#'
#' @description Generates a scalar Gaussian Random Field (GRF) in 1, 2 or 3 dimensions, with a power-law power spectrum of custom slope alpha. The field is normalized such that its mean is zero (up to floating point errors) and its expected standard deviation is one, for alpha=0. The Fourier phases are sampled in order of increasing frequency, such that the random structure is preserved when changing the output size.
#'
#' @param nside integer number of elements per dimension in the output matrix
#' @param dim integer number of dimensions
#' @param alpha spectral index, such that Fourier amplitudes scale as k^alpha. alpha=0 corresponds to uncorrelated white noise, alpha<0 is red noise, while alpha>1 is blue noise.
#' @param seed optional seed for random number generator
#'
#' @return Returns a vector, matrix or 3D-array with \code{nside} elements per dimension
#'
#' @author Danail Obreschkow
#'
#' @examples
#' # Generate the same 2D GRF in two different resolutions
#' nplot(c(0,2.1), asp=1)
#' lowres = grf(nside=30, alpha=-2, seed=1)
#' graphics::rasterImage(stretch(lowres),0,0,1,1)
#' highres = grf(nside=120, alpha=-2, seed=1)
#' graphics::rasterImage(stretch(highres),1.1,0,2.1,1)
#'
#' # Check the power spectrum of a general GRF
#' nside = 50 # change this to any integer >1
#' alpha = -1.7 # change this to any power
#' dim = 3 # change this to any positive integer (but keep nside^dim reasonable)
#' x = grf(nside=nside, dim=dim, alpha=alpha)
#' fourier.grid = dftgrid(N=nside,L=1)
#' knorm = vectornorm(expand.grid(rep(list(fourier.grid$k),dim)))
#' power = abs(dft(x))^2
#' b = bindata(knorm,power,method='equal')
#' plot(b$xmedian,b$ymedian,log='xy',pch=16,
#'      xlab='Fourier mode |k|',ylab='Power p(k)',main='Power spectrum')
#' graphics::curve(x^alpha/nside^dim,col='blue',add=TRUE) # expected power-law
#'
#' @export

grf = function(nside=100, dim=2, alpha=0, seed=NULL) {

  # input checks
  if (nside<2 | nside%%1!=0) stop('nside must be an integer larger than 1')
  if (dim<1 | dim%%1!=0) stop('dim must be a positive integer')

  # generate Fourier grid
  fourier.grid = dftgrid(N=nside,L=1)
  k = vectornorm(expand.grid(rep(list(fourier.grid$k),dim)))

  # generate random phases
  if (!is.null(seed)) set.seed(seed)
  phase = array(exp(2i*pi*runif(nside^dim))[order(order(round(k)))],dim=rep(nside,dim)) # complex phase

  # generate spectrum
  spectrum = k^(alpha/2)*phase
  spectrum[which.min(k)] = 0

  # compute inverse FT and extract real values
  return(Re(dft(spectrum,inverse=TRUE))*sqrt(2/nside^dim))

}