spatialPatternn.m

function x = spatialPatternn(DIM,BETA) 
% function x = spatialPattern(DIM, BETA)
%
% This function generates 1/f ndimensional noise, with a normal error 
% distribution (the grid must be at least 10x10x... for the errors to be normal). 
% 1/f noise is scale invariant, there is no spatial scale for which the 
% variance plateaus out, so the process is non-stationary.
%
%     DIM is a two component vector that sets the size of the pattern
%           (DIM=[10,5] is a 10x5 spatial grid)
%     BETA defines the spectral distribution. 
%          Spectral density S(f) = N f^BETA
%          (f is the frequency, N is normalisation coeff).
%               BETA = 0 is random white noise.  
%               BETA  -1 is pink noise
%               BETA = -2 is Brownian noise
%          The fractal dimension is related to BETA by, D = (6+BETA)/2
% 
% Note that the pattern will have the same class as DIM. When generating a
% big pattern, this might be useful to spare memory.
%
% Note that the spatial pattern is periodic.  If this is not wanted the
% grid size should be doubled and only the first quadrant used.
%
% Time series can be generated by setting one component of DIM to 1

% The method is briefly descirbed in Lennon, J.L. "Red-shifts and red
% herrings in geographical ecology", Ecography, Vol. 23, p101-113 (2000)
%
% Many natural systems look very similar to 1/f processes, so generating
% 1/f noise is a useful null model for natural systems.
%
% The errors are normally distributed because of the central
% limit theorem.  The phases of each frequency component are randomly
% assigned with a uniform distribution from 0 to 2*pi. By summing up the
% frequency components the error distribution approaches a normal
% distribution.

% rewritten by Maximilien Chaumon, Nov 3 2012.
%
% Based on spatialPattern.m
% Written by Jon Yearsley  1 May 2004
%     j.yearsley@macaulay.ac.uk
%
% S_f corrected to be S_f = (u.^2 + v.^2).^(BETA/2);  2/10/05
% 

if not(exist('BETA','var')) || isempty(BETA)
    BETA = 0;
end
c = class(DIM);
% Generate the grid of frequencies. u is the set of frequencies along
% each dimension
S_f=zeros(DIM,c);
for i_dim = 1:numel(DIM)
    u = [(0:floor(DIM(i_dim)/2)) -(fliplr(1:ceil(DIM(i_dim)/2)-1)) ]' / DIM(i_dim);
    r = ones(1,numel(DIM));
    r(i_dim) = DIM(i_dim);
    u = reshape(u,r);
    r = DIM;
    r(i_dim) = 1;
    u = repmat(u,r);
    u = abs(u);
    S_f = S_f + u;
end

S_f = S_f .^(BETA);

% Set any infinities to zero
S_f(S_f==inf) = 0;

% Generate a grid of random phase shifts
phi = rand(DIM)*2*pi;

% Inverse Fourier transform to obtain the the spatial pattern
x = ifftn(S_f .^ .5 .* exp(1i*phi));

% Pick just the real component
x = real(x);