Skip to content
Permalink
release
Switch branches/tags
Go to file
 
 
Cannot retrieve contributors at this time
function [cfg, artifact] = ft_artifact_threshold(cfg, data)
% FT_ARTIFACT_THRESHOLD scans data segments of interest for channels in which the
% signal exceeds a specified minimum or maximum value, or in which the peak-to-peak
% range within the trial exceeds a specified threshold.
%
% Use as
% [cfg, artifact] = ft_artifact_threshold(cfg)
% with the configuration options
% cfg.dataset = string with the filename
% or
% cfg.headerfile = string with the filename
% cfg.datafile = string with the filename
% and optionally
% cfg.headerformat
% cfg.dataformat
%
% Alternatively you can use it as
% [cfg, artifact] = ft_artifact_threshold(cfg, data)
% where the input data is a structure as obtained from FT_PREPROCESSING.
%
% In both cases the configuration should also contain
% cfg.trl = structure that defines the data segments of interest, see FT_DEFINETRIAL
% cfg.continuous = 'yes' or 'no' whether the file contains continuous data
% and
% cfg.artfctdef.threshold.channel = cell-array with channel labels
% cfg.artfctdef.threshold.bpfilter = 'no' or 'yes' (default = 'yes')
% cfg.artfctdef.threshold.bpfreq = [0.3 30]
% cfg.artfctdef.threshold.bpfiltord = 4
%
% In the same way as specifying the options for band-pass filtering, it is also
% possible to specify lpfilter, hpfilter, bsfilter, dftfilter or medianfilter, see
% FT_PREPROCESSING.
%
% The detection of artifacts is done according to the following settings,
% you should specify at least one of these thresholds
% cfg.artfctdef.threshold.min = value in uV or T, default -inf
% cfg.artfctdef.threshold.max = value in uV or T, default inf
% cfg.artfctdef.threshold.onset = value in uV or T, default inf
% cfg.artfctdef.threshold.offset = value in uV or T, default inf
%
% When cfg.artfctdef.threshold.onset and offset are used, the rising and falling
% flank are thresholded with different values. In case onset and offset are both
% positive, the data will be thresholded above their values. In case both onset and
% offset are negative, the data will be thresholded below their values.
%
% Note that this function does not support artifactpadding or filterpadding.
%
% The output argument "artifact" is a Nx2 matrix comparable to the "trl" matrix of
% FT_DEFINETRIAL. The first column of which specifying the beginsamples of an
% artifact period, the second column contains the endsamples of the artifactperiods.
%
% To facilitate data-handling and distributed computing, you can use
% cfg.inputfile = ...
% to read the input data from a *.mat file on disk. This mat files should contain
% only a single variable named 'data', corresponding to the input structure.
%
% See also FT_REJECTARTIFACT, FT_ARTIFACT_CLIP, FT_ARTIFACT_ECG, FT_ARTIFACT_EOG,
% FT_ARTIFACT_JUMP, FT_ARTIFACT_MUSCLE, FT_ARTIFACT_THRESHOLD, FT_ARTIFACT_ZVALUE
% Copyright (C) 2003-2011, Robert Oostenveld, SMI, FCDC
%
% This file is part of FieldTrip, see http://www.fieldtriptoolbox.org
% for the documentation and details.
%
% FieldTrip is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% FieldTrip is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with FieldTrip. If not, see <http://www.gnu.org/licenses/>.
%
% $Id$
% these are used by the ft_preamble/ft_postamble function and scripts
ft_revision = '$Id$';
ft_nargin = nargin;
ft_nargout = nargout;
% do the general setup of the function
ft_defaults
ft_preamble init
ft_preamble provenance
ft_preamble loadvar data
% the ft_abort variable is set to true or false in ft_preamble_init
if ft_abort
return
end
% check if the input cfg is valid for this function
cfg = ft_checkconfig(cfg, 'renamed', {'datatype', 'continuous'});
cfg = ft_checkconfig(cfg, 'renamedval', {'continuous', 'continuous', 'yes'});
% set the default options
cfg.continuous = ft_getopt(cfg, 'continuous', []);
cfg.headerformat = ft_getopt(cfg, 'headerformat', []);
cfg.dataformat = ft_getopt(cfg, 'dataformat', []);
cfg.feedback = ft_getopt(cfg, 'feedback', 'text');
cfg.representation = ft_getopt(cfg, 'representation', 'numeric'); % numeric or table
% set the default artifact detection parameters
cfg.artfctdef = ft_getopt(cfg, 'artfctdef');
cfg.artfctdef.threshold = ft_getopt(cfg.artfctdef, 'threshold');
cfg.artfctdef.threshold.channel = ft_getopt(cfg.artfctdef.threshold, 'channel', 'all');
cfg.artfctdef.threshold.bpfilter = ft_getopt(cfg.artfctdef.threshold, 'bpfilter', 'yes');
cfg.artfctdef.threshold.bpfreq = ft_getopt(cfg.artfctdef.threshold, 'bpfreq', [0.3 30]);
cfg.artfctdef.threshold.bpfiltord = ft_getopt(cfg.artfctdef.threshold, 'bpfiltord', 4);
cfg.artfctdef.threshold.range = ft_getopt(cfg.artfctdef.threshold, 'range', inf);
cfg.artfctdef.threshold.min = ft_getopt(cfg.artfctdef.threshold, 'min', -inf);
cfg.artfctdef.threshold.max = ft_getopt(cfg.artfctdef.threshold, 'max', inf);
cfg.artfctdef.threshold.onset = ft_getopt(cfg.artfctdef.threshold, 'onset', []);
cfg.artfctdef.threshold.offset = ft_getopt(cfg.artfctdef.threshold, 'offset', []);
% the data is either passed into the function by the user or read from file with cfg.inputfile
hasdata = exist('data', 'var');
% read the header, or get it from the input data
if ~hasdata
cfg = ft_checkconfig(cfg, 'dataset2files', 'yes');
cfg = ft_checkconfig(cfg, 'required', {'headerfile', 'datafile'});
hdr = ft_read_header(cfg.headerfile, 'headerformat', cfg.headerformat);
else
data = ft_checkdata(data, 'datatype', 'raw', 'hassampleinfo', 'yes');
cfg = ft_checkconfig(cfg, 'forbidden', {'dataset', 'headerfile', 'datafile'});
hdr = ft_fetch_header(data);
end
% set default cfg.continuous
if isempty(cfg.continuous)
if hdr.nTrials==1
cfg.continuous = 'yes';
else
cfg.continuous = 'no';
end
end
% get the specification of the data segments that should be scanned for artifacts
if ~isfield(cfg, 'trl') && hasdata
trl = data.sampleinfo;
for k = 1:numel(data.trial)
trl(k,3) = time2offset(data.time{k}, data.fsample);
end
elseif isfield(cfg, 'trl') && ischar(cfg.trl)
trl = loadvar(cfg.trl, 'trl');
elseif isfield(cfg, 'trl') && isnumeric(cfg.trl)
trl = cfg.trl;
else
ft_error('cannot determine which segments of data to scan for artifacts');
end
if ~isempty(cfg.artfctdef.threshold.onset) || ~isempty(cfg.artfctdef.threshold.offset)
if cfg.artfctdef.threshold.onset>0 && cfg.artfctdef.threshold.offset>0
direction = 'up';
elseif cfg.artfctdef.threshold.onset<0 && cfg.artfctdef.threshold.offset<0
direction = 'down';
else
error('incorrect specification of onset and offset');
end
else
direction = 'none';
end
% get the remaining settings
artfctdef = cfg.artfctdef.threshold;
artfctdef.trl = trl;
ntrial = size(trl,1);
label = ft_channelselection(artfctdef.channel, hdr.label);
chanindx = match_str(hdr.label, label);
nchan = length(chanindx);
artifact = table();
% only print this once
ft_notice once FieldTrip:applyingPreprocessing
ft_progress('init', cfg.feedback, ['searching for artifacts in ' num2str(nchan) ' channels']);
for trlop=1:ntrial
ft_progress(trlop/ntrial, 'searching in trial %d from %d\n', trlop, ntrial);
if hasdata
dat = ft_fetch_data(data, 'header', hdr, 'begsample', trl(trlop,1), 'endsample', trl(trlop,2), 'chanindx', chanindx, 'checkboundary', strcmp(cfg.continuous, 'no'));
else
dat = ft_read_data(cfg.datafile, 'header', hdr, 'begsample', trl(trlop,1), 'endsample', trl(trlop,2), 'chanindx', chanindx, 'checkboundary', strcmp(cfg.continuous, 'no'), 'dataformat', cfg.dataformat);
end
% determine the length of the data in this trial
nsample = trl(trlop,2)-trl(trlop,1)+1;
if size(trl,2)>2
time = offset2time(trl(trlop,3), hdr.Fs, nsample);
else
time = offset2time(0, hdr.Fs, nsample);
end
% only do the preprocessing and filtering if there is an option that suggests to have an effect
status = struct2cell(artfctdef);
status = status(cellfun(@(x) ischar(x), status));
if any(ismember(status, {'yes', 'abs', 'complex', 'real', 'imag', 'absreal', 'absimag', 'angle'}))
ft_notice('FieldTrip:applyingPreprocessing', 'applying preprocessing options')
dat = preproc(dat, label, time, artfctdef);
end
for sgnlop=1:nchan
% make a vector that indicates for each sample whether it exceeds the threshold
artval = false(1, nsample);
artval = artval | any(dat(sgnlop,:)<=artfctdef.min,1);
artval = artval | any(dat(sgnlop,:)>=artfctdef.max,1);
% compute the range as the maximum of the peak-to-peak values for each channel
ptpval = max(dat(sgnlop,:)) - min(dat(sgnlop,:));
if any(ptpval>=artfctdef.range)
artval(:) = true; % mark the whole segment as bad
end
% this is when a different onset and offset are specified
switch direction
case 'up'
onset = find(diff([0 dat(sgnlop,:)>=artfctdef.onset])>0); % find all rising flanks
offset = nan(size(onset));
for i=1:numel(onset)
rem = dat(sgnlop,onset(i)+1:end); % this is the remaining data following the artifact onset
rem = (rem<=artfctdef.offset); % threshold for the offset
if any(rem)
offset(i) = find(rem, 1, 'first'); % find the falling flank
else
offset(i) = length(rem); % take the last sample
end
offset(i) = offset(i) + onset(i);
% add it to the other artifacts in the boolean vector
artval(onset(i):offset(i)) = true;
end
case 'down'
onset = find(diff([0 dat(sgnlop,:)<=artfctdef.onset])>0); % find all rising flanks
offset = nan(size(onset));
for i=1:numel(onset)
rem = dat(sgnlop,onset(i)+1:end); % this is the remaining data following the artifact onset
rem = (rem>=artfctdef.offset);
if any(rem)
offset(i) = find(rem, 1, 'first'); % find the falling flank
else
offset(i) = length(rem); % take the last sample
end
offset(i) = offset(i) + onset(i);
% add it to the other artifacts in the boolean vector
artval(onset(i):offset(i)) = true;
end
case 'none'
% nothing to do
end
% to avoid confusion with the offset that is used further down
clear onset offset
begsample = find(diff([0 artval])>0)';
endsample = find(diff([artval 0])<0)';
offset = nan(size(begsample)); % the offset of the peak relative to the segment, just like in FT_DEFINETRIAL
channel = repmat(label(sgnlop), size(begsample));
% determine the sample at which the signal peaks
for i=1:numel(begsample)
seg = dat(sgnlop,begsample(i):endsample(i)); % get the segment of data
if all(seg>=artfctdef.max) || strcmp(direction, 'up')
[dum, indx] = max(seg);
offset(i) = 1 - indx; % relative to the start of the segment, 0 is the first sample, -1 is the 2nd, etc.
elseif all(seg<=artfctdef.min) || strcmp(direction, 'down')
[dum, indx] = min(seg);
offset(i) = 1 - indx; % relative to the start of the segment, 0 is the first sample, -1 is the 2nd, etc.
end % if up or down
end % for each artifact in this trial
% express them relative to the start of the data, not the start of the trial
begsample = begsample + trl(trlop,1) - 1;
endsample = endsample + trl(trlop,1) - 1;
% remember the parts where this channel exceeds the threshold as artifacts
if ~isempty(begsample)
artifact = vertcat(artifact, table(begsample, endsample, offset, channel));
end
end % for sgnlop
end % for trlop
ft_progress('close');
% reset it to the default
ft_notice on FieldTrip:applyingPreprocessing
if strcmp(cfg.representation, 'numeric') && istable(artifact)
if isempty(artifact)
% an empty table does not have columns
artifact = zeros(0,3);
else
% convert the table to a numeric array with the columns begsample, endsample and offset
artifact = table2array(artifact(:,1:3));
end
elseif strcmp(cfg.representation, 'table') && isnumeric(artifact)
if isempty(artifact)
% an empty table does not have columns
artifact = table();
else
% convert the numeric array to a table with the columns begsample, endsample and offset
begsample = artifact(:,1);
endsample = artifact(:,2);
offset = artifact(:,3);
artifact = table(begsample, endsample, offset);
end
end
% remember the details that were used here and store the detected artifacts
cfg.artfctdef.threshold = artfctdef;
cfg.artfctdef.threshold.artifact = artifact;
ft_notice('detected %d artifacts\n', size(artifact,1));
% do the general cleanup and bookkeeping at the end of the function
ft_postamble provenance
ft_postamble previous data
ft_postamble savevar