function [cfg, artifact] = artifact_muscle(cfg)

% ARTIFACT_MUSCLE scans the MEG channels for artifacts, due to neck muscle artifacts.
%
% [cfg, artifact] = artifact_muscle(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.
% Muscle artifacts are identified by means of thresholding the bandpass-filtered raw data.
% The artifact-definition parameters can be specified in cfg.artfctdef.muscle. If the 
% respective parameter is not specified, default values will be used.
%
% cfg.artfctdef.muscle.sgn       = 'MEG';
% cfg.artfctdef.muscle.bpfreq    = [110 140];
% cfg.artfctdef.muscle.bpfilttype= 'but';
% cfg.artfctdef.muscle.bpfiltord = 10;
% cfg.artfctdef.muscle.pretim    = 0.05;
% cfg.artfctdef.muscle.psttim    = 0.05;
% cfg.artfctdef.muscle.fltpadding= is determined based on the filter type and order (can be set to 0)
% cfg.artfctdef.muscle.padding   = 0.1;
% cfg.artfctdef.muscle.method    = 'zvalue' (default), no other options implemented.
% cfg.artfctdef.muscle.feedback  = 'no' (default), or 'yes'; displaying of raw data is 
%                                  not supported here.
% cfg.artfctdef.muscle.cutoff    = 4;
%
% if cfg.artfctdef.muscle.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_muscle_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.14  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.13  2005/06/29 12:42:00  roboos
% added version to the output configuration
%
% Revision 1.12  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.11  2005/04/08 12:08:42  roboos
% catch error when bpfilttype is not specified -> fltpadding=0
%
% Revision 1.10  2005/04/07 20:00:47  roboos
% replaced ocurrence of cfg.datatype (which might not exist) with iscontinuous
%
% Revision 1.9  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.8  2005/02/02 12:25:40  jansch
% incorporation of boundary check for trial-based data collections. inclusion
% of preproc-function for data handling.
%
% 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.2  2003/11/03 jansch
% cleaner implementation of the routine. renamed cutoff in method, and numstd in cutoff.
%
% 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 muscle artifacts if necessary.          
if ~isfield(cfg,'artfctdef'),               cfg.artfctdef                  = [];      end;
if ~isfield(cfg.artfctdef,'muscle'),        cfg.artfctdef.muscle           = [];      end;
if ~isfield(cfg.artfctdef.muscle,'sgn'),    cfg.artfctdef.muscle.sgn       = 'MEG';   end;
if ~isfield(cfg.artfctdef.muscle,'pretim'), cfg.artfctdef.muscle.pretim    = 0.05;    end;
if ~isfield(cfg.artfctdef.muscle,'psttim'), cfg.artfctdef.muscle.psttim    = 0.05;    end;
if ~isfield(cfg.artfctdef.muscle,'padding'),cfg.artfctdef.muscle.padding   = 0.1;     end;
if ~isfield(cfg.artfctdef.muscle,'method'), cfg.artfctdef.muscle.method    = 'zvalue';end;
if ~isfield(cfg.artfctdef.muscle,'feedback'),cfg.artfctdef.muscle.feedback = 'no';    end;

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

% method dependent settings
if strcmp(cfg.artfctdef.muscle.method,'zvalue'),
    if ~isfield(cfg.artfctdef.muscle,'bpfilter'),cfg.artfctdef.muscle.bpfilter = 'yes';   end;
    if ~isfield(cfg.artfctdef.muscle,'bpfreq'), cfg.artfctdef.muscle.bpfreq = [110 140];  end;
    if ~isfield(cfg.artfctdef.muscle,'bpfiltord'), cfg.artfctdef.muscle.bpfiltord = 10;         end;
    if ~isfield(cfg.artfctdef.muscle,'bpfilttype'),cfg.artfctdef.muscle.bpfilttype= 'but';      end;
    if ~isfield(cfg.artfctdef.muscle,'cutoff'), cfg.artfctdef.muscle.cutoff = 4;          end;
    if ~isfield(cfg.artfctdef.muscle,'hilbert'),cfg.artfctdef.muscle.hilbert = 'yes';     end;      
else 
    error('muscle-artifact rejection only works with method: zvalue\n');
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.muscle;
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;%compute number of samples per trial
artfctdef.trl= trl;
artfctdef.sgn= channelselection(artfctdef.sgn,hdr.label);
sgnind       = match_str(hdr.label, artfctdef.sgn);
nummusclesgn = length(sgnind);

if ~isfield(artfctdef, 'fltpadding')
  if iscontinuous
    if ~isfield(artfctdef, 'bpfilttype')
      fltpadding = 0;
    elseif strcmp(artfctdef.bpfilttype, 'but')
      fltpadding = 20*artfctdef.bpfiltord;
    elseif strcmp(artfctdef.bpfilttype, '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

fltdatavg    = zeros(nummusclesgn,1);
fltdatvar    = zeros(nummusclesgn,1);
fltdatstd    = zeros(nummusclesgn,1);
cumpernumsmp = 0;
cutfltdatstr = cell(numtrl,1);
cutdatstr = cell(numtrl,1);
artfctchn    = zeros(numtrl,1);
rejectall    = zeros(1,trl(end,2));
z_tdata      = zeros(1,trl(end,2));

%loop over all trials and calculate mean and std
for trllop = 1:numtrl  
     fprintf('scanning for muscle 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);

%compute z-values for each channel and calculate cumulative z-scores for each sample
%keep track of the channel with the highest individual z-value in each trial for displaying purposes.
for trllop = 1:numtrl
    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));
end;

clear cutfltdatstr;

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

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

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

