function [simulated] = ft_connectivitysimulation(cfg)

% FT_CONNECTIVITYSIMULATION simulates channel-level time-series data with a

% specified connectivity structure. This function returns an output data

% structure that resembles the output of FT_PREPROCESSING.

%

% Use as

% [data] = ft_connectivitysimulation(cfg)

% which will return a raw data structure that resembles the output of

% FT_PREPROCESSING.

%

% The configuration structure should contain

% cfg.method = string, can be 'linear_mix', 'mvnrnd', 'ar', 'ar_reverse' (see below)

% cfg.nsignal = scalar, number of signals

% cfg.ntrials = scalar, number of trials

% cfg.triallength = in seconds

% cfg.fsample = in Hz

%

% Method 'linear_mix' implements a linear mixing with optional time shifts

% where the number of unobserved signals can be different from the number

% of observed signals

%

% Required configuration options:

% cfg.mix = matrix, [nsignal x number of unobserved signals]

% specifying the mixing from the unobserved signals to

% the observed signals, or

% = matrix, [nsignal x number of unobserved signals x number of

% samples] specifying the mixing from the

% unobserved signals to the observed signals which

% changes as a function of time within the trial

% = cell-arry, [1 x ntrials] with each cell a matrix as

% specified above, when a trial-specific mixing is

% required

% cfg.delay = matrix, [nsignal x number of unobserved signals]

% specifying the time shift (in samples) between the

% unobserved signals and the observed signals

%

% Optional configuration options

% cfg.bpfilter = 'yes' (or 'no')

% cfg.bpfreq = [bplow bphigh] (default: [15 25])

% cfg.demean = 'yes' (or 'no')

% cfg.baselinewindow = [begin end] in seconds, the default is the complete trial

% cfg.absnoise = scalar (default: 1), specifying the standard deviation of

% white noise superimposed on top of the simulated signals

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

%

% Method 'mvnrnd' implements a linear mixing with optional timeshifts in

% where the number of unobserved signals is equal to the number of observed

% signals. This method used the MATLAB function mvnrnd. The implementation

% is a bit ad-hoc and experimental, so users are discouraged to apply it.

% The time shift occurs only after the linear mixing, so the effect of the

% parameters on the simulation is not really clear. This method will be

% disabled in the future.

%

% Required configuration options

% cfg.covmat = covariance matrix between the signals

% cfg.delay = delay vector between the signals in samples

%

% Optional configuration options

% cfg.bpfilter = 'yes' (or 'no')

% cfg.bpfreq = [bplow bphigh] (default: [15 25])

% cfg.demean = 'yes' (or 'no')

% cfg.baselinewindow = [begin end] in seconds, the default is the complete trial

% cfg.absnoise = scalar (default: 1), specifying the standard

% deviation of white noise superimposed on top

% of the simulated signals

%

% Method 'ar' implements a multivariate autoregressive model to generate

% the data.

%

% Required configuration options

% cfg.params = matrix, [nsignal x nsignal x number of lags] specifying the

% autoregressive coefficient parameters. A non-zero

% element at cfg.params(i,j,k) means a

% directional influence from signal j onto

% signal i (at lag k).

% cfg.noisecov = matrix, [nsignal x nsignal] specifying the covariance

% matrix of the innovation process

%

% Method 'ar_reverse' implements a multivariate autoregressive

% autoregressive model to generate the data, where the model coefficients

% are reverse-computed, based on the interaction pattern specified.

%

% Required configuration options

% cfg.coupling = nxn matrix, specifying coupling strength, rows causing

% column

% cfg.delay = nxn matrix, specifying the delay, in seconds, from one

% signal's spectral component to the other signal, rows

% causing column

% cfg.ampl = nxn matrix, specifying the amplitude

% cfg.bpfreq = nxnx2 matrix, specifying the lower and upper frequencies

% of the bands that are transmitted, rows causing column

%

% The generated signals will have a spectrum that is 1/f + additional

% band-limited components, as specified in the configuration.

%

% See also FT_FREQSIMULATION, FT_DIPOLESIMULATION, FT_SPIKESIMULATION,

% FT_CONNECTIVITYANALYSIS

% Copyright (C) 2009-2015, Donders Institute for Brain, Cognition and Behaviour

%

% 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$

% 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

return

end

% check input configuration for the generally applicable options

cfg = ft_checkconfig(cfg, 'required', {'nsignal' 'ntrials' 'triallength' 'fsample' 'method'});

