s01_scalpart.htb 25.2 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
%% SCRIPT DESCRIPTION
% This script features the main preprocessing and cleaning of the data before the ICA process begin. Data is flagged for high standard deveiavtion
% by channel to remove comically bad channels. The data is then rereferenced to the average site, and high/low pass filtered.
% Channel neighbor correlations are then calculated to find odd channels and bridged channels. The outermost channel is also flagged to reduce the rank of the data.
% After another rereference channel correlations are calculated by time epoch and flagged. The file is saved containing the, marks structure, then saved again 
% with all of the mentioned flags puregd from the data in order to prime it for ICA.
%
%% From Config          key_strswap         Description
%-----------------------------------------------------------------------------------------------------------------------------
%    in_path =           [in_path]           Relative path to input data files assuming cd = work_path
11
12
%    montage_info =      [montage_info]      Montage information to fit channel locations on a standard surface.
%                                            This can be a transformations matrix (e.g. [-.14,-22,-48,-.07,0,-1.57,1080,1260,1240])
13
%                                            ... or a a file containing a transformation matrix (e.g. [batch_dfn]_transmat.mat).
14
%                                            ... or a BIDS channel location tsv file (assuming that during BIDS init the channels were coregistered to the standard surface).
James Desjardins's avatar
James Desjardins committed
15
%    staging_script =    [staging_script]    Scipt file path/name that creates marks based on events in the raw data file. 
16
17
% 	 recur_sec =         [recur_sec]         Recurrence (sec) for artifact detection epoching (e.g. .5, Default 1)
%    limit_sec =         [limit_sec]         Limits (sec) for artifact detection epoching (e.g. [-.5 0], Default [0 recur_sec])
18
19
20
21
22
23
24
25
26
27
28
29
30
%    sd_t_meth =         [sd_t_meth]         Method used for flagging epochs (e.g. 'q' (quantiles), or 'na' (default))
%    sd_t_vals  =        [sd_t_vals]         Percentage trim for confidence intervals during epoch standard deviation criteria (e.g. [.3 .7])
%    sd_t_o =            [sd_t_o]            z threshold for flagging epochs during standard deviation criteria (e.g. 6)
%    sd_t_f_meth =       [sd_t_f_meth]       Fixed method used for flagging epochs (should be 'fixed')
%    sd_t_f_vals  =      [sd_t_f_vals]       Percentage trim for confidence intervals during epoch standard deviation criteria (e.g. [.3 .7], leave empty for 'fixed')
%    sd_t_f_o =          [sd_t_f_o]          z threshold for flagging epochs during fixed standard deviation criteria (e.g. .2)
%    sd_t_pad =          [sd_t_pad]          Number of windows to pad onto each side of the ch_sd time flag
%    sd_ch_meth =        [sd_ch_meth]        Method used for flagging channels (e.g. 'q' (quantiles), or 'na' (default))
%    sd_ch_vals  =       [sd_ch_vals]        Percentage trim for confidence intervals during channel standard deviation criteria (e.g. [.3 .7])
%    sd_ch_o =           [sd_ch_o]           z threshold for flagging channels during standard deviation criteria (e.g. 6)
%    sd_ch_f_meth =      [sd_ch_f_meth]      Fixed method used for flagging channels (should be 'fixed')
%    sd_ch_f_vals  =     [sd_ch_f_vals]      Percentage trim for confidence intervals during channel standard deviation criteria (e.g. [.3 .7], leave empty for 'fixed')
%    sd_ch_f_o =         [sd_ch_f_o]         z threshold for flagging channel during fixed standard deviation criteria (e.g. 6)
31
32
33
%    ref_loc_file =      [ref_loc_file]      Name of file containing the reference locations (including the relative path)    
%    low_bound_hz =      [low_bound_hz]      Lower bound of the filter bass-band
%    high_bound_hz =     [high_bound_hz]     Upper bound of the filter bass-band
34
35
36
37
38
39
40
41
%    save_f_res =        [save_f_res]        1 if you want to save the filter residuals (Default: 1)
%    n_nbr_ch =          [n_nbr_ch]          Number of channels to use in nearest neighbour r calculation (for channel criteria)
%    r_ch_meth =         [sd_ch_meth]        Method used for flagging low r channels (e.g. 'q' (quantiles), or 'na' (default))
%    r_ch_vals  =        [sd_ch_vals]        Percentage trim for confidence intervals during low r channel standard deviation criteria (e.g. [.3 .7])
%    r_ch_o =            [r_ch_o]            z threshold for flagging channel during neighbour r criteria (e.g. 3)
%    r_ch_f_meth =       [r_ch_f_meth]       Fixed method used for flagging low r channels (should be 'fixed')
%    r_ch_f_vals  =      [r_ch_f_vals]       Percentage trim for confidence intervals during low r channel standard deviation criteria (e.g. [.3 .7], leave empty for 'fixed')
%    r_ch_f_o =          [r_ch_f_o]          z threshold for flagging channel during fixed neighbour r criteria (e.g. 3)
42
43
%    bridge_trim =       [bridge_trim]       Percentage trim for z calculation of bridged channels (e.g. 40 = 20% top and bottom)
%    bridge_z =          [bridge_z]          z threshold for flagging channel during bridging criteria (e.g. 8)
44
45
46
47
48
49
50
%    n_nbr_t =           [n_nbr_t]           Number of channels to use in nearest neighbour r calculation (for epoch criteria)
%    r_t_meth =          [r_t_meth]          Method used for flagging low r epochs (e.g. 'q' (quantiles), or 'na' (default))
%    r_t_vals  =         [r_t_vals]          Percentage trim for confidence intervals during low r epoch standard deviation criteria (e.g. [.3 .7])
%    r_t_o =             [r_t_o]             z threshold for flagging epochs during neighbour r criteria (e.g. 3)
%    r_t_f_meth =        [r_t_f_meth]        Fixed method used for flagging low r epochs (should be 'fixed')
%    r_t_f_vals  =       [r_t_f_vals]        Percentage trim for confidence intervals during low r epoch standard deviation criteria (e.g. [.3 .7], leave empty for 'fixed')
%    r_t_f_o =           [r_t_f_o]           z threshold for flagging epochs during fixed neighbour r criteria (e.g. 3)
51
52
%    min_gap_ms =        [min_gap_ms]        Minimum time (ms) to allow between periods marked for rejection
%    out_path =          [out_path]          Relative path to output data files assuming cd = work_path
53
54
%    amica_param_file =  [amica_param_file]  template amicadefs.param file to modify
%    amica_threads_s02 = [amica_threads_s02] number of threads to use for running s02 amica script (Default: 8)
55

