function [hdr] = read_fcdc_header(headerfile);

% READ_FCDC_HEADER 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 by this function.
% See http://biosig.sourceforge.net for details.
%
% Use as
%   hdr = read_fcdc_header(headerfile);
%
% This returns a header structure with the following elements
%   hdr.Fs           sampling frequency
%   hdr.nChans       number of channels
%   hdr.nSamples     number of samples per trial
%   hdr.nSamplesPre  number of pre-trigger samples in each trial
%   hdr.nTrials      number of trials
%   hdr.label        cell-array with labels of each channel
%
% For continuous data, nSamplesPre=0 and nTrials=1.
%
% Depending on the dataformat, additional header information can be 
% returned in the hdr.orig subfield. 
%
% See also READ_FCDC_DATA, READ_FCDC_EVENT

% Copyright (C) 2003-2005, F.C. Donders Centre
%
% $Log: read_fcdc_header.m,v $
% Revision 1.36  2006/02/24 16:40:24  roboos
% switched back to First and LastTimeStamp that fit all channels
%
% Revision 1.35  2006/02/09 10:05:03  roboos
% cleaned up the header conversion for neuralynx
%
% Revision 1.34  2006/01/30 14:16:09  roboos
% explicitely add add a gradiometer structure for forward and inverse
% modelling in case of yokogawa data
%
% Revision 1.33  2005/12/16 13:55:46  roboos
% removed the haseeglab, openbdf and readbdf are now part of fieldtrip/private
% implemented support for bham_bdf, which does nemaing of the standard biosemi channel labels to the 5% labels
%
% Revision 1.32  2005/12/02 08:32:06  roboos
% added haseeglab 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
%
% Revision 1.31  2005/11/23 08:27:40  roboos
% added support for neuralynx_ncs
%
% Revision 1.30  2005/09/29 00:46:41  roboos
% changed from yokogawa_sqd to yokogawa_ave (as suggested by Masahiro Shimogawara)
%
% Revision 1.29  2005/10/13 08:36:15  roboos
% added FirstTimeStamp and LastTimeStamp to header for neuralynx_ds
%
% Revision 1.28  2005/10/05 06:29:44  roboos
% improved documentation
%
% Revision 1.27  2005/10/04 16:50:56  roboos
% added support for MPI Frankfurt *.dap files
%
% Revision 1.26  2005/09/09 09:22:23  roboos
% replaced round() by floor() in the computation of the number of samples for neuralynx_ds
% there seems to be some drift in the timestamps
%
% Revision 1.25  2005/09/09 09:13:07  roboos
% added support for neuralynx_dma
%
% Revision 1.24  2005/09/09 08:42:01  roboos
% added hdr.nSamples for neuralynx_ds header and moved original header information to hdr.orig
%
% Revision 1.23  2005/09/08 16:56:56  roboos
% added correct number of samples for plexon_ddt
% added original header for plexon_ddt
%
% 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:30  roboos
% added a check for the presence of the file, give error if not present
%
% Revision 1.18  2005/06/17 10:57:54  roboos
% start with an empty header and test using isempty, to prevent the use of exist function (which has compiler problems)
%
% Revision 1.17  2005/05/26 09:57:02  roboos
% always explicitely add the grad structure using one of the xxx2grad functions
% added support for getting gradiometer info from NIMH ctf inport functions
%
% Revision 1.16  2005/05/17 17:50:38  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.15  2005/05/12 07:21:49  roboos
% added support for neuralynx_ds
%
% Revision 1.14  2005/03/31 07:05:58  roboos
% implemented support for besa_avr using my own read_besa_avr
%
% Revision 1.13  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.12  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.11  2005/02/11 07:39:22  roboos
% added support for Plexon *.nex files (only continuous data)
%
% Revision 1.10  2004/12/06 13:03:32  roboos
% added support for input string with ctf dataset instead of res4/meg4 file
%
% Revision 1.9  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.8  2004/09/20 13:53:04  roboos
% moved gradiometer definition for Neuromag fif files into separate function
%
% Revision 1.7  2004/06/25 15:08:08  roberto
% fixed bug in the labels of Neuromag gradiometer definition
%
% Revision 1.6  2004/06/24 15:52:35  roberto
% added gradiometer information from Neuromag fif file (122 channel system)
%
% Revision 1.5  2004/03/29 15:11:44  roberto
% added support for EEProbe *.avg, removed original header structure from output hdr for biosig
%
% Revision 1.4  2004/03/11 10:22:12  roberto
% some changes to organization of neuroscan cnt header: its original elements are now tucked away in a separate (orig) field
%
% Revision 1.3  2004/03/10 14:04:26  roberto
% added support for Neuromag *.fif format
% changed newline from dos to unix
%
% 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.5  2003/12/05 11:05:46  roberto
% improved help, changed tmp into hdr so that it always returns as much
% header information as possible
%
% Revision 1.4  2003/12/02 12:51:55  roberto
% separated neuroscan eeg and cnt formats, both use different routine
%
% 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

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

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

