s01_scalpart.htb 22.6 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
13
14
%    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])
%                                            ... or a a file contasining a transformation matrix (e.g. [batch_dfn]_transmat.mat).
%                                            ... 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
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
%	 aref_trim =		 [aref_trim]         Trim of mean for initial average rereference (e.g. 30 = 15% off each side)
% 	 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])
%    sd_epoch_z =        [sd_epoch_z]        z threshold for flagging epochs during standard deviation criteria (e.g. 6)
%    sd_epoch_p =        [sd_epoch_p]        Percentage of flagged channels required to flag an epoch during standard deviation criteria (e.g. .1)
%    sd_epoch_trim =     [sd_epoch_trim]     Percentage trim for confidence intervals during epoch standard deviation criteria (10 = 5% top and bottom)
%    sd_epoch_pad =      [sd_epoch_pad]      Number of windows to pad onto each side of the ch_sd time flag
%    sd_chan_z =         [sd_chan_z]         z threshold for flagging channel during standard deviation criteria (e.g. 6)
%    sd_chan_p =         [sd_chan_p]         Percentage of flagged epochs required to flag a channel during standard deviation criteria (e.g. .1)
%    sd_chan_trim =      [sd_chan_trim]      Percentage trim for confidence intervals during channel standard deviation criteria (10 = 5% top and bottom)
%    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
%    n_neighb_chan =     [n_neighb_chan]     Number of channels to use in nearest neighbour r calculation (for channel criteria)
%    r_chan_z =          [r_chan_z]          z threshold for flagging channel during neighbour r criteria (e.g. 3)
%    r_chan_p =          [r_chan_p]          Percentage of flagged epochs required to flag a channel during neighbour r criteria (e.g. .1)
%    r_chan_trim =       [r_chan_trim]       Percentage trim for confidence intervals during channel neighbour r criteria (10 = 5% top and bottom)
%    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)
%    n_neighb_epoch =    [n_neighb_epoch]    Number of channels to use in nearest neighbour r calculation (for epoch criteria)
%    epoch_z =           [epoch_z]           z threshold for flagging epochs during neighbour r criteria (e.g. 3)
%    epoch_p =           [epoch_p]           Percentage of flagged channels required to flag an epoch during neighbour r criteria (e.g. .1)
%    epoch_trim =        [epoch_trim]        Percentage trim for confidence intervals during epoch neighbour r criteria (10 = 5% top and bottom)
%    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
%    amica_param_file =  [amica_param_file]  template amicadefs.param file to modify.
%    save_diff =         [savediff]          YES if you want to save the filter residuals (Default: NO)

44
%% CHECK FOR OUTPUT PATH AND CRETAE IF NECESSARY
45
46
47
if ~exist('[out_path]/[batch_dfn,/,-1]','dir');
    disp('Making directory [out_path]/[batch_dfn,/,-1]');
    mkdir [out_path]/[batch_dfn,/,-1]
48
49
end

50
51
52
53
54
55
56
57
58
59
60
%% 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

61
%Load the .set file (this will need to be BIDS data file format compliant)
62
logging_log('INFO', 'Loading set file: [batch_dfn]...');
mikec1's avatar
mikec1 committed
63
EEG = pop_loadset('filepath','[batch_dfp]','filename','[batch_dfn]');
64
65
EEG = eeg_checkset( EEG );

66
67
68
69
%add channel coregistration procedure here.
%warp locations to standard head surface: 
if ~isempty('[montage_info]');
    if isnumeric([montage_info])
James Desjardins's avatar
James Desjardins committed
70
        EEG = warp_locs( EEG,'[ref_loc_file]', ...
71
72
73
74
75
76
77
78
79
            'transform',[montage_info], ...
            'manual','off');
    else
        %if it is a BIDS channel tsv, load the tsv,
        %else read the file that is assumed to be a transformation matrix.
    end
end

%execute the event marks initiation script.
80
if exist(['[staging_script]', '.m']) || exist('[staging_script]');
James Desjardins's avatar
James Desjardins committed
81
    run [staging_script]
82
83
end

James Desjardins's avatar
James Desjardins committed
84
85
86
87
88
89
90
91
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

92
% Apply trimmed average re-reference
93
94
chan_inds = marks_label2index(EEG.marks.chan_info,{'manual'},'indexes','invert','on');
trm_m = ve_trimmean(EEG.data(chan_inds,:),[aref_trim],1);
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
trm_m_mat = repmat(trm_m,size(EEG.data,1),1);
EEG.data = EEG.data - trm_m_mat;
clear trm_m_mat;

% 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...');

%% 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');
113
[EEG,data_sd_t]=chan_variance(EEG,'data_field','data', ...
114
115
116
117
118
119
120
    '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...')
121
if strcmp('[sd_t_meth]','na');
122
123
    flag_sd_t_inds=[];
else
124
    [~,flag_sd_t_inds]=marks_array2flags(data_sd_t, ...
125
        'flag_dim','col', ...
126
        'init_method','[sd_t_meth]', ...
127
128
        'init_vals',[sd_t_vals], ...
        'init_crit',[sd_t_o], ...
129
130
131
        '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], ...
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
        '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');
158
[EEG,data_sd_ch]=chan_variance(EEG,'data_field','data', ...
159
160
161
162
163
164
165
    '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...')
166
[~,flag_sd_ch_inds]=marks_array2flags(data_sd_ch, ...
167
    'flag_dim','row', ...
168
    'init_method','[sd_ch_meth]', ...
169
170
    'init_vals',[sd_ch_vals], ...
    'init_crit',[sd_ch_o], ...
171
172
173
    'flag_method','[sd_ch_f_meth]', ... %fixed
    'flag_vals',[ [sd_ch_f_vals] ], ...
    'flag_crit',[sd_ch_f_o], ...
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
    '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...');
