ft_dipolesimulation.m

function [data] = ft_dipolesimulation(cfg)

% FT_DIPOLESIMULATION simulates channel-level time-series data that consists of the
% the spatial distribution of the the field or potential of one or multiple dipoles.
%
% Use as
%   data = ft_dipolesimulation(cfg)
% which will return a raw data structure that resembles the output of
% FT_PREPROCESSING.
%
% The dipoles position and orientation have to be specified with
%   cfg.sourcemodel.pos        = [Rx Ry Rz] (size Nx3)
%   cfg.sourcemodel.mom        = [Qx Qy Qz] (size 3xN)
%   cfg.sourcemodel.unit       = string, can be 'mm', 'cm', 'm' (default is automatic)
%
% The timecourse of the dipole activity is given as a cell-array with one
% dipole signal per trial
%   cfg.sourcemodel.signal     = cell-array with one dipole signal per trial
% or by specifying the parameters of a sine-wave signal
%   cfg.sourcemodel.frequency  =   in Hz
%   cfg.sourcemodel.phase      =   in radians
%   cfg.sourcemodel.amplitude  =   per dipole
%
% The number of trials and the time axes of the trials can be specified by
%   cfg.fsample     = simulated sample frequency (default = 1000)
%   cfg.trllen      = length of simulated trials in seconds (default = 1)
%   cfg.numtrl      = number of simulated trials (default = 10)
%   cfg.baseline    = number (default = 0.3)
% or by
%   cfg.time        = cell-array with one time axis per trial, for example obtained from an existing dataset
%
% Random white noise can be added to the data in each trial, either by
% specifying an absolute or a relative noise level
%   cfg.relnoise    = add noise with level relative to data signal
%   cfg.absnoise    = add noise with absolute level
%   cfg.randomseed  = 'yes' or a number or vector with the seed value (default = 'yes')
%
% Optional input arguments are
%   cfg.channel    = Nx1 cell-array with selection of channels (default = 'all'),
%                    see FT_CHANNELSELECTION for details
%   cfg.dipoleunit = units for dipole amplitude (default nA*m)
%   cfg.chanunit   = Nx1 cell-array with units for the channel data
%
% Optionally, you can modify the leadfields by reducing the rank, i.e. remove the weakest orientation
%   cfg.reducerank    = 'no', or number (default = 3 for EEG, 2 for MEG)
%   cfg.backproject   = 'yes' or 'no',  determines when reducerank is applied whether the
%                       lower rank leadfield is projected back onto the original linear
%                       subspace, or not (default = 'yes')
%
% The volume conduction model of the head should be specified as
%   cfg.headmodel     = structure with volume conduction model, see FT_PREPARE_HEADMODEL
%
% The EEG or MEG sensor positions should be specified as
%   cfg.elec          = structure with electrode positions or filename, see FT_READ_SENS
%   cfg.grad          = structure with gradiometer definition or filename, see FT_READ_SENS
%
% See also FT_SOURCEANALYSIS, FT_DIPOLEFITTING, FT_TIMELOCKSIMULATION,
% FT_FREQSIMULATION, FT_CONNECTIVITYSIMULATION

% Undocumented local options
% cfg.feedback
% cfg.previous
% cfg.version

% Copyright (C) 2004, 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$

% 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 if the input cfg is valid for this function
cfg = ft_checkconfig(cfg, 'forbidden',  {'channels'}); % prevent accidental typos, see issue 1729
cfg = ft_checkconfig(cfg, 'renamed',    {'elecfile', 'elec'});
cfg = ft_checkconfig(cfg, 'renamed',    {'gradfile', 'grad'});
cfg = ft_checkconfig(cfg, 'renamed',    {'optofile', 'opto'});
cfg = ft_checkconfig(cfg, 'renamed',    {'hdmfile', 'headmodel'});
cfg = ft_checkconfig(cfg, 'renamed',    {'vol',     'headmodel'});
cfg = ft_checkconfig(cfg, 'renamed',    {'dip',    'sourcemodel'});

% for consistency with FT_TIMELOCKSIMULUATION and FT_FREQSIMULATION
cfg = ft_checkconfig(cfg, 'createsubcfg', 'sourcemodel');
cfg = ft_checkconfig(cfg, 'renamed', {'ntrials', 'numtrl'});
cfg = ft_checkconfig(cfg, 'renamed', {'triallength', 'trllen'});