56
%% CHECK FOR OUTPUT PATH AND CRETAE IF NECESSARY
57
58
59
if ~exist('[out_path]/[batch_dfn,/,-1]','dir');
    disp('Making directory [out_path]/[batch_dfn,/,-1]');
    mkdir [out_path]/[batch_dfn,/,-1]
60
61
end

62
63
64
65
66
67
68
69
70
71
72
%% LOAD DATASET
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%Check if it is a .set file
logging_log('NOTICE', 'Starting scalpart');
[path name ext] = fileparts('[batch_dfn]');
if ~strcmp(ext,'.set')
    logging_log('ERROR', sprintf('Wrong file format: %s should be .set',ext));
    exit(1);
end

73
%Load the .set file (this will need to be BIDS data file format compliant)
74
logging_log('INFO', 'Loading set file: [batch_dfn]...');
mikec1's avatar
mikec1 committed
75
EEG = pop_loadset('filepath','[batch_dfp]','filename','[batch_dfn]');
76
77
EEG = eeg_checkset( EEG );

78
79
80
81
82
83
84
if ~isfield(EEG,'marks');
    if isempty(EEG.icaweights)
        EEG.marks=marks_init(size(EEG.data));
    else
        EEG.marks=marks_init(size(EEG.data),min(size(EEG.icaweights)));
    end
