cooltools (version 2.4)

grf: 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)

Value

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

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

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