function [hdr] = read_header(filename, varargin);

% READ_HEADER is a wrapper around different EEG/MEG file importers
% directly supported formats are CTF, Neuromag, EEP, BrainVision,
% Neuroscan and Neuralynx.
%
% Use as
%   hdr = read_header(filename, ...)
%
% Additional options should be specified in key-value pairs and can be
%   'headerformat'   string
%
% 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 file format, additional header information can be
% returned in the hdr.orig subfield.
%
% See also READ_DATA, READ_EVENT

% Copyright (C) 2003-2006, Robert Oostenveld, F.C. Donders Centre
%
% $Log: read_header.m,v $
% Revision 1.12  2006/12/04 10:37:10  roboos
% added support for ns_avg
% fixed bug in ns_eeg (hdr.nSamplesPre was negative)
%
% Revision 1.11  2006/10/09 15:39:51  roboos
% renamed BTi channels from 'MEGxxx' into 'Axxx'
%
% Revision 1.10  2006/09/18 21:48:33  roboos
% implemented support for fcdc_matbin, i.e. a dataset consisting of a matlab file with header and events and a seperate binary datafile
%
% Revision 1.9  2006/09/18 14:52:14  roboos
% implemented bti2grad and added it to header
%
% Revision 1.8  2006/09/18 14:22:54  roboos
% implemented support for 4D-BTi dataformat
%
% Revision 1.7  2006/09/13 11:01:01  roboos
% changed text in a error message
%
% Revision 1.6  2006/08/28 10:13:03  roboos
% use seperate filetype_check_extension function instead of subfunction, removed subfunction
%
% Revision 1.5  2006/06/22 15:07:39  roboos
% fiuxed bug in label assignment for tsl/tsh/ttl
%
% Revision 1.4  2006/06/22 13:51:47  roboos
% fixed typo in code: datatype instead of datatyppe, thanks to Thilo
%
% Revision 1.3  2006/06/22 07:53:46  roboos
% remember the original neuroscan cnt header
%
% Revision 1.2  2006/06/19 10:32:16  roboos
% added 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

% get the options
headerformat = keyval('headerformat',  varargin);

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

% start with an empty header
hdr = [];

switch headerformat
  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, headerfile)
  filename     = headerfile;                % this function will read the header
  headerformat = filetype(filename);        % update the filetype