% start with an empty header
hdr = [];

if filetype(headerfile, 'ctf_res4')
  try
    % try to read it using the default CTF importer that originates from CTF and the FCDC
    hdr = read_ctf_res4(headerfile);
    ctferror = 0;
    % add a gradiometer structure for forward and inverse modelling
    hdr.grad = ctf2grad(hdr);
  catch
    ctferror = 1;
  end
  if ctferror && haseegsf
    % try to read it using the CTF importer from the NIH and Daren Weber
    try
      tmp = ctf_read_res4(fileparts(headerfile), 0);
      ctferror = 0;
      % convert the header into a structure that FieldTrip understands
      hdr              = [];
      hdr.Fs           = tmp.setup.sample_rate;
      hdr.nChans       = length(tmp.sensor.info);
      hdr.nSamples     = tmp.setup.number_samples;
      hdr.nSamplesPre  = tmp.setup.pretrigger_samples;
      hdr.nTrials      = tmp.setup.number_trials;
      for i=1:length(tmp.sensor.info)
        hdr.label{i}   = tmp.sensor.info(i).label;
      end
      hdr.label        = hdr.label(:);
      % add a gradiometer structure for forward and inverse modelling
      hdr.grad         = nimh2grad(tmp);
      % also store the original header, so that it can be reused by read_fcdc_data
      hdr.orig         = tmp;
    catch
      ctferror = 1;
    end
  end
  if ctferror
    error('could not read CTF res4 header file');
  end

elseif filetype(headerfile, 'brainvision_vhdr')
  hdr = read_brainvision_vhdr(headerfile);

elseif filetype(headerfile, 'yokogawa_ave') || filetype(headerfile, 'yokogawa_con') || filetype(headerfile, 'yokogawa_raw')
  hdr = read_yokogawa_header(headerfile);
  % add a gradiometer structure for forward and inverse modelling
  hdr.grad = yokogawa2grad(hdr);

elseif filetype(headerfile, 'mpi_ds') || filetype(headerfile, 'mpi_dap')
  hdr = read_mpi_ds(headerfile);

elseif filetype(headerfile, 'plexon_nex')
  hdr = read_nex_header(headerfile);

elseif filetype(headerfile, 'plexon_ddt')
  orig = read_plexon_ddt(headerfile);
  hdr.nChans      = orig.NChannels;
  hdr.Fs          = orig.Freq;
  hdr.nSamples    = orig.NSamples;
  hdr.nSamplesPre = 0;      % continuous
  hdr.nTrials     = 1;      % continuous
  for i=1:hdr.nChans
    % create fake labels
    hdr.label{i} = sprintf('%03d', i);
  end
  % also remember the original header
  hdr.orig        = orig;

elseif filetype(headerfile, 'neuralynx_dma')
  hdr = read_neuralynx_dma(headerfile);

elseif filetype(headerfile, 'neuralynx_ncs')
  ncs = read_neuralynx_ncs(headerfile, 1, 0);
  [p, f, x] = fileparts(headerfile);
  hdr.Fs          = ncs.hdr.SamplingFrequency;
  hdr.label       = {f};
  hdr.nChans      = 1;
  hdr.nTrials     = 1;
  hdr.nSamplesPre = 0;
  hdr.nSamples    = ncs.NRecords * 512;
  hdr.orig        = ncs.hdr;

