function [cfg, artifact] = artifact_jump(cfg)

% ARTIFACT_JUMP scans the MEG channel for jump artifacts.
%
% [cfg, artifact] = artifact_jump(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. Jump artifacts are identified by means
% of thresholding the derivative of the detrended raw data. 
% The artifact-definition parameters can be specified in cfg.artfctdef.jump. 
% If the respective parameter is not specified, default values will be used.
%
% cfg.artfctdef.jump.sgn         = 'MEG' (default), or a cell array of channel-names/channel groups;
% cfg.artfctdef.jump.pretim      = cfg.padding, or 0.1; (pre-artifact rejection-interval in seconds) 
%                                  this depends on the filterconfigurations:
%                                  if filtering is used on the raw data, longer padding and rejection
%                                  epochs are used, to avoid filter-ringing in the data.                                  
% cfg.artfctdef.jump.psttim      = cfg.padding, or 0.1; (post-artifact rejection-interval in seconds)
% cfg.artfctdef.jump.padding     = 0.5*cfg.padding, or 0; (amount of padding at the edge of the 
%                                                           original trials, to detect artifact at 
%                                                           the edges.)
% cfg.artfctdef.jump.method,     and  
% cfg.artfctdef.jump.cutoff    
%                                  there are two methods to do jump-artifact rejection: 
%                                  -cfg.artfctdef.jump.method = 'zvalue' (default)
%                                       the differentiated, detrended data are z-transformed
%                                       and thresholded according to the number of standard
%                                       deviations in cfg.artfctdef.jump.cutoff (default value = 20)
%                                  -cfg.artfctdef.jump.method = 'abs'
%                                       the absolute difference of the minimum and maximum within each
%                                       trial is thresholded, according to the value specified in
%                                       cfg.artfctdef.jump.cutoff (in tesla: default value, in case one
%                                       chooses this method = 5e-11)
%
% if cfg.artfctdef.jump.feedback = 'yes', 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

% $Log: artifact_jump_old.m,v $
% Revision 1.2  2006/02/27 16:58:53  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.8  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.7  2005/06/29 12:42:00  roboos
% added version to the output configuration
%
% Revision 1.6  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.5  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.4  2005/02/02 12:24:36  jansch
% incorporation of boundary check for trial-based data collection. inclusion
% of the preproc function for data handling, which by default now includes
% median-filtering.
%
% Revision 1.2  2004/08/25 09:19:56  jansch
% removed carriage returns
%
% Revision 1.1  2004/08/25 09:11:57  jansch
% added artifact-routines to CVS-repository.
%
% Update!! 2004/02/06 jansch
% fixed bug in identifying maximum z-values across sensors. plus additional small fix to
% make compatible with trial-based collected datasets 
%
% Revision 1.2  2003/11/03 jansch
% cleaner implementation of the routine. renamed cutoff in method, and numstd in cutoff.
% included option to do absolute rejection
%
% Revision 1.1 2003/10/20 marbau and jansch
% included option to be given visual feedback and
% to set thresholds interactively

% set default rejection parameters for jump artifacts if necessary.          
if ~isfield(cfg,'artfctdef'),               cfg.artfctdef               = [];              end;
if ~isfield(cfg.artfctdef,'jump'),          cfg.artfctdef.jump          = [];              end;
if ~isfield(cfg.artfctdef.jump,'sgn'),      cfg.artfctdef.jump.sgn      = 'MEG';           end;
if ~isfield(cfg.artfctdef.jump,'method'),   cfg.artfctdef.jump.method   = 'zvalue';        end;
if ~isfield(cfg.artfctdef.jump,'feedback'), cfg.artfctdef.jump.feedback = 'no';            end;
if ~isfield(cfg.artfctdef.jump,'medianfilter'), cfg.artfctdef.jump.medianfilter  = 'yes';  end;
if ~isfield(cfg.artfctdef.jump,'medianfiltord'),cfg.artfctdef.jump.medianfiltord = 9;      end;

% padding dependent settings
if ~isfield(cfg,'padding'),	             cfg.padding = 0;                              end;
if cfg.padding~=0 
   if ~isfield(cfg.artfctdef.jump,'pretim'), cfg.artfctdef.jump.pretim  = cfg.padding;     end;
   if ~isfield(cfg.artfctdef.jump,'psttim'), cfg.artfctdef.jump.psttim  = cfg.padding;     end;
   if ~isfield(cfg.artfctdef.jump,'padding'),cfg.artfctdef.jump.padding = 0.5*cfg.padding; end;
else
   if ~isfield(cfg.artfctdef.jump,'pretim'), cfg.artfctdef.jump.pretim  = 0;               end;
   if ~isfield(cfg.artfctdef.jump,'psttim'), cfg.artfctdef.jump.psttim  = 0;               end;
   if ~isfield(cfg.artfctdef.jump,'padding'),cfg.artfctdef.jump.padding = 0;               end;
end;

% method dependent settings
artfctdef.fltord = 0;
if strcmp(cfg.artfctdef.jump.method,'zvalue'),
    if ~isfield(cfg.artfctdef.jump,'cutoff'), cfg.artfctdef.jump.cutoff    = 20;           end;
elseif strcmp(cfg.artfctdef.jump.method,'abs'),
    if ~isfield(cfg.artfctdef.jump,'cutoff'), cfg.artfctdef.jump.cutoff    = 5e-11;        end;
    if ~isfield(cfg.artfctdef.jump,'blc'),    cfg.artfctdef.jump.blc       = '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

fprintf('\nscanning for jump artifacts');

artfctdef     = cfg.artfctdef.jump;
cfg           = dataset2files(cfg);
hdr           = read_fcdc_header(cfg.headerfile);
trl           = cfg.trl;
numtrl        = size(trl,1);
padding       = round(artfctdef.padding*hdr.Fs);
trl(:,1)      = trl(:,1)-padding;
trl(:,2)      = trl(:,2)+padding;
if ~isempty(find(trl(:,1)<1)),
  fprintf('\nWARNING: required padding goes beyond begin of dataset');   
  trl(find(trl(:,1)<1),1) = 1;
end  
nsmp          = trl(:,2) - trl(:,1) + 1;%compute number of samples per trial
artfctdef.sgn = channelselection(artfctdef.sgn,hdr.label);
sgnind        = match_str(hdr.label, artfctdef.sgn);
artfctdef.trl = trl;
cutfltdatstr  = cell(numtrl,1);
artfctchn     = zeros(numtrl,1);
z_tdata       = zeros(1,trl(end,2));

if ~isfield(artfctdef, 'fltpadding')
  fltpadding = 0;
else
  fltpadding = round(artfctdef.fltpadding*hdr.Fs); 
end

% loop over all signals on which jump-artifact detection is required.
for sgnlop = 1:length(sgnind)
   fprintf('\nreading in signal %d from %d', sgnlop, length(sgnind)); 
   fltdatavg    = 0;
   fltdatstd    = 0;
   fltdatvar    = 0;
   cumpernumsmp = 0;
   % loop over the trials and calculate mean and std
   for trllop = 1:numtrl 
       dat = read_fcdc_data(cfg.datafile, hdr, trl(trllop,1)-fltpadding, ...
            trl(trllop,2)+fltpadding, sgnind(sgnlop), iscontinuous);%read data
       if strcmp(artfctdef.method,'zvalue'),
           fltdat = preproc(dat, artfctdef.sgn, hdr.Fs, artfctdef, [], fltpadding, fltpadding);
           fltdat = [abs(diff(fltdat,1,2)) 0];
       elseif strcmp(artfctdef.method,'abs'),
            fltdat = preproc(dat, artfctdef.sgn, hdr.Fs, artfctdef, [], fltpadding, fltpadding);
           fltdat = fltdat - (min(fltdat));
           fltdat = fltdat(1:end-1);
       end;
       cutfltdatstr{trllop} = fltdat; 
       fltdatavg = fltdatavg + mean(fltdat,2) .* nsmp(trllop); 
       fltdatvar = fltdatvar + (std(fltdat,[],2)).^2 .* nsmp(trllop); 
       cumpernumsmp = cumpernumsmp + nsmp(trllop);
   end

    fltdatavg = fltdatavg ./ cumpernumsmp;
    fltdatstd = sqrt(fltdatvar ./ cumpernumsmp);
	
    % select for each sample the highest z-value, and keep track of the highest z-value for the
    % respective channel in each trial
    for trllop = 1:numtrl
        if strcmp(artfctdef.method,'zvalue'), 
            z_tdat = (cutfltdatstr{trllop} - repmat(fltdatavg,[1 nsmp(trllop)])) ./ ...
                repmat(fltdatstd,[1 nsmp(trllop)]);
        elseif strcmp(artfctdef.method,'abs'),
            z_tdat = cutfltdatstr{trllop};    
        end;
        dummy = sort([z_tdat;z_tdata(1,trl(trllop,1):trl(trllop,2))]);
        z_tdata(1,trl(trllop,1):trl(trllop,2)) = dummy(end,:); 
        maxval(trllop,sgnlop) = max(z_tdat);
    end;
end;

% compute for each trial the channel-index with the maximum z-value
for trllop=1:numtrl
  sgnind = find(maxval(trllop,:) == max(maxval(trllop,:))); 
  artfctchn(trllop) = sgnind;
end;

clear cutfltdatstr;

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

[artfctdef, artifact] = artifact_feedback_old(cfg, z_tdata, artfctchn, 'jump');

% remember the details that were used here
cfg.artfctdef.jump          = artfctdef;
cfg.artfctdef.jump.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_jump_old.m,v 1.2 2006/02/27 16:58:53 roboos Exp $';

