gaussian_error

 PURPOSE: create a gaussian correlated random error surface 
 -------------------------------------------------------------------
 USAGE: E = gaussian_error(X,Y,r)
 where: [X] and [Y] are coordinate vectors of the target error surface
        [r] is the (correlation) range of the error surface in multiples
        of cellsize as defined in X/Y

 returns matrix [E] of size X x Y
 -------------------------------------------------------------------
 NOTES: Based on the approach of spatial moving averages according to 
 Oksanen and Sarjakoski 2005:
 http://taylorandfrancis.metapress.com/openurl.asp?genre=article&id=doi:10.1080/01431160500057947

 Felix Hebeler, Geography Dept., University Zurich, March 2006.

0001 function E = gaussian_error(X,Y,r)
0002 % PURPOSE: create a gaussian correlated random error surface
0003 % -------------------------------------------------------------------
0004 % USAGE: E = gaussian_error(X,Y,r)
0005 % where: [X] and [Y] are coordinate vectors of the target error surface
0006 %        [r] is the (correlation) range of the error surface in multiples
0007 %        of cellsize as defined in X/Y
0008 %
0009 % returns matrix [E] of size X x Y
0010 % -------------------------------------------------------------------
0011 % NOTES: Based on the approach of spatial moving averages according to
0012 % Oksanen and Sarjakoski 2005:
0013 % http://taylorandfrancis.metapress.com/openurl.asp?genre=article&id=doi:10.1080/01431160500057947
0014 %
0015 % Felix Hebeler, Geography Dept., University Zurich, March 2006.
0016 
0017 if nargin < 3
0018     error('This function needs 3 input parameters (X,Y,range)');
0019 end
0020 
0021 %% Parameters
0022 w=abs(X(1,2)-X(1,1)); %spacing of grid in the convolution (resolution of DEM, as in the paper)
0023 r=r*w;
0024 %R = Matrix of correlation coefficents
0025 %W = weight kernel
0026 
0027 %% create random gaussian error field
0028 E = randn(size(Y,2),size(X,2)); %uncorrelated gaussian (normal) matrix (white noise)
0029 
0030 %% correlation coefficients matrix R
0031 m=(2*((2*r)/w))+1; % dimension of matrix R
0032 
0033 %% Gaussian autocorrelation function
0034 % p(h)=exp(-3*(h^2/r^2)); %p(h) is the correlation coefficient of the DEM
0035 % error at lag h and range r.
0036 R = euklid_W(m,m)*w; % get euklid distance weight matrix & apply cellspacing w
0037 R = exp( -3*( (R.^2) ./ (r^2) ) ); % apply corr coeff  as in paper
0038 
0039 %% Scaling matrix R
0040 C = real(sqrtm(R.^2)); % this gives complex numbers, but close to 0, so ignore
0041 
0042 %% using s
0043 s=sum(C(:).^2);
0044 s = 1/sqrt(s);
0045 
0046 %% getting kernel W
0047 W=s.*C;
0048 
0049 %% do convolution filtering
0050 E=conv2(E,W,'same');
0051