end
85
86
87

EEG.marks = marks_add_label(EEG.marks,'time_info', ...
{'init_ind',[0,0,1],[1:EEG.pnts]});
88

89
90
91
92
%execute the staging script if specified.
if exist('[staging_script]');
    [ssp,ssn,sse]=fileparts('[staging_script]');
    addpath(ssp);
Tyler Collins's avatar
Tyler Collins committed
93
    eval(ssn);
94
95
end

96
%add channel coregistration procedure here.
97
%warp locations to standard head surface:
98
if ~isempty('[montage_info]');
99
    if isempty(str2num('[montage_info]'));
100
        %if it is a BIDS channel tsv, load the tsv,sd_t_f_vals
101
        %else read the file that is assumed to be a transformation matrix.
102
103
104
105
    else
        EEG = warp_locs( EEG,'[ref_loc_file]', ...
            'transform',[[montage_info]], ...
            'manual','off');
106
107
108
    end
end

109
110
111
112
113
114
% Apply average re-reference
%chan_inds = marks_label2index(EEG.marks.chan_info,{'manual'},'indexes','invert','on');
%chan_m = mean(EEG.data(chan_inds,:),1);
%chan_m_mat = repmat(chan_m,size(EEG.data,1),1);
%EEG.data = EEG.data - chan_m_mat;
%clear chan_m_mat;
115
116
117
118
119
120

% Window the continuous data
logging_log('INFO', 'Windowing the continous data...');
EEG = marks_continuous2epochs(EEG,'recurrence',[[recur_sec]],'limits',[[limit_sec]]);
logging_log('INFO', 'LOADED DATASET...');

121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
% rereference to selected channels..
chan_inds = marks_label2index(EEG.marks.chan_info,{'manual'},'indexes','invert','on');
epoch_inds = marks_label2index(EEG.marks.time_info,{'manual'},'indexes','invert','on');
[EEG,trim_ch_sd]=chan_variance(EEG,'data_field','data', ...
         'chan_inds',chan_inds, ...
         'epoch_inds',epoch_inds, ...
         'plot_figs','off');

chan_dist=zeros(size(trim_ch_sd));
for i=1:size(trim_ch_sd,2);
    chan_dist(:,i)=(trim_ch_sd(:,i)-median(trim_ch_sd(:,i)))/diff(quantile(trim_ch_sd(:,i),[.3,.7]));
end
mean_chan_dist=mean(chan_dist,2);
m=median(mean_chan_dist);
q=quantile(mean_chan_dist,[.3,.7]);

refchans=find(mean_chan_dist<m+6*diff(q));

EEG.data=EEG.data-repmat(mean(EEG.data(chan_inds(refchans),:,:),1),size(EEG.data,1),1);

141
142
143
144
145
146
147
148
149
%% CALCULATE DATA SD
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This flag calculates the standard deviation  of the channels. Epochs are flagged if they
% are above the SD critical value. This flag identifies comically bad epochs.

% Calculate standard deviation of activation on non-'manual' flagged channels and epochs...
logging_log('INFO', 'Calculating the data sd array for time criteria...');
chan_inds = marks_label2index(EEG.marks.chan_info,{'manual'},'indexes','invert','on');
epoch_inds = marks_label2index(EEG.marks.time_info,{'manual'},'indexes','invert','on');
150
[EEG,data_sd_t]=chan_variance(EEG,'data_field','data', ...
151
152
153
154
155
156
157
    'chan_inds',chan_inds, ...
    'epoch_inds',epoch_inds, ...
    'plot_figs','off');
logging_log('INFO', 'CALCULATED EPOCH SD...');

% Create the window criteria vector for flagging ch_sd time_info...
logging_log('INFO', 'Assessing window data sd distributions...')
158
if strcmp('[sd_t_meth]','na');
159
160
    flag_sd_t_inds=[];
