-
Notifications
You must be signed in to change notification settings - Fork 18
/
clean_channels_nolocs.m
145 lines (136 loc) · 7.25 KB
/
clean_channels_nolocs.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
function [signal,removed_channels] = clean_channels_nolocs(signal,min_corr,ignored_quantile,window_len,max_broken_time,linenoise_aware)
% Remove channels with abnormal data from a continuous data set.
% Signal = clean_channels_nolocs(Signal,MinCorrelation,IgnoredQuantile,WindowLength,MaxBrokenTime,LineNoiseAware)
%
% This is an automated artifact rejection function which ensures that the data contains no channels
% that record only noise for extended periods of time. If channels with control signals are
% contained in the data these are usually also removed. The criterion is based on correlation: if a
% channel is decorrelated from all others (pairwise correlation < a given threshold), excluding a
% given fraction of most correlated channels -- and if this holds on for a sufficiently long fraction
% of the data set -- then the channel is removed.
%
% In:
% Signal : Continuous data set, assumed to be appropriately high-passed (e.g. >0.5Hz or
% with a 0.5Hz - 2.0Hz transition band).
%
% MinCorrelation : Minimum correlation between a channel and any other channel (in a short period
% of time) below which the channel is considered abnormal for that time period.
% Reasonable range: 0.4 (very lax) to 0.6 (quite aggressive). The default is 0.45.
%
%
% The following are detail parameters that usually do not have to be tuned. If you cannot get
% the function to do what you want, you might consider adapting these to your data.
%
% IgnoredQuantile : Fraction of channels that need to have at least the given MinCorrelation value
% w.r.t. the channel under consideration. This allows to deal with channels or
% small groups of channels that measure the same noise source, e.g. if they are
% shorted. If many channels can be disconnected during an experiment and you
% have strong noise in the room, you might increase this fraction, but consider
% that this a) requires you to decrease the MinCorrelation appropriately and b)
% this can make the correlation measure more brittle. Reasonable range: 0.05 (rather
% lax) to 0.2 (very tolerant re disconnected/shorted channels).The default is
% 0.1.
%
% WindowLength : Length of the windows (in seconds) for which correlation is computed; ideally
% short enough to reasonably capture periods of global artifacts (which are
% ignored), but not shorter (for statistical reasons). Default: 2.
%
% MaxBrokenTime : Maximum time (either in seconds or as fraction of the recording) during which a
% retained channel may be broken. Reasonable range: 0.1 (very aggressive) to 0.6
% (very lax). The default is 0.5.
%
% LineNoiseAware : Whether the operation should be performed in a line-noise aware manner. If enabled,
% the correlation measure will not be affected by the presence or absence of line
% noise (using a temporary notch filter). Default: true.
%
% Out:
% Signal : data set with bad channels removed
%
% Notes:
% This function requires the Signal Processing toolbox.
%
% Christian Kothe, Swartz Center for Computational Neuroscience, UCSD
% 2010-07-06
% Copyright (C) Christian Kothe, SCCN, 2010, christian@sccn.ucsd.edu
%
% This program 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 2 of the
% License, or (at your option) any later version.
%
% This program 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 this program; if not,
% write to the Free Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
% USA
if ~exist('min_corr','var') || isempty(min_corr) min_corr = 0.45; end
if ~exist('ignored_quantile','var') || isempty(ignored_quantile) ignored_quantile = 0.1; end
if ~exist('window_len','var') || isempty(window_len) window_len = 2; end
if ~exist('max_broken_time','var') || isempty(max_broken_time) max_broken_time = 0.5; end
if ~exist('linenoise_aware','var') || isempty(linenoise_aware) linenoise_aware = true; end
% flag channels
if max_broken_time > 0 && max_broken_time < 1 %#ok<*NODEF>
max_broken_time = size(signal.data,2)*max_broken_time;
else
max_broken_time = signal.srate*max_broken_time;
end
signal.data = double(signal.data);
[C,S] = size(signal.data);
window_len = window_len*signal.srate;
wnd = 0:window_len-1;
offsets = 1:window_len:S-window_len;
W = length(offsets);
retained = 1:(C-ceil(C*ignored_quantile));
% optionally ignore both 50 and 60 Hz spectral components...
if linenoise_aware
Bwnd = design_kaiser(2*45/signal.srate,2*50/signal.srate,60,true);
if signal.srate <= 110
error('Sampling rate must be above 110 Hz');
elseif signal.srate <= 130
B = design_fir(length(Bwnd)-1,[2*[0 45 50 55]/signal.srate 1],[1 1 0 1 1],[],Bwnd);
else
B = design_fir(length(Bwnd)-1,[2*[0 45 50 55 60 65]/signal.srate 1],[1 1 0 1 0 1 1],[],Bwnd);
end
for c=signal.nbchan:-1:1
X(:,c) = filtfilt_fast(B,1,signal.data(c,:)'); end
else
X = signal.data';
end
% for each window, flag channels with too low correlation to any other channel (outside the
% ignored quantile)
flagged = zeros(C,W);
for o=1:W
sortcc = sort(abs(corrcoef(X(offsets(o)+wnd,:))));
flagged(:,o) = all(sortcc(retained,:) < min_corr);
end
% mark all channels for removal which have more flagged samples than the maximum number of
% ignored samples
removed_channels = sum(flagged,2)*window_len > max_broken_time;
% apply removal
if all(removed_channels)
disp('Warning: all channels are flagged bad according to the used criterion: not removing anything.');
elseif any(removed_channels)
disp('Now removing bad channels...');
try
signal = pop_select(signal,'nochannel',find(removed_channels));
catch e
if ~exist('pop_select','file')
disp('Apparently you do not have EEGLAB''s pop_select() on the path.');
else
disp('Could not select channels using EEGLAB''s pop_select(); details: ');
hlp_handleerror(e,1);
end
disp('Falling back to a basic substitute and dropping signal meta-data.');
if length(signal.chanlocs) == size(signal.data,1)
signal.chanlocs = signal.chanlocs(~removed_channels); end
signal.data = signal.data(~removed_channels,:);
signal.nbchan = size(signal.data,1);
[signal.icawinv,signal.icasphere,signal.icaweights,signal.icaact,signal.stats,signal.specdata,signal.specicaact] = deal([]);
end
if isfield(signal.etc,'clean_channel_mask') && sum(signal.etc.clean_channel_mask) == length(removed_channels)
signal.etc.clean_channel_mask(signal.etc.clean_channel_mask) = ~removed_channels;
else
signal.etc.clean_channel_mask = ~removed_channels;
end
end