-
Notifications
You must be signed in to change notification settings - Fork 10
/
Copy pathRELAX_ICA_subtract.m
195 lines (179 loc) · 11.5 KB
/
RELAX_ICA_subtract.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
%% RELAX EEG CLEANING PIPELINE, Copyright (C) (2022) Neil Bailey
% 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 3 of the License, or
% 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, see https://www.gnu.org/licenses/.
%% RELAX_ICA_subtract:
function [EEG] = RELAX_ICA_subtract(EEG,RELAX_cfg)
fastica_symm_Didnt_Converge=[0 0 0];
% run ICA:
if strcmp(RELAX_cfg.ICA_method,'runica')
[OUTEEG, ~] = pop_runica_nwb(EEG, 'extended',1,'interupt','on'); %runica for parametric, default extended for finding subgaussian distributions
W = OUTEEG.icaweights*OUTEEG.icasphere;
A = inv(W);
OUTEEG = eeg_checkset(OUTEEG, 'ica');
if isempty(OUTEEG.icaact)==1
OUTEEG.icaact = (OUTEEG.icaweights*OUTEEG.icasphere)*OUTEEG.data(OUTEEG.icachansind,:);
OUTEEG.icaact = reshape( OUTEEG.icaact, size(OUTEEG.icaact,1), OUTEEG.pnts, OUTEEG.trials);
end
IC=reshape(OUTEEG.icaact, size(OUTEEG.icaact,1), []);
elseif strcmp(RELAX_cfg.ICA_method,'cudaica')
[OUTEEG, ~] = pop_runica_nwb(EEG, 'cudaica', 'extended',1); %runica for parametric, default extended for finding subgaussian distributions
W = OUTEEG.icaweights*OUTEEG.icasphere;
A = inv(W);
OUTEEG = eeg_checkset(OUTEEG, 'ica');
if isempty(OUTEEG.icaact)==1
OUTEEG.icaact = (OUTEEG.icaweights*OUTEEG.icasphere)*OUTEEG.data(OUTEEG.icachansind,:);
OUTEEG.icaact = reshape( OUTEEG.icaact, size(OUTEEG.icaact,1), OUTEEG.pnts, OUTEEG.trials);
end
IC=reshape(OUTEEG.icaact, size(OUTEEG.icaact,1), []);
elseif strcmp(RELAX_cfg.ICA_method,'fastica_symm')
% NWB - the following lines repeat fastica_symm up to 3 times in
% the case of non-convergence, then switches to fastica_defl to
% ensure ICA convergence (as cleaning as adversely affected by
% non-convergence issues).
[OUTEEG, ~, NonConvergence] = pop_runica_nwb( EEG, 'icatype', 'fastica','numOfIC', EEG.nbchan, 'approach', 'symm', 'g', 'tanh', 'stabilization', 'on');
fastica_symm_Didnt_Converge(1,1)=NonConvergence;
if NonConvergence==1
[OUTEEG, ~, NonConvergence] = pop_runica_nwb( EEG, 'icatype', 'fastica','numOfIC', EEG.nbchan, 'approach', 'symm', 'g', 'tanh', 'stabilization', 'on');
fastica_symm_Didnt_Converge(1,2)=NonConvergence;
end
if NonConvergence==1
[OUTEEG, ~, NonConvergence] = pop_runica_nwb( EEG, 'icatype', 'fastica','numOfIC', EEG.nbchan, 'approach', 'symm', 'g', 'tanh', 'stabilization', 'on');
fastica_symm_Didnt_Converge(1,3)=NonConvergence;
end
if NonConvergence==1
OUTEEG = pop_runica_nwb( EEG, 'icatype', 'fastica','numOfIC', EEG.nbchan, 'approach', 'defl', 'g', 'tanh', 'stabilization', 'on');
end
W = OUTEEG.icaweights*OUTEEG.icasphere;
A = inv(W);
OUTEEG = eeg_checkset(OUTEEG, 'ica');
if isempty(OUTEEG.icaact)==1
OUTEEG.icaact = (OUTEEG.icaweights*OUTEEG.icasphere)*OUTEEG.data(OUTEEG.icachansind,:);
OUTEEG.icaact = reshape( OUTEEG.icaact, size(OUTEEG.icaact,1), OUTEEG.pnts, OUTEEG.trials);
end
IC=reshape(OUTEEG.icaact, size(OUTEEG.icaact,1), []);
elseif strcmp(RELAX_cfg.ICA_method,'fastica_defl')
OUTEEG = pop_runica_nwb( EEG, 'icatype', 'fastica','numOfIC', EEG.nbchan, 'approach', 'defl', 'g', 'tanh', 'stabilization', 'on');
W = OUTEEG.icaweights*OUTEEG.icasphere;
A = inv(W);
OUTEEG = eeg_checkset(OUTEEG, 'ica');
if isempty(OUTEEG.icaact)==1
OUTEEG.icaact = (OUTEEG.icaweights*OUTEEG.icasphere)*OUTEEG.data(OUTEEG.icachansind,:);
OUTEEG.icaact = reshape( OUTEEG.icaact, size(OUTEEG.icaact,1), OUTEEG.pnts, OUTEEG.trials);
end
IC=reshape(OUTEEG.icaact, size(OUTEEG.icaact,1), []);
elseif strcmp(RELAX_cfg.ICA_method,'amica')
OUTEEG=EEG;
% You'll need to install AMICA first, and in the folder that you
% specify in the line below (with no spaces in any part of the folder or subfolders):
% You can download AMICA via EEGLAB
amica_file = which('runamica15');
if ~exist(amica_file)
disp('!! AMICA directory not found, please ensure you have AMICA installed !!');
end
[filepath,~,~] = fileparts(amica_file);
cd(filepath);
% define parameters
numprocs = 1; % # of nodes (default = 1)
max_threads = 4; % # of threads per node
num_models = 1; % # of models of mixture ICA
max_iter = 2000; % max number of learning steps
mkdir([filepath filesep 'AMICAtmp']);
outdir = [filepath filesep 'AMICAtmp' filesep];
% Run AMICA:
[OUTEEG.icaweights, OUTEEG.icasphere, ~] = runamica15(OUTEEG.data, 'num_chans', EEG.nbchan, 'num_models',num_models,'outdir',outdir,'numprocs', numprocs, 'max_threads', max_threads, 'max_iter',max_iter,'pcakeep', EEG.nbchan, 'do_reject', 1, 'numrej', 15, 'rejsig', 3, 'rejint', 1);
W = OUTEEG.icaweights*OUTEEG.icasphere;
A = inv(W);
OUTEEG = eeg_checkset(OUTEEG, 'ica');
if isempty(OUTEEG.icaact)==1
OUTEEG.icaact = (OUTEEG.icaweights*OUTEEG.icasphere)*OUTEEG.data(OUTEEG.icachansind,:);
OUTEEG.icaact = reshape( OUTEEG.icaact, size(OUTEEG.icaact,1), OUTEEG.pnts, OUTEEG.trials);
end
IC=reshape(OUTEEG.icaact, size(OUTEEG.icaact,1), []);
elseif strcmp(RELAX_cfg.ICA_method,'picard') % RELAX 1.1.3 - NWB added to allow PICARD-O to be run using default settings
[OUTEEG, ~] = pop_runica_nwb(EEG, 'picard', 'mode','ortho','tol',1e-6,'maxiter',500); % run picard
W = OUTEEG.icaweights*OUTEEG.icasphere;
A = inv(W);
OUTEEG = eeg_checkset(OUTEEG, 'ica');
if isempty(OUTEEG.icaact)==1
OUTEEG.icaact = (OUTEEG.icaweights*OUTEEG.icasphere)*OUTEEG.data(OUTEEG.icachansind,:);
OUTEEG.icaact = reshape( OUTEEG.icaact, size(OUTEEG.icaact,1), OUTEEG.pnts, OUTEEG.trials);
end
IC=reshape(OUTEEG.icaact, size(OUTEEG.icaact,1), []);
end
% Use ICLabel to identify artifactual components, so that wICA can be
% performed on them only:
% (https://github.com/sccn/ICLabel)
OUTEEG = iclabel(OUTEEG);
IC_classifications=OUTEEG.etc.ic_classification.ICLabel.classifications; % allows user to set thresholds for classification confidence before considered an artifact
IC_classifications(IC_classifications<RELAX_cfg.ICLabel_thresholds)=0;
[~, I]=max(IC_classifications, [], 2);
if strcmp(RELAX_cfg.Clean_other_comps,'no')==1
ICsMostLikelyNotBrain=(I==2 | I ==3)';
elseif strcmp(RELAX_cfg.Clean_other_comps,'yes')==1
ICsMostLikelyNotBrain=(I>1)';
end
ArtifactICList=find(ICsMostLikelyNotBrain==1);
EEG = pop_subcomp( OUTEEG,ArtifactICList, 0); % removes artifactual components by subtraction
EEG = eeg_checkset( EEG );
EEG.RELAXProcessing_ICA.fastica_symm_Didnt_Converge=fastica_symm_Didnt_Converge; % Tracks whether fastica_symm showed convergence issues (1) or not (0), and how many non-convergences. If 3 non-convergences, then fastica_defl was implemented.
% Check if data might have been too short for effective ICA, using Makoto's rule
% of thumb that ICA requires data length of ((number of channels)^2)*30
% if data were sampled at 250 Hz (assuming that higher sampling
% rates require the same time duration of data as low sampling rates,
% so 1000Hz sampling rates require ((number of channels)^2)*120)
% (https://sccn.ucsd.edu/wiki/Makoto%27s_useful_EEGLAB_code)
ms_per_sample=(1000/EEG.srate);
if ((EEG.nbchan^2)*(120/ms_per_sample))>EEG.pnts
EEG.RELAXProcessing_ICA.DataMaybeTooShortForValidICA='yes';
else
EEG.RELAXProcessing_ICA.DataMaybeTooShortForValidICA='no';
end
if strcmp (EEG.RELAXProcessing_ICA.DataMaybeTooShortForValidICA,'yes')
warning('Data may have been shorter than recommended for effective ICA decomposition')
end
EEG.RELAXProcessing_ICA.Proportion_artifactICs_reduced_by_ICA=mean(ICsMostLikelyNotBrain);
if strcmp(RELAX_cfg.Report_all_ICA_info,'yes')
EEG.RELAXProcessing_ICA.ProportionICs_was_Brain=sum(I==1)/size(OUTEEG.etc.ic_classification.ICLabel.classifications,1);
EEG.RELAXProcessing_ICA.ProportionICs_was_Muscle=sum(I==2)/size(OUTEEG.etc.ic_classification.ICLabel.classifications,1);
EEG.RELAXProcessing_ICA.ProportionICs_was_Eye=sum(I==3)/size(OUTEEG.etc.ic_classification.ICLabel.classifications,1);
EEG.RELAXProcessing_ICA.ProportionICs_was_Heart=sum(I==4)/size(OUTEEG.etc.ic_classification.ICLabel.classifications,1);
EEG.RELAXProcessing_ICA.ProportionICs_was_LineNoise=sum(I==5)/size(OUTEEG.etc.ic_classification.ICLabel.classifications,1);
EEG.RELAXProcessing_ICA.ProportionICs_was_ChannelNoise=sum(I==6)/size(OUTEEG.etc.ic_classification.ICLabel.classifications,1);
EEG.RELAXProcessing_ICA.ProportionICs_was_Other=sum(I==7)/size(OUTEEG.etc.ic_classification.ICLabel.classifications,1);
ICsMostLikelyBrain=(I==1)';
ICsMostLikelyMuscle=(I==2)';
ICsMostLikelyEye=(I==3)';
ICsMostLikelyHeart=(I==4)';
ICsMostLikelyLineNoise=(I==5)';
ICsMostLikelyChannelNoise=(I==6)';
ICsMostLikelyOther=(I==7)';
for x=1:size(OUTEEG.etc.ic_classification.ICLabel.classifications,1)
[~, varianceWav(x)] =compvar(OUTEEG.data, OUTEEG.icaact, OUTEEG.icawinv, x);
end
BrainVariance=sum(abs(varianceWav(ICsMostLikelyBrain)));
ArtifactVariance=sum(abs(varianceWav(~ICsMostLikelyBrain)));
EEG.RELAXProcessing_wICA.ProportionVariance_was_BrainICs=(BrainVariance/(BrainVariance+ArtifactVariance));
MuscleVariance=sum(abs(varianceWav(ICsMostLikelyMuscle)));
EyeVariance=sum(abs(varianceWav(ICsMostLikelyEye)));
HeartVariance=sum(abs(varianceWav(ICsMostLikelyHeart)));
LineNoiseVariance=sum(abs(varianceWav(ICsMostLikelyLineNoise)));
ChannelNoiseVariance=sum(abs(varianceWav(ICsMostLikelyChannelNoise)));
OtherVariance=sum(abs(varianceWav(ICsMostLikelyOther)));
EEG.RELAXProcessing_ICA.ProportionVariance_was_MuscleICs=(MuscleVariance/(BrainVariance+ArtifactVariance));
EEG.RELAXProcessing_ICA.ProportionVariance_was_EyeICs=(EyeVariance/(BrainVariance+ArtifactVariance));
EEG.RELAXProcessing_ICA.ProportionVariance_was_HeartICs=(HeartVariance/(BrainVariance+ArtifactVariance));
EEG.RELAXProcessing_ICA.ProportionVariance_was_LineNoiseICs=(LineNoiseVariance/(BrainVariance+ArtifactVariance));
EEG.RELAXProcessing_ICA.ProportionVariance_was_ChannelNoiseICs=(ChannelNoiseVariance/(BrainVariance+ArtifactVariance));
EEG.RELAXProcessing_ICA.ProportionVariance_was_OtherICs=(OtherVariance/(BrainVariance+ArtifactVariance));
end
end