OSCR

Neonatal brain activity across sleep states: Evidence from resting EEG and auditory event-related potentials.

Code ↔ Paper

3 matches between paragraphs of the paper and lines of its authors' code, computed by the harvester (lexical-v1). Click a colored paragraph or line to see its counterpart.

The 3 matches · 1 of them tie a paragraph to a whole file, not to given lines: a weak match, whose lines are not tinted
  1. [1] § Methods › EEG preprocessing ↔ SA-MADE_pipeline_multiple_participants.m, lines 1–137 · score 0.91 · low pass filtered, high pass filtered, EEGLAB toolbox, voltage threshold, Developmental EEG, pipeline
  2. [2] § Methods › EEG preprocessing ↔ SA-MADE_pipeline_multiple_participants.m, lines 1–137 · score 0.61 · Missing channels, pipeline, millisecond, thresholding, voltages, segments
  3. [3] § Methods › Data acquisition ↔ Labeling_MMN_multi_participant.m, the whole file · a weak match · score 0.59 · standard tones, novel tone, blocks, deviant

Paper

Loaded from Europe PMC by your browser, not stored by OSCR: doi.org · Europe PMC

The paper is loaded when this pane is shown.

The authors' code

MATLAB · 1,672 lines · 101 KB · no license · 2 matches

  1. % ************************************************************************
  2. % The Maryland Analysis of Developmental EEG (UMADE) Pipeline
  3. % Version 1.0
  4. % Developed at the Child Development Lab, University of Maryland, College Park
  5. % Contributors to MADE pipeline:
  6. % Ranjan Debnath ([email hidden])
  7. % George A. Buzzell ([email hidden])
  8. % Santiago Morales Pamplona ([email hidden])
  9. % Stephanie Leach ([email hidden])
  10. % Maureen Elizabeth Bowers ([email hidden])
  11. % Nathan A. Fox ([email hidden])
  12. % MADE uses EEGLAB toolbox and some of its plugins. Before running the pipeline, you have to install the following:
  13. % EEGLab: https://sccn.ucsd.edu/eeglab/downloadtoolbox.php/download.php
  14. % You also need to download the following plugins/extensions from here: https://sccn.ucsd.edu/wiki/EEGLAB_Extensions
  15. % Specifically, download:
  16. % MFFMatlabIO: https://github.com/arnodelorme/mffmatlabio/blob/master/README.txt
  17. % FASTER: https://sourceforge.net/projects/faster/
  18. % ADJUST: https://www.nitrc.org/projects/adjust/
  19. % Adjusted ADJUST (included in this pipeline): https://github.com/ChildDevLab/MADE-EEG-preprocessing-pipeline
  20. % After downloading these plugins (as zip files), you need to place it in the eeglab/plugins folder.
  21. % For instance, for FASTER, you uncompress the downloaded extension file (e.g., 'FASTER.zip') and place it in the main EEGLAB "plugins" sub-directory/sub-folder.
  22. % After placing all the required plugins, add the EEGLAB folder to your path by using the following code:
  23. % addpath(genpath(('...')) % Enter the path of the EEGLAB folder in this line
  24. % Please cite the following references for in any manuscripts produced utilizing MADE pipeline:
  25. % EEGLAB: A Delorme & S Makeig (2004) EEGLAB: an open source toolbox for
  26. % analysis of single-trial EEG dynamics. Journal of Neuroscience Methods, 134, 9?21.
  27. % firfilt (filter plugin): developed by Andreas Widmann (https://home.uni-leipzig.de/biocog/content/de/mitarbeiter/widmann/eeglab-plugins/)
  28. % FASTER: Nolan, H., Whelan, R., Reilly, R.B., 2010. FASTER: Fully Automated Statistical
  29. % Thresholding for EEG artifact Rejection. Journal of Neuroscience Methods, 192, 152?162.
  30. % ADJUST: Mognon, A., Jovicich, J., Bruzzone, L., Buiatti, M., 2011. ADJUST: An automatic EEG
  31. % artifact detector based on the joint use of spatial and temporal features. Psychophysiology, 48, 229?240.
  32. % Our group has modified ADJUST plugin to improve selection of ICA components containing artifacts
  33. % This pipeline is released under the GNU General Public License version 3.
  34. % ************************************************************************
  35. % User input: user provide relevant information to be used for data processing
  36. % Preprocessing of EEG data involves using some common parameters for
  37. % every subject. This part of the script initializes the common parameters.
  38. clear % clear matlab workspace
  39. clc % clear matlab command window
  40. eeglab % restart eeglab
  41. % 1. Enter the path of the folder that has the raw data to be analyzed
  42. rawdata_location = 'lorem ipsum';
  43. % 2. Enter the path of the folder where you want to save the processed data
  44. output_location = 'lorem ipsum';
  45. scripts_location = 'lorem ipsum';
  46. % 3. Enter the path of the folder that has the individual sleep state spreadsheets
  47. event_location = 'lorem ipsum';
  48. % 3. Enter the path of the channel location.sfp file
  49. channel_locations = 'lorem ipsum';
  50. % 4. Do your data need correction for anti-aliasing filter and/or task related time offset?
  51. adjust_time_offset = 0; % 0 = NO (no correction), 1 = YES (correct time offset)
  52. % If your data need correction for time offset, initialize the offset time (in milliseconds)
  53. filter_timeoffset = 0; % anti-aliasing time offset (in milliseconds). 0 = No time offset
  54. stimulus_timeoffset = 0; % stimulus related time offset (in milliseconds). 0 = No time offset
  55. response_timeoffset = 0; % response related time offset (in milliseconds). 0 = No time offset
  56. stimulus_markers = {}; % enter the stimulus makers that need to be adjusted for time offset
  57. respose_markers = {}; % enter the response makers that need to be adjusted for time offset
  58. % 5. Do you want to down sample the data?
  59. down_sample = 1; % 0 = NO (no down sampling), 1 = YES (down sampling)
  60. sampling_rate = 500; % set sampling rate (in Hz), if you want to down sample
  61. % 6. Do you want to delete the outer layer of the channels? (Rationale has been described in MADE manuscript)
  62. % This fnction can also be used to down sample electrodes. For example, if EEG was recorded with 128 channels but you would
  63. % like to analyse only 64 channels, you can assign the list of channnels to be excluded in the 'outerlayer_channel' variable.
  64. delete_outerlayer = 1; % 0 = NO (do not delete outer layer), 1 = YES (delete outerlayer);
  65. % If you want to delete outer layer, make a list of channels to be deleted
  66. outerlayer_channel = {'E17' 'E38' 'E43' 'E44' 'E48' 'E49' 'E113' 'E114' 'E119' 'E120' 'E121' 'E125' 'E126' 'E127' 'E128' 'E56' 'E63' 'E68' 'E73' 'E81' 'E88' 'E94' 'E99' 'E107'}; % list of channels
  67. % recommended list for EGI 128 chanenl net: {'E17' 'E38' 'E43' 'E44' 'E48' 'E49' 'E113' 'E114' 'E119' 'E120' 'E121' 'E125' 'E126' 'E127' 'E128' 'E56' 'E63' 'E68' 'E73' 'E81' 'E88' 'E94' 'E99' 'E107'}
  68. % 7. Initialize the filters
  69. highpass = 0.3; % High-pass frequency
  70. lowpass = 50; % Low-pass frequency. We recommend low-pass filter at/below line noise frequency (see manuscript for detail)
  71. % 8. Are you processing task-related or resting-state EEG data?
  72. %task_eeg = 0; % 0 = resting, 1 = task
  73. task_event_markers_cell = {{'QS','AS','I','W'},{'QS','AS','I','W'},{'1','2','3'},{'DIN1'}}; % enter all the event/condition markers
  74. % 9. Do you want to epoch/segment your data?
  75. %epoch_data = 0; % 0 = NO (do not epoch), 1 = YES (epoch data)
  76. task_epoch_length_cell = {[0 2]; [0 2]; [-0.1 0.5]; [-0.2 0.4]}; % epoch length in second
  77. %rest_epoch_length = 2; % for resting EEG continuous data will be segmented into consecutive epochs of a specified length (here 2 second) by adding dummy events
  78. %overlap_epoch = 1; % 0 = NO (do not create overlapping epoch), 1 = YES (50% overlapping epoch)
  79. %dummy_events ={'rest'}; % enter dummy events name
  80. % 10. Do you want to remove/correct baseline?
  81. %remove_baseline = 0; % 0 = NO (no baseline correction), 1 = YES (baseline correction)
  82. %baseline_window = [xx xx]; % baseline period in milliseconds (MS) [] = entire epoch
  83. % 11. Do you want to remove artifact laden epoch based on voltage threshold?
  84. voltthres_rejection = 1; % 0 = NO, 1 = YES
  85. volt_threshold = [-150 150]; % lower and upper threshold (in ?V)
  86. % 12. Do you want to perform epoch level channel interpolation for artifact laden epoch? (see manuscript for detail)
  87. interp_epoch = 1; % 0 = NO, 1 = YES.
  88. frontal_channels = {'E1', 'E8', 'E14', 'E21', 'E25', 'E32', 'E17'}; % If you set interp_epoch = 1, enter the list of frontal channels to check (see manuscript for detail)
  89. % recommended list for EGI 128 channel net: {'E1', 'E8', 'E14', 'E21', 'E25', 'E32', 'E17'}
  90. %13. Do you want to interpolate the bad channels that were removed from data?
  91. interp_channels = 1; % 0 = NO (Do not interpolate), 1 = YES (interpolate missing channels)
  92. % 14. Do you want to rereference your data?
  93. %rerefer_data = 0; % 0 = NO, 1 = YES
  94. %reref=[]; % Enter electrode name/s or number/s to be used for rereferencing
  95. % For channel name/s enter, reref = {'channel_name', 'channel_name'};
  96. % For channel number/s enter, reref = [channel_number, channel_number];
  97. % For average rereference enter, reref = []; default is average rereference
  98. % 15. Do you want to save interim results?
  99. save_interim_result = 1; % 0 = NO (Do not save) 1 = YES (save interim results)
  100. % 16. How do you want to save your data? .set or .mat
  101. output_format = 1; % 1 = .set (EEGLAB data structure), 2 = .mat (Matlab data structure)
  102. % ********* no need to edit beyond this point for EGI .mff data **********
  103. % ********* for non-.mff data format edit data import function ***********
  104. % ********* below using relevant data import plugin from EEGLAB **********
  105. %% Read files to analyses
  106. datafile_names=dir(rawdata_location);
  107. datafile_names=datafile_names(~ismember({datafile_names.name},{'.', '..', '.DS_Store'}));
  108. datafile_names={datafile_names.name};
  109. ext = '.set'
  110. %% Check whether EEGLAB and all necessary plugins are in Matlab path.
  111. if exist('eeglab','file')==0
  112. error(['Please make sure EEGLAB is on your Matlab path. Please see EEGLAB' ...
  113. 'wiki page for download and instalation instructions']);
  114. end
  115. if strcmp(ext, '.mff')==1
  116. if exist('mff_import', 'file')==0
  117. error(['Please make sure "mffmatlabio" plugin is in EEGLAB plugin folder and on Matlab path.' ...
  118. ' Please see EEGLAB wiki page for download and instalation instructions of plugins.' ...
  119. ' If you are not analysing EGI .mff data, edit the data import function below.']);
  120. end
  121. else
  122. warning('Your data are not EGI .mff files. Make sure you edit data import function before using this script');
  123. end
  124. if exist('pop_firws', 'file')==0
  125. error(['Please make sure "firfilt" plugin is in EEGLAB plugin folder and on Matlab path.' ...
  126. ' Please see EEGLAB wiki page for download and instalation instructions of plugins.']);
  127. end
  128. if exist('channel_properties', 'file')==0
  129. error(['Please make sure "FASTER" plugin is in EEGLAB plugin folder and on Matlab path.' ...
  130. ' Please see EEGLAB wiki page for download and instalation instructions of plugins.']);
  131. end
  132. if exist('ADJUST', 'file')==0
  133. error(['Please make sure you download modified "ADJUST" plugin from GitHub (link is in MADE manuscript)' ...
  134. ' and ADJUST is in EEGLAB plugin folder and on Matlab path.']);
  135. end
  136. %% Create output folders to save data
  137. if save_interim_result ==1
  138. if exist([output_location filesep 'filtered_data'], 'dir') == 0
  139. mkdir([output_location filesep 'filtered_data'])
  140. end
  141. if exist([output_location filesep 'ica_data'], 'dir') == 0
  142. mkdir([output_location filesep 'ica_data'])
  143. end
  144. end
  145. if exist([output_location filesep 'processed_data'], 'dir') == 0
  146. mkdir([output_location filesep 'processed_data'])
  147. end
  148. %% Initialize output variables
  149. reference_used_for_faster=[]; % reference channel used for running faster to identify bad channel/s
  150. faster_bad_channels=[]; % number of bad channel/s identified by faster
  151. ica_preparation_bad_channels=[]; % number of bad channel/s due to channel/s exceeding xx% of artifacted epochs
  152. length_ica_data=[]; % length of data (in second) fed into ICA decomposition
  153. total_ICs=[]; % total independent components (ICs)
  154. ICs_removed=[]; % number of artifacted ICs
  155. total_epochs_before_artifact_rejection=[]; % number of total MMN trials before artifact rejection
  156. total_epochs_after_artifact_rejection=[]; % number of total MMN trials after artifact rejection
  157. total_channels_interpolated=[]; % total_channels_interpolated=faster_bad_channels+ica_preparation_bad_channels
  158. total_standard=[]; % number of MMN standard trials before artifact rejection
  159. total_standard_AS=[]; % number of MMN standard trials during AS before artifact rejection
  160. total_standard_QS=[]; % number of MMN standard trials during QS before artifact rejection
  161. total_standard_W=[]; % number of MMN standard trials during W before artifact rejection
  162. total_standard_I=[]; % number of MMN standard trials during I before artifact rejection
  163. total_standard_after_artifact_rejection=[]; % number of MMN standard trials after artifact rejection
  164. total_standard_AS_after_artifact_rejection=[]; % number of MMN standard trials during AS after artifact rejection
  165. total_standard_QS_after_artifact_rejection=[]; % number of MMN standard trials during QS after artifact rejection
  166. total_standard_W_after_artifact_rejection=[]; % number of MMN standard trials during W after artifact rejection
  167. total_standard_I_after_artifact_rejection=[]; % number of MMN standard trials during I after artifact rejection
  168. total_deviant=[]; % number of MMN deviant trials before artifact rejection
  169. total_deviant_AS=[]; % number of MMN deviant trials during AS before artifact rejection
  170. total_deviant_QS=[]; % number of MMN deviant trials during QS before artifact rejection
  171. total_deviant_W=[]; % number of MMN deviant trials during W before artifact rejection
  172. total_deviant_I=[]; % number of MMN deviant trials during I before artifact rejection
  173. total_deviant_after_artifact_rejection=[]; % number of MMN deviant trials after artifact rejection
  174. total_deviant_AS_after_artifact_rejection=[]; % number of MMN deviant trials during AS after artifact rejection
  175. total_deviant_QS_after_artifact_rejection=[]; % number of MMN deviant trials during QS after artifact rejection
  176. total_deviant_W_after_artifact_rejection=[]; % number of MMN deviant trials during W after artifact rejection
  177. total_deviant_I_after_artifact_rejection=[]; % number of MMN deviant trials during I after artifact rejection
  178. total_novel=[]; % number of MMN novel trials before artifact rejection
  179. total_novel_AS=[]; % number of MMN novel trials during AS before artifact rejection
  180. total_novel_QS=[]; % number of MMN novel trials during QS before artifact rejection
  181. total_novel_W=[]; % number of MMN novel trials during W before artifact rejection
  182. total_novel_I=[]; % number of MMN novel trials during I before artifact rejection
  183. total_novel_after_artifact_rejection=[]; % number of MMN novel trials after artifact rejection
  184. total_novel_AS_after_artifact_rejection=[]; % number of MMN novel trials during AS after artifact rejection
  185. total_novel_QS_after_artifact_rejection=[]; % number of MMN novel trials during QS after artifact rejection
  186. total_novel_W_after_artifact_rejection=[]; % number of MMN novel trials during W after artifact rejection
  187. total_novel_I_after_artifact_rejection=[]; % number of MMN novel trials during I after artifact rejection
  188. total_AS=[]; % number of Resting epochs during AS before artifact rejection
  189. total_QS=[]; % number of Resting epochs during QS before artifact rejection
  190. total_W=[]; % number of Resting epochs during W before artifact rejection
  191. total_I=[]; % number of Resting epochs during I before artifact rejection
  192. total_AS_after_artifact_rejection=[]; % number of Resting epochs during AS after artifact rejection
  193. total_QS_after_artifact_rejection=[]; % number of Resting epochs during QS after artifact rejection
  194. total_W_after_artifact_rejection=[]; % number of Resting epochs during W after artifact rejection
  195. total_I_after_artifact_rejection=[]; % number of Resting epochs during I after artifact rejection
  196. total_vep=[]; % number of vep total before artifact rejection
  197. total_AS_vep=[]; % number of vep AS trials before artifaction rejection
  198. total_QS_vep=[]; % number of vep QS trials before artifaction rejection
  199. total_I_vep=[]; % number of vep I trials before artifaction rejection
  200. total_W_vep=[]; % number of vep W trials before artifaction rejection\\
  201. total_vep_after_artifact_rejection=[]; % number of total vep after artifact rejection
  202. total_AS_vep_artifact_rejection=[]; % number of vep AS trials after artifaction rejection
  203. total_QS_vep_artifact_rejection=[]; % number of vep QS trials after artifaction rejection
  204. total_I_vep_artifact_rejection=[]; % number of vep I trials after artifaction rejection
  205. total_W_vep_artifact_rejection=[]; % number of vep W trials after artifaction rejection
  206. %% Loop over all data files
  207. for s=1:length(datafile_names)
  208. fprintf('\n\n\n*** Processing subject %d (%s) ***\n\n\n', s, datafile_names{s});
  209. ALLEEG=[];
  210. % get the current subject folder name and go into that folder
  211. subject = datafile_names{s};
  212. subject_folder = [rawdata_location subject filesep subject];
  213. cd (subject_folder);
  214. %Make a list of the EEG files in that subject's folder
  215. sub_file_list=dir('*.set');
  216. sub_file_list={sub_file_list.name};
  217. %grab the current task file name
  218. for c=1:length(sub_file_list) %1:NumberOfTasks;
  219. %Filename example: F1572_MMN_1mo_20201104_045409.mff
  220. currentTaskFilename = sub_file_list{c};
  221. subfields = strsplit(currentTaskFilename,'_');
  222. sCon = subfields{2};
  223. %check which task this file is
  224. if strcmp(sCon,'arm')
  225. task=1;
  226. elseif strcmp(sCon,'sResting')
  227. task=2;
  228. elseif strcmp(sCon, 'MMN')
  229. task=3;
  230. elseif strcmp(sCon, 'VEP')
  231. task=4;
  232. else
  233. task=0;
  234. error('Current task name does not match either task!')
  235. end
  236. %% STEP 1: Import EGI data file and relevant information
  237. cd (subject_folder);
  238. EEG = pop_loadset('filename', sub_file_list{c}, 'filepath', subject_folder);
  239. EEG = eeg_checkset(EEG);
  240. % Edit this data import function and use appropriate plugin from EEGLAB
  241. % for non-.mff data. For example, to import biosemi data, use biosig plugin.
  242. % The example codes for 64 channels biosemi data:
  243. % EEG = pop_biosig([rawdata_location, filesep, datafile_names{subject}]);
  244. % EEG = eeg_checkset(EEG);
  245. % EEG = pop_select( EEG,'nochannel', 65:72); % delete redundant channels
  246. %% STEP 2: Import channel locations
  247. EEG=pop_chanedit(EEG, 'load',{channel_locations 'filetype' 'autodetect'});
  248. EEG = eeg_checkset( EEG );
  249. % Check whether the channel locations were properly imported. The EEG signals and channel numbers should be same.
  250. if size(EEG.data, 1) ~= length(EEG.chanlocs)
  251. error('The size of the data does not match with channel numbers.');
  252. end
  253. %% STEP 3: Adjust anti-aliasing and task related time offset
  254. if adjust_time_offset==1
  255. % adjust anti-aliasing filter time offset
  256. if filter_timeoffset~=0
  257. for aafto=1:length(EEG.event)
  258. EEG.event(aafto).latency=EEG.event(aafto).latency+(filter_timeoffset/1000)*EEG.srate;
  259. end
  260. end
  261. % adjust stimulus time offset
  262. if stimulus_timeoffset~=0
  263. for sto=1:length(EEG.event)
  264. for sm=1:length(stimulus_markers)
  265. if strcmp(EEG.event(sto).type, stimulus_markers{sm})
  266. EEG.event(sto).latency=EEG.event(sto).latency+(stimulus_timeoffset/1000)*EEG.srate;
  267. end
  268. end
  269. end
  270. end
  271. % adjust response time offset
  272. if response_timeoffset~=0
  273. for rto=1:length(EEG.event)
  274. for rm=1:length(response_markers)
  275. if strcmp(EEG.event(rto).type, response_markers{rm})
  276. EEG.event(rto).latency=EEG.event(rto).latency-(response_timeoffset/1000)*EEG.srate;
  277. end
  278. end
  279. end
  280. end
  281. end
  282. %% STEP 4: Change sampling rate
  283. if down_sample==1
  284. if floor(sampling_rate) > EEG.srate
  285. error ('Sampling rate cannot be higher than recorded sampling rate');
  286. elseif floor(sampling_rate) ~= EEG.srate
  287. EEG = pop_resample( EEG, sampling_rate);
  288. EEG = eeg_checkset( EEG );
  289. end
  290. end
  291. %% STEP 5: Delete outer layer of channels
  292. chans_labels=cell(1,EEG.nbchan);
  293. for i=1:EEG.nbchan
  294. chans_labels{i}= EEG.chanlocs(i).labels;
  295. end
  296. [chans,chansidx] = ismember(outerlayer_channel, chans_labels);
  297. outerlayer_channel_idx = chansidx(chansidx ~= 0);
  298. if delete_outerlayer==1
  299. if isempty(outerlayer_channel_idx)==1
  300. error(['None of the outer layer channels present in channel locations of data.'...
  301. ' Make sure outer layer channels are present in channel labels of data (EEG.chanlocs.labels).']);
  302. else
  303. EEG = pop_select( EEG,'nochannel', outerlayer_channel_idx);
  304. EEG = eeg_checkset( EEG );
  305. end
  306. end
  307. %% STEP 5.5: Lable the task events
  308. task_epoch_length = task_epoch_length_cell{task};
  309. % Creating task, type, and latency lables
  310. if task == 1 % sResting in arm
  311. EEG = pop_editeventfield( EEG, 'indices', strcat('1:', int2str(length(EEG.event))), 'Task',task); % Creating an event field "Task" and put 1 for Resting
  312. EEG = pop_editeventfield( EEG, 'indices', strcat('1:', int2str(length(EEG.event))), 'type',NaN); % Creating an event filed "type" for inserting sleep state flags
  313. EEG = pop_editeventfield( EEG, 'indices', strcat('1:', int2str(length(EEG.event))), 'latency',1); % Creating an event filed "latency" for inserting sleep state flags
  314. EEG = eeg_checkset( EEG );
  315. % insert sleep state flags to resting
  316. EEG = pop_importevent( EEG, 'event',[event_location '/' subject '_sResting_arm.csv'] ,'fields',{'latency' 'type'},'timeunit',1,'align',0);
  317. EEG = eeg_checkset( EEG );
  318. % Add sleep state flags every second
  319. current_cond = {'QS','AS', 'W','I'};
  320. for cc=1:length(current_cond)
  321. startEventNums=[]; endEventNums=[]; startLats=[]; endLats=[];
  322. startEventNums = find(strcmp({EEG.event.type},[current_cond{cc} 'start']));
  323. endEventNums = find(strcmp({EEG.event.type},[current_cond{cc} 'end']));
  324. if isempty(startEventNums)
  325. % condition not present... can't epoch
  326. else
  327. for ll=1:length(startEventNums)
  328. startLats(ll) = (EEG.event(startEventNums(ll)).latency)/EEG.srate; % latency of the start point of each state in seconds
  329. end
  330. for ll=1:length(endEventNums)
  331. endLats(ll) = (EEG.event(endEventNums(ll)).latency)/EEG.srate; % latency of the end point of each state in seconds
  332. end
  333. for b=1:length(startEventNums)
  334. conditionleg = endLats(b)- startLats(b);
  335. for t=1:floor(conditionleg)-task_epoch_length(end)+2 % +2 because the second between two states is included in the prvious state
  336. latency = startLats(b) + t-1;
  337. EEG = pop_editeventvals(EEG,'insert',{1 [] [] []},'changefield',{1 'type' current_cond{cc}}, 'changefield',{1 'latency' latency},'changefield',{1,'Task',task});
  338. end
  339. EEG = eeg_checkset(EEG);
  340. end
  341. end
  342. end
  343. elseif task == 2 % sResting in bassinet
  344. EEG = pop_editeventfield( EEG, 'indices', strcat('1:', int2str(length(EEG.event))), 'Task',task); % Creating an event field "Task" and put 1 for Resting
  345. EEG = pop_editeventfield( EEG, 'indices', strcat('1:', int2str(length(EEG.event))), 'type',NaN); % Creating an event filed "type" for inserting sleep state flags
  346. EEG = pop_editeventfield( EEG, 'indices', strcat('1:', int2str(length(EEG.event))), 'latency',1); % Creating an event filed "latency" for inserting sleep state flags
  347. EEG = eeg_checkset( EEG );
  348. % insert sleep state flags to resting
  349. EEG = pop_importevent( EEG, 'event',[event_location '/' subject '_sResting.csv'] ,'fields',{'latency' 'type'},'timeunit',1,'align',0);
  350. EEG = eeg_checkset( EEG );
  351. % Add sleep state flags every second
  352. current_cond = {'QS','AS', 'W','I'};
  353. for cc=1:length(current_cond)
  354. startEventNums=[]; endEventNums=[]; startLats=[]; endLats=[];
  355. startEventNums = find(strcmp({EEG.event.type},[current_cond{cc} 'start']));
  356. endEventNums = find(strcmp({EEG.event.type},[current_cond{cc} 'end']));
  357. if isempty(startEventNums)
  358. % condition not present... can't epoch
  359. else
  360. for ll=1:length(startEventNums)
  361. startLats(ll) = (EEG.event(startEventNums(ll)).latency)/EEG.srate; % latency of the start point of each state in seconds
  362. end
  363. for ll=1:length(endEventNums)
  364. endLats(ll) = (EEG.event(endEventNums(ll)).latency)/EEG.srate; % latency of the end point of each state in seconds
  365. end
  366. for b=1:length(startEventNums)
  367. conditionleg = endLats(b)- startLats(b);
  368. for t=1:floor(conditionleg)-task_epoch_length(end)+2 % +2 because the second between two states is included in the prvious state
  369. latency = startLats(b) + t-1;
  370. EEG = pop_editeventvals(EEG,'insert',{1 [] [] []},'changefield',{1 'type' current_cond{cc}}, 'changefield',{1 'latency' latency},'changefield',{1,'Task',task});
  371. end
  372. EEG = eeg_checkset(EEG);
  373. end
  374. end
  375. end
  376. elseif task == 3 % MMN
  377. % Creating a task Label
  378. EEG = pop_editeventfield( EEG, 'indices', strcat('1:', int2str(length(EEG.event))), 'Task',task); % Creating an event field "Task" and put 2 for MMN
  379. EEG = pop_editeventfield( EEG, 'indices', strcat('1:', int2str(length(EEG.event))), 'State','I'); % Creating an event filed "State" for mark the sleep state for each stimulus
  380. EEG = eeg_checkset( EEG );
  381. % insert sleep state flags to MMN
  382. EEG = pop_importevent( EEG, 'event',[event_location '/' subject '_MMN.csv'] ,'fields',{'latency' 'type'},'timeunit',1,'align',0); % 0 in the sleep state file corresponds to "boundary" in EEG
  383. EEG = eeg_checkset( EEG );
  384. % Mark sleep stat flag for MMN ERP markers
  385. current_cond = {'QS','AS', 'W','I'};
  386. for cc=1:length(current_cond)
  387. startEventNums=[]; endEventNums=[]; DINSEventNums=[]; startLats=[]; endLats=[]; DIN2Lats=[];
  388. startEventNums = find(strcmp({EEG.event.type},[current_cond{cc} 'start']));
  389. endEventNums = find(strcmp({EEG.event.type},[current_cond{cc} 'end']));
  390. DINSEventNums=find(strcmp({EEG.event.type},'DIN2'));
  391. if isempty(startEventNums)
  392. % condition not present... can't epoch
  393. else
  394. for ll=1:length(startEventNums)
  395. startLats(ll) = (EEG.event(startEventNums(ll)).latency)/EEG.srate;
  396. end
  397. for ll=1:length(endEventNums)
  398. endLats(ll) = (EEG.event(endEventNums(ll)).latency)/EEG.srate;
  399. end
  400. for ll=1:length(DINSEventNums)
  401. DIN2Lats(ll)= (EEG.event(DINSEventNums(ll)).latency)/EEG.srate;
  402. end
  403. for b=1:length(startEventNums)
  404. for dl = 1: length(DINSEventNums)
  405. if DIN2Lats(dl) >= startLats(b)-task_epoch_length(1) && DIN2Lats(dl)<= endLats(b)+1-task_epoch_length(2) % select trials that fall into just one sleep state
  406. EEG.event(DINSEventNums(dl)).State=current_cond{cc};
  407. end
  408. end
  409. end
  410. end
  411. end
  412. % lable MMN
  413. cd(scripts_location);
  414. Labeling_MMN_multi();
  415. % remove data after last trps
  416. trsp = find(strcmp({EEG.event.type},'TRSP'));
  417. EEG = eeg_eegrej( EEG, [(EEG.event(trsp(end)).latency+(1.5*EEG.srate)) EEG.pnts] );
  418. EEG = eeg_checkset( EEG );
  419. cd(rawdata_location);
  420. elseif task == 4 % VEP
  421. % Creating a task Label
  422. EEG = pop_editeventfield( EEG, 'indices', strcat('1:', int2str(length(EEG.event))), 'Task',task); % Creating an event field "Task" and put 2 for MMN
  423. EEG = pop_editeventfield( EEG, 'indices', strcat('1:', int2str(length(EEG.event))), 'State','I'); % Creating an event filed "State" for mark the sleep state for each stimulus
  424. EEG = eeg_checkset( EEG );
  425. % insert sleep state flags to VEP
  426. EEG = pop_importevent( EEG, 'event',[event_location '/' subject '_VEP.csv'] ,'fields',{'latency' 'type'},'timeunit',1,'align',0);
  427. EEG = eeg_checkset( EEG );
  428. % Mark sleep stat flag for VEP ERP markers
  429. current_cond = {'QS','AS', 'W','I'};
  430. for cc=1:length(current_cond)
  431. startEventNums=[]; endEventNums=[]; DINSEventNums=[]; startLats=[]; endLats=[]; DIN1Lats=[];
  432. startEventNums = find(strcmp({EEG.event.type},[current_cond{cc} 'start']));
  433. endEventNums = find(strcmp({EEG.event.type},[current_cond{cc} 'end']));
  434. DINSEventNums=find(strcmp({EEG.event.type},'DIN1'));
  435. if isempty(startEventNums)
  436. % condition not present... can't epoch
  437. else
  438. for ll=1:length(startEventNums)
  439. startLats(ll) = (EEG.event(startEventNums(ll)).latency)/EEG.srate;
  440. end
  441. for ll=1:length(endEventNums)
  442. endLats(ll) = (EEG.event(endEventNums(ll)).latency)/EEG.srate;
  443. end
  444. for ll=1:length(DINSEventNums)
  445. DIN1Lats(ll)= (EEG.event(DINSEventNums(ll)).latency)/EEG.srate;
  446. end
  447. for b=1:length(startEventNums)
  448. for dl = 1: length(DINSEventNums)
  449. if DIN1Lats(dl) >= startLats(b)-task_epoch_length(1) && DIN1Lats(dl)<= endLats(b)+1-task_epoch_length(2) % select trials that fall into just one sleep state
  450. EEG.event(DINSEventNums(dl)).State=current_cond{cc};
  451. end
  452. end
  453. end
  454. end
  455. end
  456. end
  457. %% STEP 6: Filter data
  458. % Calculate filter order using the formula: m = dF / (df / fs), where m = filter order,
  459. % df = transition band width, dF = normalized transition width, fs = sampling rate
  460. % dF is specific for the window type. Hamming window dF = 3.3
  461. high_transband = highpass; % high pass transition band
  462. low_transband = 10; % low pass transition band
  463. hp_fl_order = 3.3 / (high_transband / EEG.srate);
  464. lp_fl_order = 3.3 / (low_transband / EEG.srate);
  465. % Round filter order to next higher even integer. Filter order is always even integer.
  466. if mod(floor(hp_fl_order),2) == 0
  467. hp_fl_order=floor(hp_fl_order);
  468. elseif mod(floor(hp_fl_order),2) == 1
  469. hp_fl_order=floor(hp_fl_order)+1;
  470. end
  471. if mod(floor(lp_fl_order),2) == 0
  472. lp_fl_order=floor(lp_fl_order)+2;
  473. elseif mod(floor(lp_fl_order),2) == 1
  474. lp_fl_order=floor(lp_fl_order)+1;
  475. end
  476. % Calculate cutoff frequency
  477. high_cutoff = highpass/2;
  478. low_cutoff = lowpass + (low_transband/2);
  479. % Performing high pass filtering
  480. EEG = eeg_checkset( EEG );
  481. EEG = pop_firws(EEG, 'fcutoff', high_cutoff, 'ftype', 'highpass', 'wtype', 'hamming', 'forder', hp_fl_order, 'minphase', 0);
  482. EEG = eeg_checkset( EEG );
  483. % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
  484. % pop_firws() - filter window type hamming ('wtype', 'hamming')
  485. % pop_firws() - applying zero-phase (non-causal) filter ('minphase', 0)
  486. % Performing low pass filtering
  487. EEG = eeg_checkset( EEG );
  488. EEG = pop_firws(EEG, 'fcutoff', low_cutoff, 'ftype', 'lowpass', 'wtype', 'hamming', 'forder', lp_fl_order, 'minphase', 0);
  489. EEG = eeg_checkset( EEG );
  490. % pop_firws() - transition band width: 10 Hz
  491. % pop_firws() - filter window type hamming ('wtype', 'hamming')
  492. % pop_firws() - applying zero-phase (non-causal) filter ('minphase', 0)
  493. % save current task in ALLEEG structure (we'll use this to merge)
  494. [ALLEEG, EEG, CURRENTSET] = eeg_store(ALLEEG, EEG, task);
  495. EEG_idx(task)=1; % save out index of current task
  496. end % end loop through subject task files
  497. %% STEP 6.5 Merge the datasets
  498. % EEG = []; % Saving over EEG
  499. EEG2Merge = find(EEG_idx==1);
  500. if length(EEG2Merge) >=2
  501. EEG = pop_mergeset(ALLEEG,EEG2Merge);
  502. EEG = eeg_checkset(EEG);
  503. EEG = pop_editset(EEG, 'setname', [datafile_names{s} '_Merged']);
  504. end
  505. %% STEP 7: Run faster to find bad channels
  506. % First check whether reference channel (i.e. zeroed channels) is present in data
  507. % reference channel is needed to run faster
  508. ref_chan=[]; FASTbadChans=[]; all_chan_bad_FAST=0;
  509. ref_chan=find(any(EEG.data, 2)==0);
  510. if numel(ref_chan)>1
  511. error(['There are more than 1 zeroed channel (i.e. zero value throughout recording) in data.'...
  512. ' Only reference channel should be zeroed channel. Delete the zeroed channel/s which is not reference channel.']);
  513. elseif numel(ref_chan)==1
  514. list_properties = channel_properties(EEG, 1:EEG.nbchan, ref_chan); % run faster
  515. FASTbadIdx=min_z(list_properties);
  516. FASTbadChans=find(FASTbadIdx==1);
  517. FASTbadChans=FASTbadChans(FASTbadChans~=ref_chan);
  518. reference_used_for_faster{s}={EEG.chanlocs(ref_chan).labels};
  519. EEG = pop_select( EEG,'nochannel', ref_chan);
  520. EEG = eeg_checkset(EEG);
  521. channels_analysed=EEG.chanlocs; % keep full channel locations to use later for interpolation of bad channels
  522. elseif numel(ref_chan)==0
  523. warning('Reference channel is not present in data. Cz channel will be used as reference channel.');
  524. ref_chan=find(strcmp({EEG.chanlocs.labels}, 'Cz')); % find Cz channel index
  525. EEG_copy=[];
  526. EEG_copy=EEG; % make a copy of the dataset
  527. EEG_copy = pop_reref( EEG_copy, ref_chan,'keepref','on'); % rerefer to Cz in copied dataset
  528. EEG_copy = eeg_checkset(EEG_copy);
  529. list_properties = channel_properties(EEG_copy, 1:EEG_copy.nbchan, ref_chan); % run faster on copied dataset
  530. FASTbadIdx=min_z(list_properties);
  531. FASTbadChans=find(FASTbadIdx==1);
  532. channels_analysed=EEG.chanlocs;
  533. reference_used_for_faster{s}={EEG.chanlocs(ref_chan).labels};
  534. end
  535. % If FASTER identifies all channels as bad channels, save the dataset
  536. % at this stage and ignore the remaining of the preprocessing.
  537. if numel(FASTbadChans)==EEG.nbchan || numel(FASTbadChans)+1==EEG.nbchan
  538. all_chan_bad_FAST=1;
  539. warning(['No usable data for datafile', datafile_names{s}]);
  540. if output_format==1
  541. EEG = eeg_checkset(EEG);
  542. EEG = pop_editset(EEG, 'setname', [datafile_names{s} '_Merged_no_usable_data_all_bad_channels']);
  543. EEG = pop_saveset(EEG, 'filename', [datafile_names{s} '_Merged_no_usable_data_all_bad_channels.set'],'filepath', [output_location filesep 'filtered_data' filesep]); % save .set format
  544. elseif output_format==2
  545. save([[output_location filesep 'processed_data' filesep] [datafile_names{s} '_Merged_no_usable_data_all_bad_channels.mat']], 'EEG'); % save .mat format
  546. end
  547. else
  548. % Reject channels that are bad as identified by Faster
  549. EEG = pop_select( EEG,'nochannel', FASTbadChans);
  550. EEG = eeg_checkset(EEG);
  551. end
  552. if numel(FASTbadChans)==0
  553. faster_bad_channels{s}='0';
  554. else
  555. faster_bad_channels{s}=num2str(FASTbadChans');
  556. end
  557. if all_chan_bad_FAST==1
  558. faster_bad_channels{s}='0';
  559. ica_preparation_bad_channels{s}='0';
  560. length_ica_data(s)=0;
  561. total_ICs(s)=0;
  562. ICs_removed{s}='0';
  563. total_epochs_before_artifact_rejection(s)=0;
  564. total_epochs_after_artifact_rejection(s)=0;
  565. total_channels_interpolated(s)=0;
  566. continue % ignore rest of the processing and go to next subject folder
  567. end
  568. %% Save data after running filter and FASTER function, if saving interim results was preferred
  569. if save_interim_result ==1
  570. if output_format==1
  571. EEG = eeg_checkset( EEG );
  572. EEG = pop_editset(EEG, 'setname', strrep(datafile_names{s}, ext, '_filtered_data'));
  573. EEG = pop_saveset( EEG,'filename',strrep(datafile_names{s}, ext, '_filtered_data.set'),'filepath', [output_location filesep 'filtered_data' filesep]); % save .set format
  574. elseif output_format==2
  575. save([[output_location filesep 'filtered_data' filesep ] strrep(datafile_names{s}, ext, '_filtered_data.mat')], 'EEG'); % save .mat format
  576. end
  577. end
  578. %% STEP 8: Prepare data for ICA
  579. EEG_copy=[];
  580. EEG_copy=EEG; % make a copy of the dataset
  581. EEG_copy = eeg_checkset(EEG_copy);
  582. % Perform 1Hz high pass filter on copied dataset
  583. transband = 1;
  584. fl_cutoff = transband/2;
  585. fl_order = 3.3 / (transband / EEG.srate);
  586. if mod(floor(fl_order),2) == 0
  587. fl_order=floor(fl_order);
  588. elseif mod(floor(fl_order),2) == 1
  589. fl_order=floor(fl_order)+1;
  590. end
  591. EEG_copy = pop_firws(EEG_copy, 'fcutoff', fl_cutoff, 'ftype', 'highpass', 'wtype', 'hamming', 'forder', fl_order, 'minphase', 0);
  592. EEG_copy = eeg_checkset(EEG_copy);
  593. % Create 1 second epoch
  594. EEG_copy=eeg_regepochs(EEG_copy,'recurrence', 1, 'limits',[0 1], 'rmbase', [NaN], 'eventtype', '999'); % insert temporary marker 1 second apart and create epochs
  595. EEG_copy = eeg_checkset(EEG_copy);
  596. % Find bad epochs and delete them from dataset
  597. vol_thrs = [-1000 1000]; % [lower upper] threshold limit(s) in mV.
  598. emg_thrs = [-100 30]; % [lower upper] threshold limit(s) in dB.
  599. emg_freqs_limit = [20 40]; % [lower upper] frequency limit(s) in Hz.
  600. % Find channel/s with xx% of artifacted 1-second epochs and delete them
  601. chanCounter = 1; ica_prep_badChans = [];
  602. numEpochs =EEG_copy.trials; % find the number of epochs
  603. all_bad_channels=0;
  604. for ch=1:EEG_copy.nbchan
  605. % Find artifaceted epochs by detecting outlier voltage
  606. EEG_copy = pop_eegthresh(EEG_copy,1, ch, vol_thrs(1), vol_thrs(2), EEG_copy.xmin, EEG_copy.xmax, 0, 0);
  607. EEG_copy = eeg_checkset( EEG_copy );
  608. % 1 : data type (1: electrode, 0: component)
  609. % 0 : display with previously marked rejections? (0: no, 1: yes)
  610. % 0 : reject marked trials? (0: no (but store the marks), 1:yes)
  611. % Find artifaceted epochs by using thresholding of frequencies in the data.
  612. % this method mainly rejects muscle movement (EMG) artifacts
  613. EEG_copy = pop_rejspec( EEG_copy, 1,'elecrange',ch ,'method','fft','threshold', emg_thrs, 'freqlimits', emg_freqs_limit, 'eegplotplotallrej', 0, 'eegplotreject', 0);
  614. % method : method to compute spectrum (fft)
  615. % threshold : [lower upper] threshold limit(s) in dB.
  616. % freqlimits : [lower upper] frequency limit(s) in Hz.
  617. % eegplotplotallrej : 0 = Do not superpose rejection marks on previous marks stored in the dataset.
  618. % eegplotreject : 0 = Do not reject marked trials (but store the marks).
  619. % Find number of artifacted epochs
  620. EEG_copy = eeg_checkset( EEG_copy );
  621. EEG_copy = eeg_rejsuperpose( EEG_copy, 1, 1, 1, 1, 1, 1, 1, 1);
  622. artifacted_epochs=EEG_copy.reject.rejglobal;
  623. % Find bad channel / channel with more than 20% artifacted epochs
  624. if sum(artifacted_epochs) > (numEpochs*20/100)
  625. ica_prep_badChans(chanCounter) = ch;
  626. chanCounter=chanCounter+1;
  627. end
  628. end
  629. % If all channels are bad, save the dataset at this stage and ignore the remaining of the preprocessing.
  630. if numel(ica_prep_badChans)==EEG.nbchan || numel(ica_prep_badChans)+1==EEG.nbchan
  631. all_bad_channels=1;
  632. warning(['No usable data for datafile', [datafile_names{s} '_Merged']]);
  633. if output_format==1
  634. EEG = eeg_checkset(EEG);
  635. EEG = pop_editset(EEG, 'setname', [datafile_names{s} '_Merged_no_usable_data_all_bad_channels']);
  636. EEG = pop_saveset(EEG, 'filename', [datafile_names{s} '_Merged_no_usable_data_all_bad_channels.set'],'filepath', [output_location filesep 'processed_data' filesep]); % save .set format
  637. elseif output_format==2
  638. save([[output_location filesep 'processed_data' filesep] [datafile_names{s} '_Merged_no_usable_data_all_bad_channels.mat']], 'EEG'); % save .mat format
  639. end
  640. else
  641. % Reject bad channel - channel with more than xx% artifacted epochs
  642. EEG_copy = pop_select( EEG_copy,'nochannel', ica_prep_badChans);
  643. EEG_copy = eeg_checkset(EEG_copy);
  644. end
  645. if numel(ica_prep_badChans)==0
  646. ica_preparation_bad_channels{s}='0';
  647. else
  648. ica_preparation_bad_channels{s}=num2str(ica_prep_badChans);
  649. end
  650. if all_bad_channels == 1
  651. length_ica_data(s)=0;
  652. total_ICs(s)=0;
  653. ICs_removed{s}='0';
  654. total_epochs_before_artifact_rejection(s)=0;
  655. total_epochs_after_artifact_rejection(s)=0;
  656. total_channels_interpolated(s)=0;
  657. continue % ignore rest of the processing and go to next datafile
  658. end
  659. % Find the artifacted epochs across all channels and reject them before doing ICA.
  660. EEG_copy = pop_eegthresh(EEG_copy,1, 1:EEG_copy.nbchan, vol_thrs(1), vol_thrs(2), EEG_copy.xmin, EEG_copy.xmax,0,0);
  661. EEG_copy = eeg_checkset(EEG_copy);
  662. % 1 : data type (1: electrode, 0: component)
  663. % 0 : display with previously marked rejections? (0: no, 1: yes)
  664. % 0 : reject marked trials? (0: no (but store the marks), 1:yes)
  665. % Find artifaceted epochs by using power threshold in 20-40Hz frequency band.
  666. % This method mainly rejects muscle movement (EMG) artifacts.
  667. EEG_copy = pop_rejspec(EEG_copy, 1,'elecrange', 1:EEG_copy.nbchan, 'method', 'fft', 'threshold', emg_thrs ,'freqlimits', emg_freqs_limit, 'eegplotplotallrej', 0, 'eegplotreject', 0);
  668. % method : method to compute spectrum (fft)
  669. % threshold : [lower upper] threshold limit(s) in dB.
  670. % freqlimits : [lower upper] frequency limit(s) in Hz.
  671. % eegplotplotallrej : 0 = Do not superpose rejection marks on previous marks stored in the dataset.
  672. % eegplotreject : 0 = Do not reject marked trials (but store the marks).
  673. % Find the number of artifacted epochs and reject them
  674. EEG_copy = eeg_checkset(EEG_copy);
  675. EEG_copy = eeg_rejsuperpose(EEG_copy, 1, 1, 1, 1, 1, 1, 1, 1);
  676. reject_artifacted_epochs=EEG_copy.reject.rejglobal;
  677. EEG_copy = pop_rejepoch(EEG_copy, reject_artifacted_epochs, 0);
  678. %% STEP 9: Run ICA
  679. length_ica_data(s)=EEG_copy.trials; % length of data (in second) fed into ICA
  680. EEG_copy = eeg_checkset(EEG_copy);
  681. EEG_copy = pop_runica(EEG_copy, 'icatype', 'runica', 'extended', 1, 'stop', 1E-7, 'interupt','off');
  682. % Find the ICA weights that would be transferred to the original dataset
  683. ICA_WINV=EEG_copy.icawinv;
  684. ICA_SPHERE=EEG_copy.icasphere;
  685. ICA_WEIGHTS=EEG_copy.icaweights;
  686. ICA_CHANSIND=EEG_copy.icachansind;
  687. % If channels were removed from copied dataset during preparation of ica, then remove
  688. % those channels from original dataset as well before transferring ica weights.
  689. EEG = eeg_checkset(EEG);
  690. EEG = pop_select(EEG,'nochannel', ica_prep_badChans);
  691. % Transfer the ICA weights of the copied dataset to the original dataset
  692. EEG.icawinv=ICA_WINV;
  693. EEG.icasphere=ICA_SPHERE;
  694. EEG.icaweights=ICA_WEIGHTS;
  695. EEG.icachansind=ICA_CHANSIND;
  696. EEG = eeg_checkset(EEG);
  697. %% STEP 10: Run adjust to find artifacted ICA components
  698. badICs=[]; EEG_copy =[];
  699. EEG_copy = EEG;
  700. EEG_copy =eeg_regepochs(EEG_copy,'recurrence', 1, 'limits',[0 1], 'rmbase', [NaN], 'eventtype', '999'); % insert temporary marker 1 second apart and create epochs
  701. EEG_copy = eeg_checkset(EEG_copy);
  702. if save_interim_result==1
  703. badICs = adjusted_ADJUST(EEG_copy, [[output_location filesep 'ica_data' filesep] [datafile_names{s} '_Merged_ica_adjust_report']]);
  704. else
  705. badICs = adjusted_ADJUST(EEG_copy, [[output_location filesep 'processed_data' filesep] [datafile_names{s} '_Merged_ica_adjust_report']]);
  706. end
  707. close all;
  708. % Mark the bad ICs found by ADJUST
  709. for ic=1:length(badICs)
  710. EEG.reject.gcompreject(1, badICs(ic))=1;
  711. EEG = eeg_checkset(EEG);
  712. end
  713. total_ICs(s)=size(EEG.icasphere, 1);
  714. if numel(badICs)==0
  715. ICs_removed{s}='0';
  716. else
  717. ICs_removed{s}=num2str(double(badICs));
  718. end
  719. %% Save dataset after ICA, if saving interim results was preferred
  720. if save_interim_result==1
  721. if output_format==1
  722. EEG = eeg_checkset(EEG);
  723. EEG = pop_editset(EEG, 'setname', [datafile_names{s} '_Merged_ica_data']);
  724. EEG = pop_saveset(EEG, 'filename', [datafile_names{s} '_Merged_ica_data.set'],'filepath', [output_location filesep 'ica_data' filesep]); % save .set format
  725. elseif output_format==2
  726. save([[output_location filesep 'ica_data' filesep] [datafile_names{s} '_Merged_ica_data.mat']], 'EEG'); % save .mat format
  727. end
  728. end
  729. %% STEP 11: Remove artifacted ICA components from data
  730. all_bad_ICs=0;
  731. ICs2remove=find(EEG.reject.gcompreject); % find ICs to remove
  732. % If all ICs and bad, save data at this stage and ignore rest of the preprocessing for this subject.
  733. if numel(ICs2remove)==total_ICs(s)
  734. all_bad_ICs=1;
  735. warning(['No usable data for datafile', [datafile_names{s} '_Merged']]);
  736. if output_format==1
  737. EEG = eeg_checkset(EEG);
  738. EEG = pop_editset(EEG, 'setname', [datafile_names{s} '_Merged_no_usable_data_all_bad_ICs']);
  739. EEG = pop_saveset(EEG, 'filename', [datafile_names{s} '_Merged_no_usable_data_all_bad_ICs.set'],'filepath', [output_location filesep 'processed_data' filesep]); % save .set format
  740. elseif output_format==2
  741. save([[output_location filesep 'processed_data' filesep] [datafile_names{s} '_Merged_no_usable_data_all_bad_ICs.mat']], 'EEG'); % save .mat format
  742. end
  743. else
  744. EEG = eeg_checkset( EEG );
  745. EEG = pop_subcomp( EEG, ICs2remove, 0); % remove ICs from dataset
  746. end
  747. if all_bad_ICs==1
  748. total_epochs_before_artifact_rejection(s)=0;
  749. total_epochs_after_artifact_rejection(s)=0;
  750. total_channels_interpolated(s)=0;
  751. continue % ignore rest of the processing and go to next datafile
  752. end
  753. %% Keep only the time locking event markers and delete all other event markers
  754. EEG = eeg_checkset(EEG);
  755. EEG3 = pop_selectevent(EEG,'type',{'QS' 'AS' 'I' 'W' '1' '2' '3' 'DIN1'},'deleteevents','on');
  756. %loop through the tasks in the merged dataset
  757. for task = unique([EEG3.event.Task])
  758. EEG=[];
  759. task_event_markers = task_event_markers_cell{task};
  760. task_epoch_length = task_epoch_length_cell{task};
  761. if task == 1 % sResting in arms
  762. % Select events specific to sResting
  763. EEG = pop_selectevent(EEG3, 'Task', [task],'deleteevents','on'); % just delete other flags/markers, but will NOT delte the EEG data
  764. %% STEP 12: Segment data into fixed length epochs
  765. EEG = eeg_checkset(EEG);
  766. EEG = pop_epoch(EEG, task_event_markers, task_epoch_length, 'epochinfo', 'yes');
  767. EEG = pop_selectevent( EEG, 'latency','-.1 <= .1','deleteevents','on');
  768. total_epochs_before_artifact_rejection(s)=EEG.trials;
  769. total_AS(s) = length(find(strcmp({EEG.event.type},'AS')));
  770. total_QS(s) = length(find(strcmp({EEG.event.type}, 'QS')));
  771. total_W(s) = length(find(strcmp({EEG.event.type}, 'W')));
  772. total_I(s) = length(find(strcmp({EEG.event.type}, 'I')));
  773. %% STEP 13: Remove baseline
  774. %if remove_baseline==1
  775. %EEG = eeg_checkset( EEG );
  776. %EEG = pop_rmbase( EEG, baseline_window);
  777. %end
  778. %% STEP 14: Artifact rejection
  779. all_bad_epochs=0;
  780. if voltthres_rejection==1 % check voltage threshold rejection
  781. if interp_epoch==1 % check epoch level channel interpolation
  782. chans=[]; chansidx=[];chans_labels2=[];
  783. chans_labels2=cell(1,EEG.nbchan);
  784. for i=1:EEG.nbchan
  785. chans_labels2{i}= EEG.chanlocs(i).labels;
  786. end
  787. [chans,chansidx] = ismember(frontal_channels, chans_labels2);
  788. frontal_channels_idx = chansidx(chansidx ~= 0);
  789. badChans = zeros(EEG.nbchan, EEG.trials);
  790. badepoch=zeros(1, EEG.trials);
  791. if isempty(frontal_channels_idx)==1 % check whether there is any frontal channel in dataset to check
  792. warning('No frontal channels from the list present in the data. Only epoch interpolation will be performed.');
  793. else
  794. % find artifaceted epochs by detecting outlier voltage in the specified channels list and remove epoch if artifacted in those channels
  795. for ch =1:length(frontal_channels_idx)
  796. EEG = pop_eegthresh(EEG,1, frontal_channels_idx(ch), volt_threshold(1), volt_threshold(2), EEG.xmin, EEG.xmax,0,0);
  797. EEG = eeg_checkset( EEG );
  798. EEG = eeg_rejsuperpose( EEG, 1, 1, 1, 1, 1, 1, 1, 1);
  799. badChans(ch,:) = EEG.reject.rejglobal;
  800. end
  801. for ii=1:size(badChans, 2)
  802. badepoch(ii)=sum(badChans(:,ii));
  803. end
  804. badepoch=logical(badepoch);
  805. end
  806. % If all epochs are artifacted, save the dataset and ignore rest of the preprocessing for this subject.
  807. if sum(badepoch)==EEG.trials || sum(badepoch)+1==EEG.trials
  808. all_bad_epochs=1;
  809. warning(['No usable data for datafile', [datafile_names{s} '_Merged']]);
  810. if output_format==1
  811. EEG = eeg_checkset(EEG);
  812. EEG = pop_editset(EEG, 'setname', [datafile_names{s} '_arm_no_usable_data_all_bad_epochs']);
  813. EEG = pop_saveset(EEG, 'filename', [datafile_names{s} '_arm_no_usable_data_all_bad_epochs.set'],'filepath', [output_location filesep 'processed_data' filesep]); % save .set format
  814. elseif output_format==2
  815. save([[output_location filesep 'processed_data' filesep] [datafile_names{s} '_arm_no_usable_data_all_bad_epochs.mat']], 'EEG'); % save .mat format
  816. end
  817. else
  818. EEG = pop_rejepoch( EEG, badepoch, 0);
  819. EEG = eeg_checkset(EEG);
  820. end
  821. if all_bad_epochs==1
  822. warning(['No usable data for datafile', [datafile_names{s} '_Merged']]);
  823. else
  824. % Interpolate artifacted data for all reaming channels
  825. badChans = zeros(EEG.nbchan, EEG.trials);
  826. % Find artifacted epochs by detecting outlier voltage but don't remove
  827. for ch=1:EEG.nbchan
  828. EEG = pop_eegthresh(EEG,1, ch, volt_threshold(1), volt_threshold(2), EEG.xmin, EEG.xmax,0,0);
  829. EEG = eeg_checkset(EEG);
  830. EEG = eeg_rejsuperpose(EEG, 1, 1, 1, 1, 1, 1, 1, 1);
  831. badChans(ch,:) = EEG.reject.rejglobal;
  832. end
  833. tmpData = zeros(EEG.nbchan, EEG.pnts, EEG.trials);
  834. for e = 1:EEG.trials
  835. % Initialize variables EEGe and EEGe_interp;
  836. EEGe = []; EEGe_interp = []; badChanNum = [];
  837. % Select only this epoch (e)
  838. EEGe = pop_selectevent( EEG, 'epoch', e, 'deleteevents', 'off', 'deleteepochs', 'on', 'invertepochs', 'off');
  839. badChanNum = find(badChans(:,e)==1); % find which channels are bad for this epoch
  840. EEGe_interp = eeg_interp(EEGe,badChanNum); %interpolate the bad channels for this epoch
  841. tmpData(:,:,e) = EEGe_interp.data; % store interpolated data into matrix
  842. end
  843. EEG.data = tmpData; % now that all of the epochs have been interpolated, write the data back to the main file
  844. % If more than 10% of channels in an epoch were interpolated, reject that epoch
  845. badepoch=zeros(1, EEG.trials);
  846. for ei=1:EEG.trials
  847. NumbadChan = badChans(:,ei); % find how many channels are bad in an epoch
  848. if sum(NumbadChan) > round((10/100)*EEG.nbchan)% check if more than 10% are bad
  849. badepoch (ei)= sum(NumbadChan);
  850. end
  851. end
  852. badepoch=logical(badepoch);
  853. end
  854. % If all epochs are artifacted, save the dataset and ignore rest of the preprocessing for this subject.
  855. if sum(badepoch)==EEG.trials || sum(badepoch)+1==EEG.trials
  856. all_bad_epochs=1;
  857. warning(['No usable data for datafile', [datafile_names{s} '_Merged']]);
  858. if output_format==1
  859. EEG = eeg_checkset(EEG);
  860. EEG = pop_editset(EEG, 'setname', [datafile_names{s} '_arm_no_usable_data_all_bad_epochs']);
  861. EEG = pop_saveset(EEG, 'filename', [datafile_names{s} '_arm_no_usable_data_all_bad_epochs.set'],'filepath', [output_location filesep 'processed_data' filesep]); % save .set format
  862. elseif output_format==2
  863. save([[output_location filesep 'processed_data' filesep] [datafile_names{s} '_arm_no_usable_data_all_bad_epochs.mat']], 'EEG'); % save .mat format
  864. end
  865. else
  866. EEG = pop_rejepoch(EEG, badepoch, 0);
  867. EEG = eeg_checkset(EEG);
  868. end
  869. else % if no epoch level channel interpolation
  870. EEG = pop_eegthresh(EEG, 1, (1:EEG.nbchan), volt_threshold(1), volt_threshold(2), EEG.xmin, EEG.xmax, 0, 0);
  871. EEG = eeg_checkset(EEG);
  872. EEG = eeg_rejsuperpose( EEG, 1, 1, 1, 1, 1, 1, 1, 1);
  873. end % end of epoch level channel interpolation if statement
  874. % If all epochs are artifacted, save the dataset and ignore rest of the preprocessing for this subject.
  875. if sum(EEG.reject.rejthresh)==EEG.trials || sum(EEG.reject.rejthresh)+1==EEG.trials
  876. all_bad_epochs=1;
  877. warning(['No usable data for datafile', [datafile_names{s} '_Merged']]);
  878. if output_format==1
  879. EEG = eeg_checkset(EEG);
  880. EEG = pop_editset(EEG, 'setname', [datafile_names{s} '_arm_no_usable_data_all_bad_epochs']);
  881. EEG = pop_saveset(EEG, 'filename', [datafile_names{s} '_arm_no_usable_data_all_bad_epochs.set'],'filepath', [output_location filesep 'processed_data' filesep]); % save .set format
  882. elseif output_format==2
  883. save([[output_location filesep 'processed_data' filesep] [datafile_names{s} '_arm_no_usable_data_all_bad_epochs.mat']], 'EEG'); % save .mat format
  884. end
  885. else
  886. EEG = pop_rejepoch(EEG,(EEG.reject.rejthresh), 0);
  887. EEG = eeg_checkset(EEG);
  888. end
  889. end % end of voltage threshold rejection if statement
  890. % if all epochs are found bad during artifact rejection
  891. if all_bad_epochs==1
  892. total_epochs_after_artifact_rejection(s)=0;
  893. total_channels_interpolated(s)=0;
  894. continue % ignore rest of the processing and go to next datafile
  895. else
  896. total_epochs_after_artifact_rejection(s)=EEG.trials;
  897. total_AS_after_artifact_rejection(s) = length(find(strcmp({EEG.event.type},'AS')));
  898. total_QS_after_artifact_rejection(s) = length(find(strcmp({EEG.event.type}, 'QS')));
  899. total_W_after_artifact_rejection(s) = length(find(strcmp({EEG.event.type}, 'W')));
  900. total_I_after_artifact_rejection(s) = length(find(strcmp({EEG.event.type}, 'I')));
  901. end
  902. %% STEP 15: Interpolate deleted channels
  903. if interp_channels==1
  904. EEG = eeg_interp(EEG, channels_analysed);
  905. EEG = eeg_checkset(EEG);
  906. end
  907. if numel(FASTbadChans)==0 && numel(ica_prep_badChans)==0
  908. total_channels_interpolated(s)=0;
  909. else
  910. total_channels_interpolated(s)=numel(FASTbadChans)+ numel(ica_prep_badChans);
  911. end
  912. %% STEP 16: Rereference data
  913. EEG = eeg_checkset(EEG);
  914. reref = [];
  915. EEG = pop_reref(EEG, reref);
  916. %% Save processed data
  917. if output_format==1
  918. EEG = eeg_checkset(EEG);
  919. EEG = pop_editset(EEG, 'setname', [datafile_names{s} '_arm_processed_data']);
  920. EEG = pop_saveset(EEG, 'filename', [datafile_names{s} '_arm_processed_data.set'],'filepath', [output_location filesep 'processed_data' filesep]); % save .set format
  921. elseif output_format==2
  922. save([[output_location filesep 'processed_data' filesep] [datafile_names{s} '_arm_processed_data.mat']], 'EEG'); % save .mat format
  923. end
  924. %% Create the report table for all the data files with relevant preprocessing outputs.
  925. report_table=table(datafile_names(s)', reference_used_for_faster(s)', faster_bad_channels(s)', ica_preparation_bad_channels(s)', length_ica_data(s)', ...
  926. total_ICs(s)', ICs_removed(s)', total_epochs_before_artifact_rejection(s)', total_AS(s)', total_QS(s)', total_W(s)', total_I(s)', ...
  927. total_epochs_after_artifact_rejection(s)', total_AS_after_artifact_rejection(s)', total_QS_after_artifact_rejection(s)', ...
  928. total_W_after_artifact_rejection(s)', total_I_after_artifact_rejection(s)', total_channels_interpolated(s)');
  929. report_table.Properties.VariableNames={'folder_list_Resting', 'reference_used_for_faster', 'faster_bad_channels', ...
  930. 'ica_preparation_bad_channels', 'length_ica_data', 'total_ICs', 'ICs_removed', 'total_epochs_before_artifact_rejection', ...
  931. 'total_AS', 'total_QS', 'total_W', 'total_I', 'total_epochs_after_artifact_rejection', 'total_AS_after_artifact_rejection', 'total_QS_after_artifact_rejection',...
  932. 'total_W_after_artifact_rejection', 'total_I_after_artifact_rejection','total_channels_interpolated'};
  933. writetable(report_table, ['/Volumes/BEHAVEBB/Amy/ADVANCES/EEG/Processed/Report/1m/', subject '_1m_arm_MADE_report_', datestr(now,'dd-mm-yyyy'),'.csv']);
  934. elseif task == 2 % sResting in bassinet
  935. % Select events specific to sResting
  936. EEG = pop_selectevent(EEG3, 'Task', [task],'deleteevents','on'); % just delete other flags/markers, but will NOT delte the EEG data
  937. %% STEP 12: Segment data into fixed length epochs
  938. EEG = eeg_checkset(EEG);
  939. EEG = pop_epoch(EEG, task_event_markers, task_epoch_length, 'epochinfo', 'yes');
  940. EEG = pop_selectevent( EEG, 'latency','-.1 <= .1','deleteevents','on');
  941. total_epochs_before_artifact_rejection(s)=EEG.trials;
  942. total_AS(s) = length(find(strcmp({EEG.event.type},'AS')));
  943. total_QS(s) = length(find(strcmp({EEG.event.type}, 'QS')));
  944. total_W(s) = length(find(strcmp({EEG.event.type}, 'W')));
  945. total_I(s) = length(find(strcmp({EEG.event.type}, 'I')));
  946. %% STEP 13: Remove baseline
  947. %if remove_baseline==1
  948. %EEG = eeg_checkset( EEG );
  949. %EEG = pop_rmbase( EEG, baseline_window);
  950. %end
  951. %% STEP 14: Artifact rejection
  952. all_bad_epochs=0;
  953. if voltthres_rejection==1 % check voltage threshold rejection
  954. if interp_epoch==1 % check epoch level channel interpolation
  955. chans=[]; chansidx=[];chans_labels2=[];
  956. chans_labels2=cell(1,EEG.nbchan);
  957. for i=1:EEG.nbchan
  958. chans_labels2{i}= EEG.chanlocs(i).labels;
  959. end
  960. [chans,chansidx] = ismember(frontal_channels, chans_labels2);
  961. frontal_channels_idx = chansidx(chansidx ~= 0);
  962. badChans = zeros(EEG.nbchan, EEG.trials);
  963. badepoch=zeros(1, EEG.trials);
  964. if isempty(frontal_channels_idx)==1 % check whether there is any frontal channel in dataset to check
  965. warning('No frontal channels from the list present in the data. Only epoch interpolation will be performed.');
  966. else
  967. % find artifaceted epochs by detecting outlier voltage in the specified channels list and remove epoch if artifacted in those channels
  968. for ch =1:length(frontal_channels_idx)
  969. EEG = pop_eegthresh(EEG,1, frontal_channels_idx(ch), volt_threshold(1), volt_threshold(2), EEG.xmin, EEG.xmax,0,0);
  970. EEG = eeg_checkset( EEG );
  971. EEG = eeg_rejsuperpose( EEG, 1, 1, 1, 1, 1, 1, 1, 1);
  972. badChans(ch,:) = EEG.reject.rejglobal;
  973. end
  974. for ii=1:size(badChans, 2)
  975. badepoch(ii)=sum(badChans(:,ii));
  976. end
  977. badepoch=logical(badepoch);
  978. end
  979. % If all epochs are artifacted, save the dataset and ignore rest of the preprocessing for this subject.
  980. if sum(badepoch)==EEG.trials || sum(badepoch)+1==EEG.trials
  981. all_bad_epochs=1;
  982. warning(['No usable data for datafile', [datafile_names{s} '_Merged']]);
  983. if output_format==1
  984. EEG = eeg_checkset(EEG);
  985. EEG = pop_editset(EEG, 'setname', [datafile_names{s} '_sResting_no_usable_data_all_bad_epochs']);
  986. EEG = pop_saveset(EEG, 'filename', [datafile_names{s} '_sResting_no_usable_data_all_bad_epochs.set'],'filepath', [output_location filesep 'processed_data' filesep]); % save .set format
  987. elseif output_format==2
  988. save([[output_location filesep 'processed_data' filesep] [datafile_names{s} '_sResting_no_usable_data_all_bad_epochs.mat']], 'EEG'); % save .mat format
  989. end
  990. else
  991. EEG = pop_rejepoch( EEG, badepoch, 0);
  992. EEG = eeg_checkset(EEG);
  993. end
  994. if all_bad_epochs==1
  995. warning(['No usable data for datafile', [datafile_names{s} '_Merged']]);
  996. else
  997. % Interpolate artifacted data for all reaming channels
  998. badChans = zeros(EEG.nbchan, EEG.trials);
  999. % Find artifacted epochs by detecting outlier voltage but don't remove
  1000. for ch=1:EEG.nbchan
  1001. EEG = pop_eegthresh(EEG,1, ch, volt_threshold(1), volt_threshold(2), EEG.xmin, EEG.xmax,0,0);
  1002. EEG = eeg_checkset(EEG);
  1003. EEG = eeg_rejsuperpose(EEG, 1, 1, 1, 1, 1, 1, 1, 1);
  1004. badChans(ch,:) = EEG.reject.rejglobal;
  1005. end
  1006. tmpData = zeros(EEG.nbchan, EEG.pnts, EEG.trials);
  1007. for e = 1:EEG.trials
  1008. % Initialize variables EEGe and EEGe_interp;
  1009. EEGe = []; EEGe_interp = []; badChanNum = [];
  1010. % Select only this epoch (e)
  1011. EEGe = pop_selectevent( EEG, 'epoch', e, 'deleteevents', 'off', 'deleteepochs', 'on', 'invertepochs', 'off');
  1012. badChanNum = find(badChans(:,e)==1); % find which channels are bad for this epoch
  1013. EEGe_interp = eeg_interp(EEGe,badChanNum); %interpolate the bad channels for this epoch
  1014. tmpData(:,:,e) = EEGe_interp.data; % store interpolated data into matrix
  1015. end
  1016. EEG.data = tmpData; % now that all of the epochs have been interpolated, write the data back to the main file
  1017. % If more than 10% of channels in an epoch were interpolated, reject that epoch
  1018. badepoch=zeros(1, EEG.trials);
  1019. for ei=1:EEG.trials
  1020. NumbadChan = badChans(:,ei); % find how many channels are bad in an epoch
  1021. if sum(NumbadChan) > round((10/100)*EEG.nbchan)% check if more than 10% are bad
  1022. badepoch (ei)= sum(NumbadChan);
  1023. end
  1024. end
  1025. badepoch=logical(badepoch);
  1026. end
  1027. % If all epochs are artifacted, save the dataset and ignore rest of the preprocessing for this subject.
  1028. if sum(badepoch)==EEG.trials || sum(badepoch)+1==EEG.trials
  1029. all_bad_epochs=1;
  1030. warning(['No usable data for datafile', [datafile_names{s} '_Merged']]);
  1031. if output_format==1
  1032. EEG = eeg_checkset(EEG);
  1033. EEG = pop_editset(EEG, 'setname', [datafile_names{s} '_sResting_no_usable_data_all_bad_epochs']);
  1034. EEG = pop_saveset(EEG, 'filename', [datafile_names{s} '_sResting_no_usable_data_all_bad_epochs.set'],'filepath', [output_location filesep 'processed_data' filesep]); % save .set format
  1035. elseif output_format==2
  1036. save([[output_location filesep 'processed_data' filesep] [datafile_names{s} '_sResting_no_usable_data_all_bad_epochs.mat']], 'EEG'); % save .mat format
  1037. end
  1038. else
  1039. EEG = pop_rejepoch(EEG, badepoch, 0);
  1040. EEG = eeg_checkset(EEG);
  1041. end
  1042. else % if no epoch level channel interpolation
  1043. EEG = pop_eegthresh(EEG, 1, (1:EEG.nbchan), volt_threshold(1), volt_threshold(2), EEG.xmin, EEG.xmax, 0, 0);
  1044. EEG = eeg_checkset(EEG);
  1045. EEG = eeg_rejsuperpose( EEG, 1, 1, 1, 1, 1, 1, 1, 1);
  1046. end % end of epoch level channel interpolation if statement
  1047. % If all epochs are artifacted, save the dataset and ignore rest of the preprocessing for this subject.
  1048. if sum(EEG.reject.rejthresh)==EEG.trials || sum(EEG.reject.rejthresh)+1==EEG.trials
  1049. all_bad_epochs=1;
  1050. warning(['No usable data for datafile', [datafile_names{s} '_Merged']]);
  1051. if output_format==1
  1052. EEG = eeg_checkset(EEG);
  1053. EEG = pop_editset(EEG, 'setname', [datafile_names{s} '_sResting_no_usable_data_all_bad_epochs']);
  1054. EEG = pop_saveset(EEG, 'filename', [datafile_names{s} '_sResting_no_usable_data_all_bad_epochs.set'],'filepath', [output_location filesep 'processed_data' filesep]); % save .set format
  1055. elseif output_format==2
  1056. save([[output_location filesep 'processed_data' filesep] [datafile_names{s} '_sResting_no_usable_data_all_bad_epochs.mat']], 'EEG'); % save .mat format
  1057. end
  1058. else
  1059. EEG = pop_rejepoch(EEG,(EEG.reject.rejthresh), 0);
  1060. EEG = eeg_checkset(EEG);
  1061. end
  1062. end % end of voltage threshold rejection if statement
  1063. % if all epochs are found bad during artifact rejection
  1064. if all_bad_epochs==1
  1065. total_epochs_after_artifact_rejection(s)=0;
  1066. total_channels_interpolated(s)=0;
  1067. continue % ignore rest of the processing and go to next datafile
  1068. else
  1069. total_epochs_after_artifact_rejection(s)=EEG.trials;
  1070. total_AS_after_artifact_rejection(s) = length(find(strcmp({EEG.event.type},'AS')));
  1071. total_QS_after_artifact_rejection(s) = length(find(strcmp({EEG.event.type}, 'QS')));
  1072. total_W_after_artifact_rejection(s) = length(find(strcmp({EEG.event.type}, 'W')));
  1073. total_I_after_artifact_rejection(s) = length(find(strcmp({EEG.event.type}, 'I')));
  1074. end
  1075. %% STEP 15: Interpolate deleted channels
  1076. if interp_channels==1
  1077. EEG = eeg_interp(EEG, channels_analysed);
  1078. EEG = eeg_checkset(EEG);
  1079. end
  1080. if numel(FASTbadChans)==0 && numel(ica_prep_badChans)==0
  1081. total_channels_interpolated(s)=0;
  1082. else
  1083. total_channels_interpolated(s)=numel(FASTbadChans)+ numel(ica_prep_badChans);
  1084. end
  1085. %% STEP 16: Rereference data
  1086. EEG = eeg_checkset(EEG);
  1087. reref = [];
  1088. EEG = pop_reref(EEG, reref);
  1089. %% Save processed data
  1090. if output_format==1
  1091. EEG = eeg_checkset(EEG);
  1092. EEG = pop_editset(EEG, 'setname', [datafile_names{s} '_sResting_processed_data']);
  1093. EEG = pop_saveset(EEG, 'filename', [datafile_names{s} '_sResting_processed_data.set'],'filepath', [output_location filesep 'processed_data' filesep]); % save .set format
  1094. elseif output_format==2
  1095. save([[output_location filesep 'processed_data' filesep] [datafile_names{s} '_sResting_processed_data.mat']], 'EEG'); % save .mat format
  1096. end
  1097. %% Create the report table for all the data files with relevant preprocessing outputs.
  1098. report_table=table(datafile_names(s)', reference_used_for_faster(s)', faster_bad_channels(s)', ica_preparation_bad_channels(s)', length_ica_data(s)', ...
  1099. total_ICs(s)', ICs_removed(s)', total_epochs_before_artifact_rejection(s)', total_AS(s)', total_QS(s)', total_W(s)', total_I(s)', ...
  1100. total_epochs_after_artifact_rejection(s)', total_AS_after_artifact_rejection(s)', total_QS_after_artifact_rejection(s)', ...
  1101. total_W_after_artifact_rejection(s)', total_I_after_artifact_rejection(s)', total_channels_interpolated(s)');
  1102. report_table.Properties.VariableNames={'folder_list_Resting', 'reference_used_for_faster', 'faster_bad_channels', ...
  1103. 'ica_preparation_bad_channels', 'length_ica_data', 'total_ICs', 'ICs_removed', 'total_epochs_before_artifact_rejection', ...
  1104. 'total_AS', 'total_QS', 'total_W', 'total_I', 'total_epochs_after_artifact_rejection', 'total_AS_after_artifact_rejection', 'total_QS_after_artifact_rejection',...
  1105. 'total_W_after_artifact_rejection', 'total_I_after_artifact_rejection','total_channels_interpolated'};
  1106. writetable(report_table, ['/Volumes/BEHAVEBB/Amy/ADVANCES/EEG/Processed/Report/1m/', subject '_1m_sResting_MADE_report_', datestr(now,'dd-mm-yyyy'),'.csv']);
  1107. elseif task == 3 % MMN
  1108. % Select events specific to MMN
  1109. EEG = pop_selectevent(EEG3, 'Task', [task],'deleteevents','on');
  1110. %% STEP 12: Segment data into fixed length epochs
  1111. EEG = eeg_checkset(EEG);
  1112. EEG = pop_epoch(EEG, task_event_markers, task_epoch_length, 'epochinfo', 'yes');
  1113. EEG = pop_selectevent( EEG, 'latency','-.1 <= .1','deleteevents','on');
  1114. total_epochs_before_artifact_rejection(s)=EEG.trials;
  1115. total_standard(s) = length(find(strcmp({EEG.event.type},'1') & [EEG.event.NumOfPrevTone]>2 )); %need at least 3 standards preceeding to use
  1116. total_deviant(s) = length(find(strcmp({EEG.event.type},'2')));
  1117. total_novel(s) = length(find(strcmp({EEG.event.type},'3')));
  1118. total_standard_AS(s) = length(find(strcmp({EEG.event.type},'1') & [EEG.event.NumOfPrevTone]>2 & strcmp({EEG.event.State},'AS'))); %need at least 3 standards preceeding to use
  1119. total_standard_QS(s) = length(find(strcmp({EEG.event.type},'1') & [EEG.event.NumOfPrevTone]>2 & strcmp({EEG.event.State}, 'QS')));
  1120. total_standard_W(s) = length(find(strcmp({EEG.event.type},'1') & [EEG.event.NumOfPrevTone]>2 & strcmp({EEG.event.State}, 'W')));
  1121. total_standard_I(s) = length(find(strcmp({EEG.event.type},'1') & [EEG.event.NumOfPrevTone]>2 & strcmp({EEG.event.State}, 'I')));
  1122. total_deviant_AS(s) = length(find(strcmp({EEG.event.type},'2') & strcmp({EEG.event.State}, 'AS')));
  1123. total_deviant_QS(s) = length(find(strcmp({EEG.event.type},'2') & strcmp({EEG.event.State}, 'QS')));
  1124. total_deviant_W(s) = length(find(strcmp({EEG.event.type},'2') & strcmp({EEG.event.State}, 'W')));
  1125. total_deviant_I(s) = length(find(strcmp({EEG.event.type},'2') & strcmp({EEG.event.State}, 'I')));
  1126. total_novel_AS(s) = length(find(strcmp({EEG.event.type},'3') & strcmp({EEG.event.State}, 'AS')));
  1127. total_novel_QS(s) = length(find(strcmp({EEG.event.type},'3') & strcmp({EEG.event.State}, 'QS')));
  1128. total_novel_W(s) = length(find(strcmp({EEG.event.type},'3') & strcmp({EEG.event.State}, 'W')));
  1129. total_novel_I(s) = length(find(strcmp({EEG.event.type},'3') & strcmp({EEG.event.State}, 'I')));
  1130. %% STEP 13: Remove baseline
  1131. EEG = eeg_checkset( EEG );
  1132. baseline_window = [-100,0];
  1133. EEG = pop_rmbase( EEG, baseline_window);
  1134. %% STEP 14: Artifact rejection
  1135. all_bad_epochs=0;
  1136. if voltthres_rejection==1 % check voltage threshold rejection
  1137. if interp_epoch==1 % check epoch level channel interpolation
  1138. chans=[]; chansidx=[];chans_labels2=[];
  1139. chans_labels2=cell(1,EEG.nbchan);
  1140. for i=1:EEG.nbchan
  1141. chans_labels2{i}= EEG.chanlocs(i).labels;
  1142. end
  1143. [chans,chansidx] = ismember(frontal_channels, chans_labels2);
  1144. frontal_channels_idx = chansidx(chansidx ~= 0);
  1145. badChans = zeros(EEG.nbchan, EEG.trials);
  1146. badepoch=zeros(1, EEG.trials);
  1147. if isempty(frontal_channels_idx)==1 % check whether there is any frontal channel in dataset to check
  1148. warning('No frontal channels from the list present in the data. Only epoch interpolation will be performed.');
  1149. else
  1150. % find artifaceted epochs by detecting outlier voltage in the specified channels list and remove epoch if artifacted in those channels
  1151. for ch =1:length(frontal_channels_idx)
  1152. EEG = pop_eegthresh(EEG,1, frontal_channels_idx(ch), volt_threshold(1), volt_threshold(2), EEG.xmin, EEG.xmax,0,0);
  1153. EEG = eeg_checkset( EEG );
  1154. EEG = eeg_rejsuperpose( EEG, 1, 1, 1, 1, 1, 1, 1, 1);
  1155. badChans(ch,:) = EEG.reject.rejglobal;
  1156. end
  1157. for ii=1:size(badChans, 2)
  1158. badepoch(ii)=sum(badChans(:,ii));
  1159. end
  1160. badepoch=logical(badepoch);
  1161. end
  1162. % If all epochs are artifacted, save the dataset and ignore rest of the preprocessing for this subject.
  1163. if sum(badepoch)==EEG.trials || sum(badepoch)+1==EEG.trials
  1164. all_bad_epochs=1;
  1165. warning(['No usable data for datafile', [datafile_names{s} '_Merged']]);
  1166. if output_format==1
  1167. EEG = eeg_checkset(EEG);
  1168. EEG = pop_editset(EEG, 'setname', [datafile_names{s} '_MMN_no_usable_data_all_bad_epoch']);
  1169. EEG = pop_saveset(EEG, 'filename', [datafile_names{s} '_MMN_no_usable_data_all_bad_epoch.set'],'filepath', [output_location filesep 'processed_data' filesep]); % save .set format
  1170. elseif output_format==2
  1171. save([[output_location filesep 'processed_data' filesep] [datafile_names{s} '_MMN_no_usable_data_all_bad_epochs.mat']], 'EEG'); % save .mat format
  1172. end
  1173. else
  1174. EEG = pop_rejepoch( EEG, badepoch, 0);
  1175. EEG = eeg_checkset(EEG);
  1176. end
  1177. if all_bad_epochs==1
  1178. warning(['No usable data for datafile', [datafile_names{s} '_Merged']]);
  1179. else
  1180. % Interpolate artifacted data for all reaming channels
  1181. badChans = zeros(EEG.nbchan, EEG.trials);
  1182. % Find artifacted epochs by detecting outlier voltage but don't remove
  1183. for ch=1:EEG.nbchan
  1184. EEG = pop_eegthresh(EEG,1, ch, volt_threshold(1), volt_threshold(2), EEG.xmin, EEG.xmax,0,0);
  1185. EEG = eeg_checkset(EEG);
  1186. EEG = eeg_rejsuperpose(EEG, 1, 1, 1, 1, 1, 1, 1, 1);
  1187. badChans(ch,:) = EEG.reject.rejglobal;
  1188. end
  1189. tmpData = zeros(EEG.nbchan, EEG.pnts, EEG.trials);
  1190. for e = 1:EEG.trials
  1191. % Initialize variables EEGe and EEGe_interp;
  1192. EEGe = []; EEGe_interp = []; badChanNum = [];
  1193. % Select only this epoch (e)
  1194. EEGe = pop_selectevent( EEG, 'epoch', e, 'deleteevents', 'off', 'deleteepochs', 'on', 'invertepochs', 'off');
  1195. badChanNum = find(badChans(:,e)==1); % find which channels are bad for this epoch
  1196. EEGe_interp = eeg_interp(EEGe,badChanNum); %interpolate the bad channels for this epoch
  1197. tmpData(:,:,e) = EEGe_interp.data; % store interpolated data into matrix
  1198. end
  1199. EEG.data = tmpData; % now that all of the epochs have been interpolated, write the data back to the main file
  1200. % If more than 10% of channels in an epoch were interpolated, reject that epoch
  1201. badepoch=zeros(1, EEG.trials);
  1202. for ei=1:EEG.trials
  1203. NumbadChan = badChans(:,ei); % find how many channels are bad in an epoch
  1204. if sum(NumbadChan) > round((10/100)*EEG.nbchan)% check if more than 10% are bad
  1205. badepoch (ei)= sum(NumbadChan);
  1206. end
  1207. end
  1208. badepoch=logical(badepoch);
  1209. end
  1210. % If all epochs are artifacted, save the dataset and ignore rest of the preprocessing for this subject.
  1211. if sum(badepoch)==EEG.trials || sum(badepoch)+1==EEG.trials
  1212. all_bad_epochs=1;
  1213. warning(['No usable data for datafile', [datafile_names{s} '_Merged']]);
  1214. if output_format==1
  1215. EEG = eeg_checkset(EEG);
  1216. EEG = pop_editset(EEG, 'setname', [datafile_names{s} '_MMN_no_usable_data_all_bad_epochs']);
  1217. EEG = pop_saveset(EEG, 'filename', [datafile_names{s} '_MMN_no_usable_data_all_bad_epochs.set'],'filepath', [output_location filesep 'processed_data' filesep]); % save .set format
  1218. elseif output_format==2
  1219. save([[output_location filesep 'processed_data' filesep] [datafile_names{s} '_MMN_no_usable_data_all_bad_epochs.mat']], 'EEG'); % save .mat format
  1220. end
  1221. else
  1222. EEG = pop_rejepoch(EEG, badepoch, 0);
  1223. EEG = eeg_checkset(EEG);
  1224. end
  1225. else % if no epoch level channel interpolation
  1226. EEG = pop_eegthresh(EEG, 1, (1:EEG.nbchan), volt_threshold(1), volt_threshold(2), EEG.xmin, EEG.xmax, 0, 0);
  1227. EEG = eeg_checkset(EEG);
  1228. EEG = eeg_rejsuperpose( EEG, 1, 1, 1, 1, 1, 1, 1, 1);
  1229. end % end of epoch level channel interpolation if statement
  1230. % If all epochs are artifacted, save the dataset and ignore rest of the preprocessing for this subject.
  1231. if sum(EEG.reject.rejthresh)==EEG.trials || sum(EEG.reject.rejthresh)+1==EEG.trials
  1232. all_bad_epochs=1;
  1233. warning(['No usable data for datafile', [datafile_names{s} '_Merged']]);
  1234. if output_format==1
  1235. EEG = eeg_checkset(EEG);
  1236. EEG = pop_editset(EEG, 'setname', [datafile_names{s} '_MMN_no_usable_data_all_bad_epochs']);
  1237. EEG = pop_saveset(EEG, 'filename', [datafile_names{s} '_MMN_no_usable_data_all_bad_epochs.set'],'filepath', [output_location filesep 'processed_data' filesep]); % save .set format
  1238. elseif output_format==2
  1239. save([[output_location filesep 'processed_data' filesep] [datafile_names{s} '_MMN_no_usable_data_all_bad_epochs.mat']], 'EEG'); % save .mat format
  1240. end
  1241. else
  1242. EEG = pop_rejepoch(EEG,(EEG.reject.rejthresh), 0);
  1243. EEG = eeg_checkset(EEG);
  1244. end
  1245. end % end of voltage threshold rejection if statement
  1246. % if all epochs are found bad during artifact rejection
  1247. if all_bad_epochs==1
  1248. total_epochs_after_artifact_rejection(s)=0;
  1249. total_channels_interpolated(s)=0;
  1250. continue % ignore rest of the processing and go to next datafile
  1251. else
  1252. total_epochs_after_artifact_rejection(s)=EEG.trials;
  1253. total_standard_after_artifact_rejection(s) = length(find(strcmp({EEG.event.type},'1') & [EEG.event.NumOfPrevTone]>2)); %need at least 3 standards preceeding to use
  1254. total_deviant_after_artifact_rejection(s) = length(find(strcmp({EEG.event.type},'2')));
  1255. total_novel_after_artifact_rejection(s) = length(find(strcmp({EEG.event.type},'3')));
  1256. total_standard_AS_after_artifact_rejection(s) = length(find(strcmp({EEG.event.type},'1') & [EEG.event.NumOfPrevTone]>2 & strcmp({EEG.event.State}, 'AS'))); %need at least 3 standards preceeding to use
  1257. total_standard_QS_after_artifact_rejection(s) = length(find(strcmp({EEG.event.type},'1') & [EEG.event.NumOfPrevTone]>2 & strcmp({EEG.event.State}, 'QS'))); %need at least 3 standards preceeding to use
  1258. total_standard_W_after_artifact_rejection(s) = length(find(strcmp({EEG.event.type},'1') & [EEG.event.NumOfPrevTone]>2 & strcmp({EEG.event.State}, 'W'))); %need at least 3 standards preceeding to use
  1259. total_standard_I_after_artifact_rejection(s) = length(find(strcmp({EEG.event.type},'1') & [EEG.event.NumOfPrevTone]>2 & strcmp({EEG.event.State}, 'I'))); %need at least 3 standards preceeding to use
  1260. total_deviant_AS_after_artifact_rejection(s) = length(find(strcmp({EEG.event.type},'2') & strcmp({EEG.event.State}, 'AS')));
  1261. total_deviant_QS_after_artifact_rejection(s) = length(find(strcmp({EEG.event.type},'2') & strcmp({EEG.event.State}, 'QS')));
  1262. total_deviant_W_after_artifact_rejection(s) = length(find(strcmp({EEG.event.type},'2') & strcmp({EEG.event.State}, 'W')));
  1263. total_deviant_I_after_artifact_rejection(s) = length(find(strcmp({EEG.event.type},'2') & strcmp({EEG.event.State}, 'I')));
  1264. total_novel_AS_after_artifact_rejection(s) = length(find(strcmp({EEG.event.type},'3') & strcmp({EEG.event.State}, 'AS')));
  1265. total_novel_QS_after_artifact_rejection(s) = length(find(strcmp({EEG.event.type},'3') & strcmp({EEG.event.State}, 'QS')));
  1266. total_novel_W_after_artifact_rejection(s) = length(find(strcmp({EEG.event.type},'3') & strcmp({EEG.event.State}, 'W')));
  1267. total_novel_I_after_artifact_rejection(s) = length(find(strcmp({EEG.event.type},'3') & strcmp({EEG.event.State}, 'I')));
  1268. end
  1269. %% STEP 15: Interpolate deleted channels
  1270. if interp_channels==1
  1271. EEG = eeg_interp(EEG, channels_analysed);
  1272. EEG = eeg_checkset(EEG);
  1273. end
  1274. if numel(FASTbadChans)==0 && numel(ica_prep_badChans)==0
  1275. total_channels_interpolated(s)=0;
  1276. else
  1277. total_channels_interpolated(s)=numel(FASTbadChans)+ numel(ica_prep_badChans);
  1278. end
  1279. %% STEP 16: Rereference data
  1280. EEG = eeg_checkset(EEG);
  1281. reref = {'E57','E100'};
  1282. reref_idx=zeros(1, length(reref));
  1283. for rr=1:length(reref)
  1284. reref_idx(rr)=find(strcmp({EEG.chanlocs.labels}, reref{rr}));
  1285. end
  1286. EEG = eeg_checkset(EEG);
  1287. EEG = pop_reref( EEG, reref_idx);
  1288. %% Save processed data
  1289. if output_format==1
  1290. EEG = eeg_checkset(EEG);
  1291. EEG = pop_editset(EEG, 'setname', [datafile_names{s} '_MMN_processed_data']);
  1292. EEG = pop_saveset(EEG, 'filename', [datafile_names{s} '_MMN_processed_data.set'],'filepath', [output_location filesep 'processed_data' filesep]); % save .set format
  1293. elseif output_format==2
  1294. save([[output_location filesep 'processed_data' filesep] [datafile_names{s} '_MMN_processed_data.mat']], 'EEG'); % save .mat format
  1295. end
  1296. %% Create the report table for all the data files with relevant preprocessing outputs.
  1297. report_table=table(datafile_names(s)', reference_used_for_faster(s)', faster_bad_channels(s)', ica_preparation_bad_channels(s)', length_ica_data(s)', ...
  1298. total_ICs(s)', ICs_removed(s)', total_epochs_before_artifact_rejection(s)', total_standard(s)', total_deviant(s)', total_novel(s)', ...
  1299. total_standard_AS(s)', total_standard_QS(s)', total_standard_W(s)', total_standard_I(s)', total_deviant_AS(s)', total_deviant_QS(s)', total_deviant_W(s)', total_deviant_I(s)', ...
  1300. total_novel_AS(s)', total_novel_QS(s)', total_novel_W(s)', total_novel_I(s)', total_epochs_after_artifact_rejection(s)', total_standard_after_artifact_rejection(s)', total_deviant_after_artifact_rejection(s)', total_novel_after_artifact_rejection(s)', ...
  1301. total_standard_AS_after_artifact_rejection(s)', total_standard_QS_after_artifact_rejection(s)', total_standard_W_after_artifact_rejection(s)', total_standard_I_after_artifact_rejection(s)', ...
  1302. total_deviant_AS_after_artifact_rejection(s)', total_deviant_QS_after_artifact_rejection(s)', total_deviant_W_after_artifact_rejection(s)', total_deviant_I_after_artifact_rejection(s)', ...
  1303. total_novel_AS_after_artifact_rejection(s)', total_novel_QS_after_artifact_rejection(s)', total_novel_W_after_artifact_rejection(s)', total_novel_I_after_artifact_rejection(s)', total_channels_interpolated(s)');
  1304. report_table.Properties.VariableNames={'folder_list_MMN', 'reference_used_for_faster', 'faster_bad_channels', ...
  1305. 'ica_preparation_bad_channels', 'length_ica_data', 'total_ICs', 'ICs_removed', 'total_epochs_before_artifact_rejection','total_standard', 'total_deviant', 'total_novel',...
  1306. 'total_standard_AS', 'total_standard_QS', 'total_standard_W', 'total_standard_I', 'total_deviant_AS', 'total_deviant_QS', 'total_deviant_W', 'total_deviant_I', ...
  1307. 'total_novel_AS', 'total_novel_QS', 'total_novel_W', 'total_novel_I', 'total_epochs_after_artifact_rejection','total_standard_after_artifact_rejection','total_deviant_after_artifact_rejection','total_novel_after_artifact_rejection', ...
  1308. 'total_standard_AS_after_artifact_rejection', 'total_standard_QS_after_artifact_rejection', 'total_standard_W_after_artifact_rejection', 'total_standard_I_after_artifact_rejection', ...
  1309. 'total_deviant_AS_after_artifact_rejection', 'total_deviant_QS_after_artifact_rejection', 'total_deviant_W_after_artifact_rejection', 'total_deviant_I_after_artifact_rejection', ...
  1310. 'total_novel_AS_after_artifact_rejection', 'total_novel_QS_after_artifact_rejection', 'total_novel_W_after_artifact_rejection', 'total_novel_I_after_artifact_rejection','total_channels_interpolated'};
  1311. writetable(report_table, ['/Volumes/BEHAVEBB/Amy/ADVANCES/EEG/Processed/Report/1m/', subject '_1m_MMN_MADE_report_', datestr(now,'dd-mm-yyyy'),'.csv']);
  1312. elseif task == 4 % VEP
  1313. % Select events specific to VEP
  1314. EEG = pop_selectevent(EEG3, 'Task', [task],'deleteevents','on');
  1315. %% STEP 12: Segment data into fixed length epochs
  1316. EEG = eeg_checkset(EEG);
  1317. EEG = pop_epoch(EEG, task_event_markers, task_epoch_length, 'epochinfo', 'yes');
  1318. EEG = pop_selectevent( EEG, 'latency','-.1 <= .1','deleteevents','on');
  1319. total_vep=EEG.trials;
  1320. total_AS_vep=length(find(strcmp({EEG.event.type},'DIN1') & strcmp({EEG.event.State}, 'AS')));
  1321. total_QS_vep=length(find(strcmp({EEG.event.type},'DIN1') & strcmp({EEG.event.State}, 'QS')));
  1322. total_I_vep=length(find(strcmp({EEG.event.type},'DIN1') & strcmp({EEG.event.State}, 'I')));
  1323. total_W_vep=length(find(strcmp({EEG.event.type},'DIN1') & strcmp({EEG.event.State}, 'W')));
  1324. %% STEP 13: Remove baseline
  1325. EEG = eeg_checkset( EEG );
  1326. baseline_window = [-200,0];
  1327. EEG = pop_rmbase( EEG, baseline_window);
  1328. %% STEP 14: Artifact rejection
  1329. all_bad_epochs=0;
  1330. if voltthres_rejection==1 % check voltage threshold rejection
  1331. if interp_epoch==1 % check epoch level channel interpolation
  1332. chans=[]; chansidx=[];chans_labels2=[];
  1333. chans_labels2=cell(1,EEG.nbchan);
  1334. for i=1:EEG.nbchan
  1335. chans_labels2{i}= EEG.chanlocs(i).labels;
  1336. end
  1337. [chans,chansidx] = ismember(frontal_channels, chans_labels2);
  1338. frontal_channels_idx = chansidx(chansidx ~= 0);
  1339. badChans = zeros(EEG.nbchan, EEG.trials);
  1340. badepoch=zeros(1, EEG.trials);
  1341. if isempty(frontal_channels_idx)==1 % check whether there is any frontal channel in dataset to check
  1342. warning('No frontal channels from the list present in the data. Only epoch interpolation will be performed.');
  1343. else
  1344. % find artifaceted epochs by detecting outlier voltage in the specified channels list and remove epoch if artifacted in those channels
  1345. for ch =1:length(frontal_channels_idx)
  1346. EEG = pop_eegthresh(EEG,1, frontal_channels_idx(ch), volt_threshold(1), volt_threshold(2), EEG.xmin, EEG.xmax,0,0);
  1347. EEG = eeg_checkset( EEG );
  1348. EEG = eeg_rejsuperpose( EEG, 1, 1, 1, 1, 1, 1, 1, 1);
  1349. badChans(ch,:) = EEG.reject.rejglobal;
  1350. end
  1351. for ii=1:size(badChans, 2)
  1352. badepoch(ii)=sum(badChans(:,ii));
  1353. end
  1354. badepoch=logical(badepoch);
  1355. end
  1356. % If all epochs are artifacted, save the dataset and ignore rest of the preprocessing for this subject.
  1357. if sum(badepoch)==EEG.trials || sum(badepoch)+1==EEG.trials
  1358. all_bad_epochs=1;
  1359. warning(['No usable data for datafile', [datafile_names{s} '_Merged']]);
  1360. if output_format==1
  1361. EEG = eeg_checkset(EEG);
  1362. EEG = pop_editset(EEG, 'setname', [datafile_names{s} '_VEP_no_usable_data_all_bad_epoch']);
  1363. EEG = pop_saveset(EEG, 'filename', [datafile_names{s} '_VEP_no_usable_data_all_bad_epoch.set'],'filepath', [output_location filesep 'processed_data' filesep]); % save .set format
  1364. elseif output_format==2
  1365. save([[output_location filesep 'processed_data' filesep] [datafile_names{s} '_VEP_no_usable_data_all_bad_epochs.mat']], 'EEG'); % save .mat format
  1366. end
  1367. else
  1368. EEG = pop_rejepoch( EEG, badepoch, 0);
  1369. EEG = eeg_checkset(EEG);
  1370. end
  1371. if all_bad_epochs==1
  1372. warning(['No usable data for datafile', [datafile_names{s} '_Merged']]);
  1373. else
  1374. % Interpolate artifacted data for all reaming channels
  1375. badChans = zeros(EEG.nbchan, EEG.trials);
  1376. % Find artifacted epochs by detecting outlier voltage but don't remove
  1377. for ch=1:EEG.nbchan
  1378. EEG = pop_eegthresh(EEG,1, ch, volt_threshold(1), volt_threshold(2), EEG.xmin, EEG.xmax,0,0);
  1379. EEG = eeg_checkset(EEG);
  1380. EEG = eeg_rejsuperpose(EEG, 1, 1, 1, 1, 1, 1, 1, 1);
  1381. badChans(ch,:) = EEG.reject.rejglobal;
  1382. end
  1383. tmpData = zeros(EEG.nbchan, EEG.pnts, EEG.trials);
  1384. for e = 1:EEG.trials
  1385. % Initialize variables EEGe and EEGe_interp;
  1386. EEGe = []; EEGe_interp = []; badChanNum = [];
  1387. % Select only this epoch (e)
  1388. EEGe = pop_selectevent( EEG, 'epoch', e, 'deleteevents', 'off', 'deleteepochs', 'on', 'invertepochs', 'off');
  1389. badChanNum = find(badChans(:,e)==1); % find which channels are bad for this epoch
  1390. EEGe_interp = eeg_interp(EEGe,badChanNum); %interpolate the bad channels for this epoch
  1391. tmpData(:,:,e) = EEGe_interp.data; % store interpolated data into matrix
  1392. end
  1393. EEG.data = tmpData; % now that all of the epochs have been interpolated, write the data back to the main file
  1394. % If more than 10% of channels in an epoch were interpolated, reject that epoch
  1395. badepoch=zeros(1, EEG.trials);
  1396. for ei=1:EEG.trials
  1397. NumbadChan = badChans(:,ei); % find how many channels are bad in an epoch
  1398. if sum(NumbadChan) > round((10/100)*EEG.nbchan)% check if more than 10% are bad
  1399. badepoch (ei)= sum(NumbadChan);
  1400. end
  1401. end
  1402. badepoch=logical(badepoch);
  1403. end
  1404. % If all epochs are artifacted, save the dataset and ignore rest of the preprocessing for this subject.
  1405. if sum(badepoch)==EEG.trials || sum(badepoch)+1==EEG.trials
  1406. all_bad_epochs=1;
  1407. warning(['No usable data for datafile', [datafile_names{s} '_Merged']]);
  1408. if output_format==1
  1409. EEG = eeg_checkset(EEG);
  1410. EEG = pop_editset(EEG, 'setname', [datafile_names{s} '_VEP_no_usable_data_all_bad_epochs']);
  1411. EEG = pop_saveset(EEG, 'filename', [datafile_names{s} '_VEP_no_usable_data_all_bad_epochs.set'],'filepath', [output_location filesep 'processed_data' filesep]); % save .set format
  1412. elseif output_format==2
  1413. save([[output_location filesep 'processed_data' filesep] [datafile_names{s} '_VEP_no_usable_data_all_bad_epochs.mat']], 'EEG'); % save .mat format
  1414. end
  1415. else
  1416. EEG = pop_rejepoch(EEG, badepoch, 0);
  1417. EEG = eeg_checkset(EEG);
  1418. end
  1419. else % if no epoch level channel interpolation
  1420. EEG = pop_eegthresh(EEG, 1, (1:EEG.nbchan), volt_threshold(1), volt_threshold(2), EEG.xmin, EEG.xmax, 0, 0);
  1421. EEG = eeg_checkset(EEG);
  1422. EEG = eeg_rejsuperpose( EEG, 1, 1, 1, 1, 1, 1, 1, 1);
  1423. end % end of epoch level channel interpolation if statement
  1424. % If all epochs are artifacted, save the dataset and ignore rest of the preprocessing for this subject.
  1425. if sum(EEG.reject.rejthresh)==EEG.trials || sum(EEG.reject.rejthresh)+1==EEG.trials
  1426. all_bad_epochs=1;
  1427. warning(['No usable data for datafile', [datafile_names{s} '_Merged']]);
  1428. if output_format==1
  1429. EEG = eeg_checkset(EEG);
  1430. EEG = pop_editset(EEG, 'setname', [datafile_names{s} '_VEP_no_usable_data_all_bad_epochs']);
  1431. EEG = pop_saveset(EEG, 'filename', [datafile_names{s} '_VEP_no_usable_data_all_bad_epochs.set'],'filepath', [output_location filesep 'processed_data' filesep]); % save .set format
  1432. elseif output_format==2
  1433. save([[output_location filesep 'processed_data' filesep] [datafile_names{s} '_VEP_no_usable_data_all_bad_epochs.mat']], 'EEG'); % save .mat format
  1434. end
  1435. else
  1436. EEG = pop_rejepoch(EEG,(EEG.reject.rejthresh), 0);
  1437. EEG = eeg_checkset(EEG);
  1438. end
  1439. end % end of voltage threshold rejection if statement
  1440. % if all epochs are found bad during artifact rejection
  1441. if all_bad_epochs==1
  1442. total_epochs_after_artifact_rejection(s)=0;
  1443. total_channels_interpolated(s)=0;
  1444. continue % ignore rest of the processing and go to next datafile
  1445. else
  1446. total_vep_after_artifact_rejection=EEG.trials;
  1447. total_AS_vep_artifact_rejection=length(find(strcmp({EEG.event.type},'DIN1') & strcmp({EEG.event.State}, 'AS')));
  1448. total_QS_vep_artifact_rejection=length(find(strcmp({EEG.event.type},'DIN1') & strcmp({EEG.event.State}, 'QS')));
  1449. total_I_vep_artifact_rejection=length(find(strcmp({EEG.event.type},'DIN1') & strcmp({EEG.event.State}, 'I')));
  1450. total_W_vep_artifact_rejection=length(find(strcmp({EEG.event.type},'DIN1') & strcmp({EEG.event.State}, 'W')));
  1451. end
  1452. %% STEP 15: Interpolate deleted channels
  1453. if interp_channels==1
  1454. EEG = eeg_interp(EEG, channels_analysed);
  1455. EEG = eeg_checkset(EEG);
  1456. end
  1457. if numel(FASTbadChans)==0 && numel(ica_prep_badChans)==0
  1458. total_channels_interpolated(s)=0;
  1459. else
  1460. total_channels_interpolated(s)=numel(FASTbadChans)+ numel(ica_prep_badChans);
  1461. end
  1462. %% STEP 16: Rereference data
  1463. EEG = eeg_checkset(EEG);
  1464. reref = [];
  1465. reref_idx=zeros(1, length(reref));
  1466. for rr=1:length(reref)
  1467. reref_idx(rr)=find(strcmp({EEG.chanlocs.labels}, reref{rr}));
  1468. end
  1469. EEG = eeg_checkset(EEG);
  1470. EEG = pop_reref( EEG, reref_idx);
  1471. %% Save processed data
  1472. if output_format==1
  1473. EEG = eeg_checkset(EEG);
  1474. EEG = pop_editset(EEG, 'setname', [datafile_names{s} '_VEP_processed_data']);
  1475. EEG = pop_saveset(EEG, 'filename', [datafile_names{s} '_VEP_processed_data.set'],'filepath', [output_location filesep 'processed_data' filesep]); % save .set format
  1476. elseif output_format==2
  1477. save([[output_location filesep 'processed_data' filesep] [datafile_names{s} '_VEP_processed_data.mat']], 'EEG'); % save .mat format
  1478. end
  1479. %% Create the report table for all the data files with relevant preprocessing outputs.
  1480. report_table=table(datafile_names(s)', reference_used_for_faster(s)', faster_bad_channels(s)', ica_preparation_bad_channels(s)', length_ica_data(s)', ...
  1481. total_ICs(s)', ICs_removed(s)', total_vep(s)', total_AS_vep(s)',total_QS_vep(s)',total_I_vep(s)',total_W_vep(s)', total_vep_after_artifact_rejection(s)', ...
  1482. total_AS_vep_artifact_rejection(s)', total_QS_vep_artifact_rejection(s)', total_I_vep_artifact_rejection(s)', total_W_vep_artifact_rejection(s)', total_channels_interpolated(s)');
  1483. report_table.Properties.VariableNames={'folder_list_VEP', 'reference_used_for_faster', 'faster_bad_channels', ...
  1484. 'ica_preparation_bad_channels', 'length_ica_data', 'total_ICs', 'ICs_removed', 'total_vep', 'total_AS_vep', 'total_QS_vep', 'total_I_vep', 'total_W_vep', 'total_vep_after_artifact_rejection', ...
  1485. 'total_AS_vep_artifact_rejection', 'total_QS_vep_artifact_rejection', 'total_I_vep_artifact_rejection', 'total_W_vep_artifact_rejection','total_channels_interpolated'};
  1486. writetable(report_table, ['/Volumes/BEHAVEBB/Amy/ADVANCES/EEG/Processed/Report/1m/', subject '_1m_VEP_MADE_report_', datestr(now,'dd-mm-yyyy'),'.csv']);
  1487. end
  1488. end
  1489. end % end of loop

SA-MADE_pipeline_multiple_participants.m at commit 0f6dc49, no license · at the source

Overview

Authors: Huiyu Yang1, Ran Liu2, Katrina R Simon3, Lissete A Gimenez4, Maureen E Bowers5, Nicolò Pini6, Stephanie C Leach7, Leilani Salas3, Lauren C Shuffrey8, William P Fifer6, Julie Herbstman3,9, Nathan A Fox10,11, Amy E Margolis12,13
13 affiliations
  1. Department of Psychology, The Ohio State University, Columbus, OH, United States
  2. Institute of Developmental Psychology, Faculty of Psychology, Beijing Normal University, Beijing, China
  3. Department of Environmental Health Sciences, Mailman School of Public Health, Columbia University Irving Medical Center, New York, NY, United States
  4. Steinhardt School of Culture, Education, and Human Development, New York University, New York, NY, United States
  5. The Pew Charitable Trusts, Washington DC, United States
  6. Division of Developmental Neuroscience, New York State Psychiatric Institute, New York, NY, United States
  7. Department of Psychological & Brain Sciences, University of Iowa, IA, United States
  8. Department of Child and Adolescent Psychiatry, NYU Grossman School of Medicine, New York, NY, United States
  9. Columbia Center for Children’s Environmental Health, Department of Environmental Health Sciences, Mailman School of Public Health, Columbia University, New York, NY, United States
  10. Neuroscience and Cognitive Science Program, University of Maryland, College Park, MD, United States
  11. Department of Human Development and Quantitative Methodology, University of Maryland, College Park, MD, United States
  12. Department of Psychiatry and Behavioral Health and the Clinical and Translational Science Institute, The Ohio State University, Columbus, OH, United States
  13. The Child Mind Institute, New York, NY, United States
Journal: Developmental cognitive neuroscience, volume 79, article 101727
Dates: received 18 May 2025; accepted 17 April 2026; published online 20 April 2026; in print June 2026
Type: Research article · Language: English
License: CC BY-NC-ND
Identifiers: DOI 10.1016/j.dcn.2026.101727 · PMID 42054975 · PMCID PMC13141766 · OpenAlex W7155001382
Open access: gold, a free copy (OpenAlex)
Status: code verified
Categories: EEG (modality), human (organism), developmental (subfield)
Methods: Preprocessing, Smoothing, state filtering, decompositions, Spectral & time-frequency, Evoked potentials, Physiology & signal measures
Keywords: Sleep EEG, Event-related potentials, EEG pipeline, Infant attention, Auditory oddball
MeSH: Brain*, Evoked Potentials, Auditory*, Sleep*, Acoustic Stimulation, Electroencephalography, Female, Humans, Infant, Infant, Newborn, Male (* major topic)
Topic: Neonatal and fetal brain pathology (Pediatrics, Perinatology and Child Health, Medicine), according to OpenAlex
Citations: cited by 1 paper (Europe PMC); 47 references in the paper

Abstract

The abstract is not reproduced here: the paper's license (CC BY-NC-ND) does not allow it. Read it in the paper, at the publisher or on Europe PMC.

Repository

Its files are read in the Code ↔ Paper reader above, with 3 matches between paragraphs and lines of code.

EBBLab/SA-MADE

License: none: the authors keep all their rights
State: the link answers, verified on 29 September 2026
Evidence: files inventoried
Commit: 0f6dc49bc1f48f7e90cb776498579865d9047cd8, 7 October 2025
Languages: MATLAB (3)
Size: 4 files, 3 scripts
Software Heritage: not archived
Found in: the text, “EEG preprocessing”
Holds: README
Not found: license file, CITATION.cff, environment file, tests, continuous integration, documentation
Tools: EEGLAB (2 files)
Availability: 1 check, the latest on 29 September 2026: the link answers
  • 29 September 2026: the link answers
4 files

Tracing map

Proposed by the machine: these links were found in the paper and verified at the source, without human review. The map will receive a Zenodo DOI once one of the paper's authors has validated it with their ORCID.

What the map holds:

  • 1 repository of the authors' code, each at its verified commit, with its license and how the link was found in the paper;
  • 3 scripts, each with its path and the digest of its content;
  • 3 matches between paragraphs of the paper and lines of the code (method lexical-v1);
  • neither the text of the paper nor the code itself.

Its JSON (tracing-map.json) is deposited on Zenodo with its DOI once the map is validated.

Data

No dataset and no data link were found in the paper.

Data availability statement

The paper has a data availability statement. Its license (CC BY-NC-ND) does not allow reproducing it here; in short, from what the harvester recognized in it:

  • it says that the data are available on request

Read it in the paper: doi.org/10.1016/j.dcn.2026.101727.

Versions

The history of this record: each version stored by the harvester or made by a correction of its authors or of the maintainers of its code, and what changed in its facts. The texts of the paper (its abstract, its availability statements) are not part of it; versions that changed only those are not listed.

Version 1, 29 September 2026: the first record

Recorded: type, language, journal, volume, pages, dates, 13 authors, 5 keywords, 10 MeSH terms, 3 funders, 44 references.

Cite

This paper

Yang, H., Liu, R., Simon, K. R., Gimenez, L. A., Bowers, M. E., Pini, N., Leach, S. C., Salas, L., Shuffrey, L. C., Fifer, W. P., Herbstman, J., Fox, N. A., & Margolis, A. E. (2026). Neonatal brain activity across sleep states: Evidence from resting EEG and auditory event-related potentials. Developmental cognitive neuroscience, 79, 101727. https://doi.org/10.1016/j.dcn.2026.101727

BibTeX

@article{yang2026neonatal,
author = {Yang, Huiyu and Liu, Ran and Simon, Katrina R and Gimenez, Lissete A and Bowers, Maureen E and Pini, Nicolò and Leach, Stephanie C and Salas, Leilani and Shuffrey, Lauren C and Fifer, William P and Herbstman, Julie and Fox, Nathan A and Margolis, Amy E},
title = {{Neonatal brain activity across sleep states: Evidence from resting EEG and auditory event-related potentials}},
journal = {Developmental cognitive neuroscience},
year = {2026},
month = apr,
volume = {79},
pages = {101727},
publisher = {Elsevier},
issn = {1878-9293},
doi = {10.1016/j.dcn.2026.101727},
url = {https://doi.org/10.1016/j.dcn.2026.101727},
pmid = {42054975},
pmcid = {PMC13141766}
}

RIS

TY - JOUR
AU - Yang, Huiyu
AU - Liu, Ran
AU - Simon, Katrina R
AU - Gimenez, Lissete A
AU - Bowers, Maureen E
AU - Pini, Nicolò
AU - Leach, Stephanie C
AU - Salas, Leilani
AU - Shuffrey, Lauren C
AU - Fifer, William P
AU - Herbstman, Julie
AU - Fox, Nathan A
AU - Margolis, Amy E
TI - Neonatal brain activity across sleep states: Evidence from resting EEG and auditory event-related potentials
T2 - Developmental cognitive neuroscience
J2 - Dev Cogn Neurosci
PY - 2026
DA - 2026/04/20
VL - 79
SP - 101727
SN - 1878-9293
PB - Elsevier
DO - 10.1016/j.dcn.2026.101727
UR - https://doi.org/10.1016/j.dcn.2026.101727
LA - en
ER -

CSL-JSON

{
"id": "10.1016/j.dcn.2026.101727",
"type": "article-journal",
"title": "Neonatal brain activity across sleep states: Evidence from resting EEG and auditory event-related potentials",
"container-title": "Developmental cognitive neuroscience",
"author": [
{
"family": "Yang",
"given": "Huiyu"
},
{
"family": "Liu",
"given": "Ran"
},
{
"family": "Simon",
"given": "Katrina R"
},
{
"family": "Gimenez",
"given": "Lissete A"
},
{
"family": "Bowers",
"given": "Maureen E"
},
{
"family": "Pini",
"given": "Nicolò"
},
{
"family": "Leach",
"given": "Stephanie C"
},
{
"family": "Salas",
"given": "Leilani"
},
{
"family": "Shuffrey",
"given": "Lauren C"
},
{
"family": "Fifer",
"given": "William P"
},
{
"family": "Herbstman",
"given": "Julie"
},
{
"family": "Fox",
"given": "Nathan A"
},
{
"family": "Margolis",
"given": "Amy E"
}
],
"container-title-short": "Dev Cogn Neurosci",
"volume": "79",
"page": "101727",
"DOI": "10.1016/j.dcn.2026.101727",
"PMID": "42054975",
"PMCID": "PMC13141766",
"ISSN": "1878-9293",
"publisher": "Elsevier",
"URL": "https://doi.org/10.1016/j.dcn.2026.101727",
"language": "en",
"issued": {
"date-parts": [
[
2026,
4,
20
]
]
}
}

The tracing map gets a citation of its own once an author has validated it and it has a DOI.

Similar papers

The papers with a page that share the most with this one: the tools found in their code, their categories, datasets, cited references and authors, the rarest counting most.

[1] doi:10.1002/dev.70128 [code]
Socioeconomic Status, the Home Language Environment, Noise Exposure, and the Mismatch Response in Infancy.
Journal: Developmental psychobiology
In common: developmental, EEG, 6 references, author Katrina R Simon
[2] doi:10.1016/j.dcn.2026.101741
Pediatric resting EEG collection, preprocessing, and analysis: A systematic review.
Journal: Developmental cognitive neuroscience
In common: developmental, EEG, 6 references
[3] doi:10.1093/cercor/bhag077 [code]
The longitudinal development of intrinsic timescales in infancy and their relation to alpha brain rhythm.
Journal: Cerebral cortex (New York, N.Y. : 1991)
In common: EEGLAB, developmental, EEG, 3 references
[4] doi:10.3389/fneur.2026.1791834 [code]
Minimum data requirements and automated preprocessing for reliable EEG biomarkers in Rett syndrome.
Journal: Frontiers in neurology
In common: EEGLAB, EEG, 3 references
[5] doi:10.1111/ejn.70670 [code]
Cross-Night Modulation of Change Detection ERPs to Foreign Speech Sound Features During N2 Sleep.
Journal: The European journal of neuroscience
In common: EEGLAB, EEG, 3 references
[6] doi:10.1093/sleep/zsaf410 [code]
Does twitch-spindle coupling differ between N2 and N3 sleep in 6-month-olds?
Journal: Sleep
In common: developmental, EEG, 3 references
[7] doi:10.7554/elife.107081 [code]
Cortical motor activity modulates respiration and reduces apnoea in neonates.
Journal: eLife
In common: EEGLAB, developmental, EEG, 1 reference
[8] doi:10.1016/j.isci.2026.115878 [code]
Twitching in sleeping premature infants provides a sensitive behavioral assay of early motor control.
Journal: iScience
In common: developmental, 2 references
[9] doi:10.1038/s41562-026-02533-1 [code]
Fluctuations in arousal reflect latent state transitions that facilitate behavioural optimization.
Journal: Nature human behaviour
In common: EEGLAB, EEG, 2 references
[10] doi:10.1038/s41598-026-47785-z [code]
Modulations of the P3b effect as a function of bilingual language experience.
Journal: Scientific reports
In common: EEGLAB, EEG, 2 references

Contribute

The authors of this paper can claim it, correct its record and validate its tracing map, and the maintainers of its code (its owner, or a public member of its organization) correct what it says of their repository; anyone signed in can ask for its removal. Every request goes to OSCR's own machine, which answers it; your account page follows them.

Sign in with ORCID to claim this paper as one of its authors, correct its record or validate its tracing map: when the paper's metadata lists your ORCID iD, you are recognized at once. Maintainers of its code: sign in with GitHub, then claim the repository on your account page.

Request its removal

To ask OSCR to remove this record, the copies of its authors' scripts or its tracing map, use the removal request page: signed in, you say who you are, what to remove and why, then review and confirm the request. Published rules decide every request (how).

Discussion, reproductions, activity

Discussion: questions and error reports about this paper and its code, from signed-in readers and its authors. It opens with sign-in.

Reproductions: reports from readers who ran the authors' code: what they reproduced, with which environment, commit and data. It opens with sign-in.

Activity: what happens around this paper: new versions of its record, its map's validation, discussions and reproductions. It opens with sign-in.