tutorial_gauss_data.m

% TUTORIAL_GAUSS_DATA
% Generates samples from three non-linearly coupled Gaussian processes
%
% data = tutorial_gauss_data(ntrials, nsamples, delay, cfunc)
%
% where
%
% NTRIALS is the number of randomly initialized repetitions of the 
% time-series. Default 50.
%
% NSAMPLES is the number of data samples to generate in each trial. Default
% 1500.
%
% DELAY are the coupling delays for the link 2<-1 and 3<-1, respectively.
% Default: [10, 15]
%
% CFUNC is the non-linearity to use in the coupling. It is given as a
% function_handle or as a string. In the latter case, a function_handle
% object will be generated by CFUNC = eval(['cf = @(x) ' CFUNC '(x);']).
% Default: @(x) sin(x).
%
% See also: tutorial_gauss_analysis, tutorial_gauss_figures

% Description: Data generation
% Documentation: tutorial_gauss.txt
% Author: German Gomez-Herrero

function data = tutorial_gauss_data(ntrials,nsamples,delay,cfunc)

OUT_FILENAME = 'gaussian_process.mat';

if nargin < 4 || isempty(cfunc),
    cf = @(x) sin(x);
elseif ischar(cfunc),
    eval(['cf = @(x) ' cfunc '(x);']);
end
if nargin < 3 || isempty(delay), delay = [10 15]; end
if nargin < 2 || isempty(nsamples), nsamples = 1500; end
if nargin < 1 || isempty(ntrials), ntrials = 50; end

% Coupling modulation parameters
f       = 0.2;      
dt      = 0.01;
pos = round((nsamples/6).*[1 3 5 6]);

taumax = max(delay);
X = randn(ntrials,nsamples+taumax);
Y = randn(ntrials,nsamples+taumax);
Z = randn(ntrials,nsamples+taumax);
for mciter = 1:ntrials    
    for i=taumax:nsamples
        if (i>=pos(1) && i <=pos(2)),
            % 2<-1 "sin" coupling factor modulation
            X(mciter,i+1)= 0.4*X(mciter,i)+ X(mciter,i+1);
            Y(mciter,i+1)= 0.5*Y(mciter,i)+ sin(2*pi*f*i*dt)*cf(X(mciter,...
                i-delay(1)+1))+ Y(mciter,i+1);
            Z(mciter,i+1)= 0.6*Z(mciter,i)+ Z(mciter,i+1);        
        elseif (i>=pos(2) && i <=pos(3)),
            % 3<-2 "cos" coupling factor modulation
            X(mciter,i+1)= 0.4*X(mciter,i)+ X(mciter,i+1);
            Y(mciter,i+1)= 0.5*Y(mciter,i)+ Y(mciter,i+1);%
            Z(mciter,i+1)= 0.6*Z(mciter,i)+ cos(2*pi*f*i*dt)*cf(Y(mciter,...
                i-delay(2)+1))+ Z(mciter,i+1);    
        else
            % The three processes are uncoupled
            X(mciter,i+1)= 0.4*X(mciter,i)+ X(mciter,i+1);
            Y(mciter,i+1)= 0.5*Y(mciter,i)+ Y(mciter,i+1);
            Z(mciter,i+1)= 0.6*Z(mciter,i)+ Z(mciter,i+1);
        end
    end    
end

data = [mat2cell(X(:,1:nsamples),ones(1,mciter),nsamples)';...
    mat2cell(Y(:,1:nsamples),ones(1,mciter),nsamples)';...
    mat2cell(Z(:,1:nsamples),ones(1,mciter),nsamples)'];

save(OUT_FILENAME, 'data');