Code covered by the BSD License

# Circulant Embedding method for generating stationary Gaussian field

### Zdravko Botev (view profile)

Implements Dietrich and Newsam's circulant embedding method for fast generation of Gaussian fields.

stationary_Gaussian_process.m
```% generated stationary Gaussian field over an m times n square grid
% for a detailed mathematical explanation of the Matlab code and further
% examples see

% Kroese, D.P. and Botev, Z.I. (2013).
% "Spatial Process Generation."
% V. Schmidt (Ed.). Lectures on Stochastic Geometry,
% Spatial Statistics and Random Fields, Volume II:
% Analysis, Modeling and Simulation of Complex Structures, Springer-Verlag, Berlin.
% http://www.maths.uq.edu.au/~kroese/ps/MCSpatial.pdf

clear all,clc, n=384; m=512; % size of grid is m*n
% size of covariance matrix is m^2*n^2
tx=[0:n-1]; ty=[0:m-1]; % create grid for field
% sample covariance function below;
% change this function to generate a different Gaussian field
rho=@(x,y)((1-x^2/50^2-x*y/(15*50)-y^2/15^2)*exp(-(x^2/50^2+y^2/15^2)));
Rows=zeros(m,n); Cols=Rows;
for i=1:n
for j=1:m
Rows(j,i)=rho(tx(i)-tx(1),ty(j)-ty(1)); % rows of blocks of cov matrix
Cols(j,i)=rho(tx(1)-tx(i),ty(j)-ty(1)); % columns of blocks of cov matrix
end
end
% create the first row of the block circulant matrix with circular blocks
% and store it as a matrix suitable for fft2;
BlkCirc_row=[Rows, Cols(:,end:-1:2);
Cols(end:-1:2,:), Rows(end:-1:2,end:-1:2)];
% compute eigen-values
lam=real(fft2(BlkCirc_row))/(2*m-1)/(2*n-1);

if abs(min(lam(lam(:)<0)))>10^-15
error('Could not find positive definite embedding!')
else
lam(lam(:)<0)=0; lam=sqrt(lam);
end

% generate field with covariance given by block circular matrix
F=fft2(lam.*complex(randn(2*m-1,2*n-1),randn(2*m-1,2*n-1)));
F=F(1:m,1:n); % extract subblock with desired covariance
field1=real(F); field2=imag(F); % two independent fields with desired covariance
imagesc(tx,ty,field1), colormap bone

```