text:
function [data] = ft_checkdata(data, varargin)

% FT_CHECKDATA checks the input data of the main FieldTrip functions, e.g. whether the
% type of data strucure corresponds with the required data. If necessary and possible,
% this function will adjust the data structure to the input requirements (e.g. change
% dimord, average over trials, convert inside from index into logical).
%
% If the input data does NOT correspond to the requirements, this function will give a
% warning message and if applicable point the user to external documentation (link to
% website).
%
% Use as
%   [data] = ft_checkdata(data, ...)
%
% Optional input arguments should be specified as key-value pairs and can include
%   feedback           = yes, no
%   datatype           = raw, freq, timelock, comp, spike, source, mesh, dip, volume, segmentation, parcellation
%   dimord             = any combination of time, freq, chan, refchan, rpt, subj, chancmb, rpttap, pos
%   senstype           = ctf151, ctf275, ctf151_planar, ctf275_planar, neuromag122, neuromag306, bti148, bti248, bti248_planar, magnetometer, electrode
%   fsample            = sampling frequency to use to go from SPIKE to RAW representation
%   ismeg              = yes, no
%   iseeg              = yes, no
%   isnirs             = yes, no
%   hasunit            = yes, no
%   hascoordsys        = yes, no
%   haschantype        = yes, no
%   haschanunit        = yes, no
%   hassampleinfo      = yes, no, ifmakessense (applies to raw and timelock data)
%   hascumtapcnt       = yes, no (only applies to freq data)
%   hasdim             = yes, no
%   hasdof             = yes, no
%   hasbrain           = yes, no (only applies to segmentation)
%   insidestyle        = logical, index, can also be empty
%   cmbstyle           = sparse, sparsewithpow, full, fullfast, fourier (applies to covariance and cross-spectral density)
%   segmentationstyle  = indexed, probabilistic (only applies to segmentation)
%   parcellationstyle  = indexed, probabilistic (only applies to parcellation)
%   trialinfostyle     = matrix, table or empty
%
% For some options you can specify multiple values, e.g.
%   [data] = ft_checkdata(data, 'senstype', {'ctf151', 'ctf275'}), e.g. in megrealign
%   [data] = ft_checkdata(data, 'datatype', {'timelock', 'freq'}), e.g. in sourceanalysis
%
% See also FT_DATATYPE_XXX for each of the respective data types.

% Copyright (C) 2007-2021, Robert Oostenveld
% Copyright (C) 2010-2012, Martin Vinck
%
% 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$

% in case of an error this function could use dbstack for more detailled
% user feedback
%
% this function should replace/encapsulate
%   fixdimord
%   fixinside
%   fixprecision
%   fixvolume
%   fixpos
%   data2raw
%   raw2data
%   grid2transform
%   transform2grid
%   fourier2crsspctrm
%   freq2cumtapcnt
%   sensortype
%   time2offset
%   offset2time
%   fixsens -> this is kept a separate function because it should also be
%              called from other modules
%
% other potential uses for this function:
%   time -> offset in freqanalysis
%   average over trials
%   csd as matrix

% FIXME the following is difficult, if not impossible, to support without knowing the parameter
% FIXME it is presently (dec 2014) not being used anywhere in FT, so can be removed
%   hastrials          = yes, no

% check whether people are using deprecated options
sel = find(strcmp(varargin(1:2:end), 'hastrialdef'));
if ~isempty(sel)
  ft_warning('the option ''hastrialdef'' is deprecated, please use ''hassampleinfo'' instead');
  varargin{2*sel-1} = 'hassampleinfo';
end
sel = find(strcmp(varargin(1:2:end), 'inside'));
if ~isempty(sel)
  ft_warning('the option ''inside'' is deprecated, please use ''insidestyle'' instead');
  varargin{2*sel-1} = 'insidestyle';
end
sel = find(strcmp(varargin(1:2:end), 'cmbrepresentation'));
if ~isempty(sel)
  ft_warning('the option ''cmbrepresentation'' is deprecated, please use ''cmbstyle'' instead');
  varargin{2*sel-1} = 'cmbstyle';
end

% get the optional input arguments
feedback             = ft_getopt(varargin, 'feedback', 'no');
dtype                = ft_getopt(varargin, 'datatype'); % should not conflict with the ft_datatype function
dimord               = ft_getopt(varargin, 'dimord');
stype                = ft_getopt(varargin, 'senstype'); % senstype is a function name which should not be masked
ismeg                = ft_getopt(varargin, 'ismeg');
iseeg                = ft_getopt(varargin, 'iseeg');
isnirs               = ft_getopt(varargin, 'isnirs');
hastrials            = ft_getopt(varargin, 'hastrials');
hasunit              = ft_getopt(varargin, 'hasunit', 'no');
hascoordsys          = ft_getopt(varargin, 'hascoordsys', 'no');
haschantype          = ft_getopt(varargin, 'haschantype', 'no');
haschanunit          = ft_getopt(varargin, 'haschanunit', 'no');
hassampleinfo        = ft_getopt(varargin, 'hassampleinfo', 'ifmakessense');
hasdim               = ft_getopt(varargin, 'hasdim');
hascumtapcnt         = ft_getopt(varargin, 'hascumtapcnt');
hasdof               = ft_getopt(varargin, 'hasdof');
hasbrain             = ft_getopt(varargin, 'hasbrain');
cmbstyle             = ft_getopt(varargin, 'cmbstyle'); % sparse, sparsewithpow, full, fullfast, fourier
channelcmb           = ft_getopt(varargin, 'channelcmb');
insidestyle          = ft_getopt(varargin, 'insidestyle'); % logical, index
segmentationstyle    = ft_getopt(varargin, 'segmentationstyle'); % this will be passed on to the corresponding ft_datatype_xxx function
parcellationstyle    = ft_getopt(varargin, 'parcellationstyle'); % this will be passed on to the corresponding ft_datatype_xxx function
trialinfostyle       = ft_getopt(varargin, 'trialinfostyle');
fsample              = ft_getopt(varargin, 'fsample');

% determine the type of input data
israw           = ft_datatype(data, 'raw');
isfreq          = ft_datatype(data, 'freq');
istimelock      = ft_datatype(data, 'timelock');
iscomp          = ft_datatype(data, 'comp');
isspike         = ft_datatype(data, 'spike');
isvolume        = ft_datatype(data, 'volume');
issegmentation  = ft_datatype(data, 'segmentation');
isparcellation  = ft_datatype(data, 'parcellation');
issource        = ft_datatype(data, 'source');
isdip           = ft_datatype(data, 'dip');
ismvar          = ft_datatype(data, 'mvar');
isfreqmvar      = ft_datatype(data, 'freqmvar');
ischan          = ft_datatype(data, 'chan');
ismesh          = ft_datatype(data, 'mesh');
% FIXME use the istrue function on ismeg and hasxxx options

if ~isequal(feedback, 'no') % can be 'yes' or 'text'
  if iscomp
    % it can be comp and raw/timelock/freq at the same time, therefore this has to go first
    nchan = size(data.topo,1);
    ncomp = size(data.topo,2);
    ft_info('the input is component data with %d components and %d original channels\n', ncomp, nchan);
  end  % if iscomp

  if ismesh
    % it can be comp and source at the same time, therefore this has to go first
    data = fixpos(data);
    npos = 0;
    ntri = 0;
    nhex = 0;
    ntet = 0;
    % the data can contain multiple surfaces
    for i=1:numel(data)
      npos = npos+size(data.pos,1);
      if isfield(data, 'tri'), ntri = ntri+size(data.tri,1); end
      if isfield(data, 'hex'), nhex = nhex+size(data.hex,1); end
      if isfield(data, 'tet'), ntet = ntet+size(data.tet,1); end
    end
    if isfield(data,'tri')
      ft_info('the input is mesh data with %d vertices and %d triangles\n', npos, ntri);
    elseif isfield(data,'hex')
      ft_info('the input is mesh data with %d vertices and %d hexahedrons\n', npos, nhex);
    elseif isfield(data,'tet')
      ft_info('the input is mesh data with %d vertices and %d tetrahedrons\n', npos, ntet);
    else
      ft_info('the input is mesh data with %d vertices', npos);
    end
  end % if ismesh

  if israw
    nchan = length(data.label);
    ntrial = length(data.trial);
    ft_info('the input is raw data with %d channels and %d trials\n', nchan, ntrial);
  elseif istimelock
    nchan = length(data.label);
    ntime = length(data.time);
    ft_info('the input is timelock data with %d channels and %d timebins\n', nchan, ntime);
  elseif isfreq
    if isfield(data, 'label')
      nchan = length(data.label);
      nfreq = length(data.freq);
      if isfield(data, 'time'), ntime = num2str(length(data.time)); else ntime = 'no'; end
      ft_info('the input is freq data with %d channels, %d frequencybins and %s timebins\n', nchan, nfreq, ntime);
    elseif isfield(data, 'labelcmb')
      nchan = length(data.labelcmb);
      nfreq = length(data.freq);
      if isfield(data, 'time'), ntime = num2str(length(data.time)); else ntime = 'no'; end
      ft_info('the input is freq data with %d channel combinations, %d frequencybins and %s timebins\n', nchan, nfreq, ntime);
    else
      ft_error('cannot infer freq dimensions');
    end
  elseif isspike
    nchan  = length(data.label);
    ft_info('the input is spike data with %d channels\n', nchan);
  elseif isvolume
    if issegmentation
      ft_info('the input is segmented volume data with dimensions [%d %d %d]\n', data.dim(1), data.dim(2), data.dim(3));
      print_voxelinfo(data)
      print_segmentationinfo(data)
    else
      ft_info('the input is volume data with dimensions [%d %d %d]\n', data.dim(1), data.dim(2), data.dim(3));
      print_voxelinfo(data)
    end
  elseif issource
    data = fixpos(data); % ensure that positions are in pos, not in pnt
    nsource = size(data.pos, 1);
    if isparcellation
      subtype = 'parcellated source';
    else
      subtype = 'source';
    end
    if isfield(data, 'dim')
      ft_info('the input is %s data with %d brainordinates on a [%d %d %d] grid\n', subtype, nsource, data.dim(1), data.dim(2), data.dim(3));
    else
      ft_info('the input is %s data with %d brainordinates\n', subtype, nsource);
    end
    clear subtype
  elseif isdip
    ft_info('the input is dipole data\n');
  elseif ismvar
    ft_info('the input is mvar data\n');
  elseif isfreqmvar
    ft_info('the input is freqmvar data\n');
  elseif ischan
    nchan = length(data.label);
    if isfield(data, 'brainordinate')
      ft_info('the input is parcellated data with %d parcels\n', nchan);
    else
      ft_info('the input is chan data with %d channels\n', nchan);
    end
  end % if israw etc.
end % give feedback

if issource && isvolume
  % it should be either one or the other: the choice here is to represent it as volume description since that is simpler to handle
  % the conversion is done by removing the grid positions
  data = rmfield(data, 'pos');
  issource = false;
end

if isfield(data, 'trialinfo')
  if strcmp(trialinfostyle, 'table')
    if ismatrix(data.trialinfo)
      data.trialinfo = array2table(data.trialinfo);
    end
  elseif strcmp(trialinfostyle, 'matrix')
    if istable(data.trialinfo)
      data.trialinfo = table2array(data.trialinfo);
    end
  else
    % no conversion is needed
  end
end

% the ft_datatype_XXX functions ensures the consistency of the XXX datatype
% and provides a detailed description of the dataformat and its history
if iscomp % this should go before israw/istimelock/isfreq
  data = ft_datatype_comp(data, 'hassampleinfo', hassampleinfo);
elseif israw
  data = ft_datatype_raw(data, 'hassampleinfo', hassampleinfo);
elseif istimelock
  data = ft_datatype_timelock(data, 'hassampleinfo', hassampleinfo);
elseif isfreq
  data = ft_datatype_freq(data);
elseif isspike
  data = ft_datatype_spike(data);
elseif issegmentation % this should go before isvolume
  data = ft_datatype_segmentation(data, 'segmentationstyle', segmentationstyle, 'hasbrain', hasbrain);
elseif isvolume
  data = ft_datatype_volume(data);
elseif isparcellation % this should go before issource
  data = ft_datatype_parcellation(data, 'parcellationstyle', parcellationstyle);
elseif issource
  data = ft_datatype_source(data);
elseif isdip
  data = ft_datatype_dip(data);
elseif ismvar || isfreqmvar
  data = ft_datatype_mvar(data);
end