cfg = ft_checkconfig(cfg, 'rename', {'blc', 'demean'});

% method specific defaults

switch cfg.method

case {'ar'}

cfg.absnoise = ft_getopt(cfg, 'absnoise', zeros(cfg.nsignal,1));

cfg = ft_checkconfig(cfg, 'required', {'params' 'noisecov'});

case {'linear_mix'}

cfg.bpfilter = ft_getopt(cfg, 'bpfilter', 'yes');

cfg.bpfreq = ft_getopt(cfg, 'bpfreq', [15 25]);

cfg.bpfilttype = ft_getopt(cfg, 'bpfilttype', 'firws');

cfg.demean = ft_getopt(cfg, 'demean', 'yes');

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

cfg = ft_checkconfig(cfg, 'required', {'mix' 'delay'});

case {'mvnrnd'}

cfg.bpfilter = ft_getopt(cfg, 'bpfilter', 'yes');

cfg.bpfreq = ft_getopt(cfg, 'bpfreq', [15 25]);

cfg.bpfilttype = ft_getopt(cfg, 'bpfilttype', 'firws');

cfg.demean = ft_getopt(cfg, 'demean', 'yes');

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

cfg = ft_checkconfig(cfg, 'required', {'covmat' 'delay'});

case {'ar_reverse'}

% reverse engineered high order ar-model

cfg = ft_checkconfig(cfg, 'required', {'coupling' 'delay' 'ampl' 'bpfreq'});

otherwise

end

trial = cell(1, cfg.ntrials);

time = cell(1, cfg.ntrials);

nsmp = round(cfg.triallength*cfg.fsample);

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

% create the labels

label = cell(cfg.nsignal,1);

for k = 1:cfg.nsignal

label{k,1} = ['signal',num2str(k, '%03d')];

end

switch cfg.method

case {'ar'}

nlag = size(cfg.params,3);

nsignal = cfg.nsignal;

params = zeros(nlag*nsignal, nsignal);

for k = 1:nlag

%params(((k-1)*nsignal+1):k*nsignal,:) = cfg.params(:,:,k);

params(((k-1)*nsignal+1):k*nsignal,:) = cfg.params(:,:,k)';

% Use the transposition to make the implementation consistent with what

% comes out of ft_mvaranalysis. The transposition is introduced on May

% 13, 2011. This swaps the directional influence for existing scripts.

end

for k = 1:cfg.ntrials

tmp = zeros(nsignal, nsmp+ceil(nlag*1.05));

noise = mvnrnd(zeros(nsignal,1), cfg.noisecov, ceil(nsmp+nlag*1.05))';

state0 = zeros(nsignal*nlag, 1);

for m = 1:nlag

indx = ((m-1)*nsignal+1):m*nsignal;

state0(indx) = params(indx,:)'*noise(:,m);

end

tmp(:,1:nlag) = flip(reshape(state0, [nsignal nlag]),2);

for m = (nlag+1):(nsmp+ceil(nlag*1.05))

state0 = reshape(flip(tmp(:,(m-nlag):(m-1)),2), [nlag*nsignal 1]);

tmp(:, m) = params'*state0 + noise(:,m);

end

trial{k} = tmp(:,(ceil(nlag*1.05)+1):end);

if any(cfg.absnoise>0)

trial{k} = trial{k} + diag(cfg.absnoise)*randn(size(trial{k}));

end

time{k} = tim;

end

% create the output data

simulated = [];

simulated.trial = trial;

simulated.time = time;

simulated.fsample = cfg.fsample;

simulated.label = label;

case {'linear_mix'}

fltpad = 50; %hard coded to avoid filtering artifacts

delay = cfg.delay;

delay = delay - min(delay(:)); %make explicitly >= 0

maxdelay = max(delay(:));

if iscell(cfg.mix)

%each trial has different mix

mix = cfg.mix;

else

%make cell-array out of mix

tmpmix = cfg.mix;

mix = cell(1,cfg.ntrials);

for tr = 1:cfg.ntrials

mix{1,tr} = tmpmix;

end

end

nmixsignal = size(mix{1}, 2); %number of "mixing signals"

nsignal = size(mix{1}, 1);

if numel(size(mix{1}))==2

%mix is static, no function of time

for tr = 1:cfg.ntrials

mix{tr} = mix{tr}(:,:,ones(1,nsmp+maxdelay));

end

elseif numel(size(mix{1}))==3 && size(mix{1},3)==nsmp

%mix changes with time

