colored_noise.m

function x = colored_noise(sz,beta,low_thresh)
% function x = colored_noise(DIM, BETA),
%
% This function generates 1/f spatial noise, with a normal error 
% distribution (the grid must be at least 10x10 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 spatial 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 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.

% 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 nargin < 3 || isempty(low_thresh)
    low_thresh = 0;
end

% Generate the grid of frequencies
freq_grid = cell(1,numel(sz));
for d=1:numel(sz)
    % 0:nyquist frequency followed by conjugate (without 0 and nyqist frequency)
    freq_grid{d} = [(0:floor(sz(d)/2)), -(ceil(sz(d)/2)-1:-1:1)]'/sz(d);
end
[freq_grid{:}] = ndgrid(freq_grid{:});

% Generate the power spectrum
S_f = zeros(sz);
for d=1:numel(sz)
    S_f = S_f + freq_grid{d}.^2;
end
if low_thresh > 0
    low_falloff_term = 1 ./ (1+exp(-10 .* (S_f.^0.5 ./ low_thresh - 1))); % sigmoid that goes through (0,~0),(low_thresh,0.5),(2*low_thresh,~1)
else
    low_falloff_term = ones(size(S_f));
end
high_falloff_term = S_f.^(beta/2);
S_f = high_falloff_term .* low_falloff_term;

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

% Generate a grid of random phase shifts
phi = rand(sz);
% phi = rand(floor(DIM./2)+1);

% Inverse Fourier transform to obtain the the spatial pattern
x = ifftn(S_f.^0.5 .* (cos(2*pi*phi)+i*sin(2*pi*phi)));
% x = ifft2(S_f .* (cos(2*pi*phi)+i*sin(2*pi*phi)));
% x = ifft2(S_f(1:floor(DIM(1)/2)+1,1:floor(DIM(2)/2)+1).^0.5 .* (cos(2*pi*phi)+i*sin(2*pi*phi)));

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