function [cfg, artifact] = artifact_eog(cfg)

% ARTIFACT_EOG scans the EOG channel(s) for eye movements.
%
% [cfg, artifact] = artifact_eog(cfg)
%
% Artifact is a matrix, the first column of which specifying the beginsamples of an artifact period. 
% The second column contains the endsamples of the artifactperiods.
% Eye movements are identified by means of thresholding the absolute of the hilbert-transformed 
% bandpass-filtered raw data. This procedure works optimally in the case of continuously recorded data. 
% The artifact-definition parameters can be specified in cfg.artfctdef.eog. If the respective
% parameter is not specified, default values will be used.
%
% cfg.artfctdef.eog.sgn          = 'EOG';
% cfg.artfctdef.eog.bpfreq       = [1 15]; passband of the bandpassfilter
% cfg.artfctdef.eog.bpfiltord    = 4;      filter-order of the bandpassfilter
% cfg.artfctdef.eog.pretim       = 0.1;    pre-artifact rejection-interval in seconds
% cfg.artfctdef.eog.psttim       = 0.1;    post-artifact rejection-interval in seconds
% cfg.artfctdef.eog.fltpadding   = is determined based on the filter type and order (can be set to 0)
% cfg.artfctdef.eog.padding      = 0.5;    amount of padding at the edge of the original
%                                          trials, to detect artifact at the edges. Note that
%                                          if the padding is reduced, the Hilbert transformation
%                                          might introduce artifacts at the edges of the data,
%                                          which in fact are not there. Note also that the padding
%                                          does not take the trial boundaries into account if you have
%                                          your data collected on a trial by trial basis.
%                                          
% cfg.artfctdef.eog.feedback     = 'no' (default), 'yes', or 'flt', in case feedback
%                                  about the filtered signal is required. 
% cfg.artfctdef.eog.method       = 'zvalue' (default), 'percnt', 'abs'
% cfg.artfctdef.eog.cutoff       = to be specified for all:
%       
%           in case of zvalue as 4 or 5 std dev.
%           in case of percnt  0.95
%           in case of abs in volts or tesla
%          
% The cutoff-criterium can be specified according to cfg.artfctdef.eog.method, and can be either in 
% 'zvalue', 'percnt', or 'abs'. The threshold is set in cfg.artfctdef.eog.cutoff, and is 
% treated as a zvalue by default.
%
% if cfg.artfctdef.eog.feedback = 'yes' or 'flt', feedback about the artifact-epochs is given, and
% the rejection threshold can be set interactively. The new value is updated in cfg.artfctdef
%
% Copyright (c) 2003, F.C. Donders Centre

% p.m. think about optimal hilbert procedure

