function [dat] = read_fcdc_data(datafile, hdr, begsample, endsample, chanindx, continuous);

% READ_FCDC_DATA is a wrapper around different EEG/MEG file importers
% directly supported formats are CTF, Neuromag, EEP, BrainVision, Neuroscan
% and Neuralynx.
%
% If you have the BIOSIG toolbox installed also EDF, BKR, CNT, BDF, GDF and
% other dataformats supported by that toolbox are supported.
% See http://biosig.sourceforge.net for details.
%
% Use as
%   dat = read_fcdc_data(datfile, hdr, begsample, endsample, chanindx, continuous)
% where hdr is obtained using READ_FCDC_HEADER.
%
% The beginsample and endsample are counted from the beginning of the raw
% data, starting at 1. This function treats all data as continuous, even if
% it is stored in disconnected data segments (i.e. trials or epochs). If
% you want to read the second trial, you should read from sample
% hdr.nSamples+1 up to 2*hdr.nSamples.
%
% The argument "chanindx" is optional, the default is to read all channels.
%
% The argument "continuous" is optional, the default is check to whether
% the data is continuous and, if not, check whether the requested data
% segment does not extend over a trial boundary. If continuous=1, this
% check is not performed (usefull for pseudo-continuous data).
%
% See also READ_FCDC_HEADER

