function [data] = besa2fieldtrip(filename);

% BESA2FIELDTRIP reads and converts various BESA datafiles into a FieldTrip
% data structure, which subsequently can be used for statistical analysis
% or other analysis methods implemented in Fieldtrip.
%
% Use as
%   [data] = besa2fieldtrip(filename)
% where the filename should point to a BESA datafile (or data that is
% exported by BESA) and where The output is a Matlab structure that is
% compatible with FieldTrip. The format of the structure depends on the
% the type of datafile:
%   *.avr is converted to a structure similar to the output of TIMELOCKANALYSIS
%   *.mul is converted to a structure similar to the output of TIMELOCKANALYSIS
%   *.tfc is converted to a structure similar to the output of FREQANALYSIS
%   *.dat is converted to a structure similar to the output of SOURCANALYSIS
%   *.pdg is converted to a structure similar to the output of ???
%
% See also EEGLAB2FIELDTRIP

% Copyright (C) 2005, Robert Oostenveld
%
% $Log: besa2fieldtrip.m,v $
% Revision 1.5  2006/02/23 10:28:16  roboos
% changed dimord strings for consistency, changed toi and foi into time and freq, added fixdimord where neccessary
%
% Revision 1.4  2006/02/01 12:26:00  roboos
% made all uses of dimord consistent with the common definition of data dimensions, see the fixdimord() function
%
% Revision 1.3  2005/09/01 09:25:33  roboos
% added beamformer source reconstruction
%
% Revision 1.2  2005/07/29 13:23:58  roboos
% various improvements and fixes
%
% Revision 1.1  2005/07/28 15:08:58  roboos
% new implementation, see documentation on website
%

if nargin<2
  fieldbox = '';
end

type = filetype(filename);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if strcmp(type, 'besa_avr')
  fprintf('reading ERP/ERF\n');
  % this should be similar to the output of TIMELOCKANALYSIS
  tmp = read_besa_avr(filename);
  % convert into a TIMELOCKANALYSIS compatible data structure
  data = [];
  data.label   = tmp.label(:);
  data.avg     = tmp.data;
  data.time    = (0:(tmp.npnt-1)) * tmp.di + tmp.tsb;
  data.time    = data.time / 1000; % convert to seconds
  data.fsample = 1000/tmp.di;
  data.dimord  = 'chan_time';

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
elseif strcmp(type, 'besa_mul')
  fprintf('reading ERP/ERF\n');
  % this should be similar to the output of TIMELOCKANALYSIS
  tmp = read_besa_mul(filename);
  % convert into a TIMELOCKANALYSIS compatible data structure
  data = [];
  data.label   = tmp.label(:);
  data.avg     = tmp.data;
  data.time    = (0:(tmp.TimePoints-1)) * tmp.SamplingInterval_ms_ + tmp.BeginSweep_ms_;
  data.time    = data.time / 1000; % convert to seconds
  data.fsample = 1000/tmp.SamplingInterval_ms_;
  data.dimord  = 'chan_time';

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
elseif strcmp(type, 'besa_tfc')
  % this should be similar to the output of FREQANALYSIS
  [ChannelLabels, Time, Frequency, Data, Info] = read_besa_tfc(filename);
  Nchan = size(ChannelLabels,1);
  % convert into a FREQANALYSIS compatible data structure
  data = [];
  data.time = Time;
  data.freq = Frequency;
  if isfield(Info, 'DataType') && strcmp(Info.DataType, 'COHERENCE_SQUARED')
    % it contains coherence between channel pairs
    fprintf('reading coherence between %d channel pairs\n', Nchan);
    for i=1:Nchan
      tmp = tokenize(deblank(ChannelLabels(i,:)), '-');
      data.labelcmb{i,1} = tmp{1};
      data.labelcmb{i,2} = tmp{2};
    end
    data.cohspctrm = permute(Data, [1 3 2]);
  else
    % it contains power on channels
    fprintf('reading power on %d channels\n', Nchan);
    for i=1:Nchan
      data.label{i} = deblank(ChannelLabels(i,:));
    end
    data.powspctrm = permute(Data, [1 3 2]);
  end
  data.dimord    = 'chan_freq_time';

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
elseif strcmp(type, 'besa_src')
  src = read_besa_src(filename);
  data.xgrid = linspace(src.X(1), src.X(2), src.X(3));
  data.ygrid = linspace(src.Y(1), src.Y(2), src.Y(3));
  data.zgrid = linspace(src.Z(1), src.Z(2), src.Z(3));
  data.avg.pow = src.vol;
  data.dim = size(src.vol);
  [X, Y, Z] = ndgrid(data.xgrid, data.ygrid, data.zgrid);
  data.pos = [X(:) Y(:) Z(:)];
  % cannot determine which voxels are inside the brain volume
  data.inside = 1:prod(data.dim);
  data.outside = [];

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
elseif strcmp(type, 'besa_pdg')
  % hmmm, I have to think about this one...
  error('sorry, pdg is not yet supported');

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
else
  error('unrecognized file format for importing BESA data');
end

% construct and add a configuration to the output
cfg = [];
cfg.filename = filename;
% add the version details of this function call to the configuration
try
  % get the full name of the function
  cfg.version.name = mfilename('fullpath');
catch
  % required for compatibility with Matlab versions prior to release 13 (6.5)
  [st, i] = dbstack;
  cfg.version.name = st(i);
end
cfg.version.id = '$Id: besa2fieldtrip.m,v 1.5 2006/02/23 10:28:16 roboos Exp $';
data.cfg = cfg;