elseif filetype(headerfile, 'neuralynx_ds')
  % read the header for all channels
  orig = read_neuralynx_header(headerfile);
  % channels can start and end at different times, therefore it is neccessary to define a consistent sample numbering over all channels
  % note that timestamps are measured in units of 1e-6 seconds
  FirstSampleTimeStamp = max(orig.FirstTimeStamp);  % sample "1" is defined as the first sample that is present in all channels
  LastSampleTimeStamp  = min(orig.LastTimeStamp);   % sample "hdr.nSamples" is defined as the last sample that is present in all channels
  if any(diff(orig.SamplingFrequency))
    warning('sampling frequencies are not consistent, assuming 32556Hz');
    hdr.Fs        = 32556;
  else
    hdr.Fs        = orig.SamplingFrequency(1);
  end
  hdr.label       = orig.label;
  hdr.nChans      = length(orig.label);           % only the continuous channels (CSC)
  hdr.nTrials     = 1;                            % it is continuous
  hdr.nSamplesPre = 0;                            % it is continuous
  hdr.nSamples    = hdr.Fs * (LastSampleTimeStamp-FirstSampleTimeStamp)/1e6;
  hdr.nSamples    = floor(hdr.nSamples/512)*512;  % it should be a multiple of 512
  % also remember the original header details
  hdr.orig        = orig;

elseif filetype(headerfile, 'ced_son')
  if ~hasneuroshare
    error('cannot find CED SON import routines, see http://www.neuroshare.org');
  end;
  orig = read_ced_son(headerfile,'readevents','no','readdata','no');
  orig = orig.header;
  % In Spike2, channels can have different sampling rates, units, length
  % etc. etc. Here, channels need to have to same properties.
  if length(unique([orig.samplerate]))>1,
    error('channels have different sampling rates');
  else,
    hdr.Fs   = orig(1).samplerate;
  end;
  hdr.nChans = length(orig);
  % nsamples of the channel with least samples
  hdr.nSamples    = min([orig.nsamples]);
  hdr.nSamplesPre = 0;
  % only continuous data supported
  if sum(strcmp(lower({orig.mode}),'continuous')) < hdr.nChans,
    error('not all channels contain continuous data');
  else,
    hdr.nTrials = 1;
  end;
  hdr.label = {orig.label};  

elseif filetype(headerfile, 'eep_cnt')
  % read the first sample from the continous data, which will also return the header
  hdr = read_eep_cnt(headerfile, 1, 1);
  hdr.Fs          = hdr.rate;
  hdr.nSamples    = hdr.nsample;
  hdr.nSamplesPre = 0;
  hdr.nChans      = hdr.nchan;
  hdr.nTrials     = 1;		% it can always be interpreted as continuous data

elseif filetype(headerfile, 'eep_avr')
  % read the whole average and keep only header info (it is a bit silly, but the easiest to do here)
  hdr = read_eep_avr(headerfile);
  hdr.Fs          = hdr.rate;
  hdr.nChans      = size(hdr.data,1);;
  hdr.nSamples    = size(hdr.data,2);
  hdr.nSamplesPre = hdr.xmin*hdr.rate/1000;
  hdr.nTrials     = 1;		% it can always be interpreted as continuous data
  % remove the data and variance if present
  hdr = rmfield(hdr, 'data');
  try, hdr = rmfield(hdr, 'variance'); end

elseif filetype(headerfile, 'ns_eeg')
  hdr = read_ns_hdr(headerfile);
  % do some reformatting/renaming of the header items
  hdr.Fs          = hdr.rate;
  hdr.nSamples    = hdr.npnt;
  hdr.nSamplesPre = hdr.rate*hdr.xmin/1000;
  hdr.nChans      = hdr.nchan;
  hdr.nTrials     = hdr.nsweeps;

