Skip to content
Permalink
release
Switch branches/tags

Name already in use

A tag already exists with the provided branch name. Many Git commands accept both tag and branch names, so creating this branch may cause unexpected behavior. Are you sure you want to create this branch?
Go to file
 
 
Cannot retrieve contributors at this time
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
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