forked from fieldtrip/fieldtrip
-
Notifications
You must be signed in to change notification settings - Fork 0
/
ft_denoise_prewhiten.m
187 lines (164 loc) · 7 KB
/
ft_denoise_prewhiten.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
function [dataout] = ft_denoise_prewhiten(cfg, datain, noise)
% FT_DENOISE_PREWHITEN applies a spatial prewhitening operation to the data using the
% inverse noise covariance matrix. The consequence is that all channels are expressed
% in singnal-to-noise units, causing different channel types to be comparable. This
% ensures equal weighting in source estimation on data with different channel types.
%
% Use as
% dataout = ft_denoise_prewhiten(cfg, datain, noise)
% where the datain is the original data from FT_PREPROCESSING and
% noise should contain the estimated noise covariance from
% FT_TIMELOCKANALYSIS.
%
% The configuration structure can contain
% cfg.channel = cell-array, see FT_CHANNELSELECTION (default = 'all')
% cfg.split = cell-array of channel types between which covariance is split, it can also be 'all' or 'no'
% cfg.lambda = scalar, or string, regularization parameter for the inverse
% cfg.kappa = scalar, truncation parameter for the inverse
%
% The channel selection relates to the channels that are pre-whitened using the same
% selection of channels in the noise covariance. All channels present in the input
% data structure will be present in the output, including trigger and other auxiliary
% channels.
%
% See also FT_DENOISE_SYNTHETIC, FT_DENOISE_PCA, FT_DENOISE_DSSP, FT_DENOISE_TSP
% Copyright (C) 2018-2019, Robert Oostenveld and Jan-Mathijs Schoffelen
%
% 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;
ft_defaults
ft_preamble init
ft_preamble debug
ft_preamble loadvar datain
ft_preamble provenance datain
ft_preamble trackconfig
% the ft_abort variable is set to true or false in ft_preamble_init
if ft_abort
% do not continue function execution in case the outputfile is present and the user indicated to keep it
return
end
% get the defaults
cfg.channel = ft_getopt(cfg, 'channel', 'all');
cfg.split = ft_getopt(cfg, 'split', 'all');
cfg.lambda = ft_getopt(cfg, 'lambda', 0);
cfg.kappa = ft_getopt(cfg, 'kappa', []);
cfg.tol = ft_getopt(cfg, 'tol', []);
cfg.realflag = ft_getopt(cfg, 'realflag', true); % for complex-valued crsspctrm
cfg.invmethod = ft_getopt(cfg, 'invmethod', 'tikhonov');
% ensure that the input data is correct, the next line is needed for a
% attempt correct detection of the data chanunit (with a hdr-field it fails
% for meggrad data)
if isfield(datain, 'hdr'), datain = rmfield(datain, 'hdr'); end
datain = ft_checkdata(datain, 'datatype', {'raw' 'timelock' 'freq'}, 'haschantype', 'yes', 'haschanunit', 'yes');
noise = ft_checkdata(noise, 'datatype', { 'timelock' 'freq'}, 'haschantype', 'yes', 'haschanunit', 'yes');
dtype_datain = ft_datatype(datain);
% check for allowed input combinations
switch dtype_datain
case 'raw'
assert(ft_datatype(noise, 'timelock'), 'noise data should be of datatype ''timelock''');
case 'timelock'
assert(ft_datatype(noise, 'timelock'), 'noise data should be of datatype ''timelock''');
case 'freq'
if ft_datatype(noise, 'freq')
% this is only allowed if both structures have the same singleton
% frequency
assert(numel(noise.freq==1) && numel(datain.freq==1) && isequal(noise.freq,datain.freq), ...
'with both datain and noise of datatype ''freq'', only singleton and equal frequency bins are allowed');
elseif ft_datatype(noise, 'timelock')
% this is OK
end
end
% select channels and trials of interest, by default this will select all channels and trials
tmpcfg = keepfields(cfg, {'trials', 'channel', 'showcallinfo'});
datain = ft_selectdata(tmpcfg, datain);
noise = ft_selectdata(tmpcfg, noise);
% restore the provenance information
[cfg, datain] = rollback_provenance(cfg, datain);
[cfg, noise] = rollback_provenance(cfg, noise);
if ft_datatype(noise, 'timelock')
if ~isfield(noise, 'cov')
ft_error('noise covariance is not present');
else
noisecov = noise.cov;
end
elseif ft_datatype(noise, 'freq')
if ~isfield(noise, 'crsspctrm')
ft_error('noise cross-spectrum is not present');
else
if istrue(cfg.realflag)
noisecov = real(noise.crsspctrm);
else
noisecov = noise.crsspctrm;
end
end
end
% determine whether it is EEG and/or MEG data
hasgrad = isfield(datain, 'grad');
haselec = isfield(datain, 'elec');
hasopto = isfield(datain, 'opto');
if isequal(cfg.split, 'no')
chantype = {};
elseif isequal(cfg.split, 'all')
chantype = unique(noise.chantype);
else
chantype = cfg.split;
end
% zero out the off-diagonal elements for the specified channel types
for i=1:numel(chantype)
sel = strcmp(noise.chantype, chantype{i});
noisecov(sel,~sel) = 0;
noisecov(~sel,sel) = 0;
end
% invert the noise covariance matrix
invnoise = ft_inv(noisecov, 'lambda', cfg.lambda, 'kappa', cfg.kappa, 'tolerance', cfg.tol, 'method', cfg.invmethod);
[U,S,V] = svd(invnoise,'econ');
% the prewhitening projection first rotates to orthogonal channels,
% then scales, and then rotates the channels back to (more or less)
% their original MEG-channel representation
prewhiten = [];
prewhiten.tra = U*sqrt(S)*U';
prewhiten.labelold = noise.label;
prewhiten.labelnew = noise.label;
prewhiten.chantypeold = noise.chantype;
prewhiten.chantypenew = noise.chantype;
prewhiten.chanunitold = noise.chanunit;
prewhiten.chanunitnew = repmat({'snr'}, size(noise.chantype));
% apply the projection to the data
dataout = ft_apply_montage(removefields(datain, {'grad', 'elec', 'opto'}), prewhiten, 'keepunused', 'yes');
if hasgrad
% the gradiometer structure needs to be updated to ensure that the forward model remains consistent with the data
dataout.grad = ft_apply_montage(datain.grad, prewhiten, 'balancename', 'prewhiten');
end
if haselec
% the electrode structure needs to be updated to ensure that the forward model remains consistent
dataout.elec = ft_apply_montage(datain.elec, prewhiten, 'balancename', 'prewhiten');
end
if hasopto
% the electrode structure needs to be updated to ensure that the forward model remains consistent
dataout.opto = ft_apply_montage(datain.opto, prewhiten, 'balancename', 'prewhiten');
end
ft_postamble debug
ft_postamble trackconfig
ft_postamble previous datain
ft_postamble provenance dataout
ft_postamble history dataout
ft_postamble savevar dataout