% Copyright (C) 2003-2005, Robert Oostenveld, F.C. Donders Centre
%
% $Log: read_fcdc_data.m,v $
% Revision 1.35  2006/02/24 16:39:09  roboos
% reorder the channels if length(chansel)==nChans but chansel~=1:nChans (only for neuralynx_dma, should be checked for other formats)
%
% Revision 1.34  2006/02/07 22:29:45  roboos
% fixed a bug in neuromag_fif, matrix should not be transposed (thanks to Jon Houck)
% reimplemented the chunck-reading for biosig, since it was completely broken
%
% Revision 1.33  2006/01/11 16:27:26  roboos
% fixed bug due to incorrect variable name for bdf (thanks to Hoi Fei Kwok)
%
% Revision 1.32  2005/12/16 13:54:22  roboos
% removed the haseeglab, openbdf and readbdf are now part of fieldtrip/private
%
% Revision 1.31  2005/12/13 10:56:43  roboos
% changed biosig reading, do not use subfunction but implement here immediately
% add "chuncked" reading for biosig
% if endsample=inf, read untill end of file
% assume that biosemi_bdf is always continuous
%
% Revision 1.30  2005/12/02 08:31:11  roboos
% added hadeeglab and hasneuroshare flags
% implemented support for ced_son fileformat, based on read_ced_son from Gijs van Elswijk
% implemented support for BDF using EEGLAB importers
% made buffer copying more efficient for neuromag (prevent matrix from growing)
%
% Revision 1.29  2005/11/23 10:42:41  roboos
% added check for begsample>=1 and endsample<=eof
%
% Revision 1.28  2005/11/23 08:29:48  roboos
% reordered the initial default settings and validity checks, added check on channel selection
% implemented support for neuralynx_ncs
%
% Revision 1.27  2005/09/29 00:46:41  roboos
% changed from yokogawa_sqd to yokogawa_ave (as suggested by Masahiro Shimogawara)
%
% Revision 1.26  2005/10/05 06:29:24  roboos
% improved documentation
%
% Revision 1.25  2005/10/04 16:50:56  roboos
% added support for MPI Frankfurt *.dap files
%
% Revision 1.24  2005/09/15 11:33:38  roboos
% Changed the neuralynx_dma reader to read in 1 second chuncks of data
% (all channels) and then select the channels of interest, instead of
% reading in all data for all channels prior to selecting the channels of
% interest. This results in much more efficient memory usage, and could
% also be applied to teh other low-level reading routines that do not
% support the reading of a selection of channels.
%
% Revision 1.23  2005/09/09 09:13:07  roboos
% added support for neuralynx_dma
%
% Revision 1.22  2005/09/08 08:43:22  roboos
% added support for Yokogawa MEG data
%
% Revision 1.21  2005/09/06 12:45:00  roboos
% renamed plextor into plexon (incorrect company name)
%
% Revision 1.20  2005/09/02 16:35:15  roboos
% added support for plexon_ddt (only version 100 sofar)
%
% Revision 1.19  2005/08/17 19:35:29  roboos
% added a check for the presence of the file, give error if not present
%
% Revision 1.18  2005/07/28 09:41:36  roboos
% fixed typo in documentation
%
% Revision 1.17  2005/05/17 17:50:37  roboos
% changed all "if" occurences of & and | into && and ||
% this makes the code more compatible with Octave and also seems to be in closer correspondence with Matlab documentation on shortcircuited evaluation of sequential boolean constructs
%
% Revision 1.16  2005/05/12 07:21:49  roboos
% added support for neuralynx_ds
%
% Revision 1.15  2005/03/31 07:05:58  roboos
% implemented support for besa_avr using my own read_besa_avr
%
% Revision 1.14  2005/03/15 17:44:29  roboos
% Implemented the option of reading CTF res4/meg4 files using the importers
% from Darren Webers EEG toolbox (also used in EEGLAB): first try to read it
% using CTF/FCDC code, then try to read it with the other code. Gradiometer
% information for forward computations is not yet converted/returned.
%
% Revision 1.13  2005/02/16 07:48:28  roboos
% replaced nex specific code by a call to read_nex_header and read_nex_data subfunctions in file_io
%
% Revision 1.12  2005/02/11 07:39:22  roboos
% added support for Plexon *.nex files (only continuous data)
%
% Revision 1.11  2005/02/07 12:46:31  jansch
% small modification for the handling of continuous neuroscan-data, so that
% a subset of channels can be specified in the configuration.
%
% Revision 1.10  2005/02/03 17:14:13  roboos
% fixed bug for neuromag data (beginepoch not defined)
%
% Revision 1.9  2005/01/25 14:03:58  roboos
% implemented check for datasegment extending beyound trial boundaries
% default is to determine continuous/epoched datatype by looking at the header
% added extra boolean argument for (pseudo)continuous data in which case boundary check is skipped
%
% Revision 1.8  2004/12/06 13:03:32  roboos
% added support for input string with ctf dataset instead of res4/meg4 file
%
% Revision 1.7  2004/10/01 10:26:24  roboos
% made detection of external toolboxes consistent with read_fcdc_mri, i.e. hasXXX at the beginning
%
% Revision 1.6  2004/03/30 07:17:57  roberto
% fixed bug in chanindx due to incorrect number of input arguments
%
% Revision 1.5  2004/03/29 15:07:25  roberto
% added support for EEProbe *.avg and Brainvision *.dat (exported from BVA)
%
% Revision 1.4  2004/03/11 10:20:45  roberto
% implemeted default chanindx argument for all channels + some code changes to reflect that
% added extra comment on loadcnt importer originating from EEGLAB
%
% Revision 1.3  2004/03/10 14:05:27  roberto
% fixed a potential bug in reading multiple epochs of a Neuroscan *.eeg file
% added support for Neuromag *.fif files
%
% Revision 1.2  2004/03/10 11:58:13  roberto
% added support for BIOSIG importers as the last (fallback) attempt
%
% Revision 1.1  2004/01/16 08:08:30  roberto
% moved from file_io toolbox to framework
%
% Revision 1.4  2003/12/04 09:01:00  roberto
% made ns_eeg consistent with other formats by changing hdr.nSample into samples per trial
%
% Revision 1.3  2003/11/28 16:36:43  roberto
% added support for Eric Maris' continuous neuroscan *.cnt files
%
% Revision 1.2  2003/10/23 11:11:25  roberto
% converted ascii text from dos to unix
%
% Revision 1.1  2003/10/23 11:10:10  roberto
% *** empty log message ***
%

% detect the external toolboxes that optionally can be used
hasbiosig     = (exist('sopen') & exist('sread'));                       % see http://biosig.sourceforge.net/
hasmegpd      = (exist('rawdata') & exist('channames'));                 % see http://www.kolumbus.fi/kuutela/programs/meg-pd/
haseegsf      = (exist('ctf_read_meg4') & exist('ctf_read_res4'));       % see http://eeg.sourceforge.net/
hasneuroshare = (exist('ns_SetLibrary') & exist('ns_GetAnalogData'));    % see http://www.ced.co.uk and http://www.neuroshare.org

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