if ~isempty(dtype)
  if ~isa(dtype, 'cell')
    dtype = {dtype};
  end

  okflag = 0;
  for i=1:length(dtype)
    % check that the data matches with one or more of the required ft_datatypes
    switch dtype{i}
      case 'raw+comp'
        okflag = okflag + (israw & iscomp);
      case 'freq+comp'
        okflag = okflag + (isfreq & iscomp);
      case 'timelock+comp'
        okflag = okflag + (istimelock & iscomp);
      case 'source+mesh'
        okflag = okflag + (issource & ismesh);
      case 'raw'
        okflag = okflag + (israw & ~iscomp);
      case 'freq'
        okflag = okflag + (isfreq & ~iscomp);
      case 'timelock'
        okflag = okflag + (istimelock & ~iscomp);
      case 'comp'
        okflag = okflag + (iscomp && ~(israw || istimelock || isfreq));
      case 'spike'
        okflag = okflag + isspike;
      case 'volume'
        okflag = okflag + isvolume;
      case 'source'
        okflag = okflag + issource;
      case 'dip'
        okflag = okflag + isdip;
      case 'mvar'
        okflag = okflag + ismvar;
      case 'freqmvar'
        okflag = okflag + isfreqmvar;
      case 'chan'
        okflag = okflag + ischan;
      case 'segmentation'
        okflag = okflag + issegmentation;
      case 'parcellation'
        okflag = okflag + isparcellation;
      case 'mesh'
        okflag = okflag + ismesh;
    end % switch dtype
  end % for dtype

  % try to convert the data if needed
  for iCell = 1:length(dtype)
    if okflag
      % the requested datatype is specified in descending order of
      % preference (if there is a preference at all), so don't bother
      % checking the rest of the requested data types if we already
      % succeeded in converting
      break;
    end
    if isequal(dtype(iCell), {'parcellation'}) && issegmentation
      data = volume2source(data); % segmentation=volume, parcellation=source
      data = ft_datatype_parcellation(data);
      issegmentation = 0;
      isvolume = 0;
      isparcellation = 1;
      issource = 1;
      okflag = 1;
    elseif isequal(dtype(iCell), {'segmentation'}) && isparcellation
      data = source2volume(data); % segmentation=volume, parcellation=source
      data = ft_datatype_segmentation(data);
      isparcellation = 0;
      issource = 0;
      issegmentation = 1;
      isvolume = 1;
      okflag = 1;
    elseif isequal(dtype(iCell), {'source'}) && isvolume
      data = volume2source(data);
      data = ft_datatype_source(data);
      isvolume = 0;
      issource = 1;
      okflag = 1;
    elseif isequal(dtype(iCell), {'volume'}) && (ischan || istimelock || isfreq)
      if isfield(data, 'brainordinate')
        data = parcellated2source(data);
        data = ft_datatype_volume(data);
      else
        ft_error('cannot convert channel-level data to volumetric representation');
      end
      ischan = 0; istimelock = 0; isfreq = 0;
      isvolume = 1;
      okflag = 1;
    elseif (isequal(dtype(iCell), {'source'}) || isequal(dtype(iCell), {'source+mesh'})) && (ischan || istimelock || isfreq)
      if isfield(data, 'brainordinate')
        data = parcellated2source(data);
        data = ft_datatype_source(data);
      else
        data = chan2source(data);
        data = ft_datatype_source(data);
      end % converting channel data
      ischan = 0; istimelock = 0; isfreq = 0;
      issource = 1;
      okflag = 1;
    elseif isequal(dtype(iCell), {'volume'}) && issource
      data = source2volume(data);
      data = ft_datatype_volume(data);
      isvolume = 1;
      issource = 0;
      okflag = 1;
    elseif isequal(dtype(iCell), {'raw'}) && issource
      data = source2raw(data);
      data = ft_datatype_raw(data, 'hassampleinfo', hassampleinfo);
      issource = 0;
      israw = 1;
      okflag = 1;
    elseif isequal(dtype(iCell), {'raw'}) && istimelock
      if iscomp
        data = removefields(data, {'topo', 'topolabel', 'topodimord', 'unmixing', 'unmixingdimord'}); % these fields are not desired
        iscomp = 0;
      end
      data = timelock2raw(data);
      data = ft_datatype_raw(data, 'hassampleinfo', hassampleinfo);
      istimelock = 0;
      israw = 1;
      okflag = 1;
    elseif isequal(dtype(iCell), {'raw+comp'}) && istimelock && iscomp
      data = timelock2raw(data);
      data = ft_datatype_raw(data, 'hassampleinfo', hassampleinfo);
      istimelock = 0;
      iscomp = 1;
      israw = 1;
      okflag = 1;
    elseif isequal(dtype(iCell), {'timelock+comp'}) && israw && iscomp
      data = raw2timelock(data);
      data = ft_datatype_timelock(data, 'hassampleinfo', hassampleinfo);
      istimelock = 1;
      iscomp = 1;
      israw = 0;
      okflag = 1;
    elseif isequal(dtype(iCell), {'comp'}) && israw  && iscomp
      data = keepfields(data, {'label', 'topo', 'topolabel', 'unmixing', 'elec', 'grad', 'cfg'}); % these are the only relevant fields
      data = ft_datatype_comp(data, 'hassampleinfo', hassampleinfo);
      israw = 0;
      iscomp = 1;
      okflag = 1;
    elseif isequal(dtype(iCell), {'comp'}) && istimelock && iscomp
      data = keepfields(data, {'label', 'topo', 'topolabel', 'unmixing', 'elec', 'grad', 'cfg'}); % these are the only relevant fields
      data = ft_datatype_comp(data, 'hassampleinfo', hassampleinfo);
      istimelock = 0;
      iscomp = 1;
      okflag = 1;
    elseif isequal(dtype(iCell), {'comp'}) && isfreq && iscomp
      data = keepfields(data, {'label', 'topo', 'topolabel', 'unmixing', 'elec', 'grad', 'cfg'}); % these are the only relevant fields
      data = ft_datatype_comp(data, 'hassampleinfo', 'no'); % freq data does not have sampleinfo
      isfreq = 0;
      iscomp = 1;
      okflag = 1;
    elseif isequal(dtype(iCell), {'raw'}) && israw
      if iscomp
        data = removefields(data, {'topo', 'topolabel', 'topodimord', 'unmixing', 'unmixingdimord'}); % these fields are not desired
        iscomp = 0;
      end
      data = ft_datatype_raw(data, 'hassampleinfo', hassampleinfo);
      okflag = 1;
    elseif isequal(dtype(iCell), {'timelock'}) && istimelock
      if iscomp
        data = removefields(data, {'topo', 'topolabel', 'topodimord', 'unmixing', 'unmixingdimord'}); % these fields are not desired
        iscomp = 0;
      end
      data = ft_datatype_timelock(data, 'hassampleinfo', hassampleinfo);
      okflag = 1;
    elseif isequal(dtype(iCell), {'freq'}) && isfreq
      if iscomp
        data = removefields(data, {'topo', 'topolabel', 'topodimord', 'unmixing', 'unmixingdimord'}); % these fields are not desired
        iscomp = 0;
      end
      data = ft_datatype_freq(data);
      okflag = 1;
    elseif isequal(dtype(iCell), {'timelock'}) && israw
      if iscomp
        data = removefields(data, {'topo', 'topolabel', 'topodimord', 'unmixing', 'unmixingdimord'}); % these fields are not desired
        iscomp = 0;
      end
      data = raw2timelock(data);
      data = ft_datatype_timelock(data, 'hassampleinfo', hassampleinfo);
      israw = 0;
      istimelock = 1;
      okflag = 1;
    elseif isequal(dtype(iCell), {'raw'}) && isfreq
      if iscomp
        data = removefields(data, {'topo', 'topolabel', 'topodimord', 'unmixing', 'unmixingdimord'}); % these fields are not desired
        iscomp = 0;
      end
      data = freq2raw(data);
      data = ft_datatype_raw(data, 'hassampleinfo', hassampleinfo);
      isfreq = 0;
      israw = 1;
      okflag = 1;

    elseif isequal(dtype(iCell), {'raw'}) && ischan
      data = chan2timelock(data);
      data = timelock2raw(data);
      data = ft_datatype_raw(data, 'hassampleinfo', hassampleinfo);
      ischan = 0;
      israw = 1;
      okflag = 1;
    elseif isequal(dtype(iCell), {'timelock'}) && ischan
      data = chan2timelock(data);
      data = ft_datatype_timelock(data, 'hassampleinfo', hassampleinfo);
      ischan = 0;
      istimelock = 1;
      okflag = 1;
    elseif isequal(dtype(iCell), {'freq'}) && ischan
      data = chan2freq(data);
      data = ft_datatype_freq(data);
      ischan = 0;
      isfreq = 1;
      okflag = 1;
    elseif isequal(dtype(iCell), {'spike'}) && israw
      data = raw2spike(data);
      data = ft_datatype_spike(data);
      israw = 0;
      isspike = 1;
      okflag = 1;
    elseif isequal(dtype(iCell), {'raw'}) && isspike
      data = spike2raw(data,fsample);
      data = ft_datatype_raw(data, 'hassampleinfo', hassampleinfo);
      isspike = 0;
      israw   = 1;
      okflag  = 1;
    end
  end % for iCell

  if ~okflag
    % construct an error message
    typestr = printor(dtype, true);
    helpfun = cell(size(dtype));
    for i=1:numel(dtype)
      helpfun{i} = sprintf('ft_datatype_%s', dtype{i});
    end
    helpfun = helpfun(cellfun(@exist, helpfun)>0);
    if ~isempty(helpfun)
      helpstr = printor(helpfun);
      ft_error('This function requires %s data as input, see %s.', typestr, helpstr);
    else
      ft_error('This function requires %s data as input.', typestr);
    end
  end % if okflag
end

if ~isempty(dimord)
  if ~isa(dimord, 'cell')
    dimord = {dimord};
  end

  if isfield(data, 'dimord')
    okflag = any(strcmp(data.dimord, dimord));
  else
    okflag = 0;
  end

  if ~okflag
    % construct an error message
    ft_error('This function requires data with a dimord of %s.', printor(dimord, true));
  end % if okflag
end

if ~isempty(stype)
  if ~isa(stype, 'cell')
    stype = {stype};
  end

  if isfield(data, 'grad') || isfield(data, 'elec') || isfield(data, 'opto')
    if any(strcmp(ft_senstype(data), stype))
      okflag = 1;
    elseif any(cellfun(@ft_senstype, repmat({data}, size(stype)), stype))
      % this is required to detect more general types, such as "meg" or "ctf" rather than "ctf275"
      okflag = 1;
    else
      okflag = 0;
    end
  else
    % the data does not contain a sensor array
    okflag = 0;
  end

  if ~okflag
    % construct an error message
    ft_error('This function requires data with an %s sensor array.', printor(stype, true));
  end % if okflag
end

if ~isempty(ismeg)
  if isequal(ismeg, 'yes')
    okflag = isfield(data, 'grad');
  elseif isequal(ismeg, 'no')
    okflag = ~isfield(data, 'grad');
  end

  if ~okflag && isequal(ismeg, 'yes')
    ft_error('This function requires MEG data with a ''grad'' field');
  elseif ~okflag && isequal(ismeg, 'no')
    ft_error('This function should not be given MEG data with a ''grad'' field');
  end % if okflag
end

if ~isempty(iseeg)
  if isequal(iseeg, 'yes')
    okflag = isfield(data, 'elec');
  elseif isequal(iseeg, 'no')
    okflag = ~isfield(data, 'elec');
  end

  if ~okflag && isequal(iseeg, 'yes')
    ft_error('This function requires EEG data with an ''elec'' field');
  elseif ~okflag && isequal(iseeg, 'no')
    ft_error('This function should not be given EEG data with an ''elec'' field');
  end % if okflag
end

if ~isempty(isnirs)
  if isequal(isnirs, 'yes')
    okflag = isfield(data, 'opto');
  elseif isequal(isnirs, 'no')
    okflag = ~isfield(data, 'opto');
  end

  if ~okflag && isequal(isnirs, 'yes')
    ft_error('This function requires NIRS data with an ''opto'' field');
  elseif ~okflag && isequal(isnirs, 'no')
    ft_error('This function should not be given NIRS data with an ''opto'' field');
  end % if okflag
end

if ~isempty(insidestyle)
  if strcmp(insidestyle, 'index')
    ft_warning('the indexed representation of inside/outside source locations is deprecated');
  end
  % TODO absorb the fixinside function into this code
  data   = fixinside(data, insidestyle);
  okflag = isfield(data, 'inside');

  if ~okflag
    % construct an error message
    ft_error('This function requires data with an ''inside'' field.');
  end % if okflag
end

if istrue(hasunit) && ~isfield(data, 'unit')
  % calling convert_units with only the input data adds the units without converting
  data = ft_determine_units(data);
end % if hasunit

if istrue(hascoordsys) && ~isfield(data, 'coordsys')
  data = ft_determine_coordsys(data);
end % if hascoordsys

if istrue(haschantype) && ~isfield(data, 'chantype')
  data.chantype = ft_chantype(data);
end % if haschantype

if istrue(haschanunit) && ~isfield(data, 'chanunit')
  data.chanunit = ft_chanunit(data);
end % if haschanunit

if isequal(hastrials, 'yes')
  hasrpt = isfield(data, 'trial');
  if ~hasrpt && isfield(data, 'dimord')
    % instead look in the dimord for rpt or subj
    hasrpt = ~isempty(strfind(data.dimord, 'rpt')) || ...
      ~isempty(strfind(data.dimord, 'rpttap')) || ...
      ~isempty(strfind(data.dimord, 'subj'));
  end
  if ~hasrpt
    ft_error('This function requires data with a ''trial'' field');
  end % if hasrpt
elseif isequal(hastrials, 'no') && istimelock
  if ~isfield(data, 'avg') && (isfield(data, 'trial') || isfield(data, 'individual'))
    % average on the fly
    tmpcfg = [];
    tmpcfg.keeptrials = 'no';
    data = ft_timelockanalysis(tmpcfg, data);
  end
end

if strcmp(hasdim, 'yes') && ~isfield(data, 'dim')
  data.dim = pos2dim(data.pos);
elseif strcmp(hasdim, 'no') && isfield(data, 'dim')
  data = rmfield(data, 'dim');
end % if hasdim

if strcmp(hascumtapcnt, 'yes') && ~isfield(data, 'cumtapcnt')
  ft_error('This function requires data with a ''cumtapcnt'' field');
elseif strcmp(hascumtapcnt, 'no') && isfield(data, 'cumtapcnt')
  data = rmfield(data, 'cumtapcnt');
end % if hascumtapcnt

if strcmp(hasdof, 'yes') && ~isfield(data, 'dof')
  ft_error('This function requires data with a ''dof'' field');
