function [data] = ft_steadystatesimulation(cfg)

% FT_STEADYSTATESIMULATION creates a simulated EEG/MEG dataset. This function

% allows to simulate the effect of several independent stimulus trains. These can

% be presented as a periodic sequence, or as single (or few) transient stimuli.

% This function creates a single block of data. You can call it repeatedly and use

% FT_APPENDDATA to combine different blocks.

%

% Use as

% data = ft_steadystatesimulation(cfg)

% where cfg is a configuration structure that should contain

% cfg.fsample = scalar, sampling frequency in Hz (default = 512)

% cfg.duration = scalar, trial length in seconds (default = 4.56)

% cfg.baseline = scalar, baseline length in seconds (default = 0)

% cfg.ntrials = integer N, number of trials (default = 320)

% cfg.iti = scalar, inter-trial interval in seconds (default = 1)

% cfg.randomseed = 'yes' or a number or vector with the seed value (default = 'yes')

%

% Each trial can contain multiple nested experimental manipulations

% cfg.level1.condition = scalar, or vector of length L1 (default = 1)

% cfg.level1.gain = scalar, or vector of length L1 (default = 1)

% cfg.level2.condition = scalar, or vector of length L2 (default = 1)

% cfg.level2.gain = scalar, or vector of length L2 (default = 1)

% cfg.level3.condition = scalar, or vector of length L3 (default = 1)

% cfg.level3.gain = scalar, or vector of length L3 (default = 1)

% If you don't need level 2 and up, specify the condition and gain as empty.

% Idem for level 3 and up.

%

% Stimuli are created at the lowest experimental level, and are modulated according to the product of the gain of all levels.

% Each trial can contain one or multiple stimuli.

% The behavior of each stimuli is specified with

% cfg.stimulus1.mode = 'periodic', 'transient' or 'off' (default = 'periodic')

% cfg.stimulus2.mode = 'periodic', 'transient' or 'off' (default = 'transient')

%

% If the stimulus is periodic (below as example for stimulus1), the following options apply

% cfg.stimulus1.number = does not apply for periodic stimuli

% cfg.stimulus1.onset = in seconds, first stimulus relative to the start of the trial (default = 0)

% cfg.stimulus1.onsetjitter = in seconds, max jitter that is added to the onset (default = 0)

% cfg.stimulus1.isi = in seconds, i.e. for 10Hz you would specify 0.1 seconds as the interstimulus interval (default = 0.1176)

% cfg.stimulus1.isijitter = in seconds, max jitter relative to the previous stimulus (default = 0)

% cfg.stimulus2.condition = does not apply for periodic stimuli

% cfg.stimulus2.gain = does not apply for periodic stimuli

% cfg.stimulus1.kernelshape = 'sine'

% cfg.stimulus1.kernelduration = in seconds (default = isi)

%

% If the stimulus is transient (below as example for stimulus2), the following options apply

% cfg.stimulus2.number = scalar M, how many transients are to be presented per trial (default = 4)

% cfg.stimulus2.onset = in seconds, first stimulus relative to the start of the trial (default = 0.7)

% cfg.stimulus2.onsetjitter = in seconds, max jitter that is added to the onset (default = 0.2)

% cfg.stimulus2.isi = in seconds as the interstimulus interval (default = 0.7)

% cfg.stimulus2.isijitter = in seconds, max jitter relative to the previous stimulus (default = 0.2)

% cfg.stimulus2.condition = 1xM vector with condition codes for each transient within a trial (default = [1 1 2 2])

% cfg.stimulus2.gain = 1xM vector with gain for each condition for each transient within a trial (default = [1 1 1 1])

% cfg.stimulus2.kernelshape = 'hanning'

% cfg.stimulus2.kernelduration = in seconds (default = 0.75*isi)

%

% RANDOMIZATIONS:

% - The onsetjitter is randomized between 0 and the value given, and is always added to the onset.

% - The isijitter is randomized between 0 and the value given, and is always added to the interstimulus interval (isi).

% - For periodic stimuli, which are constant within a trial, the condition code and gain are shuffled over all trials.

% - For transient stimuli, the condition code and gain are shuffled within each trial.

%

% Using the default settings, we model a peripherally presented flickering stimulus

% that appears at different excentricities together with a centrally presented

% transient stimulus that appears 4x per trial. To simulate the experiment described

% at , you have to call this 4 times with a different cfg.configuration and

% cfg.gain to model the task load and use FT_APPENDDATA to concatenate the trials. In

% this case cfg.condition models the factor "task load" (2 levels, low and high),

% cfg.stimulus1.condition models the factor "excentricity" (4 levels), and

