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

by

 

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

Contact us