function [dat] = read_data(filename, varargin);

% READ_DATA is a wrapper around different EEG/MEG file importers
% directly supported formats are CTF, Neuromag, EEP, BrainVision,
% Neuroscan and Neuralynx.
%
% Use as
%   dat = read_data(filename, ...)
%
% Additional options should be specified in key-value pairs and can be
%   'header'         header structure, see READ_HEADER
%   'begsample'      first sample to read
%   'endsample'      last sample to read
%   'begtrial'       first trial to read, mutually exclusive with begsample+endsample
%   'endtrial'       last trial to read, mutually exclusive with begsample+endsample
%   'chanindx'       list with channel indices to read
%   'checkboundary'  boolean, whether to check for reading segments over a trial boundary
%   'dataformat'     string
%   'headerformat'   string
%
% This function returns a 2-D matrix of size Nchans*Nsamples for
% continuous data when begevent and endevent are specified, or a 3-D
% matrix of size Nchans*Nsamples*Ntrials for epoched or trial-based
% data when begtrial and endtrial are specified.
%
% See also READ_HEADER, READ_EVENT

% Copyright (C) 2003-2006, Robert Oostenveld, F.C. Donders Centre
%
% $Log: read_data.m,v $
% Revision 1.7  2006/12/04 10:38:12  roboos
% added support for ns_avg
% fixed long-outstanding problem for reading multiple trials from ns_eeg data
%
% Revision 1.6  2006/09/18 21:47:54  roboos
% implemented support for fcdc_matbin, i.e. a dataset consisting of a matlab file with header and events and a seperate binary datafile
% added smart default for reading all data when no selection given
%
% Revision 1.5  2006/09/18 14:22:54  roboos
% implemented support for 4D-BTi dataformat
%
% Revision 1.4  2006/08/28 10:12:22  roboos
% use seperate filetype_check_extension instead of (missing) check_extension subfunction (thanks to Thilo for reporting this bug)
%
% Revision 1.3  2006/06/19 10:31:47  roboos
% translate continuous option into checkboundary (complements), added documentation
%
% Revision 1.2  2006/06/19 08:14:17  roboos
% updated documentation
%
% Revision 1.1  2006/06/07 09:32:20  roboos
% new implementation based on the read_fcdc_xxx functions, now with
% variable (key-val) input arguments, changed the control structure
% in the rpobram (switch instead of ifs), allow the user to specify
% the file format, allow the user to specify either a sample selection
% or a block selection. The reading functionality should not have
% changed compared to the read_fcdc_xxx versions.
%

% test whether the file or directory exists
if ~exist(filename)
  error(sprintf('file or directory ''%s'' does not exist', filename));
end

% for backward compatibility support with read_fcdc_data
if nargin==4 && isstruct(varargin{1})
  % input appears to be filename, hdr, begsample, endsample
  varargin = {'header', varargin{1}, 'begsample', varargin{2}, 'endsample', varargin{3}};
elseif nargin==5 && isstruct(varargin{1})
  % input appears to be filename, hdr, begsample, endsample, chanindx
  varargin = {'header', varargin{1}, 'begsample', varargin{2}, 'endsample', varargin{3}, 'chanindx', varargin{4}};
elseif nargin==6 && isstruct(varargin{1})
  % input appears to be filename, hdr, begsample, endsample, chanindx, continuous
  continuous    = varargin{5};
  checkboundary = (continuous==0);  % the default is to check in case the file is _not_ continuous
  varargin = {'header', varargin{1}, 'begsample', varargin{2}, 'endsample', varargin{3}, 'chanindx', varargin{4}, 'checkboundary', checkboundary};
end

% get the options
hdr           = keyval('header',        varargin);
begsample     = keyval('begsample',     varargin);
endsample     = keyval('endsample',     varargin);
begtrial      = keyval('begtrial',      varargin);
endtrial      = keyval('endtrial',      varargin);
chanindx      = keyval('chanindx',      varargin);
checkboundary = keyval('checkboundary', varargin);
dataformat    = keyval('dataformat',    varargin);
headerformat  = keyval('headerformat',  varargin);