% cfg.stimulation2.condition models the factor "stimulus type" (2 levels, non-target

% or target).

%

% See also FT_DIPOLESIMULATION, FT_TIMELOCKSIMULATION, FT_FREQSIMULATION,

% FT_CONNECTIVITYSIMULATION, FT_APPENDDATA

% Copyright (C) 2017, Stefan Wiens, Malina Szychowska, Robert Oostenveld

%

% This file is part of FieldTrip, see http://www.fieldtriptoolbox.org

% for the documentation and details.

%

% FieldTrip is free software: you can redistribute it and/or modify

% it under the terms of the GNU General Public License as published by

% the Free Software Foundation, either version 3 of the License, or

% (at your option) any later version.

%

% FieldTrip is distributed in the hope that it will be useful,

% but WITHOUT ANY WARRANTY; without even the implied warranty of

% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the

% GNU General Public License for more details.

%

% You should have received a copy of the GNU General Public License

% along with FieldTrip. If not, see <http://www.gnu.org/licenses/>.

%

% $Id$

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% the initial part deals with parsing the input options and data

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% these are used by the ft_preamble/ft_postamble function and scripts

ft_revision = '$Id$';

ft_nargin = nargin;

ft_nargout = nargout;

% do the general setup of the function

ft_defaults

ft_preamble init

ft_preamble debug

ft_preamble provenance

ft_preamble randomseed

% the ft_abort variable is set to true or false in ft_preamble_init

if ft_abort

% do not continue function execution in case the outputfile is present and the user indicated to keep it

return

end

% get the options or set the defaults

cfg.level1 = ft_getopt(cfg, 'level1');

cfg.level2 = ft_getopt(cfg, 'level2');

cfg.level3 = ft_getopt(cfg, 'level3');

cfg.level1.condition = ft_getopt(cfg.level1, 'condition', 1);

cfg.level2.condition = ft_getopt(cfg.level2, 'condition', []);

cfg.level3.condition = ft_getopt(cfg.level3, 'condition', []);

cfg.level1.gain = ft_getopt(cfg.level1, 'gain', 1);

cfg.level2.gain = ft_getopt(cfg.level2, 'gain', []);

cfg.level3.gain = ft_getopt(cfg.level3, 'gain', []);

cfg.ntrials = ft_getopt(cfg, 'ntrials', 320);

cfg.fsample = ft_getopt(cfg, 'fsample', 512);

cfg.duration = ft_getopt(cfg, 'duration', 4.5);

cfg.baseline = ft_getopt(cfg, 'baseline', 0);

cfg.iti = ft_getopt(cfg, 'iti', 1);

cfg.stimulus1 = ft_getopt(cfg, 'stimulus1');

cfg.stimulus2 = ft_getopt(cfg, 'stimulus2');

cfg.stimulus1.mode = ft_getopt(cfg.stimulus1, 'mode', 'periodic');

cfg.stimulus1.number = ft_getopt(cfg.stimulus1, 'number', 0);

cfg.stimulus1.onset = ft_getopt(cfg.stimulus1, 'onset', 0);

cfg.stimulus1.onsetjitter = ft_getopt(cfg.stimulus1, 'onsetjitter', 0);

cfg.stimulus1.isi = ft_getopt(cfg.stimulus1, 'isi', 0.1176);

cfg.stimulus1.isijitter = ft_getopt(cfg.stimulus1, 'isijitter', 0);

cfg.stimulus1.condition = ft_getopt(cfg.stimulus1, 'condition');

cfg.stimulus1.gain = ft_getopt(cfg.stimulus1, 'gain');

cfg.stimulus1.kernelshape = ft_getopt(cfg.stimulus1, 'kernelshape', 'sine');

cfg.stimulus1.kernelduration = ft_getopt(cfg.stimulus1, 'kernelduration', cfg.stimulus1.isi);

cfg.stimulus2.mode = ft_getopt(cfg.stimulus2, 'mode', 'transient');

cfg.stimulus2.number = ft_getopt(cfg.stimulus2, 'number', 4);

cfg.stimulus2.onset = ft_getopt(cfg.stimulus2, 'onset', 0.7);

cfg.stimulus2.onsetjitter = ft_getopt(cfg.stimulus2, 'onsetjitter', 0);

cfg.stimulus2.isi = ft_getopt(cfg.stimulus2, 'isi', 0.7);

cfg.stimulus2.isijitter = ft_getopt(cfg.stimulus2, 'isijitter', 0);

cfg.stimulus2.condition = ft_getopt(cfg.stimulus2, 'condition', [1 2 3 4]);