for tr = 1:cfg.ntrials

mix{tr} = cat(3,mix{tr},mix{tr}(:,:,nsmp*ones(1,maxdelay)));

end

%FIXME think about this

%due to the delay the mix cannot be defined instantaneously with respect to all signals

end

for tr = 1:cfg.ntrials

mixsignal = randn(nmixsignal, nsmp + 2*fltpad + maxdelay);

if nmixsignal==size(cfg.bpfreq,1)

for sg = 1:nmixsignal

tmpcfg = cfg;

tmpcfg.bpfreq = cfg.bpfreq(sg,:);

newmixsignal(sg,:) = preproc(mixsignal(sg,:), label, offset2time(-fltpad, cfg.fsample, size(mixsignal,2)), tmpcfg, fltpad, fltpad);

end

else

% it can be done with a single set of cfg settings to preproc

newmixsignal = preproc(mixsignal, label, offset2time(-fltpad, cfg.fsample, size(mixsignal,2)), cfg, fltpad, fltpad);

end

tmp = zeros(cfg.nsignal, nsmp);

for i=1:cfg.nsignal

for j=1:nmixsignal

begsmp = 1 + delay(i,j);

endsmp = nsmp + delay(i,j);

tmpmix = reshape(mix{tr}(i,j,:),[1 nsmp+maxdelay]) .* newmixsignal(j,:);

tmp(i,:) = tmp(i,:) + tmpmix(begsmp:endsmp);

end

end

trial{tr} = tmp;

% add some noise

trial{tr} = ft_preproc_baselinecorrect(trial{tr} + cfg.absnoise*randn(size(trial{tr})));

% define time axis for this trial

time{tr} = tim;

end

% create the output data

simulated = [];

simulated.trial = trial;

simulated.time = time;

simulated.fsample = cfg.fsample;

simulated.label = label;

case {'mvnrnd'}

fltpad = 100; % hard coded

shift = max(cfg.delay(:,1)) - cfg.delay(:,1);

for k = 1:cfg.ntrials

% create the multivariate time series plus some padding

tmp = mvnrnd(zeros(1,cfg.nsignal), cfg.covmat, nsmp+2*fltpad+max(shift))';

% add the delays

newtmp = zeros(cfg.nsignal, nsmp+2*fltpad);

for kk = 1:cfg.nsignal

begsmp = + shift(kk) + 1;

endsmp = nsmp + 2*fltpad + shift(kk);

newtmp(kk,:) = ft_preproc_baselinecorrect(tmp(kk,begsmp:endsmp));

end

% apply preproc

newtmp = preproc(newtmp, label, offset2time(-fltpad, cfg.fsample, size(newtmp,2)), cfg, fltpad, fltpad);

trial{k} = newtmp;

% add some noise

trial{k} = ft_preproc_baselinecorrect(trial{k} + cfg.absnoise*randn(size(trial{k})));

% define time axis for this trial

time{k} = tim;

end

% create the output data

simulated = [];

simulated.trial = trial;

simulated.time = time;

simulated.fsample = cfg.fsample;

simulated.label = label;

case 'ar_reverse'

% generate a spectral transfer matrix, and a cross-spectral matrix

% according to the specifications

% predefine some variables

fstep = 1/5;

fs = cfg.fsample;

Nyq = fs./2;

foi = (0:fstep:Nyq);

omega = foi./fs;

n = numel(foi);

% local renaming

nsignal = cfg.nsignal;

fband = cfg.bpfreq;

coupling = cfg.coupling;

ampl = cfg.ampl;

delay = cfg.delay;

% create a 1/f spectrum

slope = 0.5;

oneoverf = sqrt(max(omega(2)./10,omega).^-slope); % takes sqrt for amplitude

oneoverf = oneoverf./oneoverf(1);

%oneoverf(1) = 0;

%z = firws_filter(5.*fs, fs, Nyq./1.01);

%z = z(1:numel(foi)); %.*exp(-1i.*pi.*foi.*rand(1)./100);

%oneoverf = z.*oneoverf;

% convert into indices

findx = fband;

for k = 1:numel(fband)

if isfinite(fband(k))

findx(k) = nearest(foi, fband(k));

end

end

% allocate some memory

mask = false(nsignal, nsignal, n);

krn = zeros(size(mask));

phi = zeros(size(krn));

dat = zeros(size(krn));

coupling_ampl = zeros(size(krn));

for k = 1:nsignal

for m = 1:nsignal

if all(isfinite(squeeze(findx(k,m,:))))