% set the defaults
cfg.sourcemodel         = ft_getopt(cfg, 'sourcemodel', []);
cfg.sourcemodel.pos     = ft_getopt(cfg.sourcemodel, 'pos', [-5 0 15]);
cfg.sourcemodel.mom     = ft_getopt(cfg.sourcemodel, 'mom', [1 0 0]');
cfg.sourcemodel.time    = ft_getopt(cfg.sourcemodel, 'time', {});
cfg.sourcemodel.signal  = ft_getopt(cfg.sourcemodel, 'signal', {});
cfg.fsample             = ft_getopt(cfg, 'fsample', 250);
cfg.relnoise            = ft_getopt(cfg, 'relnoise', 0);
cfg.absnoise            = ft_getopt(cfg, 'absnoise', 0);
cfg.feedback            = ft_getopt(cfg, 'feedback', 'text');
cfg.channel             = ft_getopt(cfg, 'channel',  'all');
cfg.dipoleunit          = ft_getopt(cfg, 'dipoleunit', 'nA*m');
cfg.chanunit            = ft_getopt(cfg, 'chanunit', {});

% collect and preprocess the electrodes/gradiometer and head model
% this will also update cfg.channel to match the electrodes/gradiometers
[headmodel, sens, cfg] = prepare_headmodel(cfg, []);

% construct the low-level options for the leadfield computation as key-value pairs, these are passed to FT_COMPUTE_LEADFIELD
leadfieldopt = {};
leadfieldopt = ft_setopt(leadfieldopt, 'reducerank',     ft_getopt(cfg, 'reducerank'));
leadfieldopt = ft_setopt(leadfieldopt, 'backproject',    ft_getopt(cfg, 'backproject'));
leadfieldopt = ft_setopt(leadfieldopt, 'normalize',      ft_getopt(cfg, 'normalize'));
leadfieldopt = ft_setopt(leadfieldopt, 'normalizeparam', ft_getopt(cfg, 'normalizeparam'));
leadfieldopt = ft_setopt(leadfieldopt, 'weight',         ft_getopt(cfg, 'weight'));

cfg.sourcemodel = fixdipole(cfg.sourcemodel);
Ndipoles = size(cfg.sourcemodel.pos,1);

% in case no time or signal was given, set some additional defaults
if ~isempty(cfg.sourcemodel.time) && ~isempty(cfg.sourcemodel.signal)
  assert(length(cfg.sourcemodel.signal)==length(cfg.sourcemodel.time)); % these must match
  cfg.numtrl    = length(cfg.sourcemodel.time);
  cfg.fsample   = 1/mean(diff(cfg.sourcemodel.time{1}));  % determine from time-axis
  cfg.trllen    = length(cfg.sourcemodel.time{1})/cfg.fsample;
  cfg.baseline  = -cfg.sourcemodel.time{1}(1);
elseif ~isempty(cfg.sourcemodel.time)
  cfg.numtrl    = length(cfg.sourcemodel.time);
  cfg.fsample   = 1/mean(diff(cfg.sourcemodel.time{1}));  % determine from time-axis
  cfg.trllen    = length(cfg.sourcemodel.time{1})/cfg.fsample;
  cfg.baseline  = -cfg.sourcemodel.time{1}(1);
elseif ~isempty(cfg.sourcemodel.signal)
  cfg.numtrl    = length(cfg.sourcemodel.signal);
  cfg.fsample   = ft_getopt(cfg, 'fsample', 1000);
  cfg.trllen    = length(cfg.sourcemodel.signal{1})/cfg.fsample;
  cfg.baseline  = ft_getopt(cfg, 'baseline', 0);
else
  cfg.numtrl    = ft_getopt(cfg, 'numtrl', 10);
  cfg.fsample   = ft_getopt(cfg, 'fsample', 1000);
  cfg.trllen    = ft_getopt(cfg, 'trllen', 1);
  cfg.baseline  = ft_getopt(cfg, 'baseline', 0);
end

% no signal was given, set some additional defaults
if isempty(cfg.sourcemodel.signal)
  cfg.sourcemodel.frequency   = ft_getopt(cfg, 'frequency', ones(Ndipoles,1)*10);
  cfg.sourcemodel.phase       = ft_getopt(cfg, 'phase', zeros(Ndipoles,1));
  cfg.sourcemodel.amplitude   = ft_getopt(cfg, 'amplitude', ones(Ndipoles,1));
end

if isfield(cfg.sourcemodel, 'frequency')
  % this should be a column vector
  cfg.sourcemodel.frequency = cfg.sourcemodel.frequency(:);
end

if isfield(cfg.sourcemodel, 'phase')
  % this should be a column vector
  cfg.sourcemodel.phase = cfg.sourcemodel.phase(:);
end

if ~isempty(cfg.sourcemodel.time)
  % use the user-supplied time vectors
  diptime = cfg.sourcemodel.time;
else
  % construct a time axis for every trial
  nsample = round(cfg.trllen*cfg.fsample);
  diptime = cell(1, cfg.numtrl);
  for iTr = 1:cfg.numtrl
    diptime{iTr} = (((1:nsample)-1)/cfg.fsample) - cfg.baseline;
  end
end

if ~isempty(cfg.sourcemodel.signal)
  % use the user-supplied signal for the dipoles
  dipsignal = cfg.sourcemodel.signal;
else
  dipsignal = cell(1, cfg.numtrl);
  for iTr = 1:cfg.numtrl
    % compute a cosine signal with the desired frequency, phase and amplitude for each dipole
    for i=1:Ndipoles
      dipsignal{iTr}(i,:) = cos(cfg.sourcemodel.frequency(i)*diptime{iTr}*2*pi + cfg.sourcemodel.phase(i)) * cfg.sourcemodel.amplitude(i);
    end
  end
end

dippos = cfg.sourcemodel.pos;
dipmom = cfg.sourcemodel.mom;

if ~iscell(dipmom)
  dipmom = {dipmom};
end

if ~iscell(dippos)
  dippos = {dippos};
end

if length(dippos)==1
  dippos = repmat(dippos, 1, cfg.numtrl);
elseif length(dippos)~=cfg.numtrl
  ft_error('incorrect number of trials specified in the dipole position');
end

if length(dipmom)==1
  dipmom = repmat(dipmom, 1, cfg.numtrl);
elseif length(dipmom)~=cfg.numtrl
  ft_error('incorrect number of trials specified in the dipole moment');
end

data.time  = diptime;
data.trial = {};

ft_progress('init', cfg.feedback, 'computing data data');
for trial=1:cfg.numtrl
  ft_progress(trial/cfg.numtrl, 'computing data data for trial %d\n', trial);
  if numel(cfg.chanunit) == numel(cfg.channel)
    lf = ft_compute_leadfield(dippos{trial}, sens, headmodel, 'dipoleunit', cfg.dipoleunit, 'chanunit', cfg.chanunit, leadfieldopt{:});
  else
    lf = ft_compute_leadfield(dippos{trial}, sens, headmodel, leadfieldopt{:});
  end
  nsamples = size(dipsignal{trial},2);
  nchannels = size(lf,1);
  data.trial{trial} = zeros(nchannels,nsamples);
  for i = 1:3
    data.trial{trial} = data.trial{trial} + ...
      lf(:,i:3:end) * (repmat(dipmom{trial}(i:3:end),1,nsamples) .* dipsignal{trial});
  end
end
ft_progress('close');

if ft_senstype(sens, 'meg')
  data.grad = sens;
elseif ft_senstype(sens, 'eeg')
  data.elec = sens;
end

% determine RMS value of data data
ss = 0;
sc = 0;
for trial=1:cfg.numtrl
  ss = ss + sum(data.trial{trial}(:).^2);
  sc = sc + length(data.trial{trial}(:));
end
rms = sqrt(ss/sc);
ft_info('RMS value of data data is %g\n', rms);

% add noise to the data data
for trial=1:cfg.numtrl
  relnoise = randn(size(data.trial{trial})) * cfg.relnoise * rms;
  absnoise = randn(size(data.trial{trial})) * cfg.absnoise;
  data.trial{trial} = data.trial{trial} + relnoise + absnoise;
end

data.fsample = cfg.fsample;
data.label   = sens.label;

% 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