cfg.stimulus2.gain = ft_getopt(cfg.stimulus2, 'gain', [1 1 1 1]);

cfg.stimulus2.kernelshape = ft_getopt(cfg.stimulus2, 'kernelshape', 'hanning');

cfg.stimulus2.kernelduration = ft_getopt(cfg.stimulus2, 'kernelduration', 0.75*cfg.stimulus2.isi);

ncondition1 = length(cfg.level1.condition);

ncondition2 = length(cfg.level2.condition);

ncondition3 = length(cfg.level3.condition);

switch cfg.stimulus1.mode

case 'off'

cfg.stimulus1.number = 0;

cfg.stimulus1.condition = [];

cfg.stimulus1.gain = [];

case 'transient'

assert(numel(cfg.stimulus1.condition)==cfg.stimulus1.number);

assert(numel(cfg.stimulus1.gain)==cfg.stimulus1.number);

end

switch cfg.stimulus2.mode

case 'off'

cfg.stimulus2.number = 0;

cfg.stimulus2.condition = [];

cfg.stimulus2.gain = [];

case 'transient'

assert(numel(cfg.stimulus2.condition)==cfg.stimulus2.number);

assert(numel(cfg.stimulus2.gain)==cfg.stimulus2.number);

end

% they cannot be the same, as that would confuse the trialinfo

assert(~isequal(cfg.stimulus1.mode, cfg.stimulus2.mode));

% determine which levels are being used

level1 = ncondition1>0;

level2 = ncondition2>0;

level3 = ncondition3>0;

% count the number of blocks and number of trials per block

nblock = max(ncondition1,1) * max(ncondition2,1) * max(ncondition3,1);

ntrial = cfg.ntrials / nblock;

assert(ntrial==round(ntrial), 'number of trials inconsistent');

% the time (followinf onset) and baseline are the same for each trial

baseline = (-cfg.baseline):(1/cfg.fsample):0;

baseline(end) = []; % remove the last sample

time = 0:(1/cfg.fsample):cfg.duration;

% construct a sine-wave kernel, feval fails on this one

sine = @(n) sin(2*pi*((1:n)-1)/n);

% construct the kernels for the response to each stimulus

if strcmp(cfg.stimulus1.kernelshape, 'sine')

kernel1 = sine(round(cfg.stimulus1.kernelduration*cfg.fsample));

else

kernel1 = feval(cfg.stimulus1.kernelshape, round(cfg.stimulus1.kernelduration*cfg.fsample));

end

if strcmp(cfg.stimulus2.kernelshape, 'sine')

kernel2 = sine(round(cfg.stimulus2.kernelduration*cfg.fsample));

else

kernel2 = feval(cfg.stimulus2.kernelshape, round(cfg.stimulus2.kernelduration*cfg.fsample));

end

% these should be row-vectors

kernel1 = kernel1(:)';

kernel2 = kernel2(:)';

% start with empty data

data = [];

% this keeps track of the trials

k = 1;

shuffle = randperm(ncondition3);

tmpcondition3 = cfg.level3.condition(shuffle);

tmpgain3 = cfg.level3.gain(shuffle);

for cond3=1:max(ncondition3, 1)

if level3

trialgain = tmpgain3(cond3);

else

trialgain = 1;

end

shuffle = randperm(ncondition2);

tmpcondition2 = cfg.level2.condition(shuffle);

tmpgain2 = cfg.level2.gain(shuffle);

for cond2=1:max(ncondition2, 1)

if level2

trialgain = trialgain * tmpgain2(cond2);

end

shuffle = randperm(ncondition1);

tmpcondition1 = cfg.level1.condition(shuffle);

tmpgain1 = cfg.level1.gain(shuffle);

for cond1=1:max(ncondition1, 1)

if level1

trialgain = trialgain * tmpgain1(cond1);

end

for trial=1:ntrial

trialinfo = k; % start with the trial number

varname = {'trial'};

if level1

trialinfo = [trialinfo tmpcondition1(cond1)];

varname(end+1) = {'l1_cond'};

end

if level2>0

trialinfo = [trialinfo tmpcondition2(cond2)];

varname(end+1) = {'l2_cond'};

end

if level3>0

trialinfo = [trialinfo tmpcondition3(cond3)];

varname(end+1) = {'l3_cond'};

end

dat = [];

stimcfg = cfg.stimulus1;

signal = zeros(size(time)); % the baseline will be added later

switch stimcfg.mode

case 'periodic'

sample = nearest(time, stimcfg.onset + rand(1)*stimcfg.onsetjitter);

for i=1:ceil((cfg.duration-stimcfg.onset)/stimcfg.isi)