else
161
    [~,flag_sd_t_inds]=marks_array2flags(data_sd_t, ...
162
        'flag_dim','col', ...
163
        'init_method','[sd_t_meth]', ...
164
165
        'init_vals',[sd_t_vals], ...
        'init_crit',[sd_t_o], ...
166
167
168
        'flag_method','[sd_t_f_meth]',... % 'fixed', ...
        'flag_vals',[ [sd_t_f_vals] ], ... % NEW
        'flag_crit',[sd_t_f_o], ... % 'flag_val',[sd_t_p], ...
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
        'plot_figs','off');
end
logging_log('INFO', 'CREATED EPOCH CRITERIA VECTOR...');

% Edit the time flag info structure
logging_log('INFO', 'Updating epflaginfo structure...');
chsd_epoch_flags = zeros(size(EEG.data(1,:,:)));
chsd_epoch_flags(1,:,epoch_inds(flag_sd_t_inds))=1;
chsd_epoch_flags=padflags(EEG,chsd_epoch_flags,[sd_t_pad],'value',.5);
EEG.marks = marks_add_label(EEG.marks,'time_info', ...
{'ch_sd',[0,0,1],chsd_epoch_flags});
logging_log('INFO', 'EDITED TIMEFLAGINFO STRUCT...');

% Combine ch_sd time_info flags into 'manual'...
EEG = pop_marks_merge_labels(EEG,'time_info',{'ch_sd'},'target_label','manual');
logging_log('INFO', 'COMBINED MARKS STRUCTURE INTO MANUAL FLAGS...');

%% CALCULATE DATA SD
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This flag calculates the standard deviation  of the channels. Channels are flagged if they
% are above the SD critical value. This flag identifies comically bad channels.

% Calculate standard deviation of activation on non-'manual' flagged channels and epochs...
logging_log('INFO', 'Calculating the data sd array for channel criteria...');
chan_inds = marks_label2index(EEG.marks.chan_info,{'manual'},'indexes','invert','on');
epoch_inds = marks_label2index(EEG.marks.time_info,{'manual'},'indexes','invert','on');
195
[EEG,data_sd_ch]=chan_variance(EEG,'data_field','data', ...
196
197
198
199
200
201
202
    'chan_inds',chan_inds, ...
    'epoch_inds',epoch_inds, ...
    'plot_figs','off');
logging_log('INFO', 'CALCULATED CHAN SD...');

% Create the window criteria vector for flagging ch_sd chan_info...
logging_log('INFO', 'Assessing window data sd distributions...')
203
[~,flag_sd_ch_inds]=marks_array2flags(data_sd_ch, ...
204
    'flag_dim','row', ...
205
    'init_method','[sd_ch_meth]', ...
206
207
    'init_vals',[sd_ch_vals], ...
    'init_crit',[sd_ch_o], ...
208
209
210
    'flag_method','[sd_ch_f_meth]', ... %fixed
    'flag_vals',[ [sd_ch_f_vals] ], ...
    'flag_crit',[sd_ch_f_o], ...
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
    'plot_figs','off');
logging_log('INFO', 'CREATED CHANNEL CRITERIA VECTOR...');

% Edit the channel flag info structure
logging_log('INFO', 'Updating chflaginfo structure...');
chsd_chan_flags = zeros(EEG.nbchan,1);
chsd_chan_flags(chan_inds(flag_sd_ch_inds)) = 1;
EEG.marks = marks_add_label(EEG.marks,'chan_info', ...
{'ch_sd',[.7,.7,1],[.2,.2,1],-1,chsd_chan_flags});
logging_log('INFO', 'EDITED CHANFLAGINFO STRUCT...');

% Combine ch_sd chan_info flags into 'manual'...
EEG = pop_marks_merge_labels(EEG,'chan_info',{'ch_sd'},'target_label','manual');
logging_log('INFO', 'COMBINED MARKS STRUCTURE INTO MANUAL FLAGS...');