% determine the filetype
if isempty(dataformat)
  dataformat = filetype(filename);
end

switch dataformat
  case '4d_pdf'
    datafile   = filename;
    headerfile = [datafile '.m4d'];
    sensorfile = [datafile '.xyz'];
  case {'4d_m4d', '4d_xyz'}
    datafile   = filename(1:(end-4)); % remove the extension
    headerfile = [datafile '.m4d'];
    sensorfile = [datafile '.xyz'];
  case 'ctf_ds'
    % convert CTF filename into filenames
    [path, file, ext] = fileparts(filename);
    headerfile = fullfile(filename, [file '.res4']);
    datafile   = fullfile(filename, [file '.meg4']);
  case 'ctf_meg4'
    [path, file, ext] = fileparts(filename);
    headerfile = fullfile(path, [file '.res4']);
    datafile   = fullfile(path, [file '.meg4']);
  case 'ctf_res4'
    [path, file, ext] = fileparts(filename);
    headerfile = fullfile(path, [file '.res4']);
    datafile   = fullfile(path, [file '.meg4']);
  case 'brainvision_vhdr'
    [path, file, ext] = fileparts(filename);
    headerfile = fullfile(path, [file '.vhdr']);
    if exist(fullfile(path, [file '.eeg']))
      datafile   = fullfile(path, [file '.eeg']);
    elseif exist(fullfile(path, [file '.seg']))
      datafile   = fullfile(path, [file '.seg']);
    elseif exist(fullfile(path, [file '.dat']))
      datafile   = fullfile(path, [file '.dat']);
    end
  case 'brainvision_eeg'
    [path, file, ext] = fileparts(filename);
    headerfile = fullfile(path, [file '.vhdr']);
    datafile   = fullfile(path, [file '.eeg']);
  case 'brainvision_seg'
    [path, file, ext] = fileparts(filename);
    headerfile = fullfile(path, [file '.vhdr']);
    datafile   = fullfile(path, [file '.seg']);
  case 'brainvision_dat'
    [path, file, ext] = fileparts(filename);
    headerfile = fullfile(path, [file '.vhdr']);
    datafile   = fullfile(path, [file '.dat']);
  case 'fcdc_matbin'
    [path, file, ext] = fileparts(filename);
    headerfile = fullfile(path, [file '.mat']);
    datafile   = fullfile(path, [file '.bin']);
  otherwise
    % convert filename into filenames, assume that the header and data are the same
    datafile   = filename;
    headerfile = filename;
end

if ~strcmp(filename, datafile)
  filename   = datafile;                % this function will read the data
  dataformat = filetype(filename);      % update the filetype
end

% for backward compatibility, default is to check when it is not continous
if isempty(checkboundary)
  checkboundary = ~keyval('continuous', varargin);
end

% read the header if it is not provided
if isempty(hdr)
  hdr = read_header(filename, 'headerformat', headerformat);
end

% set the default channel selection, which is all channels
if isempty(chanindx)
  chanindx = 1:hdr.nChans;
end

% read untill the end of the file if the endsample is "inf"
if any(isinf(endsample)) && any(endsample>0)
  endsample = hdr.nSamples*hdr.nTrials;
end

% test whether the file exists
if ~exist(filename)
  error(sprintf('file ''%s'' does not exist', filename));
end

% test whether the requested channels can be accomodated
if min(chanindx)<1 || max(chanindx)>hdr.nChans
  error('selected channels are not present in the data');
end

% test whether the requested data segment is not outside the file
if any(begsample<1) || any(endsample>(hdr.nSamples*hdr.nTrials))
  error('cannot read data before the start or after the end of the file');
end