elseif strcmp(hasdof, 'no') && isfield(data, 'dof')
  data = rmfield(data, 'dof');
end % if hasdof

if ~isempty(cmbstyle)
  if istimelock
    data = fixcov(data, cmbstyle);
  elseif isfreq
    data = fixcsd(data, cmbstyle, channelcmb);
  elseif isfreqmvar
    data = fixcsd(data, cmbstyle, channelcmb);
  else
    ft_error('this function requires data with a covariance, coherence or cross-spectrum');
  end
end % cmbstyle

if isfield(data, 'grad')
  % ensure that the gradiometer structure is up to date
  data.grad = ft_datatype_sens(data.grad);
end

if isfield(data, 'elec')
  % ensure that the electrode structure is up to date
  data.elec = ft_datatype_sens(data.elec);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% represent the covariance matrix in a particular manner
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [data] = fixcov(data, desired)
if any(isfield(data, {'cov', 'corr'}))
  if ~isfield(data, 'labelcmb')
    current = 'full';
  else
    current = 'sparse';
  end
else
  ft_error('Could not determine the current representation of the covariance matrix');
end
if isequal(current, desired)
  % nothing to do
elseif strcmp(current, 'full') && strcmp(desired, 'sparse')
  % FIXME should be implemented
  ft_error('not yet implemented');
elseif strcmp(current, 'sparse') && strcmp(desired, 'full')
  % FIXME should be implemented
  ft_error('not yet implemented');
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% represent the cross-spectral density matrix in a particular manner
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [data] = fixcsd(data, desired, channelcmb)

% FIXCSD converts univariate frequency domain data (fourierspctrm) into a bivariate
% representation (crsspctrm), or changes the representation of bivariate frequency
% domain data (sparse/full/sparsewithpow, sparsewithpow only works for crsspctrm or
% fourierspctrm)

% Copyright (C) 2010, Jan-Mathijs Schoffelen, Robert Oostenveld

if isfield(data, 'crsspctrm') && isfield(data, 'powspctrm')
  current = 'sparsewithpow';
elseif isfield(data, 'powspctrm')
  current = 'sparsewithpow';
elseif isfield(data, 'fourierspctrm') && ~isfield(data, 'labelcmb')
  current = 'fourier';
elseif ~isfield(data, 'labelcmb')
  current = 'full';
elseif isfield(data, 'labelcmb')
  current = 'sparse';
else
  ft_error('Could not determine the current representation of the %s matrix', param);
end

% first go from univariate fourier to the required bivariate representation
if isequal(current, desired)
  % nothing to do

elseif strcmp(current, 'fourier') && strcmp(desired, 'sparsewithpow')
  if startsWith(data.dimord, 'rpttap')
    nrpt = size(data.cumtapcnt,1);
  else
    nrpt = 1;
  end
  if contains(data.dimord, 'freq'), nfrq = length(data.freq); else nfrq = 1; end
  if contains(data.dimord, 'time'), ntim = length(data.time); else ntim = 1; end

  fastflag = all(data.cumtapcnt(:)==data.cumtapcnt(1));
  flag     = nrpt==1; % needed to truncate the singleton dimension upfront

  % create auto-spectra
  nchan     = length(data.label);
  if fastflag
    % all trials have the same amount of tapers
    powspctrm = zeros(nrpt,nchan,nfrq,ntim);
    ntap      = data.cumtapcnt(1);
    for p = 1:ntap
      powspctrm = powspctrm + abs(data.fourierspctrm(p:ntap:end,:,:,:,:)).^2;
    end
    powspctrm = powspctrm./ntap;
  else
    % different amount of tapers
    powspctrm = zeros(nrpt,nchan,nfrq,ntim) + zeros(nrpt,nchan,nfrq,ntim)*1i;
    sumtapcnt = [0;cumsum(data.cumtapcnt(:))];
    for p = 1:nrpt
      indx   = (sumtapcnt(p)+1):sumtapcnt(p+1);
      tmpdat = data.fourierspctrm(indx,:,:,:);
      powspctrm(p,:,:,:) = (sum(tmpdat.*conj(tmpdat),1))./data.cumtapcnt(p);
    end
  end

  % create cross-spectra
  if ~isempty(channelcmb)
    ncmb      = size(channelcmb,1);
    cmbindx   = zeros(ncmb,2);
    labelcmb  = cell(ncmb,2);
    for k = 1:ncmb
      ch1 = find(strcmp(data.label, channelcmb(k,1)));
      ch2 = find(strcmp(data.label, channelcmb(k,2)));
      if ~isempty(ch1) && ~isempty(ch2)
        cmbindx(k,:)  = [ch1 ch2];
        labelcmb(k,:) = data.label([ch1 ch2])';
      end
    end

    crsspctrm = zeros(nrpt,ncmb,nfrq,ntim) + 1i.*zeros(nrpt,ncmb,nfrq,ntim);
    if fastflag
      for p = 1:ntap
        tmpdat1   = data.fourierspctrm(p:ntap:end,cmbindx(:,1),:,:,:);
        tmpdat2   = data.fourierspctrm(p:ntap:end,cmbindx(:,2),:,:,:);
        crsspctrm = crsspctrm + tmpdat1.*conj(tmpdat2);
      end
      crsspctrm = crsspctrm./ntap;
    else
      for p = 1:nrpt
        indx    = (sumtapcnt(p)+1):sumtapcnt(p+1);
        tmpdat1 = data.fourierspctrm(indx,cmbindx(:,1),:,:);
        tmpdat2 = data.fourierspctrm(indx,cmbindx(:,2),:,:);
        crsspctrm(p,:,:,:) = (sum(tmpdat1.*conj(tmpdat2),1))./data.cumtapcnt(p);
      end
    end
    data.crsspctrm = crsspctrm;
    data.labelcmb  = labelcmb;
  end
  data.powspctrm = powspctrm;
  data           = rmfield(data, 'fourierspctrm');
  if ntim>1
    data.dimord = 'chan_freq_time';
  else
    data.dimord = 'chan_freq';
  end

  if nrpt>1
    data.dimord = ['rpt_' data.dimord];
  end

  if flag
    siz = size(data.powspctrm);
    data.powspctrm = reshape(data.powspctrm, [siz(2:end) 1]);
    if isfield(data, 'crsspctrm')
      siz = size(data.crsspctrm);
      data.crsspctrm = reshape(data.crsspctrm, [siz(2:end) 1]);
    end
  end

elseif strcmp(current, 'fourier') && strcmp(desired, 'sparse')
  if isempty(channelcmb), ft_error('no channel combinations are specified'); end
  if startsWith(data.dimord, 'rpttap')
    nrpt = size(data.cumtapcnt,1);
  else
    nrpt = 1;
  end
  if contains(data.dimord, 'freq'), nfrq = length(data.freq); else nfrq = 1; end
  if contains(data.dimord, 'time'), ntim = length(data.time); else ntim = 1; end

  flag      = nrpt==1; % flag needed to squeeze first dimension if singleton
  ncmb      = size(channelcmb,1);
  cmbindx   = zeros(ncmb,2);
  labelcmb  = cell(ncmb,2);
  for k = 1:ncmb
    ch1 = find(strcmp(data.label, channelcmb(k,1)));
    ch2 = find(strcmp(data.label, channelcmb(k,2)));
    if ~isempty(ch1) && ~isempty(ch2)
      cmbindx(k,:)  = [ch1 ch2];
      labelcmb(k,:) = data.label([ch1 ch2])';
    end
  end

  sumtapcnt = [0;cumsum(data.cumtapcnt(:))];
  fastflag  = all(data.cumtapcnt(:)==data.cumtapcnt(1));

  if fastflag && nrpt>1
    ntap = data.cumtapcnt(1);

    % compute running sum across tapers
    siz = [size(data.fourierspctrm) 1];

    for p = 1:ntap
      indx      = p:ntap:nrpt*ntap;

      if p==1
        tmpc = zeros(numel(indx), size(cmbindx,1), siz(3), siz(4)) + ...
          1i.*zeros(numel(indx), size(cmbindx,1), siz(3), siz(4));
      end

      for k = 1:size(cmbindx,1)
        tmpc(:,k,:,:) = data.fourierspctrm(indx,cmbindx(k,1),:,:).*  ...
          conj(data.fourierspctrm(indx,cmbindx(k,2),:,:));
      end

      if p==1
        crsspctrm = tmpc;
      else
        crsspctrm = tmpc + crsspctrm;
      end
    end
    crsspctrm = crsspctrm./ntap;
  else
    crsspctrm = zeros(nrpt, ncmb, nfrq, ntim);
    for p = 1:nrpt
      indx    = (sumtapcnt(p)+1):sumtapcnt(p+1);
      tmpdat1 = data.fourierspctrm(indx,cmbindx(:,1),:,:);
      tmpdat2 = data.fourierspctrm(indx,cmbindx(:,2),:,:);
      crsspctrm(p,:,:,:) = (sum(tmpdat1.*conj(tmpdat2),1))./data.cumtapcnt(p);
    end
  end
  data.crsspctrm = crsspctrm;
  data.labelcmb  = labelcmb;
  data           = rmfield(data, 'fourierspctrm');
  data           = rmfield(data, 'label');
  if ntim>1
    data.dimord = 'chancmb_freq_time';
  else
    data.dimord = 'chancmb_freq';
  end

  if nrpt>1
    data.dimord = ['rpt_' data.dimord];
  end

  if flag
    if isfield(data,'powspctrm')
      % deal with the singleton 'rpt', i.e. remove it
      siz = size(data.powspctrm);
      data.powspctrm = reshape(data.powspctrm, [siz(2:end) 1]);
    end
    if isfield(data,'crsspctrm')
      % this conditional statement is needed in case there's a single channel
      siz            = size(data.crsspctrm);
      data.crsspctrm = reshape(data.crsspctrm, [siz(2:end) 1]);
    end
  end