% Concatenate epoched data back to continuous data
logging_log('INFO', 'Concatenating windowed data...');
EEG = marks_epochs2continuous(EEG);
EEG = eeg_checkset(EEG,'eventconsistency');
logging_log('INFO', 'CONCATENATED THE WINDOWED DATA INTO CONTINUOUS DATA...');

%% REREFERENCE TO INTERPOLATED AVERAGE SITE
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Rereference the data to an average interpolated site containing 19 channels
% ... excluding 'manual' flagged channels from the calculation

logging_log('INFO', 'Rereferencing to averaged interpolated site...');
chan_inds = marks_label2index(EEG.marks.chan_info,{'manual'},'indexes','invert','on');
EEG = interp_ref(EEG,'[ref_loc_file]','chans',chan_inds);
EEG = eeg_checkset(EEG);
logging_log('INFO', 'REREFERENCED TO INTERPOLATED AVERAGE SITE...');

%% FILTER HIGH PASS
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Filters the data to remove frequencies lower than the selected value. The residuals that
% were removed can be saved for further analysis if needed.
logging_log('INFO', 'High pass filtering the data...');
248
if([[low_bound_hz]]);
249
    pre_hp_data = EEG.data;
250
    EEG = pop_eegfiltnew(EEG,[],[[low_bound_hz]],[],true,[],0);
251
252
253
254
255
256
257
258
259
260
261
262
    EEG.setname = 'filtHP';
    EEG = eeg_checkset( EEG );
    if [save_f_res]
        % Save high pass filter residual data array
        logging_log('DEBUG', 'Saving file: [out_path]/[batch_dfn,_,-1]_hpd.mat...');
        data = pre_hp_data - EEG.data;
        try OCTAVE_VERSION;
            save('-mat7-binary','[out_path]/[batch_dfn,_,-1]_hpd.mat','data');
        catch
            save('[out_path]/[batch_dfn,_,-1]_hpd.mat','data');
        end
        logging_log('DEBUG', 'TIME TO: SAVE [out_path]/[batch_dfn,_,-1]_hpd.mat FILE...');
263
    end
264
265
266
267
    logging_log('INFO', 'FILTERED HIGH PASS...');
    clear pre_hp_data
else
    logging_log('INFO', 'High pass step skipped...');
268
end
269
270


271
272
273
274
275
276

%% FILTER LOW PASS
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Filters the data to remove frequencies higher the selected value. The residuals that
% were removed can be saved for further analysis if needed.
logging_log('INFO', 'Low pass filtering the data...');
277
if([[high_bound_hz]]);
278
    pre_lp_data = EEG.data;
279
    EEG = pop_eegfiltnew(EEG,[],[[high_bound_hz]],[],0,[],0);
280
281
282
283
284
285
286
287
288
289
290
291
    EEG.setname = 'filtLP';
    EEG = eeg_checkset( EEG );
    if [save_f_res]
        % Save high pass filter residual data array
        logging_log('DEBUG', 'Saving file: [out_path]/[batch_dfn,_,-1]_lpd.mat...');
        data = pre_lp_data - EEG.data;
        try OCTAVE_VERSION;
            save('-mat7-binary','[out_path]/[batch_dfn,_,-1]_lpd.mat','data');
        catch
            save('[out_path]/[batch_dfn,_,-1]_lpd.mat','data');
        end
        logging_log('DEBUG', 'TIME TO: SAVE [out_path]/[batch_dfn,_,-1]_lpd.mat FILE...');
292
    end
293
294
295
296
    logging_log('INFO', 'FILTERED LOW PASS...');
    clear pre_lp_data
else
    logging_log('INFO', 'Low pass filtering skipped...');
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
end

%% CALCULATE NEAREST NEIGHBOUR R VALUES
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Checks neighboring channels for too high or low of a correlation.

% Window the continuous data
logging_log('INFO', 'Windowing the continous data...');
EEG = marks_continuous2epochs(EEG,'recurrence',[[recur_sec]],'limits',[[limit_sec]]);
logging_log('INFO', 'WINDOWED THE CONTINUOUS DATA...');