elseif filetype(headerfile, 'ns_cnt')
  % read_ns_cnt originates from the EEGLAB package (loadcnt.m) but is
  % an old version since the new version is not compatible any more
  orig = read_ns_cnt(headerfile, 'ldheaderonly', 1);
  % do some reformatting/renaming of the header items
  hdr.Fs          = orig.rate;
  hdr.nChans      = orig.nchannels;
  hdr.nSamples    = orig.nsamples;
  hdr.nSamplesPre = 0;
  hdr.nTrials     = 1;
  for i=1:hdr.nChans
    hdr.label{i} = deblank(orig.chan.names(i,:));
  end

elseif filetype(headerfile, 'besa_avr')
  orig = read_besa_avr(headerfile);
  hdr.Fs          = 1000/orig.di;
  hdr.nChans      = size(orig.data,1);
  hdr.nSamples    = size(orig.data,2);
  hdr.nSamplesPre = -(hdr.Fs * orig.tsb/1000);   % convert from ms to samples
  hdr.nTrials     = 1;
  if isfield(orig, 'label') && iscell(orig.label)
    hdr.label = orig.label;
  elseif isfield(orig, 'label') && ischar(orig.label)
    hdr.label = tokenize(orig.label, ' ');
  else
    for i=1:hdr.nChans
      warning('creating fake channel names for besa_avr');
      hdr.label{i} = sprintf('%03d', i);
    end
  end

elseif filetype(headerfile, 'neuromag_fif')
  % Neuromag MEG data
  if ~hasmegpd
    error('cannot find Neuromag import routines, see http://boojum.hut.fi/~kuutela/meg-pd/');
  end
  % add a gradiometer structure for forward and inverse modelling
  hdr.grad = fif2grad(headerfile);
  [hdr.label, type, number] = channames(headerfile);
  rawdata('any',headerfile);
  hdr.Fs = rawdata('sf');
  hdr.nTrials = 0;
  rawdata('goto', 0);
  [buf, status] = rawdata('next');
  hdr.nChans   = size(buf,1);
  hdr.nSamples = size(buf,2);
  while strcmp(status, 'ok')
    hdr.nTrials = hdr.nTrials + 1;
    [buf, status] = rawdata('next');
  end
  % I don't know how to get this out of the file
  hdr.nSamplesPre  = 0;
  % this would give some information that I don't know how to use
  % [range, cal] = rawdata('range');
  % this would give the current latency
  % rawdata('t');
  rawdata('close');

