forked from fieldtrip/fieldtrip
-
Notifications
You must be signed in to change notification settings - Fork 0
/
ft_artifact_threshold.m
205 lines (181 loc) · 8.39 KB
/
ft_artifact_threshold.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
function [cfg, artifact] = ft_artifact_threshold(cfg, data)
% FT_ARTIFACT_THRESHOLD scans for trials in which the range, i.e. the minimum,
% the maximum or the range (min-max difference) of the signal in any
% channel 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
%
% Alternatively you can use it as
% [cfg, artifact] = ft_artifact_threshold(cfg, data)
%
% 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
%
% The following configuration options can be specified
% cfg.artfctdef.threshold.channel = cell-array with channel labels
% cfg.artfctdef.threshold.bpfilter = 'no' or 'yes'
% cfg.artfctdef.threshold.bpfreq = [0.3 30]
% cfg.artfctdef.threshold.bpfiltord = 4
%
% The detection of artifacts is done according to the following settings,
% you should specify at least one of these thresholds
% cfg.artfctdef.threshold.range = value in uV/T, default inf
% cfg.artfctdef.threshold.min = value in uV/T, default -inf
% cfg.artfctdef.threshold.max = value in uV/T, default inf
%
% When cfg.artfctdef.threshold.range is used, the within-channel
% peak-to-peak range is checked against the specified maximum range (so not
% the overall range across channels).
%
% Contrary to the other artifact detection functions, this function
% will mark the whole trial as an artifact if the threshold is exceeded.
% Furthermore, this function does not support artifact- or filterpadding.
%
% To facilitate data-handling and distributed computing you can use
% cfg.inputfile = ...
% If you specify one of these (or both) the input data will be read from a *.mat
% file on disk and/or the output data will be written to a *.mat file. These mat
% files should contain only a single variable, corresponding with the
% input/output 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 default rejection parameters for clip artifacts if necessary
if ~isfield(cfg, 'artfctdef'), cfg.artfctdef = []; end
if ~isfield(cfg.artfctdef,'threshold'), cfg.artfctdef.threshold = []; end
if ~isfield(cfg, 'headerformat'), cfg.headerformat = []; end
if ~isfield(cfg, 'dataformat'), cfg.dataformat = []; end
% copy the specific configuration for this function out of the master cfg
artfctdef = cfg.artfctdef.threshold;
% rename some cfg fields for backward compatibility
if isfield(artfctdef, 'sgn') && ~isfield(artfctdef, 'channel')
artfctdef.channel = artfctdef.sgn;
artfctdef = rmfield(artfctdef, 'sgn');
end
if isfield(artfctdef, 'cutoff') && ~isfield(artfctdef, 'range')
artfctdef.range = artfctdef.cutoff;
artfctdef = rmfield(artfctdef, 'cutoff');
end
% set default preprocessing parameters if necessary
if ~isfield(artfctdef, 'channel'), artfctdef.channel = 'all'; end
if ~isfield(artfctdef, 'bpfilter'), artfctdef.bpfilter = 'yes'; end
if ~isfield(artfctdef, 'bpfreq'), artfctdef.bpfreq = [0.3 30]; end
if ~isfield(artfctdef, 'bpfiltord'), artfctdef.bpfiltord = 4; end
% set the default artifact detection parameters
if ~isfield(artfctdef, 'range'), artfctdef.range = inf; end
if ~isfield(artfctdef, 'min'), artfctdef.min = -inf; end
if ~isfield(artfctdef, 'max'), artfctdef.max = inf; end
% 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 given as input
data = ft_checkdata(data, 'hassampleinfo', 'yes');
cfg = ft_checkconfig(cfg, 'forbidden', {'dataset', 'headerfile', 'datafile'});
hdr = ft_fetch_header(data);
end
% set default cfg.continuous
if ~isfield(cfg, 'continuous')
if hdr.nTrials==1
cfg.continuous = 'yes';
else
cfg.continuous = 'no';
end
end
if ~isfield(cfg, 'trl')
% get it from the data itself
cfg.trl = data.trialinfo;
cfg.trl(:,3) = 0;
end
% get the remaining settings
numtrl = size(cfg.trl,1);
channel = ft_channelselection(artfctdef.channel, hdr.label);
channelindx = match_str(hdr.label,channel);
artifact = [];
for trlop = 1:numtrl
if hasdata
dat = ft_fetch_data(data, 'header', hdr, 'begsample', cfg.trl(trlop,1), 'endsample', cfg.trl(trlop,2), 'chanindx', channelindx, 'checkboundary', strcmp(cfg.continuous, 'no'));
else
dat = ft_read_data(cfg.datafile, 'header', hdr, 'begsample', cfg.trl(trlop,1), 'endsample', cfg.trl(trlop,2), 'chanindx', channelindx, 'checkboundary', strcmp(cfg.continuous, 'no'), 'dataformat', cfg.dataformat);
end
dat = preproc(dat, channel, offset2time(cfg.trl(trlop,3), hdr.Fs, size(dat,2)), artfctdef);
% compute the min, max and range over all channels and samples
minval = min(dat(:));
maxval = max(dat(:));
% compute the range as the maximum of the peak-to-peak values for each
% channel
ptpval = max(dat, [], 2) - min(dat, [], 2);
% track for bad trials for each channel
badChnInd = find(ptpval > artfctdef.range);
% determine range and index of 'worst' channel
worstChanRange = max(ptpval);
worstChanInd = find(worstChanRange == ptpval);
% test the min, max and range against the specified thresholds
if ~isempty(artfctdef.min) && minval<artfctdef.min
fprintf('threshold artifact scanning: trial %d from %d exceeds min-threshold\n', trlop, numtrl);
artifact(end+1,1:2) = cfg.trl(trlop,1:2);
elseif ~isempty(artfctdef.max) && maxval>artfctdef.max
fprintf('threshold artifact scanning: trial %d from %d exceeds max-threshold\n', trlop, numtrl);
artifact(end+1,1:2) = cfg.trl(trlop,1:2);
elseif ~isempty(artfctdef.range) && worstChanRange>artfctdef.range
fprintf('threshold artifact scanning: trial %d from %d exceeds range-threshold; max-range channel = %s\n', trlop, numtrl, hdr.label{channelindx(worstChanInd)});
artifact(end+1,1:2) = cfg.trl(trlop,1:2);
else
fprintf('threshold artifact scanning: trial %d from %d is ok\n', trlop, numtrl);
end
end
fprintf('detected %d artifacts\n', size(artifact,1));
% remember the details that were used here
cfg.artfctdef.threshold = artfctdef;
cfg.artfctdef.threshold.trl = cfg.trl; % trialdefinition prior to rejection
cfg.artfctdef.threshold.channel = channel; % exact channels used for detection
cfg.artfctdef.threshold.artifact = artifact; % detected artifacts
% do the general cleanup and bookkeeping at the end of the function
ft_postamble provenance
ft_postamble previous data