requesttrials  = isempty(begsample) && isempty(endsample);
requestsamples = isempty(begtrial)  && isempty(endtrial);

if isempty(begsample) && isempty(endsample) && isempty(begtrial) && isempty(endtrial)
  % neither samples nor trials are specified, set the defaults to read the complete data trial-wise (also works for continuous)
  requestsamples = 0;
  requesttrials  = 1;
  begtrial       = 1;
  endtrial       = hdr.nTrials;
end

% set the default, which is to assume that it is should check for boundaries when samples are requested
if isempty(checkboundary) && requesttrials
  checkboundary = 0;
elseif isempty(checkboundary) && requestsamples
  checkboundary = 1;
end

if requesttrials && requestsamples
  error('you cannot select both trials and samples at the same time');
elseif requesttrials
  % this allows for support using a continuous reader
  begsample = (begtrial-1)*hdr.nSamples + 1;
  endsample = (endtrial  )*hdr.nSamples;
elseif requestsamples
  % this allows for support using a trial-based reader
  begtrial = floor((begsample-1)/hdr.nSamples)+1;
  endtrial = floor((endsample-1)/hdr.nSamples)+1;
else
  error('you should either specify begin/end trial or begin/end sample');
end

% test whether the requested data segment does not extend over a discontinuous trial boundary
if checkboundary && hdr.nTrials>1
  if begtrial~=endtrial
    error('requested data segment extends over a discontinuous trial boundary');
  end
end

