% Generate 1-dimensional Gaussian random field with given covariance function

% this version does it directly
% -- all helper functions are included in the file
% -- we use sins and cosines instead of exponentials

clear;

whichC = 1;  % change to try different covariance functions

ifcheckC = 1;  % 1 to verify covariance function is correct; 0 to skip
ifplotall = 0;  % 1 to plot all realizations in check; 0 to skip
nr = 1000;  % how many realizations to check covariance function

% Ck = spectral density function
% this is the fourier transform of covariance function.
% It can be any function of k that is nonnegative and symmetric about k=0.
switch whichC
    case 1
        Ckf = @(k) ( (abs(k)<5) );
    case 2
        Ckf = @(k)(  exp(-k.^2/2)/sqrt(2*pi) );
    case 3
        Ckf = @(k) (  1./(k.^2+1)/pi  );
end

% domain discretization parameters
Lx = 10;     % length of spatial domain (note boundary conditions are periodic)
Nx = 2^7;    % number of points used to represent function
dx = Lx/Nx;
x = linspace(0,Lx-dx,Nx);  % x-grid points
k = k_of_x(x);             % k (wavenumber) grid points
dk = k(2)-k(1);

% spatial covariance function
Ck = Ckf(k);
[~,C] = mhifft(k,Ck);  % spatial covariance function


% generate random Fourier coefficients
h = grf(k,Ck);

%ht = grf1(k,2*Ck(k),1);   % generates random Fourier coefficients;            
%[x,h] = mhifft(k,ht,1);     % generates GRF by taking inverse Fourier transform


% set up plots so they're nicer
fs = 18;  % default font size
fsl = 12; % font size for legend
set(groot,'defaultAxesFontSize',fs);
set(groot,'DefaultLineLineWidth',2);
set(groot,'defaultAxesLineWidth',2);
set(groot,'defaultAxesFontName','Times');
set(0,'defaulttextInterpreter','latex')

% Plot the Gaussian random field
figure(1)
plot(x,h);
xlabel('x');
title('realization of a Gaussian random field');

% plot the spatial covariance function
figure(2)
tiledlayout(1,2);
nexttile;
xpl = [-fliplr(x(2:end)),x];
Cpl = [fliplr(C(2:end)),C];
plot(xpl,Cpl);
xlabel('x');
title('covariance function $C(t)$');

% plot the spectral density function
nexttile;
plot(k,Ck);
xlabel('k');
title('spectral density function $\hat{C}(k)$');

figure(1)  % puts realization on top




%%% ----------------------------%%%
%%%%%%%%%   Code check   %%%%%%%%%%
%%% ----------------------------%%%
% Now let's check the code is doing the right thing.
% We'll do this by generating many GRFs, and computing the empirical
% covariance function and comparing to C(t).

if(ifcheckC)

    % generate many realizations of a GRF
    h2 = zeros(nr,Nx);
    c2 = zeros(nr,Nx);
    for ir=1:nr
        hh = grf(k,Ck);
        h2(ir,:) = hh;
        c2(ir,:) = mhccov(hh);  % autocovariance of this realization
    end

    % empirical covariance function
    Cest = mean(c2,1);
    Cestpl = [fliplr(Cest(2:end)),Cest];

    % compare empirical and desired covariance function
    figure(4)
    plot(xpl,Cpl,xpl,Cestpl);
    legend({'desired covariance';'estimated covariance'},'location','best');
    %plot(c2');  % plot all estimated covariance functions

    if(ifplotall)
        figure(3)
        plot(x,h2)
    end

end



% remove changes to default settings
set(groot,'defaultAxesFontSize','remove');
set(groot,'defaultAxesLineWidth','remove');
set(groot,'defaultAxesFontName','remove');
set(groot,'DefaultLineLineWidth','remove');




% create a Gaussian random field with given spectral density function
function h = grf(k,Ck)
Nk = length(k);
dk = k(2)-k(1);
x = x_of_k(k);
Nx = length(x);
Ak = sqrt(Ck/2).*randn([1,Nk]);   % sin coefficients
Bk = sqrt(Ck/2).*randn([1,Nk]);   % cos coefficients

% construct random function as 
%     \sum_k Ak*sin(kx)/sqrt(2) +  Bk*cos(kx)/sqrt(2)
% factor of sqrt(2) is a normalization factor for sin/cos
h = dk*(sum(repmat(Ak',1,Nx).*sin(k'*x),1)/sqrt(2)+ sum(repmat(Bk',1,Nx).*cos(k'*x),1)/sqrt(2));

end


% compute autocovariance function of a sequence y by hand
% (since Matlab's function doesn't seem to work)
% assumes y is periodic
function cc = mhccov(y)
ny = length(y);
cc = zeros(1,ny);
for ss = 0:ny-1   % ss is the shift
    ys = circshift(y,ss);
    cc(ss+1) = mean(ys.*y);
end
end


% Fast Fourier INVERSE Transform of the pair (k,ft) into (x,f).  The length of 
% x and f must be an even number, preferably a power of two.  
% x starts with zero.
% If rflag is set then only the real part of f is returned.
function [x , f] = mhifft(k,ft,rflag)
N         = max(size(k));
x         = x_of_k(k);
Period    = N*(x(2)-x(1));
inull     = N/2;
f         = (N/Period)*ifft(circshift(ft,-[0 inull-1]));
if (nargin == 3)
    if (rflag == '1')
        f = real(f);
    end
end
end


% Fast Fourier Transform of the pair (x,f) into (k,ft).  The length of 
% x and f must be an even number, preferably a power of two.  The index of
% the zero mode is inull=N/2.
function [k , ft] = mhfft(x,f)
N         = max(size(x));
k         = k_of_x(x);
Period    = N*(x(2)-x(1));
inull     = N/2;
ft        = (Period/N)*circshift(fft(f),[0 inull-1]);
end

% x-vector with leading zero from k-vector 
function x = x_of_k(k)
N      = max(size(k));
dk     = k(2)-k(1);
dx     = 2*pi/(N*dk);
x      = dx*(linspace(1,N,N)-1);
end


% Wave vector from x-vector
function k = k_of_x(x)
N      = max(size(x));
dx     = x(2)-x(1);
dk     = (2*pi)/(N*dx);
inull  = N/2;
k      = dk*(linspace(1,N,N)-inull);
end

