-
Notifications
You must be signed in to change notification settings - Fork 2
/
Copy pathDefineArtificialICAs_AdjustAndSASICA.m
143 lines (117 loc) · 6.64 KB
/
DefineArtificialICAs_AdjustAndSASICA.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
function [icarejcomps EEG_ICA_edited] = DefineArtificialICAs_AdjustAndSASICA(cfg,EEG)
% This program aims to isolate artificial ICAs, e.g., VEOG, HEOG, and focal component, from the EEG data.
% However, "No automated method can accurately isolate artifacts without supervision (Chaumon et al., 2015),"
% which I totally agree; Therefore, if you don't trust the following automatic pipeline please use SASICA or
% Adjust or other toolboxes to visually inspect all the ICAs, e.g., "PlotAndCheckICAsWithAdjust.m";
if ~isfield(cfg, 'focalcomp'), cfg.focalcomp.enable = 1; end
if ~isfield(cfg, 'EOGcorr'), cfg.EOGcorr.enable = 1; end
if ~isfield(cfg, 'autocorr'), cfg.autocorr.enable = 1; end
if ~isfield(cfg, 'opts'), cfg.opts.noplot = 1; end
if ~isfield(cfg, 'plotarg'), cfg.plotarg = 0; end; plotarg = cfg.plotarg;
%% Load the EEG data with ICA weights;
%EEG1 = pop_eegchanoperator( EEG, { 'ch125=(ch21+ch14+ch15-ch12-ch5-ch6)/3 label VEOG' 'ch126=(ch25+ch32+ch26-ch8-ch2-ch1)/3 label HEOG' });
%% run the ICAs detection with SASICA and Adjust with the program modified by WX;
% EOG correlation test;
cfg.EOGcorr.corthreshH = 'auto 2.5'; % M +/- 2.5*SD (99%); M +/- 2*SD (95%);
cfg.EOGcorr.Veogchannames = 'VEOG';
cfg.EOGcorr.corthreshV = 'auto 2.5';
cfg.EOGcorr.Heogchannames = 'HEOG';
% Focal component test;
cfg.focalcomp.focalICAout = 'auto 2.5';
% Muscle movement (auto correlation);
cfg.autocorr.dropautocorr = 'auto 2.5';
% ADJUST
cfg.ADJUST.enable = 1;
%[EEG cfg] = WX_eeg_SASICA(EEG,cfg);
[EEG1 output] = WX_eeg_SASICA(EEG,cfg);
%% Remove the bad components;
% Results from SASICA
%focalcomps = find(EEG.reject.SASICA.icarejfocalcomp);
%eogcomps = find(EEG.reject.SASICA.icarejchancorr);
focalcomps = find(EEG1.reject.SASICA.icarejfocalcomp);
eogcomps = find(EEG1.reject.SASICA.icarejchancorr);
if cfg.autocorr.enable == 0;
mclcomps = [];
else
mclcomps = find(EEG1.reject.SASICA.icarejautocorr);
gdsf_adjust = find(EEG1.reject.SASICA.icaADJUST.GDSF > EEG1.reject.SASICA.icaADJUST.soglia_GDSF);
mclcomps = intersect(mclcomps,gdsf_adjust);
end
% all bad components in SASICA
rejcomps_sasica = union(focalcomps, eogcomps);
% V Eye movements from ADJSUT:
% 1. Maximum Epoch Variance (MEV)
% 2. Spatial Average Difference (SAD)
mevcomps = find(EEG1.reject.SASICA.icaADJUST.nuovaV > EEG1.reject.SASICA.icaADJUST.soglia_V);
sadcomps = find(EEG1.reject.SASICA.icaADJUST.SAD > EEG1.reject.SASICA.icaADJUST.soglia_SAD);
vecomps_adjust = intersect(mevcomps,sadcomps);
% H Eye movements from ADJSUT:
% 1. MEV 2.Spatial Eye Difference (SED);
sedcomps = find(EEG1.reject.SASICA.icaADJUST.SED > EEG1.reject.SASICA.icaADJUST.soglia_SED);
hecomps_adjust = intersect(mevcomps,sedcomps);
% Blink
% 1. SAD 2. Spatial Variance Difference (SVD)
blkcomps = EEG1.reject.SASICA.icaADJUST.blink; % This blink detection with infant data in ADJSUT is really shitty;
% over all eye related components by ADJUST;
eyecomps_adjust = unique([vecomps_adjust,hecomps_adjust,blkcomps]);
% Focal (Discontinuities) components from ADJUST;
gdsfcomps_adjust = EEG1.reject.SASICA.icaADJUST.disc;
% find out the overlapping problematic componencts by sasica and ADJUST;
icarejcomps_eye = intersect(eogcomps,eyecomps_adjust);
icarejcomps_focal = intersect(focalcomps,gdsfcomps_adjust); % The focal comps detection is faily consistent between SASICA and ADJSUT;
icarejcomps = union(icarejcomps_eye,icarejcomps_focal);
% add the mucle movement components to the icarejcomps array;
if cfg.autocorr.enable == 1 % it means auto correlation has been run and we get a 1 x 124 matrix for the results;
icarejcomps = union(icarejcomps,mclcomps);
end
% Both SASICA and ADJUST save their results in the gcompreject field, which is the default place for bad comps in EEGLAB;
EEG1.reject.gcompreject=EEG1.reject.gcompreject.*0; % clear ADJUST's results;
EEG1.reject.gcompreject(icarejcomps) = 1; % Substitute it with the final rejected comps for the dataset;
disp(['Remove ' num2str(length(icarejcomps)) ' components: ' num2str(icarejcomps)]);
EEG_ICA_edited = pop_subcomp( EEG1, icarejcomps, 0);
if plotarg;
% [ALLEEG EEG CURRENTSET] = pop_newset(ALLEEG, EEG, 1,'overwrite','on','gui','on');
% eeglab redraw;
pop_eegplot(EEG1,1,0,0,'','dispchans',10,'spacing',100);
set (gcf,'name','Plot of EEG data before ICA rejection');
EEG_ICA_edited.reject = EEG1.reject;
pop_eegplot(EEG_ICA_edited,1,0,0,'','dispchans',10,'spacing',100,'children',gcf);
set (gcf,'name','Plot of EEG data after ICA rejection');
disp(['These components are removed: ' num2str(icarejcomps)]);
SASICA(EEG);
end
clear EEG1
return
% plot
if plotarg;
% [ALLEEG EEG CURRENTSET] = pop_newset(ALLEEG, EEG, 1,'overwrite','on','gui','on');
% eeglab redraw;
pop_eegplot(EEG1,1,0,0,'','dispchans',10,'spacing',100);
set (gcf,'name','Plot of EEG data before ICA rejection');
EEG_ICA_edited.reject = EEG1.reject;
pop_eegplot(EEG_ICA_edited,1,0,0,'','dispchans',10,'spacing',100,'children',gcf);
set (gcf,'name','Plot of EEG data after ICA rejection');
disp(['These components are removed: ' num2str(icarejcomps)]);
SASICA(EEG);
else % if not for plot then save the data;
% remove channel
EEG_ICA_edited = pop_select( EEG_ICA_edited,'nochannel',{'VEOG' 'HEOG'});
% EEGLAB creates a new reject structure after removing the two EOG channels with pop_select;
%% Save the ICA_edited dataset;
pop_saveset(EEG_ICA_edited, 'filename',EEG1filename,'filepath',SegmentAverageFiles);
% Save the data into Fieldtrip format;
EEG_FT = eeglab2fieldtrip(EEG_ICA_edited,'preprocessing','none');
if typeofdata
ftfilename = ['Experiment 2 Subject ' num2str(participantnumber) ' Fieldtrip ' num2str(hz) 'Hz Highpass_ICA_Edited_sERP.mat']; %sERP means single trials ERP;
rejname = ['Experiment 2 Subject ' num2str(participantnumber) ' ' num2str(hz) 'Hz Highpass_ICA_RejInfo_sERP.mat'];
else
ftfilename = ['Experiment 2 Subject ' num2str(participantnumber) ' Age ' num2str(age) ' Fieldtrip ' num2str(hz) 'Hz Highpass_ICA_Edited.mat'];
rejname = ['Experiment 2 Subject ' num2str(participantnumber) ' Age ' num2str(age) ' ' num2str(hz) 'Hz Highpass_ICA_RejInfo.mat'];
end
disp(['Save ' ftfilename]);
save([ftpath ftfilename],'EEG_FT');
% Save the rejection structure;
rejstrct = EEG.reject;
save([ftpath rejname], 'rejstrct');
end
disp('*********************** Finished ***********************');