% Calculate nearest neighbout correlation on non-'manual' flagged channels and epochs...
logging_log('INFO', 'Calculating nearst neighbour r array for channel criteria...');
chan_inds = marks_label2index(EEG.marks.chan_info,{'manual'},'indexes','invert','on');
epoch_inds = marks_label2index(EEG.marks.time_info,{'manual'},'indexes','invert','on');
312
[EEG,data_r_ch,~,~,~] = chan_neighbour_r(EEG, ...
313
314
315
316
317
318
319
    [n_nbr_ch],'max', ...
    'chan_inds',chan_inds, ...
    'epoch_inds',epoch_inds, ...
    'plot_figs','off');

% Create the window criteria vector for flagging low_r chan_info...
logging_log('INFO', 'Assessing channel r distributions criteria...')
320
[~,flag_r_ch_inds]=marks_array2flags(data_r_ch, ...
321
    'flag_dim','row', ...
322
    'init_method','[r_ch_meth]', ...
323
    'init_vals',[r_ch_vals], ...
324
    'init_dir','neg', ...
325
    'init_crit',[r_ch_o], ...
326
327
328
    'flag_method','[r_ch_f_meth]', ...
    'flag_vals',[ [r_ch_f_vals] ], ...
    'flag_crit',[r_ch_f_o], ...
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
    'plot_figs','off');

% Edit the channel flag info structure
logging_log('INFO', 'Updating chflaginfo structure...');
lowr_chan_flags = zeros(EEG.nbchan,1);
lowr_chan_flags(chan_inds(flag_r_ch_inds))=1;
EEG.marks = marks_add_label(EEG.marks,'chan_info', ...
    {'low_r',[1,.7,.7],[1,.2,.2],-1,lowr_chan_flags});

logging_log('INFO', 'CALCULATED NEAREST NEIGHBOUR R VALUES...');

%% IDENTIFY BRIDGED CHANNELS
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Uses the correlation of neighboors calculated to flag bridged channels.

logging_log('INFO', 'Examing nearest neighbour r array for linked channels...');
345
346
mr = median(data_r_ch,2);
iqrr = iqr(data_r_ch,2);
347
msr = mr./iqrr;
348
flag_b_chan_inds = find(msr>ve_trimmean(msr,[bridge_trim],1)+ve_trimstd(msr,[bridge_trim],1)*[bridge_z]);
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
logging_log('INFO', 'IDENTIFIED BRIDGED CHANNELS...');

% Edit the channel flag info structure
logging_log('INFO', 'Updating chflaginfo structure...');
lnkflags = zeros(EEG.nbchan,1);
lnkflags(chan_inds(flag_b_chan_inds))=1;
EEG.marks = marks_add_label(EEG.marks,'chan_info', ...
    {'bridge',[.7,1,.7],[.2,1,.2],-1,lnkflags});
logging_log('INFO', 'EDITED CHANFLAGINFO STRUCT...');

% Combine low_rand bridge chan_info flags into 'manual'...
EEG = pop_marks_merge_labels(EEG,'chan_info',{'low_r','bridge'},'target_label','manual');
logging_log('INFO', 'COMBINED MARKS STRUCTURE INTO MANUAL FLAGS...');