% $Log: artifact_eog_old.m,v $
% Revision 1.2  2006/02/27 16:58:52  roboos
% added dataset2files prior to reading the headerfile
%
% Revision 1.1  2006/01/12 13:35:58  roboos
% renamed JMs artifact_eog/jump/muscle functions to xxx_old
% also internally the dependency on private/artifact_threshold has been renamed
% these functions only remain present for people who insist on using the old detection routines
%
% Revision 1.19  2005/12/20 08:36:47  roboos
% add the artifact Nx2 matrix to the output configuration
% changed some indentation and white space, renamed a few variables
%
% Revision 1.18  2005/09/19 08:31:11  jansch
% fixed small bug in filterpadding
%
% Revision 1.17  2005/06/29 12:42:00  roboos
% added version to the output configuration
%
% Revision 1.16  2005/05/17 17:50:36  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/04/08 12:08:42  roboos
% catch error when bpfilttype is not specified -> fltpadding=0
%
% Revision 1.14  2005/04/07 20:00:47  roboos
% replaced ocurrence of cfg.datatype (which might not exist) with iscontinuous
%
% Revision 1.13  2005/04/07 19:53:00  roboos
% fixed another typo, added check+errormessage if no channels selected
%
% Revision 1.12  2005/04/07 19:47:15  roboos
% changed bpfltpadding into fltpadding (did not work due to recent change)
%
% Revision 1.11  2005/04/06 13:12:24  roboos
% fixed incorrect default cfgs for flttype and fltord (should be bpfilttype and bpfiltord)
% changed default settings for fltpadding, to be on the safe side for trial based data
%
% Revision 1.10  2005/02/04 14:58:26  jansch
% *** empty log message ***
%
% Revision 1.9  2005/02/04 14:57:32  jansch
% fixed bug in setting of defaults for 'abs' and 'percent' methods.
%
% Revision 1.8  2005/02/02 12:23:07  jansch
% incorporation of boundary check for trial-based data collections +
% inclusion of the preproc-function for handling the data.
%
% Revision 1.6  2004/08/25 09:19:56  jansch
% removed carriage returns
%
% Revision 1.5  2004/08/25 09:11:57  jansch
% added artifact-routines to CVS-repository.
%
% Revision 1.3  2003/11/03 jansch
% cleaner implementation of the routine. renamed cutoff in method, and numstd in cutoff.
% included option to set an absolute rejection-threshold
% 
% Revision 1.2  2003/10/20 marbau and jansch
% included option to be given visual feedback and
% to set thresholds interactively
%
% Revision 1.1  2003/05/21 10:14:33  roberto
% renamed all artefact into artifact

% set default rejection parameters for eog artifacts if necessary.
if ~isfield(cfg,'artfctdef'),            cfg.artfctdef               = [] ;      end;
if ~isfield(cfg.artfctdef,'eog'),        cfg.artfctdef.eog           = [];       end;
if ~isfield(cfg.artfctdef.eog,'sgn'),    cfg.artfctdef.eog.sgn       = {'EOG'};  end;
if ~isfield(cfg.artfctdef.eog,'method'), cfg.artfctdef.eog.method    = 'zvalue'; end;

% for backward-compatibility
if isfield(cfg.artfctdef.eog,'pssbnd'),  
   cfg.artfctdef.eog.bpfreq   = cfg.artfctdef.eog.pssbnd; 
   cfg.artfctdef.eog          = rmfield(cfg.artfctdef.eog,'pssbnd');
end;

% method independent settings
if ~isfield(cfg.artfctdef.eog,'pretim'), cfg.artfctdef.eog.pretim    = 0.1;      end;
if ~isfield(cfg.artfctdef.eog,'psttim'), cfg.artfctdef.eog.psttim    = 0.1;      end;
if ~isfield(cfg.artfctdef.eog,'padding'),cfg.artfctdef.eog.padding   = 0.5;      end;
if ~isfield(cfg.artfctdef.eog,'feedback'),cfg.artfctdef.eog.feedback = 'no';     end;

% method dependent settings
if strcmp(cfg.artfctdef.eog.method,'zvalue'),
    if ~isfield(cfg.artfctdef.eog,'cutoff'),   cfg.artfctdef.eog.cutoff   = 4;      end;
    if ~isfield(cfg.artfctdef.eog,'bpfilter'), cfg.artfctdef.eog.bpfilter = 'yes';  end;
    if ~isfield(cfg.artfctdef.eog,'bpfreq'),   cfg.artfctdef.eog.bpfreq   = [1 15]; end;
    if ~isfield(cfg.artfctdef.eog,'bpfiltord'),cfg.artfctdef.eog.bpfiltord= 4;      end;
    if ~isfield(cfg.artfctdef.eog,'bpfilttype'),cfg.artfctdef.eog.bpfilttype= 'but';  end;
    if ~isfield(cfg.artfctdef.eog,'hilbert'),  cfg.artfctdef.eog.hilbert  = 'yes';  end;