elseif strcmp(current, 'fourier') && strcmp(desired, 'full')
  % this is how it is currently and the desired functionality of prepare_freq_matrices
  if startsWith(data.dimord, 'rpttap')
    nrpt = size(data.cumtapcnt, 1);
    flag = 0;
  else
    nrpt = 1;
    flag = 1;
  end
  if contains(data.dimord, 'freq'), nfrq = length(data.freq); else nfrq = 1; end
  if contains(data.dimord, 'time'), ntim = length(data.time); else ntim = 1; end
  if any(data.cumtapcnt(1,:) ~= data.cumtapcnt(1,1)), ft_error('this only works when all frequencies have the same number of tapers'); end
  nchan     = length(data.label);
  crsspctrm = zeros(nrpt,nchan,nchan,nfrq,ntim);
  sumtapcnt = [0;cumsum(data.cumtapcnt(:,1))];
  for k = 1:ntim
    for m = 1:nfrq
      for p = 1:nrpt
        % FIXME speed this up in the case that all trials have equal number of tapers
        indx   = (sumtapcnt(p)+1):sumtapcnt(p+1);
        tmpdat = transpose(data.fourierspctrm(indx,:,m,k));
        crsspctrm(p,:,:,m,k) = (tmpdat*tmpdat')./data.cumtapcnt(p);
        clear tmpdat;
      end
    end
  end
  data.crsspctrm = crsspctrm;
  data           = rmfield(data, 'fourierspctrm');

  if ntim>1
    data.dimord = 'chan_chan_freq_time';
  else
    data.dimord = 'chan_chan_freq';
  end

  if nrpt>1
    data.dimord = ['rpt_' data.dimord];
  end

  % remove first singleton dimension
  if flag || nrpt==1, siz = size(data.crsspctrm); data.crsspctrm = reshape(data.crsspctrm, siz(2:end)); end

elseif strcmp(current, 'fourier') && strcmp(desired, 'fullfast')
  nrpt = size(data.fourierspctrm, 1);
  nchn = numel(data.label);
  nfrq = length(data.freq);
  if contains(data.dimord, 'time'), ntim = length(data.time); else ntim = 1; end

  data.fourierspctrm = reshape(data.fourierspctrm, [nrpt nchn nfrq*ntim]);
  %data.fourierspctrm(~isfinite(data.fourierspctrm)) = 0;
  crsspctrm = complex(zeros(nchn,nchn,nfrq*ntim));
  for k = 1:nfrq*ntim
    tmp = transpose(data.fourierspctrm(:,:,k));
    tmp(~isfinite(tmp)) = 0;
    n   = sum(tmp~=0,2);
    crsspctrm(:,:,k) = tmp*tmp'./n(1);
  end
  data           = rmfield(data, 'fourierspctrm');
  data.crsspctrm = reshape(crsspctrm, [nchn nchn nfrq ntim]);
  if isfield(data, 'time')
    data.dimord = 'chan_chan_freq_time';
  else
    data.dimord = 'chan_chan_freq';
  end

  if isfield(data, 'trialinfo'),  data = rmfield(data, 'trialinfo'); end
  if isfield(data, 'sampleinfo'), data = rmfield(data, 'sampleinfo'); end
  if isfield(data, 'cumsumcnt'),  data = rmfield(data, 'cumsumcnt');  end
  if isfield(data, 'cumtapcnt'),  data = rmfield(data, 'cumtapcnt');  end

end % convert to the requested bivariate representation

% from one bivariate representation to another
if isequal(current, desired)
  % nothing to do

elseif (strcmp(current, 'full')       && strcmp(desired, 'fourier')) || ...
    (strcmp(current, 'sparse')        && strcmp(desired, 'fourier')) || ...
    (strcmp(current, 'sparsewithpow') && strcmp(desired, 'fourier'))
  % this is not possible
  ft_error('converting the cross-spectrum into a Fourier representation is not possible');

elseif strcmp(current, 'full') && strcmp(desired, 'sparsewithpow')
  ft_error('not yet implemented');

elseif strcmp(current, 'sparse') && strcmp(desired, 'sparsewithpow')
  % convert back to crsspctrm/powspctrm representation: useful for plotting functions etc
  indx     = labelcmb2indx(data.labelcmb);
  autoindx = indx(indx(:,1)==indx(:,2), 1);
  cmbindx  = setdiff(1:size(indx,1), autoindx);

  if strcmp(data.dimord(1:3), 'rpt')
    data.powspctrm = data.crsspctrm(:, autoindx, :, :);
    data.crsspctrm = data.crsspctrm(:, cmbindx,  :, :);
  else
    data.powspctrm = data.crsspctrm(autoindx, :, :);
    data.crsspctrm = data.crsspctrm(cmbindx,  :, :);
  end
  data.label    = data.labelcmb(autoindx,1);
  data.labelcmb = data.labelcmb(cmbindx, :);

  if isempty(cmbindx)
    data = rmfield(data, 'crsspctrm');
    data = rmfield(data, 'labelcmb');
  end

elseif strcmp(current, 'full') && strcmp(desired, 'sparse')
  if contains(data.dimord, 'rpt'),  nrpt = size(data.cumtapcnt,1); else nrpt = 1; end
  if contains(data.dimord, 'freq'), nfrq = length(data.freq); else nfrq = 1; end
  if contains(data.dimord, 'time'), ntim = length(data.time); else ntim = 1; end
  nchan    = length(data.label);
  ncmb     = nchan*nchan;
  labelcmb = cell(ncmb, 2);
  cmbindx  = zeros(nchan, nchan);
  k = 1;
  for j=1:nchan
    for m=1:nchan
      labelcmb{k, 1} = data.label{m};
      labelcmb{k, 2} = data.label{j};
      cmbindx(m,j)   = k;
      k = k+1;
    end
  end

  % reshape all possible fields
  fn = fieldnames(data);
  for ii=1:numel(fn)
    if numel(data.(fn{ii})) == nrpt*ncmb*nfrq*ntim
      if nrpt>1
        data.(fn{ii}) = reshape(data.(fn{ii}), nrpt, ncmb, nfrq, ntim);
      else
        data.(fn{ii}) = reshape(data.(fn{ii}), ncmb, nfrq, ntim);
      end
    end
  end
  % remove obsolete fields
  data = removefields(data, {'label', 'dof'});

  % replace updated fields
  data.labelcmb  = labelcmb;
  if ntim>1
    data.dimord = 'chancmb_freq_time';
  else
    data.dimord = 'chancmb_freq';
  end

  if nrpt>1
    data.dimord = ['rpt_' data.dimord];
  end

elseif strcmp(current, 'sparsewithpow') && strcmp(desired, 'sparse')
  % this representation for sparse data contains autospectra as e.g. {'A' 'A'} in labelcmb
  if isfield(data, 'crsspctrm')
    dimtok         = tokenize(data.dimord, '_');
    catdim         = match_str(dimtok, {'chan' 'chancmb'});
    data.crsspctrm = cat(catdim, data.powspctrm, data.crsspctrm);
    data.labelcmb  = [data.label(:) data.label(:); data.labelcmb];
    data           = rmfield(data, 'powspctrm');
    data.dimord    = strrep(data.dimord, 'chan_', 'chancmb_');
  else
    data.crsspctrm = data.powspctrm;
    data.labelcmb  = [data.label(:) data.label(:)];
    data           = rmfield(data, 'powspctrm');
    data.dimord    = strrep(data.dimord, 'chan_', 'chancmb_');
  end
  data = rmfield(data, 'label');

elseif strcmp(current, 'sparse') && strcmp(desired, 'full')
  if contains(data.dimord, 'rpt'),  nrpt = size(data.cumtapcnt,1); else nrpt = 1; end
  if contains(data.dimord, 'freq'), nfrq = length(data.freq); else nfrq = 1; end
  if contains(data.dimord, 'time'), ntim = length(data.time); else ntim = 1; end

  if ~isfield(data, 'label')
    % ensure that the bivariate spectral factorization results can be
    % processed. FIXME this is experimental and will not work if the user
    % did something weird before
    for k = 1:numel(data.labelcmb)
      tmp = tokenize(data.labelcmb{k}, '[');
      data.labelcmb{k} = tmp{1};
    end
    data.label = unique(data.labelcmb(:));
  end

  nchan     = length(data.label);
  ncmb      = size(data.labelcmb,1);
  cmbindx   = zeros(nchan,nchan);

  for k = 1:size(data.labelcmb,1)
    ch1 = find(strcmp(data.label, data.labelcmb(k,1)));
    ch2 = find(strcmp(data.label, data.labelcmb(k,2)));
    if ~isempty(ch1) && ~isempty(ch2)
      cmbindx(ch1,ch2) = k;
    end
  end

  complete = all(cmbindx(:)~=0);

  % remove obsolete fields
  data = removefields(data, {'powspctrm', 'labelcmb', 'dof'});

  fn = setdiff(fieldnames(data), {'time' 'freq' 'dimord', 'label', 'cfg'});
  for ii=1:numel(fn)
    if numel(data.(fn{ii})) == nrpt*ncmb*nfrq*ntim
      if nrpt==1
        data.(fn{ii}) = reshape(data.(fn{ii}), [nrpt ncmb nfrq ntim]);
      end

      tmpall = nan(nrpt,nchan,nchan,nfrq,ntim);

      for j = 1:nrpt
        for k = 1:ntim
          for m = 1:nfrq
            tmpdat = nan(nchan,nchan);
            indx   = find(cmbindx);
            if ~complete
              % this realizes the missing combinations to be represented as the
              % conjugate of the corresponding combination across the diagonal
              tmpdat(indx) = reshape(data.(fn{ii})(j,cmbindx(indx),m,k),[numel(indx) 1]);
              tmpdat       = ctranspose(tmpdat);
            end
            tmpdat(indx)    = reshape(data.(fn{ii})(j,cmbindx(indx),m,k),[numel(indx) 1]);
            tmpall(j,:,:,m,k) = tmpdat;
          end % for m
        end % for k
      end % for j

      % replace the data in the old representation with the new representation
      if nrpt>1
        data.(fn{ii}) = tmpall;
      else
        data.(fn{ii}) = reshape(tmpall, [nchan nchan nfrq ntim]);
      end
    end % if numel
  end % for ii

  if ntim>1
    data.dimord = 'chan_chan_freq_time';
  else
    data.dimord = 'chan_chan_freq';
  end

  if nrpt>1
    data.dimord = ['rpt_' data.dimord];
  end

elseif strcmp(current, 'sparse') && strcmp(desired, 'fullfast')
  if contains(data.dimord, 'rpt'),  nrpt = size(data.cumtapcnt,1); else nrpt = 1; end
  if contains(data.dimord, 'freq'), nfrq = length(data.freq); else nfrq = 1; end
  if contains(data.dimord, 'time'), ntim = length(data.time); else ntim = 1; end

  if ~isfield(data, 'label')
    data.label = unique(data.labelcmb(:));
  end

  nchan     = length(data.label);
  ncmb      = size(data.labelcmb,1);
  cmbindx   = zeros(nchan,nchan);

  for k = 1:size(data.labelcmb,1)
    ch1 = find(strcmp(data.label, data.labelcmb(k,1)));
    ch2 = find(strcmp(data.label, data.labelcmb(k,2)));
    if ~isempty(ch1) && ~isempty(ch2)
      cmbindx(ch1,ch2) = k;
    end
  end

  complete = all(cmbindx(:)~=0);

  fn = setdiff(fieldnames(data), {'time' 'freq' 'dimord' 'label' 'cfg'});
  for ii=1:numel(fn)
    if numel(data.(fn{ii})) == nrpt*ncmb*nfrq*ntim
      if nrpt==1
        data.(fn{ii}) = reshape(data.(fn{ii}), [nrpt ncmb nfrq ntim]);
      end

      tmpall = nan(nchan,nchan,nfrq,ntim);

      for k = 1:ntim
        for m = 1:nfrq
          tmpdat = nan(nchan,nchan);
          indx   = find(cmbindx);
          if ~complete
            % this realizes the missing combinations to be represented as the
            % conjugate of the corresponding combination across the diagonal
            tmpdat(indx) = reshape(nanmean(data.(fn{ii})(:,cmbindx(indx),m,k)),[numel(indx) 1]);
            tmpdat       = ctranspose(tmpdat);
          end
          tmpdat(indx)    = reshape(nanmean(data.(fn{ii})(:,cmbindx(indx),m,k)),[numel(indx) 1]);
          tmpall(:,:,m,k) = tmpdat;
        end % for m
      end % for k

      % replace the data in the old representation with the new representation
      if nrpt>1
        data.(fn{ii}) = tmpall;
      else
        data.(fn{ii}) = reshape(tmpall, [nchan nchan nfrq ntim]);
      end
    end % if numel
  end % for ii

  % remove obsolete fields
  data = removefields(data, {'powspctrm', 'labelcmb', 'dof'});

  if ntim>1
    data.dimord = 'chan_chan_freq_time';
  else
    data.dimord = 'chan_chan_freq';
  end

elseif strcmp(current, 'sparsewithpow') && any(strcmp(desired, {'full', 'fullfast'}))
  % recursively call ft_checkdata, but ensure channel order to be the same as the original input.
  origlabelorder = data.label; % keep track of the original order of the channels
  data       = ft_checkdata(data, 'cmbstyle', 'sparse');
  data.label = origlabelorder; % this avoids the labels to be alphabetized in the next call
  data       = ft_checkdata(data, 'cmbstyle', 'full');

end % convert from one to another bivariate representation

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% convert between datatypes
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [source] = chan2source(data)
chanpos = zeros(0,3);
chanlab = cell(0,1);
posunit = [];
if isfield(data, 'elec')
  chanpos = cat(1, chanpos, data.elec.chanpos);
  chanlab = cat(1, chanlab, data.elec.label);
  if isfield(data.elec, 'unit')
    posunit = data.elec.unit;
  end
end
if isfield(data, 'grad')
  chanpos = cat(1, chanpos, data.grad.chanpos);
  chanlab = cat(1, chanlab, data.grad.label);
  if isfield(data.grad, 'unit')
    posunit = data.grad.unit;
  end
end
if isfield(data, 'opto')
  chanpos = cat(1, chanpos, data.opto.chanpos);
  chanlab = cat(1, chanlab, data.opto.label);
  if isfield(data.opto, 'unit')
    posunit = data.opto.unit;
  end
end

fn = fieldnames(data);
fn = setdiff(fn, {'label', 'time', 'freq', 'hdr', 'cfg', 'grad', 'elec', 'dimord', 'unit'}); % remove irrelevant fields
fn(~cellfun(@isempty, regexp(fn, 'dimord$'))) = []; % remove irrelevant (dimord) fields
sel = false(size(fn));
for i=1:numel(fn)
  try
    sel(i) = ismember(getdimord(data, fn{i}), {'chan', 'chan_time', 'chan_freq', 'chan_freq_time', 'chan_chan'});
  end
end
parameter = fn(sel);

% determine the channel indices for which the position is known
[datsel, possel] = match_str(data.label, chanlab);

source = [];
source.pos = chanpos(possel, :);
if ~isempty(posunit)
  source.unit = posunit;
end
for i=1:numel(parameter)
  dat = data.(parameter{i});
  dimord = getdimord(data, parameter{i});
  dimtok = tokenize(dimord, '_');
  for dim=1:numel(dimtok)
    if strcmp(dimtok{dim}, 'chan')
      dat = dimindex(dat, dim, {datsel});
      dimtok{dim} = 'pos';
    end
  end
  dimord = sprintf('%s_', dimtok{:});
  dimord = dimord(1:end-1); % remove the last '_'
  % copy the data to the source representation
  source.(parameter{i})            = dat;
  source.([parameter{i} 'dimord']) = dimord;
end
% copy the descriptive fields, these are necessary for visualising the data in ft_sourceplot
source = copyfields(data, source, {'time', 'freq'});

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% convert between datatypes
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [source] = parcellated2source(data)
if ~isfield(data, 'brainordinate')
  ft_error('projecting parcellated data onto the full brain model geometry requires the specification of brainordinates');
end
% the main structure contains the functional data on the parcels
% the brainordinate sub-structure contains the original geometrical model
source = ft_checkdata(data.brainordinate, 'datatype', 'source');
data   = rmfield(data, 'brainordinate');
if isfield(data, 'cfg')
  source.cfg = data.cfg;
end

fn = fieldnames(data);
fn = setdiff(fn, {'label', 'time', 'freq', 'hdr', 'cfg', 'grad', 'elec', 'dimord', 'unit'}); % remove irrelevant fields
fn(~cellfun(@isempty, regexp(fn, 'dimord$'))) = []; % remove irrelevant (dimord) fields
sel = false(size(fn));
for i=1:numel(fn)
  try
    sel(i) = ismember(getdimord(data, fn{i}), {'chan', 'chan_time', 'chan_freq', 'chan_freq_time', 'chan_chan'});
  end
end
parameter = fn(sel);

fn = fieldnames(source);
sel = false(size(fn));
for i=1:numel(fn)
  tmp = source.(fn{i});
  sel(i) = iscell(tmp) && isequal(sort(tmp(:)), sort(data.label(:)));
end
parcelparam = fn(sel);
if numel(parcelparam)~=1
  ft_error('cannot determine which parcellation to use');
else
  parcelparam = parcelparam{1}(1:(end-5)); % minus the 'label'
end

for i=1:numel(parameter)
  source.(parameter{i}) = unparcellate(data, source, parameter{i}, parcelparam);
end

% copy the descriptive fields, these are necessary for visualising the data in ft_sourceplot
source = copyfields(data, source, {'time', 'freq'});

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% convert between datatypes
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function data = volume2source(data)
if ~isfield(data, 'pos')
  xgrid = 1:data.dim(1);
  ygrid = 1:data.dim(2);
  zgrid = 1:data.dim(3);
  [x, y, z] = ndgrid(xgrid, ygrid, zgrid);
  data.pos = ft_warp_apply(data.transform, [x(:) y(:) z(:)]);
  data.dim = data.dim(1:3); % remove the 4th and further dimensions
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% convert between datatypes
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function data = source2volume(data)

fn = fieldnames(data);
fd = nan(size(fn));
for i=1:numel(fn)
  fd(i) = ndims(data.(fn{i}));
end

if ~isfield(data, 'dim')
  % this part depends on the assumption that the list of positions is describing a full 3D volume in
  % an ordered way which allows for the extraction of a transformation matrix, i.e. slice by slice
  data.dim = pos2dim(data.pos);
  try
    % if the dim is correct, it should be possible to obtain the transform
    ws = warning('off', 'MATLAB:rankDeficientMatrix');
    pos2transform(data.pos, data.dim);
    warning(ws);
  catch
    % remove the incorrect dim
    data = rmfield(data, 'dim');
  end
end

if isfield(data, 'dim')
  data.transform = pos2transform(data.pos, data.dim);
end

% make inside a volume
data = fixinside(data, 'logical');

% only process the fields for which the dimord starts with 'pos_'
fn = fieldnames(data);
fn = setdiff(fn, {'cfg', 'pos', 'dim', 'transform', 'unit', 'xgrid', 'ygrid', 'zgrid', 'tri', 'tet', 'hex'});
keep = false(size(fn));
for k=1:numel(fn)
  keep(k) = ~endsWith(fn{k}, 'dimord') && (startsWith(getdimord(data, fn{k}), 'pos_') || startsWith(getdimord(data, fn{k}), '{pos}_'));
end
fn = fn(keep);

% reshape everything with a 'pos' in the dimord
for k = 1:numel(fn)
  dimord = getdimord(data, fn{k});
  dimsiz = getdimsiz(data, fn{k});
  if startsWith(dimord, 'pos_pos')
    % reshape the first two dimensions
    tmp = getsubfield(data, fn{k});
    data.(fn{k}) = reshape(tmp, [data.dim data.dim dimsiz(3:end)]);
  elseif startsWith(dimord, 'pos')
    % reshape the first dimension
    tmp = getsubfield(data, fn{k});
    data.(fn{k}) = reshape(tmp, [data.dim dimsiz(2:end)]);
  elseif startsWith(dimord, '{pos}')
    % this is a cell array
    tmp = getsubfield(data, fn{k});
    data.(fn{k}) = reshape(tmp, data.dim);
  elseif startsWith(dimord, '{pos_pos}')
    % this is a cell array
    tmp = getsubfield(data, fn{k});
    data.(fn{k}) = reshape(tmp, [data.dim data.dim]);
  elseif contains(dimord, 'pos')
    % the position should always come as the first
    ft_error('unsupported data representation');
  end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% convert between datatypes
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function data = freq2raw(freq)

if isfield(freq, 'powspctrm')
  param = 'powspctrm';
elseif isfield(freq, 'fourierspctrm')
  param = 'fourierspctrm';
else
  ft_error('not supported for this data representation');
end

if strcmp(freq.dimord, 'rpt_chan_freq_time') || strcmp(freq.dimord, 'rpttap_chan_freq_time')
  dat = freq.(param);
elseif strcmp(freq.dimord, 'chan_freq_time')
  dat = freq.(param);
  dat = reshape(dat, [1 size(dat)]); % add a singleton dimension
else
  ft_error('not supported for dimord %s', freq.dimord);
end

nrpt  = size(dat,1);
nchan = size(dat,2);
nfreq = size(dat,3);
ntime = size(dat,4);
data = [];
% create the channel labels like "MLP11@12Hz""
k = 0;
for i=1:nfreq
  for j=1:nchan
    k = k+1;
    data.label{k} = sprintf('%s@%dHz', freq.label{j}, freq.freq(i));
  end
end
% reshape and copy the data as if it were timecourses only
for i=1:nrpt
  data.time{i}  = freq.time;
  data.trial{i} = reshape(dat(i,:,:,:), nchan*nfreq, ntime);
  if any(sum(isnan(data.trial{i}),1)==size(data.trial{i},1))
    tmp = sum(~isfinite(data.trial{i}),1)==size(data.trial{i},1);
    begsmp = find(~tmp,1, 'first');
    endsmp = find(~tmp,1, 'last' );
    data.trial{i} = data.trial{i}(:, begsmp:endsmp);
    data.time{i}  = data.time{i}(begsmp:endsmp);
  end
end

if isfield(freq, 'trialinfo'), data.trialinfo = freq.trialinfo; end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% convert between datatypes
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [tlck] = raw2timelock(data)

data   = ft_checkdata(data, 'hassampleinfo', 'yes');
ntrial = numel(data.trial);
nchan  = numel(data.label);

if ntrial==1
  tlck.time   = data.time{1};
  tlck.avg    = data.trial{1};
  tlck.label  = data.label;
  tlck.dimord = 'chan_time';
  tlck        = copyfields(data, tlck, {'grad', 'elec', 'opto', 'cfg', 'trialinfo', 'topo', 'topodimord', 'topolabel', 'unmixing', 'unmixingdimord'});

else
  % the code below tries to construct a general time-axis where samples of all trials can fall on
  % find the earliest beginning and latest ending
  begtime = min(cellfun(@min, data.time));
  endtime = max(cellfun(@max, data.time));
  % find 'common' sampling rate
  fsample = 1./nanmean(cellfun(@mean, cellfun(@diff,data.time, 'uniformoutput', false)));
  % estimate number of samples
  nsmp = round((endtime-begtime)*fsample) + 1; % numerical round-off issues should be dealt with by this round, as they will/should never cause an extra sample to appear
  % construct general time-axis
  time = linspace(begtime,endtime,nsmp);

  % concatenate all trials
  tmptrial = nan(ntrial, nchan, length(time));

  begsmp = nan(ntrial, 1);
  endsmp = nan(ntrial, 1);
  for i=1:ntrial
    begsmp(i) = nearest(time, data.time{i}(1));
    endsmp(i) = nearest(time, data.time{i}(end));
    tmptrial(i,:,begsmp(i):endsmp(i)) = data.trial{i};
  end

  % construct the output timelocked data
  tlck.trial   = tmptrial;
  tlck.time    = time;
  tlck.dimord  = 'rpt_chan_time';
  tlck.label   = data.label;
  tlck         = copyfields(data, tlck, {'grad', 'elec', 'opto', 'cfg', 'trialinfo', 'sampleinfo', 'topo', 'topodimord', 'topolabel', 'unmixing', 'unmixingdimord'});
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% convert between datatypes
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [data] = timelock2raw(data)
fn = getdatfield(data);
if any(ismember(fn, {'trial', 'individual', 'avg'}))
  % trial, individual and avg (in that order) should be preferred over all other data fields
  % see http://bugzilla.fieldtriptoolbox.org/show_bug.cgi?id=2965#c12
  fn = fn(ismember(fn, {'trial', 'individual', 'avg'}));
end
dimord = cell(size(fn));
for i=1:numel(fn)
  % determine the dimensions of each of the data fields
  dimord{i} = getdimord(data, fn{i});
end
% the fields trial, individual and avg (with their corresponding default dimord) are preferred
if sum(strcmp(dimord, 'rpt_chan_time'))==1
  fn = fn{strcmp(dimord, 'rpt_chan_time')};
  ft_info('constructing trials from "%s"\n', fn);
  dimsiz = getdimsiz(data, fn);
  ntrial = dimsiz(1);
  nchan  = dimsiz(2);
  ntime  = dimsiz(3);
  tmptrial = {};
  tmptime  = {};
  for j=1:ntrial
    tmptrial{j} = reshape(data.(fn)(j,:,:), [nchan, ntime]);
    tmptime{j}  = data.time;
  end
  data       = rmfield(data, fn);
  data.trial = tmptrial;
  data.time  = tmptime;
elseif sum(strcmp(dimord, 'subj_chan_time'))==1
  fn = fn{strcmp(dimord, 'subj_chan_time')};
  ft_info('constructing trials from "%s"\n', fn);
  dimsiz = getdimsiz(data, fn);
  nsubj = dimsiz(1);
  nchan  = dimsiz(2);
  ntime  = dimsiz(3);
  tmptrial = {};
  tmptime  = {};
  for j=1:nsubj
    tmptrial{j} = reshape(data.(fn)(j,:,:), [nchan, ntime]);
    tmptime{j}  = data.time;
  end
  data       = rmfield(data, fn);
  data.trial = tmptrial;
  data.time  = tmptime;
elseif sum(strcmp(dimord, 'chan_time'))==1
  fn = fn{strcmp(dimord, 'chan_time')};
  ft_info('constructing single trial from "%s"\n', fn);
  tmptime  = {data.time};
  tmptrial = {data.(fn)};
  data = rmfield(data, fn);
  data.trial = tmptrial;
  data.time  = tmptime;
else
  ft_error('unsupported data structure');
end
% remove unwanted fields
data = removefields(data, {'avg', 'var', 'cov', 'dimord', 'numsamples' ,'dof'});

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% convert between datatypes
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [data] = chan2freq(data)
data.dimord = [data.dimord '_freq'];
data.freq   = 0;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% convert between datatypes
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [data] = chan2timelock(data)
data.dimord = [data.dimord '_time'];
data.time   = 0;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% convert between datatypes
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [spike] = raw2spike(data)
ft_info('converting raw data into spike data\n');
nTrials 	 = length(data.trial);
[spikelabel] = detectspikechan(data);
spikesel     = match_str(data.label, spikelabel);
nUnits       = length(spikesel);
if nUnits==0
  ft_error('cannot convert raw data to spike format since the raw data structure does not contain spike channels');
end

trialTimes  = zeros(nTrials,2);
for iUnit = 1:nUnits
  unitIndx = spikesel(iUnit);
  spikeTimes  = []; % we dont know how large it will be, so use concatenation inside loop
  trialInds   = [];
  for iTrial = 1:nTrials

    % read in the spike times
    [spikeTimesTrial]    = getspiketimes(data, iTrial, unitIndx);
    nSpikes              = length(spikeTimesTrial);
    spikeTimes           = [spikeTimes; spikeTimesTrial(:)];
    trialInds            = [trialInds; ones(nSpikes,1)*iTrial];

    % get the begs and ends of trials
    hasNum = find(~isnan(data.time{iTrial}));
    if iUnit==1, trialTimes(iTrial,:) = data.time{iTrial}([hasNum(1) hasNum(end)]); end
  end

  spike.label{iUnit}     = data.label{unitIndx};
  spike.waveform{iUnit}  = [];
  spike.time{iUnit}      = spikeTimes(:)';
  spike.trial{iUnit}     = trialInds(:)';

  if iUnit==1, spike.trialtime             = trialTimes; end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% convert between datatypes
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [data] = spike2raw(spike, fsample)

if nargin<2 || isempty(fsample)
  timeDiff = abs(diff(sort([spike.time{:}])));
  fsample  = 1/min(timeDiff(timeDiff>0));
  ft_warning('Desired sampling rate for spike data not specified, automatically resampled to %f', fsample);
end

% get some sizes
nUnits  = length(spike.label);
nTrials = size(spike.trialtime,1);

% preallocate
data.trial(1:nTrials) = {[]};
data.time(1:nTrials)  = {[]};
for iTrial = 1:nTrials

  % make bins: note that the spike.time is already within spike.trialtime
  x = [spike.trialtime(iTrial,1):(1/fsample):spike.trialtime(iTrial,2)];
  timeBins   = [x x(end)+1/fsample] - (0.5/fsample);
  time       = (spike.trialtime(iTrial,1):(1/fsample):spike.trialtime(iTrial,2));

  % convert to continuous
  trialData = zeros(nUnits,length(time));
  for iUnit = 1:nUnits

    % get the timestamps and only select those timestamps that are in the trial
    ts       = spike.time{iUnit};
    hasTrial = spike.trial{iUnit}==iTrial;
    ts       = ts(hasTrial);

    N = histc(ts,timeBins);
    if isempty(N)
      N = zeros(1,length(timeBins)-1);
    else
      N(end) = [];
    end

    % store it in a matrix
    trialData(iUnit,:) = N;
  end

  data.trial{iTrial} = trialData;
  data.time{iTrial}  = time;

end % for all trials

% create the associated labels and other aspects of data such as the header
data.label = spike.label;
data.fsample = fsample;
if isfield(spike,'hdr'), data.hdr = spike.hdr; end
if isfield(spike,'cfg'), data.cfg = spike.cfg; end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% convert between datatypes
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [data] = source2raw(source)

fn = fieldnames(source);
fn = setdiff(fn, {'pos', 'dim', 'transform', 'time', 'freq', 'cfg'});
for i=1:length(fn)
  dimord{i} = getdimord(source, fn{i});
end
sel = strcmp(dimord, 'pos_time');
assert(sum(sel)>0, 'the source structure does not contain a suitable field to represent as raw channel-level data');
assert(sum(sel)<2, 'the source structure contains multiple fields that can be represented as raw channel-level data');
fn     = fn{sel};
dimord = dimord{sel};

switch dimord
  case 'pos_time'
    % add fake raw channel data to the original data structure
    data.trial{1} = source.(fn);
    data.time{1}  = source.time;
    % add fake channel labels
    data.label = {};
    for i=1:size(source.pos,1)
      data.label{i} = sprintf('source%d', i);
    end
    data.label = data.label(:);
    data.cfg = source.cfg;
  otherwise
    % FIXME other formats could be implemented as well
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SUBFUNCTION for detection of channels
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [spikelabel, eeglabel] = detectspikechan(data)

maxRate = 2000; % default on what we still consider a neuronal signal: this firing rate should never be exceeded

% autodetect the spike channels
ntrial = length(data.trial);
nchans  = length(data.label);
spikechan = zeros(nchans,1);
for i=1:ntrial
  for j=1:nchans
    hasAllInts    = all(isnan(data.trial{i}(j,:)) | data.trial{i}(j,:) == round(data.trial{i}(j,:)));
    hasAllPosInts = all(isnan(data.trial{i}(j,:)) | data.trial{i}(j,:)>=0);
    T = nansum(diff(data.time{i}),2); % total time
    fr            = nansum(data.trial{i}(j,:),2) ./ T;
    spikechan(j)  = spikechan(j) + double(hasAllInts & hasAllPosInts & fr<=maxRate);
  end
end
spikechan = (spikechan==ntrial);

spikelabel = data.label(spikechan);
eeglabel   = data.label(~spikechan);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SUBFUNCTION
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [spikeTimes] = getspiketimes(data, trial, unit)
spikeIndx       = logical(data.trial{trial}(unit,:));
spikeCount      = data.trial{trial}(unit,spikeIndx);
spikeTimes      = data.time{trial}(spikeIndx);
if isempty(spikeTimes), return; end
multiSpikes     = find(spikeCount>1);
% get the additional samples and spike times, we need only loop through the bins
[addSamples, addTimes]   = deal([]);
for iBin = multiSpikes(:)' % looping over row vector
  addTimes     = [addTimes ones(1,spikeCount(iBin))*spikeTimes(iBin)];
  addSamples   = [addSamples ones(1,spikeCount(iBin))*spikeIndx(iBin)];
end
% before adding these times, first remove the old ones
spikeTimes(multiSpikes) = [];
spikeTimes              = sort([spikeTimes(:); addTimes(:)]);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SUBFUNCTION
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function print_segmentationinfo(segmentation)
% give feedback about the volume of tissue compartments

fn = fieldnames(segmentation);
fn = setdiff(fn, 'inside');
[indexed, probabilistic] = determine_segmentationstyle(segmentation, fn, segmentation.dim);

% ignore the fields that do not contain a segmentation
sel           = indexed | probabilistic;
fn            = fn(sel);
indexed       = indexed(sel);
probabilistic = probabilistic(sel);

% get the volume of a cubic element
if isfield(segmentation, 'unit')
  voxelvolume = abs(det(segmentation.transform(1:3,1:3)));
  voxelunit   = sprintf('%s^3', segmentation.unit);
  % convert to cubic centimeter, which corresponds to milliliter
  voxelvolume = voxelvolume*ft_scalingfactor(voxelunit, 'cm^3');
  voxelunit   = 'ml';
else
  voxelvolume = 1;
  voxelunit = 'voxels';
end

if all(indexed)

  % give feedback about each of the tissues in each of the volumnes
  totalvolume = prod(segmentation.dim)*voxelvolume;
  for k = 1:numel(fn)
    ft_info('the volume of each of the segmented compartments in "%s" is', fn{k});
    if ~isfield(segmentation, [fn{k} 'label'])
      % this will add the xxxlabel field with default labels
      segmentation = fixsegmentation(segmentation, fn(k), 'indexed');
    end
    tissuelabel = segmentation.([fn{k} 'label']);
    tissueindex = segmentation.(fn{k});
    width = max(cellfun(@length, tissuelabel)); width = max(width, 15);
    summedvolume = 0;
    for m = 1:numel(tissuelabel)
      volume = sum(tissueindex(:)==m)*voxelvolume;
      summedvolume = summedvolume + volume;
      ft_info('%s : %8.0f %s (%6.2f %%)', pad(tissuelabel{m}, width), volume, voxelunit, 100*volume/totalvolume);
    end
    % print the summary of the totals
    ft_info('%s : %8.0f %s (%6.2f %%)', pad('total segmented', width), summedvolume, voxelunit, 100*summedvolume/totalvolume);
    ft_info('%s : %8.0f %s (%6.2f %%)', pad('total volume',    width), totalvolume, voxelunit, 100*totalvolume/totalvolume);

  end

elseif all(probabilistic)

  % give feedback about each of the tissues in each of the volumnes
  width = max(cellfun(@length, fn)); width = max(width, 15);
  totalvolume  = prod(segmentation.dim)*voxelvolume;
  summedvolume = 0;
  ft_info('the volume of each of the segmented compartments is');
  for k = 1:numel(fn)
    tissuelabel = fn{k};
    tissueprobability = segmentation.(tissuelabel);
    volume = sum(tissueprobability(:)*voxelvolume);
    summedvolume = summedvolume + volume;
    ft_info('%s : %8.0f %s (%6.2f %%)', pad(tissuelabel, width), volume, voxelunit, 100*volume/totalvolume);
  end
  % print the summary of the totals
  ft_info('%s : %8.0f %s (%6.2f %%)', pad('total segmented', width), summedvolume, voxelunit, 100*summedvolume/totalvolume);
  ft_info('%s : %8.0f %s (%6.2f %%)', pad('total volume',    width), totalvolume, voxelunit, 100*totalvolume/totalvolume);

end % if all inxexed or probabilistic

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SUBFUNCTION
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function print_voxelinfo(mri)
% give feedback about the size and volume of voxels

if isfield(mri, 'transform') && isfield(mri, 'unit')
  voxelpos = [
    1 1 1
    2 1 1 % shifted by one voxel along 1st dimension
    1 2 1 % shifted by one voxel along 2nd dimension
    1 1 2 % shifted by one voxel along 3rd dimension
    ];

  headpos = ft_warp_apply(mri.transform, voxelpos);

  fprintf('voxel size along 1st dimension (i) : %f %s\n', norm(headpos(2,:)-headpos(1,:)), mri.unit);
  fprintf('voxel size along 2nd dimension (j) : %f %s\n', norm(headpos(3,:)-headpos(1,:)), mri.unit);
  fprintf('voxel size along 3rd dimension (k) : %f %s\n', norm(headpos(4,:)-headpos(1,:)), mri.unit);
  fprintf('volume per voxel                   : %f %s^3\n', abs(det(mri.transform(1:3,1:3))), mri.unit);
end % if hasfield

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function makeGifAndJpg(fig_num)

    saveas(figure(4), 'postion','jpg');
    pause(0.1);
    saveas(figure(5), 'velocity','jpg');
    pause(0.1);
    filename = 'testAnimated6.gif';

    h = figure(fig_num);
    for i = 1:90
        view(4*i+15, 15);
        drawnow

        % Capture the plot as an image
        frame = getframe(h);
        im = frame2im(frame);
        [imind,cm] = rgb2ind(im,256);

        % Write to the GIF File
        if i == 1
          imwrite(imind,cm,filename,'gif', 'Loopcount',inf);
        else
          imwrite(imind,cm,filename,'gif','WriteMode','append');
        end
    end

end

author: LenaShengzhen

repo: AerialRobotics


text:
function quality_test ( )

%*****************************************************************************80
%
%% QUALITY_TEST tests the QUALITY library.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    28 August 2005
%
%  Author:
%
%    John Burkardt
%
  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'QUALITY_TEST\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '  Test the QUALITY library.\n' );

  quality_test_cvt ( );
  quality_test_halton ( );
  quality_test_sphere ( );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'QUALITY_TEST\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
% Entropy: Returns entropy (in bits) of each column of 'X'
% by Will Dwinnell
%
% H = Entropy(X)
%
% H = row vector of calculated entropies (in bits)
% X = data to be analyzed
%
% Example: Measure sample entropy of observations of variables with
%   1, 2, 3 and 4 bits of entropy.
%
% Note: Estimated entropy values are slightly less than true, due to
% finite sample size.
%
% X = ceil(repmat([2 4 8 16],[1e3,1]) .* rand(1e3,4));
% Entropy(X)
%
% Last modified: Nov-12-2006

function H = Entropy(X)

% Establish size of data
[n m] = size(X);

% Housekeeping
H = zeros(1,m);

for Column = 1:m,
    % Assemble observed alphabet
    Alphabet = unique(X(:,Column));

    % Housekeeping
    Frequency = zeros(size(Alphabet));

    % Calculate sample frequencies
    for symbol = 1:length(Alphabet)
        Frequency(symbol) = sum(X(:,Column) == Alphabet(symbol));
    end

    % Calculate sample class probabilities
    P = Frequency / sum(Frequency);

    % Calculate entropy in bits
    % Note: floating point underflow is never an issue since we are
    %   dealing only with the observed alphabet
    H(Column) = -sum(P .* log2(P));
end

% God bless Claude Shannon.

% EOF

author: Sable

repo: mcbench-benchmarks


text:
function visualizeCellVectorsRadial2D(phi_cell)
%VISUALIZECELLS plots the values of cell variable phi
%
% SYNOPSIS:
%  visualizeGradsRadial2D(phi)
%
% PARAMETERS:
%  phi_cell: A CellVector variable
%
% RETURNS:
%  None
%
% EXAMPLE:
%
% SEE ALSO:
%

% Written by Ali A. Eftekhari
% See the license file

% phi_cell=gradientCellTerm(phi);
L = phi_cell.domain.cellcenters.x(end);
x = phi_cell.domain.cellcenters.x;
y = phi_cell.domain.cellcenters.y;
[TH,R] = meshgrid(y, x);
[X,Y] = pol2cart(TH,R);
h = polar([0 2*pi], [0 L]);
delete(h);
quiver(X,Y,phi_cell.xvalue.*cos(TH)-phi_cell.yvalue.*sin(TH), phi_cell.xvalue.*sin(TH)+phi_cell.yvalue.*cos(TH))

author: simulkade

repo: FVTool

save_path: github-repos/MATLAB/simulkade-FVTool


text:
function y = MVL(Ai, Bi, Aj, Bj, t)

% MVL: compute the integral Mij(t) of equation (51)
%
% Usage:    y = MVL(Ai, Bi, Aj, Bj, t)
%
% INPUTS:
%       Ai, Aj, Bi, Bj: matrices of the operators Phi_i(s) and Phi_j(s)
%       t: upper limit of the integral
%
% OUTPUT:
%       y = MVL(t) = int_0^t expm(tau*Ai)*Bi*Bj'*expm(tau*Aj')dtau
%
% SEE ALSO: reference [19]

