grf: Gaussian Random Field generator

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

Gaussian Random Field generator

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.

Usage

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

Arguments

nside	integer number of elements per dimension in the output matrix
dim	integer number of dimensions
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.
seed	optional seed for random number generator

Value

Returns a vector, matrix or 3D-array with nside elements per dimension

Author(s)

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

obreschkow/cooltools documentation built on Nov. 16, 2024, 2:46 a.m.