elseif strcmp(cfg.artfctdef.eog.method,'percnt'),
    if ~isfield(cfg.artfctdef.eog,'cutoff'),   cfg.artfctdef.eog.cutoff   = 0.995;  end;
    if ~isfield(cfg.artfctdef.eog,'bpfilter'), cfg.artfctdef.eog.bpfilter = 'yes';  end;
    if ~isfield(cfg.artfctdef.eog,'bpfreq'),   cfg.artfctdef.eog.bpfreq   = [1 15]; end;
    if ~isfield(cfg.artfctdef.eog,'bpfiltord'),cfg.artfctdef.eog.bpfiltord= 4;      end;
    if ~isfield(cfg.artfctdef.eog,'bpfilttype'),cfg.artfctdef.eog.bpfilttype= 'but';  end;
    if ~isfield(cfg.artfctdef.eog,'hilbert'),  cfg.artfctdef.eog.hilbert  = 'yes';  end;
elseif strcmp(cfg.artfctdef.eog.method,'abs'),
    if ~isfield(cfg.artfctdef.eog,'cutoff'),   cfg.artfctdef.eog.cutoff   = 50e-6;  end;
    if ~isfield(cfg.artfctdef.eog,'hpfilter'), cfg.artfctdef.eog.bpfilter = 'yes';  end;
    if ~isfield(cfg.artfctdef.eog,'hpfreq'),   cfg.artfctdef.eog.bpfreq   = 1;      end;
    if ~isfield(cfg.artfctdef.eog,'bpfiltord'),cfg.artfctdef.eog.bpfiltord= 4;      end;
    if ~isfield(cfg.artfctdef.eog,'bpfilttype'),cfg.artfctdef.eog.bpfilttype= 'but';  end;
    if ~isfield(cfg.artfctdef.eog,'hilbert'),  cfg.artfctdef.eog.hilbert  = 'yes';  end;
end    

if ~isfield(cfg, 'datatype') || ~strcmp(cfg.datatype, 'continuous')
  % datatype is unknown or not continuous, perform epoch boundary check
  iscontinuous = 0;
else
  % do not perform epoch boundary check, usefull for pseudo-continuous data
  iscontinuous = strcmp(cfg.datatype, 'continuous');
end

artfctdef    = cfg.artfctdef.eog;
cfg          = dataset2files(cfg);
hdr          = read_fcdc_header(cfg.headerfile); 
padsmp       = round(artfctdef.padding*hdr.Fs);
trl          = cfg.trl;
trl(:,1)     = trl(:,1) - padsmp; % pad the trial with padsmp-samples, in order to detect
trl(:,2)     = trl(:,2) + padsmp; % artifacts at the edges of the relevant trials.
numtrl       = size(trl,1);
trllength    = trl(:,2) - trl(:,1) + 1;
artfctdef.trl= trl;
artfctdef.sgn= channelselection(artfctdef.sgn,hdr.label);
sgnind       = match_str(hdr.label, artfctdef.sgn);
numeogsgn    = length(sgnind);

if numeogsgn<1
  error('no EOG channels selected');
end

if ~isfield(artfctdef, 'fltpadding')
  if iscontinuous
    if ~isfield(artfctdef, 'bpfilttype') || ~isfield(artfctdef, 'hpfilttype'),
       fltpadding = 0;
    elseif strcmp(artfctdef.bpfilttype, 'but') || strcmp(artfctdef.hpfilttype, 'but'),
      fltpadding = 20*artfctdef.bpfiltord;
    elseif strcmp(artfctdef.bpfilttype, 'fir') || strcmp(artfctdef.hpfilttype, 'fir'),
      fltpadding = artfctdef.bpfiltord;
    else
      warning('unknown filter type, cannot determine filter padding');
      fltpadding = 0;
    end
  else
    warning('default is not to apply filter padding on trial based data');
    fltpadding = 0;
  end
else
  fltpadding = round(artfctdef.fltpadding*hdr.Fs);
end

% initialize some variables
fltdatavg    = zeros(numeogsgn,1);
fltdatvar    = zeros(numeogsgn,1);
fltdatstd    = zeros(numeogsgn,1);
cumpernumsmp = 0;
cutfltdatstr = cell(numtrl,1);
artfctchn    = zeros(numtrl,1);
rejectall    = zeros(1,trl(end,2));
z_tdata      = zeros(1,trl(end,2));