mask(k,m,findx(k,m,1):findx(k,m,2)) = true;

end

krn(k,m,mask(k,m,:)) = hanning(sum(mask(k,m,:)))';

phi(k,m,:) = 2.*pi.*delay(k,m).*foi;

%phi(k,m,:) = phi(k,m,:).*mask(k,m,:);

%phi(k,m,mask(k,m,:)) = phi(k,m,mask(k,m,:))-mean(phi(k,m,mask(k,m,:)));

if all(isfinite(squeeze(findx(k,m,:))))

phi(k,m,1:findx(k,m,1)) = phi(k,m,findx(k,m,1));

phi(k,m,findx(k,m,2):end) = phi(k,m,findx(k,m,2));

phi(k,m,:) = phi(k,m,:)-mean(phi(k,m,:));

end

coupling_ampl(k,m,:) = coupling(k,m).*krn(k,m,:);

end

end

% this matrix contains the intrinsic amplitude spectra on the diagonal

for k = 1:nsignal

if all(isfinite(squeeze(fband(k,k,:))))

z = firws_filter((1/fstep).*fs, fs, [fband(k,k,1) fband(k,k,2)]);

z = z(1:numel(foi)); %.*exp(-1i.*pi.*foi.*rand(1)./100);

z = z.*ampl(k,k);

plateau = nearest(foi,fband(k,k,1)):nearest(foi,fband(k,k,2));

oneoverf(plateau) = mean(abs(oneoverf(plateau)));

dat(k,k,:) = -(abs(oneoverf)+abs(z)).*exp(1i.*(angle(z)+angle(oneoverf)));

else

dat(k,k,:) = oneoverf;

end

end

% now we can create a spectral transfer matrix

tf = zeros(nsignal,nsignal,n)+1i.*zeros(nsignal,nsignal,n);

for k = 1:nsignal

for m = 1:nsignal

if k~=m && all(isfinite(squeeze(fband(k,m,:))))

z = firws_filter((1/fstep).*fs, fs, [fband(k,m,1) fband(k,m,2)]);

z = z(1:numel(foi));

tf(m,k,:) = coupling(k,m).*exp(-1i.*phi(k,m,:)).*shiftdim(z,-1); % deliberate index swap!

elseif k==m

tf(k,m,:) = dat(k,m,:);

end

end

end

% create the cross spectral matrix

c = zeros(size(tf));

for k = 1:n

c(:,:,k) = tf(:,:,k)*tf(:,:,k)'; % assume noise to be I, i.e. the tf to swallow the amplitudes

end

% scale the Nyquist and DC bins

c(:,:,1) = real(c(:,:,1)./2);

c(:,:,end) = real(c(:,:,end)./2);

% create a freq-structure

freq = [];

freq.crsspctrm = c;

freq.label = label;

freq.freq = foi;

freq.dimord = 'chan_chan_freq';

% estimate the transfer-matrix non-parametrically

tmpcfg = [];

tmpcfg.method = 'transfer';

tmpcfg.granger.stabilityfix = true;

t = ft_connectivityanalysis(tmpcfg, freq);

% estimate the ar-model coefficients

a = transfer2coeffs(t.transfer,t.freq);

% recursively call this function to generate the data, this is

% somewhate tricky with respect to keeping the provenance info. Here,

% it is solved by removing from the cfg the original user-specified

% fields

cfgorig = cfg;

cfg = removefields(cfgorig, {'coupling' 'ampl' 'delay' 'bpfreq'});

cfg.method = 'ar';

cfg.params = a;

cfg.noisecov = diag(diag(t.noisecov.*cfg.fsample./2));

simulated = ft_connectivitysimulation(cfg);

cfg.previous = keepfields(cfgorig, {'coupling' 'ampl' 'delay' 'bpfreq'});

otherwise

ft_error('unknown method');

end

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

ft_postamble debug

ft_postamble randomseed

ft_postamble provenance simulated

ft_postamble history simulated

ft_postamble savevar simulated

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

% SUBFUNCTION

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

function z = firws_filter(N, Fs, Fbp)

switch numel(Fbp)

case 1

[dum, B] = ft_preproc_lowpassfilter(randn(1,N), Fs, Fbp, [], 'firws', 'onepass-minphase');

z = fft(B, N);

case 2

[dum, B] = ft_preproc_bandpassfilter(randn(1,N), Fs, Fbp, [], 'firws', 'onepass-minphase');

z = fft(B, N);

end