% set the default continuous flag, which is to assume that it is not continuous
if nargin<6 || isempty(continuous)
  continuous = 0;
end

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

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

% test whether the file exists
if ~exist(datafile)
  error(sprintf('file ''%s'' does not exist', datafile));
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 begsample<1 || endsample>(hdr.nSamples*hdr.nTrials)
  error('cannot read data before the start or after the end of the file');
end

% test whether the requested data segment does not extend over a discontinuous trial boundary
if ~continuous && hdr.nTrials>1
  begtrl = floor((begsample-1)/hdr.nSamples) + 1;
  endtrl = floor((endsample-1)/hdr.nSamples) + 1;
  if begtrl~=endtrl
    error('requested data segment extends over a discontinuous trial boundary');
  end
end

% detect the filetype and reat the requested channels and samples
if filetype(datafile, 'ctf_meg4')
  try
    % try to read it using the default CTF importer that originates from CTF and the FCDC
    dat = read_ctf_meg4(datafile, hdr, begsample, endsample, chanindx);
    ctferror = 0;
  catch
    ctferror = 1;
  end
  if ctferror && haseegsf
    % try to read it using the CTF importer from the NIH and Daren Weber
    try
      % determine the trial containing the begin and end sample
      begtrial = floor((begsample-1)/hdr.nSamples)+1;
      endtrial = floor((endsample-1)/hdr.nSamples)+1;
      tmp = ctf_read_meg4(fileparts(datafile), hdr.orig, chanindx, 'all', begtrial:endtrial);
      % convert the data into the format that FieldTrip expects, i.e. continuous Nchans X Nsamples
      dat = permute(tmp.data, [2 1 3]);
      dat = reshape(dat, length(chanindx), length(begtrial:endtrial)*hdr.nSamples);
      begsample = begsample - (begtrial-1)*hdr.nSamples;
      endsample = endsample - (begtrial-1)*hdr.nSamples;
      dat = dat(:, begsample:endsample);
      ctferror = 0;
    catch
      ctferror = 1;
    end
  end
  if ctferror
    disp(lasterr);
    error('could not read CTF meg4 data file');
  end

elseif filetype(datafile, 'plexon_nex')
  dat = read_nex_data(datafile, hdr, begsample, endsample, chanindx);

elseif filetype(datafile, 'plexon_ddt')
  dat = read_plexon_ddt(datafile, begsample, endsample);
  dat = dat.data(chanindx,:);