%% FLAG RANK CHAN
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Flags the channel that is the least unique (the best channel to remove prior 
% to ICA in order to account for the rereference rank deficiency.

logging_log('INFO', 'Updating chflaginfo structure...');
chan_inds = marks_label2index(EEG.marks.chan_info,{'manual'},'indexes','invert','on');
370
[r_max,rank_ind] = max(data_r_ch(chan_inds));
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
rankflags = zeros(EEG.nbchan,1);
rankflags(chan_inds(rank_ind))=1;
EEG.marks = marks_add_label(EEG.marks,'chan_info', ...
    {'rank',[.1,.1,.24],[.1,.1,.24],-1,rankflags});
logging_log('INFO', 'EDITED CHANFLAGINFO STRUCT...');

%% REREFERENCE TO INTERPOLATED AVERAGE SITE
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Rereference the data to an average interpolated site containined 19 channels
logging_log('INFO', 'Rereferencing to averaged interpolated site...');
chan_inds = marks_label2index(EEG.marks.chan_info,{'manual'},'indexes','invert','on');
EEG = interp_ref(EEG,'[ref_loc_file]','chans',chan_inds);
EEG = eeg_checkset(EEG);
logging_log('INFO', 'REREFERENCED TO INTERPOLATED AVERAGE SITE...');

%% CALCULATE NEAREST NEIGHBOUR R VALUES
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Similarly to the neighbor r calculation done between channels this section looks
% at the correlation, but between all channels and for epochs of time. Time segmenents 
% are flagged for removal.

logging_log('INFO', 'Calculating nearest neighbour r array for window criteria...');
chan_inds = marks_label2index(EEG.marks.chan_info,{'manual'},'indexes','invert','on');
epoch_inds = marks_label2index(EEG.marks.time_info,{'manual'},'indexes','invert','on');
395
[EEG,data_r_t,~,~,~] = chan_neighbour_r(EEG, ...
396
397
398
399
400
401
402
403
    [n_nbr_t],'max', ...
    'chan_inds',chan_inds, ...
    'epoch_inds',epoch_inds, ...
    'plotfigs','off');
logging_log('INFO', 'CALCULATED NEAREST NEIGHBOUR R VALUES...');

% Create the window criteria vector for flagging low_r time_info...
logging_log('INFO', 'Assessing epoch r distributions criteria...')
404
[~,flag_r_t_inds]=marks_array2flags(data_r_t, ...
405
    'flag_dim','col', ...
406
    'init_method','[r_t_meth]', ...
407
    'init_vals',[r_t_vals], ...
408
    'init_dir','neg', ...
409
    'init_crit',[r_t_o], ...
410
411
412
    'flag_method','[r_t_f_meth]', ...
    'flag_vals',[ [r_t_f_vals] ], ...
    'flag_crit',[r_t_f_o], ...
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
    'plot_figs','off');
logging_log('INFO', 'CREATED EPOCH CRITERIA VECTOR...');

% Edit the time flag info structure
logging_log('INFO', 'Updating latflaginfo structure...');
lowr_epoch_flags = zeros(size(EEG.data(1,:,:)));
lowr_epoch_flags(1,:,epoch_inds(flag_r_t_inds))=1;
EEG.marks = marks_add_label(EEG.marks,'time_info', ...
    {'low_r',[0,1,0],lowr_epoch_flags});
clear lowr_epoch_flags;
logging_log('INFO', 'TIME TO: UPDATE REJECTION STRUCTURE...');

% Combine low_r time_info flags into 'manual'...
EEG = pop_marks_merge_labels(EEG,'time_info',{'low_r'},'target_label','manual');
logging_log('INFO', 'COMBINED MARKS STRUCTURE INTO MANUAL FLAGS...');

% Concatenate epoched data back to continuous data
logging_log('INFO', 'Concatenating windowed data...');
EEG = marks_epochs2continuous(EEG);
EEG = eeg_checkset(EEG,'eventconsistency');
logging_log('INFO', 'CONCATENATED THE WINDOWED DATA INTO CONTINUOUS DATA...');

EEG=pop_marks_flag_gap(EEG,{'manual'},[min_gap_ms],'mark_gap',[.8,.8,.8],'offsets',[0 0],'ref_point','both');

% Combine mark_gap time_info flags into 'manual'...
EEG = pop_marks_merge_labels(EEG,'time_info',{'mark_gap'},'target_label','manual');
logging_log('INFO', 'COMBINED MARKS STRUCTURE INTO MANUAL FLAGS...');

%%%-----------------------------------------------------------------------
%% SAVE sa.set FILE
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
logging_log('INFO', 'Saving file: [batch_dfn,_,-1]_sa.set...');
445
446
if ~exist('[out_path]/[batch_dfn,/,-1]','dir');
    mkdir('[out_path]/[batch_dfn,/,-1]');
mikec1's avatar
mikec1 committed
447
end
448
EEG = pop_saveset( EEG, 'filename','[out_path]/[batch_dfn,_,-1]_sa.set');
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
logging_log('INFO', 'SAVED sa FILE...');

%%%-----------------------------------------------------------------------
%% PREPARE FILE FOR AMICA
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Purges the flagged data to make sure it is primed for ICA and creates the amica paramaeter file.

% Removed flagged channels and time segments
sprintf('%s','Purging flagged channels...\n');
EEG = pop_marks_select_data(EEG,'channel marks',[],'labels',{'manual','rank'},'remove','on');
EEG = pop_marks_select_data(EEG,'time marks',[],'labels',{'manual'},'remove','on');
EEG = eeg_checkset(EEG);
logging_log('INFO', 'TIME TO: PURGE DATA...');

% Make sure EEG.data is doubles for use in amica
EEG.data = double(EEG.data);

% Save the SA purge file
467
468
logging_log('INFO',sprintf('%s','Saving file: [out_path]/[batch_dfn,_,-1]_sa_purge.set...'));
EEG = pop_saveset( EEG, 'filename','[out_path]/[batch_dfn,_,-1]_sa_purge.set');
469

470
471
% Save diagnostic arrays
try OCTAVE_VERSION;
472
473
474
475
    save('-mat7-binary','[out_path]/[batch_dfn,_,-1]_data_sd_t.mat','data_sd_t');
    save('-mat7-binary','[out_path]/[batch_dfn,_,-1]_data_sd_ch.mat','data_sd_ch');
    save('-mat7-binary','[out_path]/[batch_dfn,_,-1]_data_r_ch.mat','data_r_ch');
    save('-mat7-binary','[out_path]/[batch_dfn,_,-1]_data_r_t.mat','data_r_t');
476
catch
477
478
479
480
    save('[out_path]/[batch_dfn,_,-1]_data_sd_t.mat','data_sd_t');
    save('[out_path]/[batch_dfn,_,-1]_data_sd_ch.mat','data_sd_ch');
    save('[out_path]/[batch_dfn,_,-1]_data_r_ch.mat','data_r_ch');
    save('[out_path]/[batch_dfn,_,-1]_data_r_t.mat','data_r_t');
481
482
end

483
484
% Save the AMICA parameter file
swapstr = file_strswap('[amica_param_file]', ...
485
486
    '[repstr_fdt_fname]','[out_path]/[batch_dfn,_,-1]_sa_purge.fdt', ...
    '[repstr_outpath]','[out_path]/[batch_dfn,_,-1]_amicaout_init', ...
487
488
489
490
    '[repstr_nbchan]', num2str(EEG.nbchan), ...
    '[repstr_pnts]', sprintf('%12.0f',EEG.pnts), ...
    '[repstr_tpp]', '[amica_threads_s02]', ...
    '[repstr_pca]', num2str(EEG.nbchan));
491
fidparam = fopen('[out_path]/[batch_dfn,_,-1]_init.param','w');
492
493
494
495
496
497
498
fprintf(fidparam,'%s',swapstr);
fclose(fidparam);

% CREATE THE AMICA OUTPUT FOLDER
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Needed for MPI jobs on cluster only in order to prevent a race effect bug.

499
500
[~] = rmdir('[out_path]/[batch_dfn,_,-1]_amicaout_init','s');
mkdir [out_path]/[batch_dfn,_,-1]_amicaout_init;
501
502
503
504
505
logging_log('INFO', 'Created Amica Output Folder');

logging_log('INFO', 'Scalpart complete');
logging_log('INFO', 'Scheduler: [scheduler]');
print_chan_sample(EEG);