for trllop = 1:numtrl  
     fprintf('scanning for eog artifacts in trial %d from %d\n', trllop, numtrl);
     dat = read_fcdc_data(cfg.datafile, hdr, trl(trllop,1)-fltpadding, ...
         trl(trllop,2)+fltpadding, sgnind, iscontinuous);%read data
     fltdat = preproc(dat, artfctdef.sgn, hdr.Fs, artfctdef, [], fltpadding, fltpadding);
     cutfltdatstr{trllop} = fltdat; 
     fltdatavg = fltdatavg + mean(fltdat,2) .* trllength(trllop); 
     fltdatvar = fltdatvar + (std(fltdat,[],2)).^2 .* trllength(trllop); 
     cumpernumsmp = cumpernumsmp + trllength(trllop);
end

fltdatavg = fltdatavg ./ cumpernumsmp;
fltdatstd = sqrt(fltdatvar ./ cumpernumsmp);

for trllop = 1:numtrl
    if strcmp(artfctdef.method,'zvalue'),
        dummy = ((cutfltdatstr{trllop} - repmat(fltdatavg,1,trllength(trllop))) ./ ...
                repmat(fltdatstd,1,trllength(trllop)));
        [dummy, index] = sort(dummy,1);
        [dum indx] = max(dummy(end,:));
        artfctchn(trllop)=index(end,indx);
        z_tdata(1,trl(trllop,1):trl(trllop,2)) = sum(dummy,1)./sqrt(size(dummy,1));
    elseif strcmp(artfctdef.method,'abs'),
        dummy = cutfltdatstr{trllop};
        [dummy,index] = sort(dummy);
        [dum indx] = max(dummy(end,:));
        artfctchn(trllop)=index(end,indx);
        z_tdata(1,trl(trllop,1):trl(trllop,2)) = sum(dummy,1)./sqrt(size(dummy,1));
    end;    
end;

if strcmp(artfctdef.method,'percnt'),
    dummy = zeros(numeogsgn,trl(end,2));
    for trllop = 1:numtrl
        dummy(:,trl(trllop,1):trl(trllop,2)) = cutfltdatstr{trllop};
        clear cutfltdatstr{trllop};    
    end;
    
    for sgnlop = 1:numeogsgn
            [sd,id]         = sort(dummy(sgnlop,:),2);
            sel             = find(sd);
            dummy(sgnlop,id(sel)) = [1:length(sel)]./length(sel);                
    end
    [dummy, index] = sort(dummy,1);
    
    for trllop = 1:numtrl
        dum = dummy(:,trl(trllop,1):trl(trllop,2));
        [dum,index] = sort(dum,1);
        artfctchn(trllop)=index(end,find(dum(end,:)==max(dum(end,:))));
        z_tdata(1,trl(trllop,1):trl(trllop,2)) = dum(end,:);
    end
end

clear cutfltdatstr;

% update the configuration prior to calling artifact_feedback
cfg.artfctdef.eog = artfctdef;

if (strcmp(cfg.artfctdef.eog.feedback,'yes'))
    [artfctdef,artifact] = artifact_feedback_old(cfg, z_tdata, artfctchn,'EOGraw');
elseif (strcmp(cfg.artfctdef.eog.feedback,'flt'))
    [artfctdef,artifact] = artifact_feedback_old(cfg, z_tdata, artfctchn,'EOGflt');
else    
    [artfctdef,artifact] = artifact_feedback_old(cfg, z_tdata, artfctchn,'EOG');   
end;

% remember the details that were used here
cfg.artfctdef.eog          = artfctdef;
cfg.artfctdef.eog.artifact = artifact;

% add version information 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: artifact_eog_old.m,v 1.2 2006/02/27 16:58:52 roboos Exp $';