% translate the name of the dataset directory into the name of the data file
if filetype(filename, 'ctf_ds')
  [path, file, ext] = fileparts(filename);
  filename = fullfile(filename, [file '.meg4']);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% read the data with the low-level reading function
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
switch dataformat

  case {'4d_pdf', '4d_m4d', '4d_xyz'}
    % 4D/BTi MEG data is multiplexed, can be epoched/discontinuous
    offset        = begsample-1;
    numsamples    = endsample-begsample+1;
    [fid,message] = fopen(datafile,'rb','ieee-be');
    % determine the type and size of the samples
    switch hdr.orig.fileformat
      case 1
        samplesize = 2;
        sampletype = 'short';
        % FIXME the data probably should be calibrated
        warning('integer data is not yet calilbrated');
      case 2
        samplesize = 4;
        sampletype = 'long';
        % FIXME the data probably should be calibrated
        warning('integer data is not yet calilbrated');
      case 3
        samplesize = 4;
        sampletype = 'float';
        % assume that the data is calibrated
      case 4
        samplesize = 8;
        sampletype = 'double';
        % assume that the data is calibrated
      otherwise
        error('unsupported data format');
    end
    % jump to the desired data
    fseek(fid, offset*samplesize, 'cof');
    % read the desired data
    if length(chanindx)==1
      % read only one channel
      fseek(fid, (chanindx-1)*samplesize, 'cof');                            % seek to begin of channel
      dat = fread(fid, numsamples, ['1*' sampletype], (hdr.nChans-1)*samplesize)'; % read one channel, skip the rest
    else
      % read all channels
      dat = fread(fid, [hdr.nChans, numsamples], sampletype);
    end
    fclose(fid);
    if length(chanindx)==1
      % only one channel was selected, which is managed by the code above
      % nothing to do
    elseif length(chanindx)==hdr.nChans
      % all channels have been selected
      % nothing to do
    else
      % select the desired channel(s)
      dat = dat(chanindx,:);
    end

  case 'besa_avr'
    % BESA average data
    orig = read_besa_avr(filename);
    dat  = orig.data(chanindx, begsample:endsample);

  case 'besa_swf'
    % BESA source waveform
    orig = read_besa_swf(filename);
    dat  = orig.data(chanindx, begsample:endsample);

  case {'biosemi_bdf', 'bham_bdf'}
    % this uses the openbdf and readbdf functions that I copied from the EEGLAB toolbox
    epochlength = hdr.orig.Head.SampleRate(1);
    % it has already been checked in read_header that all channels have the same sampling rate
    begepoch = floor((begsample-1)/epochlength) + 1;
    endepoch = floor((endsample-1)/epochlength) + 1;
    nepochs  = endepoch - begepoch + 1;
    orig = openbdf(filename);
    dat  = zeros(length(chanindx),nepochs*epochlength);
    for i=begepoch:endepoch
      % read and concatenate all required data epochs
      [orig, buf] = readbdf(orig, i, 0);
      if size(buf,2)~=hdr.nChans || size(buf,1)~=epochlength
        error('error reading selected data from bdf-file');
      else
        dat(:,((i-begepoch)*epochlength+1):((i-begepoch+1)*epochlength)) = buf(:,chanindx)';
      end
    end
    begsample = begsample - (begepoch-1)*epochlength;  % correct for the number of bytes that were skipped
    endsample = endsample - (begepoch-1)*epochlength;  % correct for the number of bytes that were skipped
    dat = dat(:, begsample:endsample);
    % close the file between seperate read operations
    fclose(orig.Head.FILE.FID);

  case {'brainvision_eeg', 'brainvision_dat'}
    dat = read_brainvision_eeg(filename, hdr, begsample, endsample);
    dat = dat(chanindx,:);	% select the desired channels

  case 'brainvision_seg'
    dat = read_brainvision_seg(filename, hdr, begsample, endsample);
    dat = dat(chanindx,:);	% select the desired channels

  case 'ced_son'
    % chek the availability of the required low-level toolbox
    hastoolbox('neuroshare', 1);
    % use the reading function supplied by Gijs van Elswijk
    %
    % CED ADC is done in sequence, thus the analog channels
    % do not share the same time axis. This is _ignored_ here.
    %
    % Set the READ_CED_SON parameter 'readtimestamps' to
    % 'yes' to get time axis for each data channel returned.
    % This time information can be used to do
    % a temporal realignment of the data.
    tmp = read_ced_son(filename,'readdata','yes',...
      'begsample',begsample,...
      'endsample',endsample,...
      'channels',chanindx);
    dat = cell2mat(tmp.data');

  case {'ctf_ds', 'ctf_meg4', 'ctf_res4', 'read_ctf_meg4'} % the default reader for CTF is read_ctf_meg4
    if strcmp(dataformat, 'ctf_ds')
      [p, f, x] = fileparts(filename);
      filename = fullfile(p, [f '.meg4']);
    elseif strcmp(dataformat, 'ctf_res4')
      [p, f, x] = fileparts(filename);
      filename = fullfile(p, [f '.meg4']);
    end
    % read it using the default CTF importer that originates from CTF and the FCDC
    dat = read_ctf_meg4(filename, hdr.orig, begsample, endsample, chanindx);

  case 'ctf_read_meg4'
    % check that the required low-level toolbos ix available
    hastoolbox('eegsf', 1);
    if strcmp(dataformat, 'ctf_res4')
      [p, f, x] = fileparts(filename);
      filename = p;
    elseif strcmp(dataformat, 'ctf_meg4')
      [p, f, x] = fileparts(filename);
      filename = p;
    end
    % read it using the CTF importer from the NIH and Daren Weber
    tmp = ctf_read_meg4(filename, hdr.orig, chanindx, 'all', begtrial:endtrial);
    dat = cat(3, tmp.data{:});
    % the data is shaped in a 3-D array
    dimord = 'samples_chans_trials';

  case 'eep_avr'
    % check that the required low-level toolbos ix available
    hastoolbox('eeprobe', 1);
    dat = read_eep_avr(filename);
    dat = dat.data(chanindx,begsample:endsample);       % select the desired channels and samples

  case 'eep_cnt'
    % check that the required low-level toolbos ix available
    hastoolbox('eeprobe', 1);
    dat = read_eep_cnt(filename, begsample, endsample);
    dat = dat.data(chanindx,:);                         % select the desired channels

  case 'fcdc_matbin'
    % multiplexed data in a *.bin file, accompanied by a matlab file containing the header
    offset        = begsample-1;
    numsamples    = endsample-begsample+1;
    samplesize    = 8;
    sampletype    = 'double';
    [fid,message] = fopen(datafile,'rb','ieee-le');
    % jump to the desired data
    fseek(fid, offset*samplesize*hdr.nChans, 'cof');
    % read the desired data
    if length(chanindx)==1
      % read only one channel
      fseek(fid, (chanindx-1)*samplesize, 'cof');                            % seek to begin of channel
      dat = fread(fid, numsamples, ['1*' sampletype], (hdr.nChans-1)*samplesize)'; % read one channel, skip the rest
    else
      % read all channels
      dat = fread(fid, [hdr.nChans, numsamples], sampletype);
    end
    fclose(fid);
    if length(chanindx)==1
      % only one channel was selected, which is managed by the code above
      % nothing to do
    elseif length(chanindx)==hdr.nChans
      % all channels have been selected
      % nothing to do
    else
      % select the desired channel(s)
      dat = dat(chanindx,:);
    end

  case {'mpi_ds', 'mpi_dap'}
    [hdr, dat] = read_mpi_ds(filename);
    dat = dat(chanindx, begsample:endsample); % select the desired channels and samples

  case 'neuralynx_dma'
    if length(chanindx)==hdr.nChans && all(chanindx==1:length(chanindx))
      % read the selected data for all channels, no reordering needed
      dat = read_neuralynx_dma(filename, begsample, endsample);
      dat = dat.Data;
    elseif length(chanindx)==hdr.nChans && ~all(chanindx==1:length(chanindx))
      % read the selected data for all channels, channels have to be reordered
      dat = read_neuralynx_dma(filename, begsample, endsample);
      dat = dat.Data(chanindx,:);
    else
      % read one-second chuncks of data with all channels, subsequently select the channels of interest
      chunck = begsample:hdr.Fs:endsample;
      if chunck(end)~=endsample
        % the selected data does not consist of an integer number of seconds
        chunck(end+1) = endsample;
      end
      % preallocate memory to hold the result
      dat = zeros(length(chanindx), endsample-begsample+1);
      for i=1:(length(chunck)-1)
        % read each chunck and select the channels of interest
        chunck_begsample = chunck(i);
        chunck_endsample = chunck(i+1);
        buf = read_neuralynx_dma(filename, chunck_begsample, chunck_endsample);
        dat(:,(chunck_begsample-begsample+1):(chunck_endsample-begsample+1)) = buf.Data(chanindx,:);
      end
    end

  case 'neuralynx_sdma'
    dat = read_neuralynx_sdma(filename, begsample, endsample, chanindx);

  case 'neuralynx_ncs'
    NRecords  = hdr.nSamples/512;
    begrecord = floor(begsample/512)+1;
    endrecord = ceil(endsample/512);
    % read the records that contain the desired samples
    ncs = read_neuralynx_ncs(filename, begrecord, endrecord);
    dat = zeros(1,(endrecord-begrecord+1)*512);
    for i=1:(endrecord-begrecord+1)
      dat(((i-1)*512+1):(i*512)) = ncs.dat{i}(:);
    end
    % cut out the desired samples
    begsample = begsample - (begrecord-1)*512;
    endsample = endsample - (begrecord-1)*512;
    dat = dat(begsample:endsample);

  case 'neuralynx_nse'
    % read all records
    nse = read_neuralynx_nse(filename);
    % convert timestamps to samples
    sample = round((nse.TimeStamp - hdr.FirstTimeStamp)./hdr.TimeStampPerSample + 1);
    % select the timestamps that are between begin and endsample
    sample = sample(sample>=begsample & sample<=endsample) - begsample + 1;
    dat = zeros(1,endsample-begsample+1);
    dat(sample) = 1;

  case {'neuralynx_ttl', 'neuralynx_tsl', 'neuralynx_tsh'}
    dat = read_neuralynx_ttl(filename, begsample, endsample);

  case 'neuralynx_ds'
    dat = read_neuralynx_ds(filename, hdr, begsample, endsample, chanindx);

  case 'neuralynx_cds'
    % this is a combined Neuralynx dataset with seperate subdirectories for the LFP, MUA and spike channels
    dirlist = dir(filename);
    haslfp   = any(filetype_check_extension({dirlist.name}, 'lfp'));
    hasmua   = any(filetype_check_extension({dirlist.name}, 'mua'));
    hasspike = any(filetype_check_extension({dirlist.name}, 'spike'));

    lfpdir   = '';
    muadir   = '';
    spikedir = '';
    if haslfp
      lfpdir = fullfile(filename, dirlist(find(filetype_check_extension({dirlist.name}, 'lfp'))).name);
    end
    if hasmua
      muadir = fullfile(filename, dirlist(find(filetype_check_extension({dirlist.name}, 'mua'))).name);
    end
    if hasspike
      spikedir = fullfile(filename, dirlist(find(filetype_check_extension({dirlist.name}, 'spike'))).name);
    end

    chandir    = {lfpdir, muadir, spikedir};
    chanprefix = {'lfp_', 'mua_', 'spike_'};
    chanext    = {'.ncs', '.ncs', '.nse'};
    chantype   = zeros(hdr.nChans,1);
    chantype(strncmp(hdr.label, 'lfp_'  , 4)) = 1;  % these channels in the combined dataset contain LFP
    chantype(strncmp(hdr.label, 'mua_'  , 4)) = 2;  % these channels in the combined dataset contain MUA
    chantype(strncmp(hdr.label, 'spike_', 6)) = 3;  % these channels in the combined dataset contain spike
    chanfile = cell(hdr.nChans, 1);
    orglabel = cell(hdr.nChans, 1);
    for i=1:hdr.nChans
      % determine the original name of the channel, i.e. without the prefix
      orglabel{i} = hdr.label{i}((1+length(chanprefix{chantype(i)})):end);
      % determine the original filename of the channel including the directory and file extension
      chanfile{i} = fullfile(chandir{chantype(i)}, [orglabel{i} chanext{chantype(i)}]);
    end

    dat = zeros(length(chanindx), endsample-begsample+1);
    for i=1:length(chanindx)
      chan     = chanindx(i);                                     % channel number in the combined dataset
      orgfile  = chanfile{chan};                                  % original file name, including path
      orghdr   = hdr.orig{chantype(chan)};                        % header of the original dataset containing this channel
      orgindx  = strmatch(orglabel{chan}, orghdr.label, 'exact'); % index of the channel in the original dataset
      dat(i,:) = read_data(orgfile, orghdr, begsample, endsample, orgindx);
    end

  case 'ns_avg'
    % NeuroScan average data
    orig = read_ns_avg(filename);
    dat  = orig.data(chanindx, begsample:endsample);

  case 'ns_cnt'
    % Neuroscan continuous data
    sample1    = begsample-1;
    ldnsamples = endsample-begsample+1; % number of samples to read
    ldchan     = 1:hdr.nChans;          % must be row vector
    chanoi     = chanindx(:)';          % channels of interest
    if sample1<0
      error('begin sample cannot be for the beginning of the file');
    end
    % read_ns_cnt originates from the EEGLAB package (loadcnt.m) but is
    % an old version since the new version is not compatible any more
    % all data is read, and only the relevant data is kept.
    tmp = read_ns_cnt(filename, 'sample1', sample1, 'ldnsamples', ldnsamples, 'ldchan', ldchan);
    dat = tmp.dat(chanoi,:);

  case 'ns_eeg'
    % Neuroscan epoched file
    tmp       = read_ns_eeg(filename, begtrial:endtrial);
    siz       = [(endtrial-begtrial+1) hdr.nChans hdr.nSamples];
    dat       = reshape(tmp.data, siz); % ensure 3-D array
    dat       = dat(:,chanindx,:);      % select channels
    dimord    = 'trials_chans_samples'; % selection using begsample and endsample will be done later

  case 'neuromag_fif'
    % check that the required low-level toolbox is available
    hastoolbox('meg-pd', 1);
    begtime = (begsample-1)/hdr.Fs;
    begepoch = floor((begsample-1)/hdr.nSamples) + 1;
    endepoch = floor((endsample-1)/hdr.nSamples) + 1;
    rawdata('any',filename);
    rawdata('goto', begtime);
    dat = [];
    for i=begepoch:endepoch
      [buf, status] = rawdata('next');
      if ~strcmp(status, 'ok')
        error('error reading selected data from fif-file');
      else
        dat(:,((i-begepoch)*hdr.nSamples+1):((i-begepoch+1)*hdr.nSamples)) = buf(chanindx,:);
      end
    end
    rawdata('close');
    begsample = begsample - (begepoch-1)*hdr.nSamples;  % correct for the number of bytes that were skipped
    endsample = endsample - (begepoch-1)*hdr.nSamples;  % correct for the number of bytes that were skipped
    dat = dat(:, begsample:endsample);

  case 'plexon_ddt'
    dat = read_plexon_ddt(filename, begsample, endsample);
    dat = dat.data(chanindx,:);

  case 'plexon_nex'
    dat = read_nex_data(filename, hdr, begsample, endsample, chanindx);

  case {'yokogawa_ave', 'yokogawa_con', 'yokogawa_raw'}
    % check that the required low-level toolbox is available
    hastoolbox('yokogawa', 1);
    dat = read_yokogawa_data(filename, hdr, begsample, endsample, chanindx);

  otherwise
    error('unsupported data format');
end

if ~exist('dimord', 'var')
  dimord = 'chans_samples';  % almost all low-level readers return the data as 2D array
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% reshape the 2-D or 3-D matrix to a common order of the dimensions
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
switch dimord
  case {'chans_samples', 'chans_samples_trials'}
    % nothing to do
  case 'samples_chans'
    dat = permute(dat, [2 1]);
    dimord = 'chans_samples';
  case 'samples_chans_trials'
    dat = permute(dat, [2 1 3]);
    dimord = 'chans_samples_trials';
  case 'trials_samples_chans'
    dat = permute(dat, [3 2 1]);
    dimord = 'chans_samples_trials';
  case 'trials_chans_samples'
    dat = permute(dat, [2 3 1]);
    dimord = 'chans_samples_trials';
  otherwise
    error('unexpected dimord');
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% convert between 3-D trial based and 2-D continuous output
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if requesttrials  && strcmp(dimord, 'chans_samples')
  % reformat the continuous representation into trials
  nchans   = size(dat,1);
  nsamples = hdr.nSamples;
  ntrials  = size(dat,2)/hdr.nSamples;
  dat = reshape(dat, [nchans nsamples ntrials]); % convert into a 3-D array

elseif requestsamples && strcmp(dimord, 'chans_samples_trials')
  % reformat the trials into a continuous representation
  nchans   = size(dat,1);
  nsamples = size(dat,2);
  ntrials  = size(dat,3);
  dat = reshape(dat, [nchans nsamples*ntrials]); % convert into a 2-D array
  % determine the selection w.r.t. the data as it is on disk
  begselection = (begtrial-1)*hdr.nSamples + 1;
  endselection = (endtrial  )*hdr.nSamples;
  % determine the selection w.r.t. the data that has been read in
  begselection = begsample - begselection + 1;
  endselection = endsample - begselection + 1;
  dat = dat(:,begselection:endselection);
end