212
pre_hp_data = EEG.data;
213
214
215
EEG = pop_eegfiltnew(EEG,[],[low_bound_hz],[],true,[],0);
EEG.setname = 'filtHP';
EEG = eeg_checkset( EEG );
216
if [save_f_res]
217
    % Save high pass filter residual data array
218
    logging_log('DEBUG', 'Saving file: [out_path]/[batch_dfn,_,-1]_hpd.mat...');
219
220
    data = pre_hp_data - EEG.data;
    try OCTAVE_VERSION;
221
        save('-mat7-binary','[out_path]/[batch_dfn,_,-1]_hpd.mat','data');
222
    catch
223
        save('[out_path]/[batch_dfn,_,-1]_hpd.mat','data');
224
    end
225
    logging_log('DEBUG', 'TIME TO: SAVE [out_path]/[batch_dfn,_,-1]_hpd.mat FILE...');
226
227
end
logging_log('INFO', 'FILTERED HIGH PASS...');
228
clear pre_hp_data
229
230
231
232
233
234
235

%% 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...');
236
pre_lp_data = EEG.data;
237
238
239
EEG = pop_eegfiltnew(EEG,[],[high_bound_hz],[],0,[],0);
EEG.setname = 'filtLP';
EEG = eeg_checkset( EEG );
240
if [save_f_res]
241
    % Save high pass filter residual data array
242
    logging_log('DEBUG', 'Saving file: [out_path]/[batch_dfn,_,-1]_lpd.mat...');
243
244
    data = pre_lp_data - EEG.data;
    try OCTAVE_VERSION;
245
        save('-mat7-binary','[out_path]/[batch_dfn,_,-1]_lpd.mat','data');
246
    catch
247
        save('[out_path]/[batch_dfn,_,-1]_lpd.mat','data');
248
    end
249
    logging_log('DEBUG', 'TIME TO: SAVE [out_path]/[batch_dfn,_,-1]_lpd.mat FILE...');
250
251
end
logging_log('INFO', 'FILTERED LOW PASS...');
252
clear pre_lp_data
253
254
255
256
257
258
259
260
261
262
263
264
265
266

%% 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');
267
[EEG,data_r_ch,~,~,~] = chan_neighbour_r(EEG, ...
268
269
270
271
272
273
274
    [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...')
275
[~,flag_r_ch_inds]=marks_array2flags(data_r_ch, ...
276
    'flag_dim','row', ...
277
    'init_method','[r_ch_meth]', ...
278
    'init_vals',[r_ch_vals], ...
279
    'init_dir','neg', ...
280
    'init_crit',[r_ch_o], ...
281
282
283
    'flag_method','[r_ch_f_meth]', ...
    'flag_vals',[ [r_ch_f_vals] ], ...
    'flag_crit',[r_ch_f_o], ...
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
    '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...');
300
301
mr = median(data_r_ch,2);
iqrr = iqr(data_r_ch,2);
302
msr = mr./iqrr;
303
flag_b_chan_inds = find(msr>ve_trimmean(msr,[bridge_trim],1)+ve_trimstd(msr,[bridge_trim],1)*[bridge_z]);
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
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');
325
[r_max,rank_ind] = max(data_r_ch(chan_inds));
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
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');
350
[EEG,data_r_t,~,~,~] = chan_neighbour_r(EEG, ...
351
352
353
354
355
356
357
358
    [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...')
359
[~,flag_r_t_inds]=marks_array2flags(data_r_t, ...
360
    'flag_dim','col', ...
361
    'init_method','[r_t_meth]', ...
362
    'init_vals',[r_t_vals], ...
363
    'init_dir','neg', ...
364
    'init_crit',[r_t_o], ...
365
366
367
    'flag_method','[r_t_f_meth]', ...
    'flag_vals',[ [r_t_f_vals] ], ...
    'flag_crit',[r_t_f_o], ...
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
    '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...');
400
401
if ~exist('[out_path]/[batch_dfn,/,-1]','dir');
    mkdir('[out_path]/[batch_dfn,/,-1]');
mikec1's avatar
mikec1 committed
402
end
403
EEG = pop_saveset( EEG, 'filename','[out_path]/[batch_dfn,_,-1]_sa.set');
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
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
422
423
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');
424

425
426
% Save diagnostic arrays
try OCTAVE_VERSION;
427
428
429
430
    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');
431
catch
432
433
434
435
    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');
436
437
end

438
439
% Save the AMICA parameter file
swapstr = file_strswap('[amica_param_file]', ...
440
441
    '[repstr_fdt_fname]','[out_path]/[batch_dfn,_,-1]_sa_purge.fdt', ...
    '[repstr_outpath]','[out_path]/[batch_dfn,_,-1]_amicaout_init', ...
442
443
444
445
    '[repstr_nbchan]', num2str(EEG.nbchan), ...
    '[repstr_pnts]', sprintf('%12.0f',EEG.pnts), ...
    '[repstr_tpp]', '[amica_threads_s02]', ...
    '[repstr_pca]', num2str(EEG.nbchan));
446
fidparam = fopen('[out_path]/[batch_dfn,_,-1]_init.param','w');
447
448
449
450
451
452
453
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.

454
455
[~] = rmdir('[out_path]/[batch_dfn,_,-1]_amicaout_init','s');
mkdir [out_path]/[batch_dfn,_,-1]_amicaout_init;
456
457
458
459
460
logging_log('INFO', 'Created Amica Output Folder');

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