elseif filetype(headerfile, 'biosemi_bdf') || filetype(headerfile, 'bham_bdf')
  % this uses the openbdf and readbdf functions that I copied from the EEGLAB toolbox
  orig = openbdf(headerfile);
  if any(orig.Head.SampleRate~=orig.Head.SampleRate(1))
    error('channels with different sampling rate not supported');
  end
  hdr.Fs          = orig.Head.SampleRate(1);
  hdr.nChans      = orig.Head.NS;
  hdr.label       = cellstr(orig.Head.Label);
  % it is continuous data, therefore append all records in one trial
  hdr.nSamples    = orig.Head.NRec * orig.Head.Dur * orig.Head.SampleRate(1);
  hdr.nSamplesPre = 0;
  hdr.nTrials     = 1;
  hdr.orig        = orig;
  % close the file between seperate read operations
  fclose(orig.Head.FILE.FID); 
  if filetype(headerfile, 'bham_bdf')
    % this is for the Biosemi system used at the University of Birmingham
    labelold = { 'A1' 'A2' 'A3' 'A4' 'A5' 'A6' 'A7' 'A8' 'A9' 'A10' 'A11' 'A12' 'A13' 'A14' 'A15' 'A16' 'A17' 'A18' 'A19' 'A20' 'A21' 'A22' 'A23' 'A24' 'A25' 'A26' 'A27' 'A28' 'A29' 'A30' 'A31' 'A32' 'B1' 'B2' 'B3' 'B4' 'B5' 'B6' 'B7' 'B8' 'B9' 'B10' 'B11' 'B12' 'B13' 'B14' 'B15' 'B16' 'B17' 'B18' 'B19' 'B20' 'B21' 'B22' 'B23' 'B24' 'B25' 'B26' 'B27' 'B28' 'B29' 'B30' 'B31' 'B32' 'C1' 'C2' 'C3' 'C4' 'C5' 'C6' 'C7' 'C8' 'C9' 'C10' 'C11' 'C12' 'C13' 'C14' 'C15' 'C16' 'C17' 'C18' 'C19' 'C20' 'C21' 'C22' 'C23' 'C24' 'C25' 'C26' 'C27' 'C28' 'C29' 'C30' 'C31' 'C32' 'D1' 'D2' 'D3' 'D4' 'D5' 'D6' 'D7' 'D8' 'D9' 'D10' 'D11' 'D12' 'D13' 'D14' 'D15' 'D16' 'D17' 'D18' 'D19' 'D20' 'D21' 'D22' 'D23' 'D24' 'D25' 'D26' 'D27' 'D28' 'D29' 'D30' 'D31' 'D32' 'EXG1' 'EXG2' 'EXG3' 'EXG4' 'EXG5' 'EXG6' 'EXG7' 'EXG8' 'Status'};
    labelnew = { 'P9' 'PPO9h' 'PO7' 'PPO5h' 'PPO3h' 'PO5h' 'POO9h' 'PO9' 'I1' 'OI1h' 'O1' 'POO1' 'PO3h' 'PPO1h' 'PPO2h' 'POz' 'Oz' 'Iz' 'I2' 'OI2h' 'O2' 'POO2' 'PO4h' 'PPO4h' 'PO6h' 'POO10h' 'PO10' 'PO8' 'PPO6h' 'PPO10h' 'P10' 'P8' 'TPP9h' 'TP7' 'TTP7h' 'CP5' 'TPP7h' 'P7' 'P5' 'CPP5h' 'CCP5h' 'CP3' 'P3' 'CPP3h' 'CCP3h' 'CP1' 'P1' 'Pz' 'CPP1h' 'CPz' 'CPP2h' 'P2' 'CPP4h' 'CP2' 'CCP4h' 'CP4' 'P4' 'P6' 'CPP6h' 'CCP6h' 'CP6' 'TPP8h' 'TP8' 'TPP10h' 'T7' 'FTT7h' 'FT7' 'FC5' 'FCC5h' 'C5' 'C3' 'FCC3h' 'FC3' 'FC1' 'C1' 'CCP1h' 'Cz' 'FCC1h' 'FCz' 'FFC1h' 'Fz' 'FFC2h' 'FC2' 'FCC2h' 'CCP2h' 'C2' 'C4' 'FCC4h' 'FC4' 'FC6' 'FCC6h' 'C6' 'TTP8h' 'T8' 'FTT8h' 'FT8' 'FT9' 'FFT9h' 'F7' 'FFT7h' 'FFC5h' 'F5' 'AFF7h' 'AF7' 'AF5h' 'AFF5h' 'F3' 'FFC3h' 'F1' 'AF3h' 'Fp1' 'Fpz' 'Fp2' 'AFz' 'AF4h' 'F2' 'FFC4h' 'F4' 'AFF6h' 'AF6h' 'AF8' 'AFF8h' 'F6' 'FFC6h' 'FFT8h' 'F8' 'FFT10h' 'FT10'};
    % rename the channel labels
    for i=1:length(labelnew)
      chan = strmatch(labelold(i), hdr.label, 'exact');
      hdr.label(chan) = labelnew(chan);
    end
  end

elseif hasbiosig
  % as a final resort try reading the file using the BIOSIG toolbox
  try
    hdr = read_biosig_header(headerfile);
    if isfield(hdr, 'orig')
      % remove the original header information, since fieldtrip is not interested in it
      hdr = rmfield(hdr, 'orig');
    end
  catch
    error('error reading data using BIOSIG toolbox, see http://biosig.sourceforge.net/');
  end
end

if isempty(hdr)
  % none of the importers above was able to read the header
  error('unsupported fileformat for EEG/MEG header file');
end