sample = sample + round(cfg.fsample*(stimcfg.isi + stimcfg.isijitter*rand(1)));

if sample>length(signal)

% this periodic stimulus falls outside the trial

sample = nan;

else

signal(sample) = 1;

end

trialinfo = [trialinfo sample+length(baseline)];

varname(end+1) = {sprintf('p%d_sample', i)};

end % while

case 'transient'

shuffle = randperm(numel(stimcfg.gain));

sample = nearest(time, stimcfg.onset + rand(1)*stimcfg.onsetjitter);

transientgain = stimcfg.gain(shuffle);

transientcondition = stimcfg.condition(shuffle);

for i=1:stimcfg.number

sample = sample + round(cfg.fsample*(stimcfg.isi + stimcfg.isijitter*rand(1)));

if sample>length(signal)

warning('transient stimulus falls outside of trial');

sample = nan;

else

signal(sample) = transientgain(i);

end

trialinfo = [trialinfo transientcondition(i)];

varname(end+1) = {sprintf('t%d_cond', i)};

trialinfo = [trialinfo sample+length(baseline)];

varname(end+1) = {sprintf('t%d_sample', i)};

end % while

end % switch mode

dat(2,:) = [zeros(size(baseline)) signal*trialgain ];

dat(3,:) = [zeros(size(baseline)) erpconvolve(signal*trialgain, kernel1)];

stimcfg = cfg.stimulus2;

signal = zeros(size(time)); % the baseline will be added later

switch stimcfg.mode

case 'periodic'

sample = nearest(time, stimcfg.onset + rand(1)*stimcfg.onsetjitter);

for i=1:ceil((cfg.duration-stimcfg.onset)/stimcfg.isi)

sample = sample + round(cfg.fsample*(stimcfg.isi + stimcfg.isijitter*rand(1)));

if sample>length(signal)

% this periodic stimulus falls outside the trial

sample = nan;

else

signal(sample) = 1;

end

trialinfo = [trialinfo sample+length(baseline)];

varname(end+1) = {sprintf('p%d_sample', i)};

end % while

case 'transient'

shuffle = randperm(numel(stimcfg.gain));

sample = nearest(time, stimcfg.onset + rand(1)*stimcfg.onsetjitter);

transientgain = stimcfg.gain(shuffle);

transientcondition = stimcfg.condition(shuffle);

for i=1:stimcfg.number

sample = sample + round(cfg.fsample*(stimcfg.isi + stimcfg.isijitter*rand(1)));

if sample>length(signal)

warning('transient stimulus falls outside of trial');

sample = nan;

else

signal(sample) = transientgain(i);

end

trialinfo = [trialinfo transientcondition(i)];

varname(end+1) = {sprintf('t%d_cond', i)};

trialinfo = [trialinfo sample+length(baseline)];

varname(end+1) = {sprintf('t%d_sample', i)};

end % while

end % switch mode

dat(4,:) = [zeros(size(baseline)) signal*trialgain ];

dat(5,:) = [zeros(size(baseline)) erpconvolve(signal*trialgain, kernel2)];

% the first channel contains the sum of the two convolved signals

dat(1,:) = sum(dat([3 5],:), 1);

data.trial{k} = dat;

data.time{k} = [baseline time];

data.trialinfo(k,:) = trialinfo;

k = k+1;

end % for each trial in each block

end % for all conditions in level 1

end % for all conditions in level 2

end % for all conditions in level 3

data.trialinfo = array2table(data.trialinfo, 'VariableNames', varname);

data.fsample = cfg.fsample;

data.label = {

'combined'

'stickfunction1'

'convolved1'

'stickfunction2'

'convolved2'

};

% construct the sample information, start with consecutive segments

nsample = length(baseline) + length(time);

data.sampleinfo(:,1) = ((1:ntrial)-1)*nsample + 1;

data.sampleinfo(:,2) = ((1:ntrial) )*nsample;

% add the inter-trial-interval

data.sampleinfo(:,1) = data.sampleinfo(:,1) + transpose((1:ntrial)-1)*round(cfg.iti*cfg.fsample);

data.sampleinfo(:,2) = data.sampleinfo(:,2) + transpose((1:ntrial)-1)*round(cfg.iti*cfg.fsample);

% do the general cleanup and bookkeeping at the end of the function

ft_postamble debug

ft_postamble randomseed

ft_postamble provenance data

ft_postamble history data

ft_postamble savevar data

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% SUBFUNCTION

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function dat = erpconvolve(stickfunction, kernel)

dat = convn(stickfunction, kernel, 'full');

dat = dat(1:length(stickfunction));