end

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

  case {'4d_pdf', '4d_m4d', '4d_xyz'}
    % check that the required low-level toolbox is available
    hastoolbox('4d-version', 1);
    orig = read4dhdr(headerfile);
    hdr.Fs          = orig.samp_rate;
    hdr.nChans      = orig.TotalChannels;
    hdr.nSamples    = orig.TotalSlices;
    hdr.nSamplesPre = -round(orig.FirstLatency*orig.samp_rate);
    hdr.nTrials     = orig.TotalEpochs;
    hdr.label = {};
    % start with intrinsic channel names
    for i=1:hdr.nChans
      hdr.label{i} = sprintf('chan%03d', i);
    end
    % rename known MEG channels
    for i=1:length(orig.MegIndexArray)
      indx = orig.MegIndexArray(i);
      % hdr.label{indx} = sprintf('MEG %03d', i);
      hdr.label{indx} = sprintf('A%d', i); % according to BTi convention
    end
    % rename known EEG channels
    for i=1:length(orig.EegIndexArray)
      indx = orig.EegIndexArray(i);
      hdr.label{indx} = sprintf('EEG%03d', i);
    end
    % rename known trigger channels
    for i=1:length(orig.TriggerIndex)
      indx = orig.TriggerIndex(i);
      hdr.label{indx} = sprintf('TRIG%03d', i);
    end
    hdr.label       = hdr.label(:);
    % add a gradiometer structure for forward and inverse modelling
    hdr.grad = bti2grad(orig);
    % remember original header details
    hdr.orig = orig;

  case 'besa_avr'
    orig = read_besa_avr(filename);
    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

  case 'besa_swf'
    orig = read_besa_swf(filename);
    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;
    hdr.label       = orig.label;

  case {'biosemi_bdf', 'bham_bdf'}
    % this uses the openbdf and readbdf functions that I copied from the EEGLAB toolbox
    orig = openbdf(filename);
    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(filename, '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

  case {'brainvision_vhdr', 'brainvision_seg', 'brainvision_eeg', 'brainvision_dat'}
    hdr = read_brainvision_vhdr(filename);

  case 'ced_son'
    % check that the required low-level toolbox is available
    hastoolbox('neuroshare', 1);
    % use the reading function supplied by Gijs van Elswijk
    orig = read_ced_son(filename,'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 with different sampling rates are not supported');
    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};

  case {'ctf_ds', 'ctf_meg4', 'ctf_res4', 'read_ctf_res4'} % the default reader for CTF is read_ctf_res4
    if strcmp(headerformat, 'ctf_ds')
      [p, f, x] = fileparts(filename);
      filename = fullfile(filename, [f '.res4']);
    elseif strcmp(headerformat, 'ctf_meg4')
      [p, f, x] = fileparts(filename);
      filename = fullfile(p, [f '.res4']);
    end
    % read it using the default CTF importer that originates from CTF and the FCDC
    orig = read_ctf_res4(filename);
    hdr.Fs           = orig.Fs;
    hdr.nChans       = orig.nChans;
    hdr.nSamples     = orig.nSamples;
    hdr.nSamplesPre  = orig.nSamplesPre;
    hdr.nTrials      = orig.nTrials;
    hdr.label        = orig.label;
    % add a gradiometer structure for forward and inverse modelling
    hdr.grad = ctf2grad(orig);
    % add the original header details
    hdr.orig = orig;

  case 'ctf_read_res4'
    % check that the required low-level toolbos ix available
    hastoolbox('eegsf', 1);
    if strcmp(headerformat, 'ctf_res4')
      [p, f, x] = fileparts(filename);
      filename = p;
    elseif strcmp(headerformat, 'ctf_meg4')
      [p, f, x] = fileparts(filename);
      filename = p;
    end
    % read it using the CTF importer from the NIH and Daren Weber
    orig = ctf_read_res4(filename, 0);
    % convert the header into a structure that FieldTrip understands
    hdr              = [];
    hdr.Fs           = orig.setup.sample_rate;
    hdr.nChans       = length(orig.sensor.info);
    hdr.nSamples     = orig.setup.number_samples;
    hdr.nSamplesPre  = orig.setup.pretrigger_samples;
    hdr.nTrials      = orig.setup.number_trials;
    for i=1:length(orig.sensor.info)
      hdr.label{i}   = orig.sensor.info(i).label;
    end
    hdr.label        = hdr.label(:);
    % add a gradiometer structure for forward and inverse modelling
    hdr.grad         = nimh2grad(orig);
    % add the original header details
    hdr.orig = orig;

  case 'eep_avr'
    % check that the required low-level toolbox is available
    hastoolbox('eeprobe', 1);
    % read the whole average and keep only header info (it is a bit silly, but the easiest to do here)
    hdr = read_eep_avr(filename);
    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

  case 'eep_cnt'
    % check that the required low-level toolbox is available
    hastoolbox('eeprobe', 1);
    % read the first sample from the continous data, which will also return the header
    hdr = read_eep_cnt(filename, 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

  case 'fcdc_matbin'
    % this is multiplexed data in a *.bin file, accompanied by a matlab file containing the header
    load(headerfile, 'hdr');

  case {'mpi_ds', 'mpi_dap'}
    hdr = read_mpi_ds(filename);

  case 'neuralynx_dma'
    hdr = read_neuralynx_dma(filename);
    hdr.TimeStampPerSample = (hdr.LastTimeStamp-hdr.FirstTimeStamp)/(hdr.nSamples-1);

  case 'neuralynx_sdma'
    hdr = read_neuralynx_sdma(filename);
    hdr.TimeStampPerSample = (hdr.LastTimeStamp-hdr.FirstTimeStamp)/(hdr.nSamples-1);

  case 'neuralynx_ncs'
    ncs = read_neuralynx_ncs(filename, 1, 0);
    [p, f, x]       = fileparts(filename);
    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;
    % FIXME add hdr.FirstTimeStamp and hdr.LastTimeStamp

  case 'neuralynx_nse'
    nse = read_neuralynx_nse(filename, 1, 0);
    [p, f, x]       = fileparts(filename);
    hdr.Fs          = nse.hdr.SamplingFrequency;
    hdr.label       = {f};
    hdr.nChans      = 1;
    hdr.nTrials     = nse.NRecords;  % each record contains one waveform
    hdr.nSamples    = 32;            % there are 32 samples in each waveform
    hdr.nSamplesPre = 0;
    hdr.orig        = nse.hdr;
    % FIXME add hdr.FirstTimeStamp and hdr.LastTimeStamp

  case {'neuralynx_ttl', 'neuralynx_tsl', 'neuralynx_tsh'}
    hdr             = [];
    hdr.Fs          = 32556;
    hdr.nChans      = 1;
    hdr.nSamples    = (filesize(filename)-8)/4;
    hdr.nSamplesPre = 1;
    hdr.nTrials     = 1;
    switch headerformat
      case 'neuralynx_ttl'
        hdr.label = {'Parallel_in'};
      case 'neuralynx_tsl'
        hdr.label = {'TimeStampLow'};
      case 'neuralynx_tsh'
        hdr.label = {'TimeStampHi'};
      otherwise
        error('unknown header format');
    end

  case 'neuralynx_ds'
    % read the header for all channels
    orig  = read_neuralynx_ds(filename);
    if any(diff(orig.SamplingFrequency))
      % the header of the ncs files sometimes specifies a sampling frequency of zero
      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
    % do a sanity check on the data
    if any(orig.NumRecords~=orig.NumRecords(1))
      % although channels can in principle have different length, that is not supported here
      warning('the different channels do not have the same number of records');
    end
    if any(orig.FirstTimeStamp~=orig.FirstTimeStamp(1))
      % although channels can in principle start at different times, that is not supported here
      warning('the different channels do not start at the same time');
    end
    if any(orig.LastTimeStamp~=orig.LastTimeStamp(1))
      % although channels can in principle end at different times, that is not supported here
      warning('the different channels do not end at the same time');
    end
    % each record contains 512 samples
    hdr.nSamples = orig.NumRecords(1)*512;
    % timestamps are measured in units of approximately 1us
    % in case of 32556 Hz sampling rate, there are approximately 30.7 timestamps per sample
    % in case of 1000 Hz sampling rate,  there are approximately 1000 timestamps per sample
    % note that the last timestamp in the original header corresponds with the
    % first sample of the last record, and not with the last sample
    hdr.FirstTimeStamp     = orig.FirstTimeStamp(1);
    hdr.TimeStampPerSample = (orig.LastTimeStamp(1)-orig.FirstTimeStamp(1))/((orig.NumRecords(1)-1)*512);
    hdr.LastTimeStamp      = orig.LastTimeStamp(1) + round(hdr.TimeStampPerSample*512);
    % remember the original header details
    hdr.orig = orig;

  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'));
    %hasttl  = any(filetype_check_extension({dirlist.name}, 'ttl'));   % seperate file with original Parallel_in
    %hastsl  = any(filetype_check_extension({dirlist.name}, 'tsl'));   % seperate file with original TimeStampLow
    %hastsh  = any(filetype_check_extension({dirlist.name}, 'tsh'));   % seperate file with original TimeStampHi

    if haslfp
      lfpdir = fullfile(filename, dirlist(find(filetype_check_extension({dirlist.name}, 'lfp'))).name);
      lfphdr = read_header(lfpdir);
      for i=1:lfphdr.nChans
        lfplab{i} = sprintf('lfp_%s', lfphdr.label{i});
      end
    else
      lfphdr = [];
      lfplab = {};
    end

    if hasmua
      muadir = fullfile(filename, dirlist(find(filetype_check_extension({dirlist.name}, 'mua'))).name);
      muahdr = read_header(muadir);
      for i=1:muahdr.nChans
        mualab{i} = sprintf('mua_%s', muahdr.label{i});
      end
    else
      muahdr = [];
      mualab = {};
    end

    if hasspike
      spikedir = fullfile(filename, dirlist(find(filetype_check_extension({dirlist.name}, 'spike'))).name);
      spikehdr = read_header(spikedir);
      for i=1:spikehdr.nChans
        spikelab{i} = sprintf('spike_%s', spikehdr.label{i});
      end
    else
      spikehdr = [];
      spikelab = {};
    end

    if haslfp
      % assume that the LFP header describes the MUA and the spikes as well
      hdr = lfphdr;
    else hasmua
      % assume that the MUA header describes the spikes as well
      hdr = muahdr;
    end

    % concatenate the lfp/mua/spike channels, they have already been renamed
    hdr.label  = cat(1, lfplab(:), mualab(:), spikelab(:));
    hdr.nChans = length(hdr.label);
    % remember the original headers
    hdr.orig    = {};
    hdr.orig{1} = lfphdr;
    hdr.orig{2} = muahdr;
    hdr.orig{3} = spikehdr;

  case 'neuromag_fif'
    % check that the required low-level toolbox is available
    hastoolbox('meg-pd', 1);
    % add a gradiometer structure for forward and inverse modelling
    hdr.grad = fif2grad(filename);
    [hdr.label, type, number] = channames(filename);
    rawdata('any',filename);
    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');

  case 'ns_avg'
    orig = read_ns_hdr(filename);
    % do some reformatting/renaming of the header items
    hdr.Fs          = orig.rate;
    hdr.nSamples    = orig.npnt;
    hdr.nSamplesPre = round(-orig.rate*orig.xmin/1000);
    hdr.nChans      = orig.nchan;
    hdr.nTrials     = 1; % the number of trials in this datafile is only one, i.e. the average
    % remember the original header details
    hdr.orig = orig;

  case '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(filename, '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
    % remember the original header details
    hdr.orig = orig;

  case 'ns_eeg'
    orig = read_ns_hdr(filename);
    % do some reformatting/renaming of the header items
    hdr.Fs          = orig.rate;
    hdr.nSamples    = orig.npnt;
    hdr.nSamplesPre = round(-orig.rate*orig.xmin/1000);
    hdr.nChans      = orig.nchan;
    hdr.nTrials     = orig.nsweeps;
    % remember the original header details
    hdr.orig = orig;

  case 'plexon_ddt'
    orig = read_plexon_ddt(filename);
    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;

  case 'plexon_nex'
    hdr = read_nex_header(filename);

  case {'yokogawa_ave', 'yokogawa_con', 'yokogawa_raw'}
    % chek that the required low-level toolbox is available
    hastoolbox('yokogawa', 1);
    hdr = read_yokogawa_header(filename);
    % add a gradiometer structure for forward and inverse modelling
    hdr.grad = yokogawa2grad(hdr);

  otherwise
    error('unsupported data format');
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SUBFUNCTION to determine the file size in bytes
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [siz] = filesize(filename);
l = dir(filename);
if l.isdir
  error(sprintf('"%s" is not a file', filename));
end
siz = l.bytes;