Code covered by the BSD License  

Highlights from
Circulant Embedding method for generating stationary Gaussian field

image thumbnail

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.

% 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.
% weblink:

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
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
% 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

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

% generate field with covariance given by block circular matrix
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

Contact us