elseif filetype(datafile, '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(datafile, 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(datafile, 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(datafile, chunck_begsample, chunck_endsample);
      dat(:,(chunck_begsample-begsample+1):(chunck_endsample-begsample+1)) = buf.Data(chanindx,:);
    end
  end

elseif filetype(datafile, '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(datafile, 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);

elseif filetype(datafile, 'neuralynx_ds')
  dat = read_neuralynx_data(datafile, hdr, begsample, endsample, chanindx);

elseif filetype(datafile,'ced_son')
  % Read data using read_ced_son routines. These routines need the
  % NeuroShare library (http://www.neuroshare.org). As of now
  % there is only a Windows version available of the NeuroShare lib.
  if ~hasneuroshare
    error('cannot find CED SON import routines, see http://www.neuroshare.org');
  end;
  if ~continuous,
    error('only continuous data supported for SON files')
  else,
    % 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(datafile,'readdata','yes',...
      'begsample',begsample,...
      'endsample',endsample,...
      'channels',chanindx);
    dat = cell2mat(tmp.data');
  end;

elseif filetype(datafile, 'brainvision_eeg')
  dat = read_brainvision_eeg(datafile, hdr, begsample, endsample);
  dat = dat(chanindx,:);	% select the desired channels

elseif filetype(datafile, 'brainvision_dat')
  dat = read_brainvision_eeg(datafile, hdr, begsample, endsample);
  dat = dat(chanindx,:);	% select the desired channels

elseif filetype(datafile, 'brainvision_seg')
  dat = read_brainvision_seg(datafile, hdr, begsample, endsample);
  dat = dat(chanindx,:);	% select the desired channels

elseif filetype(datafile, 'yokogawa_ave') || filetype(datafile, 'yokogawa_con') || filetype(datafile, 'yokogawa_raw')
  dat = read_yokogawa_data(datafile, hdr, begsample, endsample, chanindx);

elseif filetype(datafile, 'mpi_ds') || filetype(datafile, 'mpi_dap')
  [hdr, dat] = read_mpi_ds(datafile);
  dat = dat(chanindx, begsample:endsample); % select the desired channels and samples

elseif filetype(datafile, 'eep_cnt')
  dat = read_eep_cnt(datafile, begsample, endsample);
  dat = dat.data(chanindx,:);	% select the desired channels

elseif filetype(datafile, 'eep_avr')
  dat = read_eep_avr(datafile);
  dat = dat.data(chanindx,begsample:endsample);	% select the desired channels and samples

elseif filetype(datafile, 'ns_eeg')
  % Neuroscan epoched file
  begepoch  = floor((begsample-1)/hdr.npnt) + 1;  % determine epoch with the requested samples
  tmp       = read_ns_eeg(datafile, begepoch);
  begsample = begsample - (begepoch-1)*hdr.npnt;  % correct for the number of bytes that were skipped
  endsample = endsample - (begepoch-1)*hdr.npnt;  % correct for the number of bytes that were skipped
  dat = tmp.data(chanindx, begsample:endsample);

elseif filetype(datafile, '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(datafile, 'sample1', sample1, 'ldnsamples', ldnsamples, 'ldchan', ldchan);
  dat = tmp.dat(chanoi,:);

elseif filetype(datafile, 'besa_avr')
  % BESA average data
  orig = read_besa_avr(datafile);
  dat  = orig.data(chanindx, begsample:endsample);

elseif filetype(datafile, 'neuromag_fif')
  % Neuromag MEG data
  if ~hasmegpd
    error('cannot find Neuromag import routines, see http://boojum.hut.fi/~kuutela/meg-pd/');
  end

  begtime = (begsample-1)/hdr.Fs;
  begepoch = floor((begsample-1)/hdr.nSamples) + 1;
  endepoch = floor((endsample-1)/hdr.nSamples) + 1;
  rawdata('any',datafile);
  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);

elseif filetype(datafile, 'biosemi_bdf') || filetype(datafile, '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_fcdc_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(datafile);
  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);

elseif hasbiosig && ~isempty(hdr)
  % Try reading the data with the BIOSIG import filters if they are available.
  % This only needs to be attempted if the READ_FCDC_HEADER file was able to
  % extract header information out of the file, otherwise the file cannot
  % be recognized anyway.
  try
    if length(chanindx)==hdr.nChans
      % read the selected data for all channels at once
      biosig = sopen(datafile,'r');
      begtime  = (begsample-1) / hdr.Fs;
      endtime  = (endsample  ) / hdr.Fs;
      duration = endtime-begtime;
      dat = sread(biosig, duration, begtime)';
      if ~all(chanindx(:)'==1:hdr.nChans)
        % the desired ordering of channels is different from that in the file
        dat = dat(chanindx, :);
      end
    else
      % it could be that the file is continuous or epoched, but to save (intermediate) memory it is read 
      % in segemts of one second with all channels, subsequently the channels of interest are selected
      % preallocate memory to hold the result
      dat = zeros(length(chanindx), endsample-begsample+1);
      biosig = sopen(datafile,'r');
      numsegment = ceil((endsample-begsample+1)/hdr.Fs);
      for i=1:numsegment
        % read each segment and select the channels of interest
        segment_begsample = (i-1)*hdr.Fs + begsample;
        segment_endsample = (i  )*hdr.Fs + begsample - 1;
        if segment_endsample>endsample
          % the last segment can be less than one second
          segment_endsample = endsample;
        end
        begtime  = (segment_begsample-1) / hdr.Fs;
        endtime  = (segment_endsample  ) / hdr.Fs;
        duration = endtime-begtime;
        buf = sread(biosig, duration, begtime)';
        dat(:,(segment_begsample:segment_endsample)-begsample+1) = buf(chanindx,:);
      end
    end
    sclose(biosig);
  catch
    error('error reading data using BIOSIG toolbox, see http://biosig.sourceforge.net/');
  end
end

if ~exist('dat', 'var')
  % none of the importers above has been able to read the data
  error('unsupported fileformat for EEG/MEG data file');
end