F = expm([-Ai                                   Bi*Bj';
          zeros( size(Aj',1), size(Ai,2) )      Aj'] * t );

y = expm(Ai*t) * F( 1:size(Ai,1), (1+size(Ai,2)):end );


text:
function [V,E,I] = spline_to_poly(P,C,tol)
  % SPLINE_TO_POLY Evaluate a cubic Bezier spline as a polyline where each
  % segment corresponds to a locally flat segment of the curve up to given
  % tolerance
  %
  % [V,E] = spline_to_poly(P,C,tol)
  %
  % Inputs:
  %   P  #P by dim list of control point locations
  %   C  #C by 4 list of indices into P of cubic Bezier curves
  %   tol  tolerance
  % Outputs:
  %   V  #V by dim list of vertex locations
  %   E  #E by dim list of edge indices into V
  %   I  #E list of indices into 1:#C
  %

  % use transpose for amortized concatenation...
  V = P';
  E = [];
  I = [];
  % consider each cubic
  for c = 1:size(C,1)
    Pc = cubic_flat_eval(P(C(c,:),:),tol);
    J = [C(c,1) size(V,2)+(1:size(Pc,1)-2) C(c,4)];
    Ec = [J(1:end-1);J(2:end)]';
    % use transpose for amortized concatenation...
    E = [E,                    Ec'];
    V = [V,         Pc(2:end-1,:)'];
    I = [I;repmat(c,size(Ec,1),1)];
  end
  V = V';
  E = E';
  [V,~,~,E] = remove_unreferenced(V,E);

end

author: alecjacobson

repo: gptoolbox


text:
function SER = scfde(SP)

numSymbols = SP.FFTsize;
H_channel = fft(SP.channel,SP.FFTsize);

for n = 1:length(SP.SNR),
    tic;
    errCount = 0;

    for k = 1:SP.numRun,
        tmp = round(rand(2,numSymbols));
        tmp = tmp*2 - 1;
        inputSymbols = (tmp(1,:) + i*tmp(2,:))/sqrt(2);

        TxSymbols = [inputSymbols(numSymbols-SP.CPsize+1:numSymbols) inputSymbols];

        RxSymbols = filter(SP.channel, 1, TxSymbols); % Multipath Channel

        tmp = randn(2, numSymbols+SP.CPsize);
        complexNoise = (tmp(1,:) + i*tmp(2,:))/sqrt(2);
        noisePower = 10^(-SP.SNR(n)/10);
        RxSymbols = RxSymbols + sqrt(noisePower)*complexNoise;

        EstSymbols = RxSymbols(SP.CPsize+1:numSymbols+SP.CPsize);
        Y = fft(EstSymbols, SP.FFTsize);

        if SP.equalizerType == 'ZERO'
            Y = Y./H_channel;
        elseif SP.equalizerType == 'MMSE'
            C = conj(H_channel)./(conj(H_channel).*H_channel + 10^(-SP.SNR(n)/10));
            Y = Y.*C;
        end

        EstSymbols = ifft(Y);

        EstSymbols = sign(real(EstSymbols)) + i*sign(imag(EstSymbols));
        EstSymbols = EstSymbols/sqrt(2);

        I = find((inputSymbols-EstSymbols) == 0);
        errCount = errCount + (numSymbols-length(I));
    end
    SER(n,:) = errCount / (numSymbols*SP.numRun);
    [SP.SNR(n) SER(n,:)]
    toc
end


text:
function stla_node_xyz_print ( node_num, node_xyz )

%*****************************************************************************80
%
%% STLA_NODE_XYZ_PRINT prints the node coordinates.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    24 September 2005
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    3D Systems, Inc,
%    Stereolithography Interface Specification,
%    October 1989.
%
%  Parameters:
%
%    Input, integer NODE_NUM, the number of nodes.
%
%    Input, real NODE_XYZ(3,NODE_NUM), the coordinates of the nodes.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '    Node         Coordinates\n' );
  fprintf ( 1, '\n' );

  for node = 1 : node_num

    fprintf ( 1, '  %6d  %14f  %14f  %14f\n', node, node_xyz(1:3,node) );

  end

  return
end


function varths = Get_Varths_swing(zetas,nms,Cont_dt)
%UNTITLED Summary of this function goes here
%   Detailed explanation goes here

zetas2     = zetas.^2; zetas3 = zetas.*zetas2; zetas4 = zetas.*zetas3;

varths = zeros(1,4);
gamms = zeros(1,3);

gamms(1) = -1/24 +zetas4(1)/8 -zetas3(1)/3 +zetas2(1)/4;
gamms(2) = 5/12 - zetas4(1)/4 +zetas3(1)/3 +zetas2(1)/2 -zetas(1);
gamms(3) = (1 + zetas4(1))/8 - zetas2(1)/4;
varths(1) = Cont_dt(nms(1)-1)*gamms(1) + Cont_dt(nms(1))*gamms(2) + Cont_dt(nms(1)+1)*gamms(3);

gamms(1)  = 1/12 - zetas4(1)/8 +(zetas3(1) - zetas2(1))/2;
gamms(2)  = 5/6  +zetas4(1)/4 - 2*zetas3(1)/3;
gamms(3)  = 1/12 -zetas4(1)/8 + zetas3(1)/6;
varths(2) = Cont_dt(nms(1))*gamms(1) + Cont_dt(nms(1)+1)*gamms(2) + Cont_dt(nms(1)+2)*gamms(3);

gamms(1)  = (1 - zetas4(2))/8 +zetas3(2)/3 - zetas2(2)/4;
gamms(2)  = 5/12 + zetas4(2)/4 - zetas3(2)/3 - zetas2(2)/2 + zetas(2);
gamms(3)  = 1/12 - (1 + zetas4(2))/8 + zetas2(2)/4;
varths(3) = Cont_dt(nms(2)-1)*gamms(1) + Cont_dt(nms(2))*gamms(2) + Cont_dt(nms(2)+1)*gamms(3);

gamms(1)  = zetas4(2)/8 - zetas3(2)/2 + zetas2(2)/2;
gamms(2)  = -zetas4(2)/4 + 2*zetas3(2)/3;
gamms(3)  = zetas4(2)/8 - zetas3(2)/6;
varths(4) = Cont_dt(nms(2))*gamms(1) + Cont_dt(nms(2)+1)*gamms(2) + Cont_dt(nms(2)+2)*gamms(3);

end


text:
function copyDataFromCluster(localLoc, remoteLoc, clusterHost)
%COPYDATAFROMCLUSTER Copies files or directories from a location
% on a remote host to the local machine.

% Copyright 2006 The MathWorks, Inc.

% Use scp to copy files.
copyCmd = sprintf('scp -r "%s:%s" %s', clusterHost, remoteLoc, localLoc);

[s, r] = system(copyCmd);
if s ~= 0
    error('distcomp:scheduler:FailedRemoteOperation', ...
        ['Failed to copy files from "%s" on the host "%s"\n' ...
        'to "%s".\n' ...
        'Command Output:\n' ...
        '"%s"\n' ...
        ], remoteLoc, clusterHost, localLoc, r);
end


text:
function [fx,dfdx,dfdp] = f_AR(x,theta,u,in)
% AR(1) evolution function
fx = x;
dfdx = eye(length(x));
dfdp = [];

author: MBB-team

repo: VBA-toolbox


text:
function test_failed=test_realout
%TEST_REALOUT  Test if functions produce real-valued output

  test_failed=0;

  disp(' ===============  TEST_REALOUT ================');

  a = 7;
  M = 19;
  W = 3;
  L = a*M*4;
  Nwil = L/(2*M);
  Nmd  = L/M;
  Ngab = L/a;

  Nfft = 19;

  test_failed=realhelper(test_failed,'dwilt',randn(L,W),randn(L,1),M);
  test_failed=realhelper(test_failed,'dwilt',randn(L,W),randn(2*M,1),M);

  test_failed=realhelper(test_failed,'wmdct',randn(L,W),randn(L,1),M);
  test_failed=realhelper(test_failed,'wmdct',randn(L,W),randn(2*M,1),M);

  test_failed=realhelper(test_failed,'idwilt',randn(2*M,Nwil),randn(L,1));
  test_failed=realhelper(test_failed,'idwilt',randn(2*M,Nwil),randn(2*M,1));

  test_failed=realhelper(test_failed,'iwmdct',randn(M,Nmd),randn(L,1));
  test_failed=realhelper(test_failed,'iwmdct',randn(M,Nmd),randn(2*M,1));

  test_failed=realhelper(test_failed,'gabdual',randn(L,1),a,M);
  test_failed=realhelper(test_failed,'gabdual',randn(L,1),a,M);
  test_failed=realhelper(test_failed,'gabtight',randn(M,1),a,M);
  test_failed=realhelper(test_failed,'gabtight',randn(M,1),a,M);

  test_failed=realhelper(test_failed,'dcti',randn(Nfft,1));
  test_failed=realhelper(test_failed,'dctii',randn(Nfft,1));
  test_failed=realhelper(test_failed,'dctiii',randn(Nfft,1));
  test_failed=realhelper(test_failed,'dctiv',randn(Nfft,1));
  test_failed=realhelper(test_failed,'dsti',randn(Nfft,1));
  test_failed=realhelper(test_failed,'dstii',randn(Nfft,1));
  test_failed=realhelper(test_failed,'dstiii',randn(Nfft,1));
  test_failed=realhelper(test_failed,'dstiv',randn(Nfft,1));

  test_failed=realhelper(test_failed,'pfilt',randn(L,1),randn(L,1));

  c=dgtreal(randn(L,W),randn(L,1),a,M);
  test_failed=realhelper(test_failed,'idgtreal',c,randn(L,1),a,M);
  test_failed=realhelper(test_failed,'idgtreal',c,randn(M,1),a,M);

  c=fftreal(randn(Nfft,1));
  test_failed=realhelper(test_failed,'ifftreal',c,Nfft);

function test_failed=realhelper(test_failed,funname,varargin)

  outres=feval(funname,varargin{:});
  res=~isreal(outres);

  if res>0
      outres
  end;
  [test_failed,fail]=ltfatdiditfail(res,test_failed);

  fprintf('REAL %s %i %s\n',funname,res,fail);

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
clear; close all; clc;

fx = 1;
x = 0;

for i = 1:5
    x = [x, x(i)+0.5];
    fx = [fx, fx(i)*1.5];
end

stem(x, fx,'linewidth',2)
grid on;
set(gca,'ytick',round(fx,2))
xlabel('$$x$$','interpreter','latex');
ylabel('$$f(x)$$','interpreter','latex');
set(gca,'fontsize',11)

author: angeloyeo

repo: gongdols

save_path: github-repos/MATLAB/angeloyeo-gongdols

path: github-repos/MATLAB/angeloyeo-gongdols/gongdols-7be9fbd988dec6edab1dc881cb22d63e6f69398d/기초수학/e^x는 왜 미분해도 e^x일까/pics.m


text:
function [ t_start, y_start ] = p26_start ( neqn )

%*****************************************************************************80
%
%% P26_START returns the starting point for problem p26.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    23 February 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer NEQN, the number of equations.
%
%    Output, real T_START, Y_START(NEQN,1), the initial data.
%
  y_start = zeros ( neqn, 1 );

  t_start = 0.0;
  y_start = [ 0.0; 0.0 ];

  return
end


text:
% This subroutine "overlay.m" is called from varios
% program (view_*.m, subcata.m). It plots an overlay
% of coastlines, faults, earthquakes etc on a map.
% This file should be customized for each region

%  Stefan Wiemer   11/94

report_this_filefun(mfilename('fullpath'));

global main mainfault faults coastline vo

figure_w_normalized_uicontrolunits(map)
ax = axis;
if ax(1) < min(coastline(:,1)) | ax(2) > min(coastline(:,1))
    disp([' Regetting coastline ... ']);
    struc = gshhs2('gshhs_h.b',[s4_south s3_north],[s2_west s1_east] );
    s2_west = ax(1);s1_east = ax(2);s3_north = ax(4);s4_south = ax(3);
end
%delete(mapplot);delete(plovo);delete(epimax);

[i,j ] = size(struc);
lat = [];
lon = [];
lat2 = [];
lon2 = [];
figure_w_normalized_uicontrolunits(map)
hold on
for i = 1:j
    if struc(i).tag(1:4) == 'land' | struc(i).tag(1:4) == 'isla'
        lat = [lat ; struc(i).lat ];
        lon = [lon ; struc(i).long];
        l1 = [  struc(i).lat ];
        l2 =   [  struc(i).long];
        l = isnan(l1); l1(l) = [];
        l = isnan(l2); l2(l) = [];
        fill(l2,l1,[0.9 0.9 0.9]);
        hold on
    end

    if struc(i).tag(1:4) == 'lake'
        lat2 = [struc(i).lat ];
        lon2 = [struc(i).long];
        l1 = [  struc(i).lat ];
        l2 =   [  struc(i).long];
        l = isnan(l1); l1(l) = [];
        l = isnan(l2); l2(l) = [];
        fill(l2,l1,[0.85 0.85 1]);
        hold on
    end
end

coastline = [lon lat];

hold on
if exist('coastline', 'var')
    if isempty(coastline) ==  0
        mapplot = plot(coastline(:,1),coastline(:,2));
        set(mapplot,'LineWidth',1.0,'Color',[0  0      0 ])
    end
end
if exist('vo', 'var')
    if isempty(vo) ==  0
        plovo = plot(vo.Longitude,vo.Latitude,'^r');
        set(plovo,'LineWidth',1.5,'MarkerSize',6,...
            'MarkerFaceColor','w','MarkerEdgeColor','r');
    end
end

% plot the well location
if exist('well', 'var')
    if isempty(well) ==  0
        i = find(well(:,1) == inf);
        plowe = plot(well(i+1,1),well(i+1,2),'d');
        set(plowe,'LineWidth',1.5,'MarkerSize',6,...
            'MarkerFaceColor','k','MarkerEdgeColor','k');
    end
end

%plot main faultline

if exist('mainfault', 'var')
    if isempty(mainfault) == 0
        plo3 = plot(mainfault(:,1),mainfault(:,2),'b');
        set(plo3,'LineWidth',3.0)
    end  % if exist mainfault
end

%
% plot big earthquake epicenters with a 'x' and the data/magnitude
%
if exist('maepi', 'var')
    if isempty(maepi) == 0
        epimax = plot(maepi(:,1),maepi(:,2),'hm');
        set(epimax,'LineWidth',1.5,'MarkerSize',12,...
            'MarkerFaceColor','y','MarkerEdgeColor','k')

        stri2 = [];
        for i = 1:length(maepi(:,1));
            s = sprintf('   %3.1f M=%3.1f',maepi(i,3),maepi(i,6));
            stri2 = [stri2 ; s];
        end   % for i
        te1 = text(maepi(:,1),maepi(:,2),stri2);
        set(te1,'FontWeight','normal','Color','k','FontSize',9,'Clipping','on')
    end  %  if length(maepi)
end  %  if length(maepi)

%plot mainshock(s)
%
if exist('main', 'var')
    if isempty(main) == 0
        plo1 = plot(main(:,1),main(:,2),'*k');
        set(plo1,'MarkerSize',12,'LineWidth',2.0)
    end  % if main
end

if exist('faults', 'var')
    if isempty(faults) == 0
        plo4 = plot(faults(:,1),faults(:,2),'k');
        set(plo4,'LineWidth',0.2)
    end  % if exist faults
end
%axis([ s2_west s1_east s4_south s3_north])

author: CelsoReyes


text:
function n = dim(A)
%DIM          Dimension of a square matrix
%
%    n = dim(A)
%

% written  05/21/09     S.M. Rump
%

  n = A.size(1);
  if ( length(A.size)~=2 ) | ( n ~= A.size(2) )
    error('function dim called with non-square matrix')
  end;

repo: OpenMAS


text:
function toeplitz_cholesky_test ( )

%*****************************************************************************80
%
%% TOEPLITZ_CHOLESKY_TEST tests the TOEPLITZ_CHOLESKY library.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    13 November 2012
%
%  Author:
%
%    John Burkardt
%
  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TOEPLITZ_CHOLESKY_TEST:\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '  Test the TOEPLITZ_CHOLESKY library.\n' );

  toeplitz_cholesky_test01 ( );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TOEPLITZ_CHOLESKY_TEST:\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  return
end


text:
% using MDL and AIC to find source number, and MUSIC to find DOAs
clear(); close all;

wavelength = 1; % normalized
d = wavelength / 2;
design_ula = design_array_1d('ula', 12, d);
doas = linspace(-pi/8, pi/3, 8);
power_source = 1;
power_noise = 1;
snapshot_count = 1000;
source_count = length(doas);

% stochastic (unconditional) model
[~, R] = snapshot_gen_sto(design_ula, doas, wavelength, snapshot_count, power_noise, power_source);

% source number detection
[~, l] = eig(0.5*(R+R'), 'vector');
l = flipud(l);
n_mdl = sn_mdl(l, design_ula.element_count, snapshot_count);
n_aic = sn_aic(l, design_ula.element_count, snapshot_count);
fprintf('There are %d sources.\n', source_count);
fprintf('# of sources estimated by MDL = %d\n', n_mdl);
fprintf('# of sources estimated by AIC = %d\n', n_aic);

% normal MUSIC
tic;
sp_normal = music_1d(R, source_count, design_ula, wavelength, 180, 'RefineEstimates', true);
toc;
sp_normal.true_positions = doas;
fprintf('[MUSIC] Estimated DOAs:\n');
disp(sp_normal.x_est);
plot_sp(sp_normal, 'title', 'Normal MUSIC', 'PlotType', 'Polar');

% root MUSIC
tic;
sp_root = rmusic_1d(R, source_count, 2*pi*design_ula.element_spacing/wavelength);
toc;
sp_root.true_positions = doas;
fprintf('[Root-MUSIC] Estimated DOAs:\n');
disp(sp_root.x_est);
plot_sp(sp_root, 'Title', 'Root-MUSIC');

author: morriswmz

repo: doa-tools

save_path: github-repos/MATLAB/morriswmz-doa-tools


text:
% eegplugin_ascinstep() - EEGLAB plugin for importing ASC INStep files.
%
% Usage:
%   >> eegplugin_ascinstep(menu);
%   >> eegplugin_ascinstep(menu, trystrs, catchstrs);
%
% Inputs:
%   menu       - [float]  EEGLAB menu handle
%   trystrs    - [struct] "try" strings for menu callbacks. See notes on EEGLab plugins.
%                (http://www.sccn.ucsd.edu/eeglab/contrib.html)
%   catchstrs  - [struct] "catch" strings for menu callbacks. See notes on EEGLab plugins.
%                (http://www.sccn.ucsd.edu/eeglab/contrib.html)
%
%
% Notes:
%   This plugins consist of the following Matlab files:
%   pop_loadeep.m           pop_loadeep_avg.m
%   loadeep.m               loadeep_avg.m
%   read_eep_cnt.m          read_eep_avr.m
%   read_eep_cnt.mexglx     read_eep_avr.mexglx
%   read_eep_cnt.dll        read_eep_avr.dll
%
% Authors: Arnaud Delorme, SCCN, INC, UCSD, 2006

% Copyright (C) 2006 Arnaud Delorme
%
% This program 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 2 of the License, or
% (at your option) any later version.
%
% This program 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 this program; if not, write to the Free Software
% Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

function vers = eegplugin_ascinstep(fig, trystrs, catchstrs)

    vers = 'eepimport1.02';

    if nargin < 3
        error('eegplugin_ascinstep requires 3 arguments');
    end;

    % add folder to path
    % ------------------
    if ~exist('pop_loadeep')
        p = which('eegplugin_ascinstep.m');
        p = p(1:findstr(p,'eegplugin_ascinstep.m')-1);
        addpath( p );
    end;

    % find import data menu
    % ---------------------
    menu = findobj(fig, 'tag', 'import data');

    % menu callbacks
    % --------------
    comcnt = [ trystrs.no_check '[EEG LASTCOM] = pop_loadascinstep;' catchstrs.new_and_hist ];

    uimenu( menu, 'label', 'From INStep .ASC file', 'callback', comcnt, 'separator', 'on' );

author: PatternRecognition

repo: OpenBMI


text:
function [ fx, fy ] = f13_f1 ( n, x, y )

%*****************************************************************************80
%
%% F13_F1 returns first derivatives of function f13.
%
%  Discussion:
%
%    Currently, the derivative information is of no interest to me.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    01 August 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of evaluation points.
%
%    Input, real X(N,1), Y(N,1), the evalution points.
%
%    Output, real FX(N,1), FY(N,1), the derivative values.
%
  fx(1:n,1) = zeros ( n, 1 );
  fy(1:n,1) = zeros ( n, 1 );

  return
end


text:
%% This script is used to generate the Brain Network for user to perform further analysis by themselves.

clear all
addpath(genpath(pwd));
input_folder='./data/'; % this is the directory of the time course data for all subjects, each file for each subject, similar as you prepared for GUI
output_dir='./Generated_BrainNet/'; %this is the output directory;
meth_Net='aHOFC'; %Here you can set the brain network construction method;
%["PC","aHOFC","tHOFC","SR","WSR","SLR","SGR","WSGR","GSR","SSGSR","dHOFC"];
switch meth_Net
    case {'SR','WSR','GSR'}
        lambda=0.01:0.01:0.1; %User can change the parameter range by themselves;
    case {'SLR','SGR','WSGR','SSGSR'}
        lambda_1=0.01:0.01:0.1; %User can change the parameter range by themselves;
        lambda_2=0.01:0.01:0.1;
    case 'dHOFC'
        C=100:100:800;
        W=20:10:70;
        s=1; %step size is usually set to 1 or 2;
    case {'PC','tHOFC','aHOFC'}

end

dirOutput=dir(fullfile(input_folder,'*.txt'));
fileName={dirOutput.name}';
%folder={dirOutput.folder}';
for i=1:length(fileName)
    BOLD{i,1}=load([input_folder,'/',fileName{i}]);
end
%label=importdata(label_input);

[~,nROI]=size(BOLD{1});
nSubj=length(BOLD);

fprintf('Begin network construction\n');

switch meth_Net
    case 'PC'       % Pearson's correlation
        BrainNet{1}=PC(BOLD);
    case 'tHOFC'    % Topographical high-order FC
        BrainNet{1}=tHOFC(BOLD);
    case 'aHOFC'    % Associated high-order FC
        BrainNet{1}=aHOFC(BOLD);
    case 'SR'       % Sparse representation
        parfor i=1:length(lambda)
            BrainNet{i}=SR(BOLD,lambda(i));
        end

    case 'WSR'      % PC weighted SR
        parfor i=1:length(lambda)
            BrainNet{i}=WSR(BOLD,lambda(i));
        end

    case 'SLR'      % Sparse low-rank representation
        lambda1=lambda_1;
        lambda2=lambda_2;
        num_lambda1=length(lambda1);
        num_lambda2=length(lambda2);
        parfor i=1:num_lambda1
            for j=1:num_lambda2
                BrainNet{i,j}=SLR(BOLD,lambda1(i),lambda2(j));
            end
        end
        BrainNet=reshape(BrainNet,1,num_lambda1*num_lambda2);

    case 'SGR'      % Sparse group representation
        lambda1=lambda_1;
        lambda2=lambda_2;
        num_lambda1=length(lambda1);
        num_lambda2=length(lambda2);
        parfor i=1:num_lambda1
            for j=1:num_lambda2
                BrainNet{i,j}=SGR(BOLD,lambda1(i),lambda2(j));
            end
        end

        BrainNet=reshape(BrainNet,1,num_lambda1*num_lambda2);

    case 'WSGR'     % PC weighted SGR
        lambda1=lambda_1; % parameter for sparsity
        lambda2=lambda_2; % parameter for group sparsity
        num_lambda1=length(lambda1);
        num_lambda2=length(lambda2);

        parfor i=1:num_lambda1
            for j=1:num_lambda2
                BrainNet{i,j}=WSGR(BOLD,lambda1(i),lambda2(j));
            end
        end
        BrainNet=reshape(BrainNet,1,num_lambda1*num_lambda2);

    case 'GSR'      % Group sparse representation
        parfor i=1:length(lambda)
            BrainNet{i}=GSR(BOLD,lambda(i));
        end

    case 'SSGSR'    % Strength and Similarity guided GSR
        %lambda1=lambda_1(1:6); % parameter for group sparsity
        lambda1=lambda_1;
        lambda2=lambda_2; % parameter for inter-subject LOFC-pattern similarity
        num_lambda1=length(lambda1);
        num_lambda2=length(lambda2);
        parfor i=1:num_lambda1
            for j=1:num_lambda2
                BrainNet{i,j}=SSGSR(BOLD,lambda1(i),lambda2(j));
            end
        end
        BrainNet=reshape(BrainNet,1,num_lambda1*num_lambda2);

    case 'dHOFC'    % Dynamic high-order FC
        num_W=length(W);
        num_C=length(C);
        parfor i=1:num_W % number of clusters
            for j=1:num_C
                [BrainNet{i,j},IDX{i,j}]=dHOFC(BOLD,W(i),s,C(j));
            end
        end
        BrainNet=reshape(BrainNet,1,num_W*num_C);
end
save(char(strcat(output_dir,meth_Net,'.mat')),'BrainNet','-v7.3');
% All the generated networks (and those resulted from different parameters) are saved as a cell array, nROIxnROIxnSubject in each cell, different cells for results with different combinations of parameters
fprintf('Network construction finished\n');

author: zzstefan

repo: BrainNetClass


text:
function [ y, m, d, ierror ] = ymd_check_julian ( y, m, d )

%*****************************************************************************80
%
%% YMD_CHECK_JULIAN checks a Julian YMD date.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    09 December 2008
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer Y, integer M, integer D, the YMD date.
%
%    Output, integer Y, integer M, integer D, the YMD date, which may
%    be corrected if necessary and possible.
%
%    Output, integer IERROR, is 0 if the date is legal.
%

%
%  Check the month.
%
  [ y, m, ierror ] = ym_check_julian ( y, m );

  if ( ierror ~= 0 )
    return
  end
%
%  Check the day.
%
  [ y, m, d ] = day_borrow_julian ( y, m, d );

  [ y, m, d ] = day_carry_julian ( y, m, d );

  return
end


text:
% Author: Housam Binous

% Dynamic and control of a tank using the genetic algorithm toolbox

% National Institute of Applied Sciences and Technology, Tunis, TUNISIA

% Email: binoushousam@yahoo.com

function f=obj2(par)

global sys

% real PID controller's transfer function

num=[par(2) 1];
den=[par(2) 0];

pid1=tf(num,den);

num=[par(3) 1];
den=[0.1*par(3) 1];

pid2=tf(num,den);

pid=par(1)*pid1*pid2;

% response to change in set point
% set point for tank's height is equal to 2

sys_series=series(pid,sys);

sys_controlled=feedback(sys_series,1);

u=2*ones(501,1);

[y1 t3]=lsim(sys_controlled,u,0:0.01:5);

% objective function to be minimized by
% genetic algorithm in order to obtain
% real PID controller's parameters

f=0;

for i=1:201
    f=f+abs(2-y1(i))*t3(i);
end


text:
function g = ceil(f)
%CEIL   Pointwise ceiling of a CHEBFUN.
%   G = CEIL(F) returns the CHEBFUN G such that G(x) = CEIL(F(x)) for each x in
%   F.domain.
%
%   If F is complex, then the G = CEIL(REAL(F)) + 1i*CEIL(IMAG(F)).
%
% See also FLOOR, ROUND, FIX.

% Copyright 2017 by The University of Oxford and The Chebfun Developers.
% See http://www.chebfun.org/ for Chebfun information.

% Deal with unbounded functions:
if ( ~isfinite(f) )
    error('CHEBFUN:CHEBFUN:ceil:inf', ...
        'Ceil is not defined for functions which diverge to infinity.');
end

g = -floor(-f);

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
clear all, close all, clc

load ../DATA/allFaces.mat

% We use the first 36 people for training data
trainingFaces = faces(:,1:sum(nfaces(1:36)));
avgFace = mean(trainingFaces,2);  % size n*m by 1;

% Compute eigenfaces on mean-subtracted training data
X = trainingFaces-avgFace*ones(1,size(trainingFaces,2));
[U,S,V] = svd(X,'econ');

imagesc(reshape(avgFace,n,m)) % Plot avg face
imagesc(reshape(U(:,1),n,m))  % Plot first eigenface

%% Now show eigenface reconstruction of image that was omitted from test set

testFace = faces(:,1+sum(nfaces(1:36))); % First face of person 37
imagesc(reshape(testFace,n,m))

testFaceMS = testFace - avgFace;
for r=[25 50 100 200 400 800 1600]
    reconFace = avgFace + (U(:,1:r)*(U(:,1:r)'*testFaceMS));
    imagesc(reshape(reconFace,n,m))
end

%% Project person 2 and 7 onto PC5 and PC6

P1num = 2;  % Person number 2
P2num = 7;  % Person number 7

P1 = faces(:,1+sum(nfaces(1:P1num-1)):sum(nfaces(1:P1num)));
P2 = faces(:,1+sum(nfaces(1:P2num-1)):sum(nfaces(1:P2num)));

P1 = P1 - avgFace*ones(1,size(P1,2));
P2 = P2 - avgFace*ones(1,size(P2,2));

PCAmodes = [5 6];    % Project onto PCA modes 5 and 6
PCACoordsP1 = U(:,PCAmodes)'*P1;
PCACoordsP2 = U(:,PCAmodes)'*P2;

plot(PCACoordsP1(1,:),PCACoordsP1(2,:),'kd')
plot(PCACoordsP2(1,:),PCACoordsP2(2,:),'r^')

author: dynamicslab


text: %ISNAN Datafile overload

author: marianux

repo: ecg-kit

save_path: github-repos/MATLAB/marianux-ecg-kit