/
ft_rejectcomponent.m
252 lines (215 loc) · 8.97 KB
/
ft_rejectcomponent.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
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
function [data] = ft_rejectcomponent(cfg, comp, data)
% FT_REJECTCOMPONENT backprojects an ICA (or similar) decomposition to the
% channel level after removing the independent components that contain
% the artifacts. This function does not automatically detect the artifact
% components, you will have to do that yourself.
%
% Use as
% [data] = ft_rejectcomponent(cfg, comp)
% or as
% [data] = ft_rejectcomponent(cfg, comp, data)
%
% where the input comp is the result of FT_COMPONENTANALYSIS. The output
% data will have the same format as the output of FT_PREPROCESSING.
%
% An optional input argument data can be provided. In that case
% componentanalysis will do a subspace projection of the input data
% onto the space which is spanned by the topographies in the unmixing
% matrix in comp, after removal of the artifact components. Please use
% this option of including data as input, if you wish to use the output
% data.grad in further computation, for example for leadfield computation.
%
% The configuration structure can contain
% cfg.component = list of components to remove, e.g. [1 4 7] or see FT_CHANNELSELECTION
% cfg.demean = 'no' or 'yes', whether to demean the input data (default = 'yes')
% cfg.updatesens = 'no' or 'yes' (default = 'yes')
%
% To facilitate data-handling and distributed computing you can use
% cfg.inputfile = ...
% cfg.outputfile = ...
% 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_COMPONENTANALYSIS, FT_PREPROCESSING
% Copyright (C) 2005-2014, Robert Oostenveld
%
% 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 debug
ft_preamble loadvar comp data
ft_preamble provenance comp data
% the ft_abort variable is set to true or false in ft_preamble_init
if ft_abort
return
end
% set the defaults
cfg.component = ft_getopt(cfg, 'component', []);
cfg.demean = ft_getopt(cfg, 'demean', 'yes');
cfg.feedback = ft_getopt(cfg, 'feedback', 'text');
cfg.updatesens = ft_getopt(cfg, 'updatesens', 'yes');
% the data can be passed as input arguments or can be read from disk
hascomp = exist('comp', 'var');
hasdata = exist('data', 'var');
if hascomp && hasdata
% check if the input data is valid for this function
istlck = ft_datatype(data, 'timelock'); % this will be temporary converted into raw
data = ft_checkdata(data, 'datatype', 'raw');
comp = ft_checkdata(comp, 'datatype', 'comp');
label = data.label;
nchans = length(data.label);
ncomps = length(comp.label);
elseif hascomp
% check if the input data is valid for this function
istlck = ft_datatype(comp, 'timelock'); % this will be temporary converted into raw
comp = ft_checkdata(comp, 'datatype', 'raw+comp');
label = comp.topolabel;
ncomps = length(comp.label);
else
ft_error('incorrect number of input arguments');
end
% cfg.component can be indicated by number or by label
cfg.component = ft_channelselection(cfg.component, comp.label);
reject = match_str(comp.label, cfg.component);
if isempty(reject)
ft_warning('no components were selected for rejection');
end
if min(reject)<1
ft_error('you cannot remove components that are not present in the data');
end
if max(reject)>ncomps
ft_error('you cannot remove components that are not present in the data');
end
if hasdata && strcmp(cfg.demean, 'yes')
% optionally perform baseline correction on each trial
ft_info('baseline correcting data \n');
for trial=1:numel(data.trial)
data.trial{trial} = ft_preproc_baselinecorrect(data.trial{trial});
end
end
% set the rejected component amplitudes to zero
ft_info('removing %d components\n', length(reject));
ft_info('keeping %d components\n', ncomps-length(reject));
% create a projection matrix by subtracting the subspace spanned by the
% topographies of the to-be-removed components from identity
[seldat, selcomp] = match_str(label, comp.topolabel);
if hasdata && length(seldat)~=length(label)
ft_warning('the subspace projection is not guaranteed to be correct for non-orthogonal components');
end
if hasdata
mixing = comp.topo(selcomp,:);
unmixing = comp.unmixing(:,selcomp);
% I am not sure about this, but it gives comparable results to the ~hasdata case
% when comp contains non-orthogonal (=ica) topographies, and contains a complete decomposition
montage = [];
montage.tra = eye(length(selcomp)) - mixing(:, reject)*unmixing(reject, :);
% we are going from data to components, and back again
montage.labelold = comp.topolabel(selcomp);
montage.labelnew = comp.topolabel(selcomp);
keepunused = 'yes'; % keep the original data which are not present in the mixing provided
bname = 'reject';
else
mixing = comp.topo(selcomp, :);
mixing(:, reject) = 0;
montage = [];
montage.tra = mixing;
% we are going from components to data
montage.labelold = comp.label;
montage.labelnew = comp.topolabel(selcomp);
keepunused = 'no'; % don't need to keep the original rejected components
bname = 'invcomp';
% create the initial data structure, remove all component details
data = keepfields(comp, {'trial', 'time', 'label', 'fsample', 'grad', 'elec', 'opto', 'trialinfo', 'sampleinfo'});
end % if hasdata
% apply the linear projection to the data
data = ft_apply_montage(data, montage, 'keepunused', keepunused, 'feedback', cfg.feedback, 'showcallinfo', cfg.showcallinfo);
sensfield = cell(0,1);
if isfield(data, 'grad')
sensfield{end+1} = 'grad';
end
if isfield(data, 'elec')
sensfield{end+1} = 'elec';
end
if isfield(data, 'opto')
sensfield{end+1} = 'opto';
end
% apply the linear projection also to the sensor description
if ~isempty(sensfield)
if strcmp(cfg.updatesens, 'yes')
for m = 1:numel(sensfield)
ft_info('also applying the backprojection matrix to the %s structure\n', sensfield{m});
% the balance field is needed to keep the sequence of linear projections
if ~isfield(data.(sensfield{m}), 'balance')
data.(sensfield{m}).balance.current = 'none';
end
% keepunused = 'yes' is required to get back e.g. reference or otherwise
% unused sensors in the sensor description. the unused components need to
% be removed in a second step
sens = ft_apply_montage(data.(sensfield{m}), montage, 'keepunused', 'yes', 'balancename', bname, 'feedback', cfg.feedback);
% remove the unused channels from the grad/elec/opto
[junk, remove] = match_str(comp.label, sens.label);
sens.tra(remove,:) = [];
sens.label(remove) = [];
sens.chanpos(remove,:) = [];
if isfield(sens, 'chanori')
sens.chanori(remove,:) = [];
end
% there could have been sequential subspace projections, so the
% invcomp-field may have been renamed into invcompX. If this it the case,
% take the one with the highest suffix
invcompfield = bname;
if ~isfield(sens.balance, invcompfield)
for k = 10:-1:1
if isfield(sens.balance, [bname num2str(k)])
invcompfield = [invcompfield num2str(k)];
break;
end
end
end
% remove the unused components from the balancing
[junk, remove] = match_str(comp.label, sens.balance.(invcompfield).labelnew);
sens.balance.(invcompfield).tra(remove, :) = [];
sens.balance.(invcompfield).labelnew(remove) = [];
data.(sensfield{m}) = sens;
end
else
for m = 1:numel(sensfield)
ft_info('not applying the backprojection matrix to the %s structure\n', sensfield{m});
% simply copy it over
comp.(sensfield{m}) = data.(sensfield{m});
end
end
end % if sensfield
if istlck
% convert the raw structure back into a timelock structure
data = ft_checkdata(data, 'datatype', 'timelock');
end
% do the general cleanup and bookkeeping at the end of the function
ft_postamble debug
ft_postamble previous comp data
ft_postamble provenance data
ft_postamble history data
ft_postamble savevar data