Tutorial 28: Scripting

[TUTORIAL UNDER DEVELOPMENT: NOT READY FOR PUBLIC USE]

Authors: Francois Tadel, Elizabeth Bock, Sylvain Baillet

The previous tutorials explained how to use Brainstorm in an interactive way to process one subject with two acquisition runs. In the context of a typical neuroimaging study, you may have tens or hundreds of subjects to process in the same way, it is unrealistic to do everything manually. Some parts of the analysis can be processed in batches with no direct supervision, others require more attention. This tutorial introduces tools and tricks that will help you assemble an efficient analysis pipeline.

How to process many subjects

This section proposes a standard workflow for processing a full group study with Brainstorm. It contains the same steps of analysis as the introduction tutorials, but separating what can be done automatically from what should be done manually. This workflow can be adapted to most ERP studies (stimulus-based).

Script generation

http://neuroimage.usc.edu/brainstorm/Tutorials/PipelineEditor

Script edition

Add loops, load files, ...

Loops: http://neuroimage.usc.edu/forums/showthread.php?2429-Problem-using-tags

File manipulation

Final script

The following script from the Brainstorm distribution reproduces all the introduction tutorials: brainstorm3/toolbox/script/tutorial_introduction.m

1 function tutorial_introduction(tutorial_dir, reports_dir) 2 % TUTORIAL_INTRODUCTION: Script that runs all the Brainstorm introduction tutorials. 3 % 4 % INPUTS: 5 % - tutorial_dir : Directory where the sample_introduction.zip file has been unzipped 6 % - reports_dir : Directory where to save the execution report (instead of displaying it) 7 8 % @============================================================================= 9 % This function is part of the Brainstorm software: 10 % https://neuroimage.usc.edu/brainstorm 11 % 12 % Copyright (c) University of Southern California & McGill University 13 % This software is distributed under the terms of the GNU General Public License 14 % as published by the Free Software Foundation. Further details on the GPLv3 15 % license can be found at http://www.gnu.org/copyleft/gpl.html. 16 % 17 % FOR RESEARCH PURPOSES ONLY. THE SOFTWARE IS PROVIDED "AS IS," AND THE 18 % UNIVERSITY OF SOUTHERN CALIFORNIA AND ITS COLLABORATORS DO NOT MAKE ANY 19 % WARRANTY, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO WARRANTIES OF 20 % MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, NOR DO THEY ASSUME ANY 21 % LIABILITY OR RESPONSIBILITY FOR THE USE OF THIS SOFTWARE. 22 % 23 % For more information type "brainstorm license" at command prompt. 24 % =============================================================================@ 25 % 26 % Author: Francois Tadel, 2016-2017 27 28 29 % ===== FILES TO IMPORT ===== 30 % Output folder for reports 31 if (nargin < 2) || isempty(reports_dir) || ~isdir(reports_dir) 32 reports_dir = []; 33 end 34 % You have to specify the folder in which the tutorial dataset is unzipped 35 if (nargin == 0) || isempty(tutorial_dir) || ~file_exist(tutorial_dir) 36 error('The first argument must be the full path to the dataset folder.'); 37 end 38 % Subject name 39 SubjectName = 'Subject01'; 40 % Build the path of the files to import 41 AnatDir = fullfile(tutorial_dir, 'sample_introduction', 'anatomy'); 42 Run1File = fullfile(tutorial_dir, 'sample_introduction', 'data', 'S01_AEF_20131218_01_600Hz.ds'); 43 Run2File = fullfile(tutorial_dir, 'sample_introduction', 'data', 'S01_AEF_20131218_02_600Hz.ds'); 44 NoiseFile = fullfile(tutorial_dir, 'sample_introduction', 'data', 'S01_Noise_20131218_02_600Hz.ds'); 45 % Check if the folder contains the required files 46 if ~file_exist(Run1File) 47 error(['The folder ' tutorial_dir ' does not contain the folder from the file sample_introduction.zip.']); 48 end 49 % Re-inialize random number generator 50 if (bst_get('MatlabVersion') >= 712) 51 rng('default'); 52 end 53 54 55 %% ===== TUTORIAL #1: CREATE PROTOCOL ================================================ 56 % =================================================================================== 57 disp([10 'DEMO> Tutorial #1: Create protocol' 10]); 58 % The protocol name has to be a valid folder name (no spaces, no weird characters...) 59 ProtocolName = 'TutorialIntroduction'; 60 % Start brainstorm without the GUI 61 if ~brainstorm('status') 62 brainstorm nogui 63 end 64 % Delete existing protocol 65 gui_brainstorm('DeleteProtocol', ProtocolName); 66 % Create new protocol 67 gui_brainstorm('CreateProtocol', ProtocolName, 0, 0); 68 % Start a new report 69 bst_report('Start'); 70 % Reset colormaps 71 bst_colormaps('RestoreDefaults', 'meg'); 72 % Set the current display mode to 'butterfly' 73 bst_set('TSDisplayMode', 'butterfly'); 74 75 76 %% ===== TUTORIAL #2: IMPORT ANATOMY ================================================= 77 % =================================================================================== 78 disp([10 'DEMO> Tutorial #2: Import anatomy' 10]); 79 % Process: Import FreeSurfer folder 80 bst_process('CallProcess', 'process_import_anatomy', [], [], ... 81 'subjectname', SubjectName, ... 82 'mrifile', {AnatDir, 'FreeSurfer'}, ... 83 'nvertices', 15000, ... 84 'nas', [127, 213, 139], ... 85 'lpa', [ 52, 113, 96], ... 86 'rpa', [202, 113, 91]); 87 % This automatically calls the SPM registration procedure because the AC/PC/IH points are not defined 88 89 90 91 %% ===== TUTORIAL #3: EXPLORE ANATOMY ================================================ 92 % =================================================================================== 93 disp([10 'DEMO> Tutorial #3: Explore anatomy' 10]); 94 % Get subject definition 95 sSubject = bst_get('Subject', SubjectName); 96 % Get MRI file and surface files 97 MriFile = sSubject.Anatomy(sSubject.iAnatomy).FileName; 98 CortexFile = sSubject.Surface(sSubject.iCortex).FileName; 99 HeadFile = sSubject.Surface(sSubject.iScalp).FileName; 100 % Display MRI 101 hFigMri1 = view_mri(MriFile); 102 hFigMri3 = view_mri_3d(MriFile, [], [], 'NewFigure'); 103 hFigMri2 = view_mri_slices(MriFile, 'x', 20); 104 pause(0.5); 105 % Close figures 106 close([hFigMri1 hFigMri2 hFigMri3]); 107 % Display scalp and cortex 108 hFigSurf = view_surface(HeadFile); 109 hFigSurf = view_surface(CortexFile, [], [], hFigSurf); 110 hFigMriSurf = view_mri(MriFile, CortexFile); 111 % Figure configuration 112 iTess = 2; 113 panel_surface('SetShowSulci', hFigSurf, iTess, 1); 114 panel_surface('SetSurfaceColor', hFigSurf, iTess, [1 0 0]); 115 panel_surface('SetSurfaceSmooth', hFigSurf, iTess, 0.5, 0); 116 panel_surface('SetSurfaceTransparency', hFigSurf, iTess, 0.8); 117 figure_3d('SetStandardView', hFigSurf, 'left'); 118 pause(0.5); 119 % Close figures 120 close([hFigSurf hFigMriSurf]); 121 122 123 124 %% ===== TUTORIAL #4: CHANNEL FILE =================================================== 125 % =================================================================================== 126 disp([10 'DEMO> Tutorial #4: Channel file' 10]); 127 % Process: Create link to raw files 128 sFilesRun1 = bst_process('CallProcess', 'process_import_data_raw', [], [], ... 129 'subjectname', SubjectName, ... 130 'datafile', {Run1File, 'CTF'}, ... 131 'channelalign', 1); 132 sFilesRun2 = bst_process('CallProcess', 'process_import_data_raw', [], [], ... 133 'subjectname', SubjectName, ... 134 'datafile', {Run2File, 'CTF'}, ... 135 'channelalign', 1); 136 sFilesNoise = bst_process('CallProcess', 'process_import_data_raw', [], [], ... 137 'subjectname', SubjectName, ... 138 'datafile', {NoiseFile, 'CTF'}, ... 139 'channelalign', 0); 140 sFilesRaw = [sFilesRun1, sFilesRun2, sFilesNoise]; 141 % Process: Snapshot: Sensors/MRI registration 142 bst_process('CallProcess', 'process_snapshot', [sFilesRun1, sFilesRun2], [], ... 143 'target', 1, ... % Sensors/MRI registration 144 'modality', 1, ... % MEG (All) 145 'orient', 1, ... % left 146 'Comment', 'MEG/MRI Registration'); 147 148 % View sensors 149 hFig = view_surface(HeadFile); 150 hFig = view_channels(sFilesRun1.ChannelFile, 'MEG', 1, 1, hFig); 151 % Hide sensors 152 pause(0.5); 153 hFig = view_channels(sFilesRun1.ChannelFile, 'MEG', 0, 0, hFig); 154 % View coils 155 hFig = view_channels(sFilesRun1.ChannelFile, 'CTF', 1, 1, hFig); 156 % View helmet 157 pause(0.5); 158 hFig = view_helmet(sFilesRun1.ChannelFile, hFig); 159 pause(0.5); 160 close(hFig); 161 % Edit good/bad channel for current file 162 gui_edit_channel(sFilesRun1.ChannelFile); 163 pause(0.5); 164 % Unload everything 165 bst_memory('UnloadAll', 'Forced'); 166 167 168 169 %% ===== TUTORIAL #5: REVIEW RAW ===================================================== 170 % =================================================================================== 171 disp([10 'DEMO> Tutorial #5: Review raw' 10]); 172 % Process: Convert to continuous (CTF): Continuous 173 bst_process('CallProcess', 'process_ctf_convert', sFilesRaw, [], ... 174 'rectype', 2); % Continuous 175 176 % View recordings 177 hFigMeg = view_timeseries(sFilesRun1.FileName, 'MEG'); 178 hFigEeg = view_timeseries(sFilesRun1.FileName, 'Misc'); 179 hFigSel = view_timeseries(sFilesRun1.FileName, 'MEG', {'MLT11','MLT12','MLT13'}); 180 % Figure configuration 181 pause(0.5); 182 panel_record('SetTimeLength', 3); 183 panel_record('SetStartTime', 100); 184 panel_record('SetDisplayMode', hFigMeg, 'column'); 185 panel_montage('SetCurrentMontage', hFigMeg, 'CTF LT'); 186 % Set filters: panel_filter('SetFilters', LowPassEnabled, LowPassValue, HighPassEnabled, HighPassValue, SinRemovalEnabled, SinRemovalValue, MirrorEnabled, FullSourcesEnabled) 187 panel_filter('SetFilters', 1, 100, 1, 1, 0, [], 0, 0); 188 pause(0.5); 189 panel_record('SetDisplayMode', hFigMeg, 'butterfly'); 190 panel_montage('SetCurrentMontage', hFigMeg, ''); 191 % Close figures 192 close([hFigMeg hFigEeg hFigSel]); 193 194 195 196 %% ===== TUTORIAL #8: STIM DELAYS ==================================================== 197 % =================================================================================== 198 disp([10 'DEMO> Tutorial #8: Stim delays' 10]); 199 % Process: Detect: standard_fix 200 bst_process('CallProcess', 'process_evt_detect_analog', [sFilesRun1, sFilesRun2], [], ... 201 'eventname', 'standard_fix', ... 202 'channelname', 'UADC001', ... 203 'timewindow', [], ... 204 'threshold', 1.2, ... 205 'blanking', 0.2, ... 206 'highpass', 0, ... 207 'lowpass', 0, ... 208 'refevent', 'standard', ... 209 'isfalling', 0, ... 210 'ispullup', 0, ... 211 'isclassify', 0); 212 % Process: Detect: deviant_fix 213 bst_process('CallProcess', 'process_evt_detect_analog', [sFilesRun1, sFilesRun2], [], ... 214 'eventname', 'deviant_fix', ... 215 'channelname', 'UADC001', ... 216 'timewindow', [], ... 217 'threshold', 1.2, ... 218 'blanking', 0.2, ... 219 'highpass', 0, ... 220 'lowpass', 0, ... 221 'refevent', 'deviant', ... 222 'isfalling', 0, ... 223 'ispullup', 0, ... 224 'isclassify', 0); 225 % Process: Read from channel 226 bst_process('CallProcess', 'process_evt_read', [sFilesRun1, sFilesRun2], [], ... 227 'stimchan', 'UDIO001', ... 228 'trackmode', 1, ... % Value: detect the changes of channel value 229 'zero', 0); 230 231 % Process: Delete events 232 bst_process('CallProcess', 'process_evt_delete', [sFilesRun1, sFilesRun2], [], ... 233 'eventname', 'standard, deviant, button'); 234 % Process: Rename event (standard_fix>standard) 235 bst_process('CallProcess', 'process_evt_rename', [sFilesRun1, sFilesRun2], [], ... 236 'src', 'standard_fix', ... 237 'dest', 'standard'); 238 % Process: Rename event (deviant_fix>deviant) 239 bst_process('CallProcess', 'process_evt_rename', [sFilesRun1, sFilesRun2], [], ... 240 'src', 'deviant_fix', ... 241 'dest', 'deviant'); 242 % Process: Rename event (64>button) 243 bst_process('CallProcess', 'process_evt_rename', [sFilesRun1, sFilesRun2], [], ... 244 'src', '64', ... 245 'dest', 'button'); 246 247 248 249 %% ===== TUTORIAL #10: FREQUENCY FILTERS ============================================= 250 % =================================================================================== 251 disp([10 'DEMO> Tutorial #10: Frequency filters' 10]); 252 % Process: Sinusoid removal: 60Hz 120Hz 180Hz 300Hz 253 sFilesNotch = bst_process('CallProcess', 'process_notch', sFilesRaw, [], ... 254 'freqlist', [60, 120, 180], ... 255 'sensortypes', 'MEG', ... 256 'read_all', 0); 257 % Process: Power spectrum density (Welch) 258 sFilesPsd = bst_process('CallProcess', 'process_psd', [sFilesRaw, sFilesNotch], [], ... 259 'timewindow', [], ... 260 'win_length', 4, ... 261 'win_overlap', 50, ... 262 'clusters', {}, ... 263 'sensortypes', 'MEG', ... 264 'edit', struct(... 265 'Comment', 'Power', ... 266 'TimeBands', [], ... 267 'Freqs', [], ... 268 'ClusterFuncTime', 'none', ... 269 'Measure', 'power', ... 270 'Output', 'all', ... 271 'SaveKernel', 0)); 272 % Process: Snapshot: Frequency spectrum 273 bst_process('CallProcess', 'process_snapshot', sFilesPsd, [], ... 274 'target', 10, ... % Frequency spectrum 275 'modality', 1, ... % MEG (All) 276 'Comment', 'Power spectrum density'); 277 % Process: Delete folders 278 bst_process('CallProcess', 'process_delete', sFilesRaw, [], ... 279 'target', 2); % Delete folders 280 % Separate the three outputs 281 sFilesRun1 = {sFilesNotch(1).FileName}; 282 sFilesRun2 = {sFilesNotch(2).FileName}; 283 sFilesNoise = {sFilesNotch(3).FileName}; 284 285 286 287 %% ===== TUTORIAL #11: BAD CHANNELS ================================================== 288 % =================================================================================== 289 % % Process: Set bad channels 290 % sFiles = bst_process('CallProcess', 'process_channel_setbad', sFilesRun2, [], ... 291 % 'sensortypes', 'MRT51, MLO52, MLO42, MLO43'); 292 293 294 295 %% ===== TUTORIAL #12: ARTIFACTS DETECTION =========================================== 296 % =================================================================================== 297 disp([10 'DEMO> Tutorial #12: Artifacts detection' 10]); 298 % Process: Detect heartbeats 299 bst_process('CallProcess', 'process_evt_detect_ecg', [sFilesRun1, sFilesRun2], [], ... 300 'channelname', 'ECG', ... 301 'timewindow', [], ... 302 'eventname', 'cardiac'); 303 % Process: Detect eye blinks 304 bst_process('CallProcess', 'process_evt_detect_eog', [sFilesRun1, sFilesRun2], [], ... 305 'channelname', 'VEOG', ... 306 'timewindow', [], ... 307 'eventname', 'blink'); 308 % Process: Remove simultaneous 309 bst_process('CallProcess', 'process_evt_remove_simult', [sFilesRun1, sFilesRun2], [], ... 310 'remove', 'cardiac', ... 311 'target', 'blink', ... 312 'dt', 0.25, ... 313 'rename', 0); 314 315 316 317 %% ===== TUTORIAL #13: SSP =========================================================== 318 % =================================================================================== 319 disp([10 'DEMO> Tutorial #13: SSP' 10]); 320 % Process: SSP ECG: cardiac 321 bst_process('CallProcess', 'process_ssp_ecg', [sFilesRun1, sFilesRun2], [], ... 322 'eventname', 'cardiac', ... 323 'sensortypes', 'MEG', ... 324 'usessp', 0, ... 325 'select', 1); 326 % Process: SSP EOG: blink 327 bst_process('CallProcess', 'process_ssp_eog', sFilesRun1, [], ... 328 'eventname', 'blink', ... 329 'sensortypes', 'MEG', ... 330 'usessp', 0, ... 331 'select', [1 2]); 332 bst_process('CallProcess', 'process_ssp_eog', sFilesRun2, [], ... 333 'eventname', 'blink', ... 334 'sensortypes', 'MEG', ... 335 'usessp', 0, ... 336 'select', 1); 337 338 339 %% ===== TUTORIAL #14: BAD SEGMENTS ================================================== 340 % =================================================================================== 341 disp([10 'DEMO> Tutorial #14: Bad segments' 10]); 342 % Process: Detect other artifacts 343 bst_process('CallProcess', 'process_evt_detect_badsegment', [sFilesRun1, sFilesRun2], [], ... 344 'timewindow', [], ... 345 'sensortypes', 'MEG', ... 346 'threshold', 3, ... % 3 347 'isLowFreq', 1, ... 348 'isHighFreq', 1); 349 350 % Process: Rename event (1-7Hz > saccade) (Run02 only) 351 bst_process('CallProcess', 'process_evt_rename', sFilesRun2, [], ... 352 'src', '1-7Hz', ... 353 'dest', 'saccade'); 354 355 % Manual selection of saccades (cannot be done from the pipeline editor: manual edition of the structures) 356 sMatRun2 = in_bst_data(sFilesRun2{1}, 'F'); 357 iEvtSaccade = find(strcmpi({sMatRun2.F.events.label}, 'saccade')); 358 sMatRun2.F.events(iEvtSaccade).times = [30, 81.5, 104, 142.5, 167, 187.5, 246.5, 319; 31, 83, 105, 144, 168, 188.5, 248, 320]; 359 sMatRun2.F.events(iEvtSaccade).epochs = ones(1, size(sMatRun2.F.events(iEvtSaccade).times, 2)); 360 sMatRun2.F.events(iEvtSaccade).channels = []; 361 sMatRun2.F.events(iEvtSaccade).notes = []; 362 bst_save(file_fullpath(sFilesRun2{1}), sMatRun2, 'v6', 1); 363 364 % Process: SSP: saccade (Run02 only) 365 bst_process('CallProcess', 'process_ssp', sFilesRun2, [], ... 366 'timewindow', [], ... 367 'eventname', 'saccade', ... 368 'eventtime', [-0.2, 0.2], ... 369 'bandpass', [1.5, 7], ... 370 'sensortypes', 'MEG', ... 371 'usessp', 1, ... 372 'saveerp', 0, ... 373 'method', 1, ... 374 'select', 1); 375 % Process: Detect other artifacts (Run02 only) 376 bst_process('CallProcess', 'process_evt_detect_badsegment', sFilesRun2, [], ... 377 'timewindow', [], ... 378 'sensortypes', 'MEG', ... 379 'threshold', 3, ... % 3 380 'isLowFreq', 1, ... 381 'isHighFreq', 1); 382 383 % Process: Rename event (1-7Hz > bad_1-7Hz) 384 bst_process('CallProcess', 'process_evt_rename', [sFilesRun1, sFilesRun2], [], ... 385 'src', '1-7Hz', ... 386 'dest', 'bad_1-7Hz'); 387 % Process: Rename event (40-240Hz > bad_40-240Hz) 388 bst_process('CallProcess', 'process_evt_rename', [sFilesRun1, sFilesRun2], [], ... 389 'src', '40-240Hz', ... 390 'dest', 'bad_40-240Hz'); 391 % Process: Snapshot: SSP projectors 392 bst_process('CallProcess', 'process_snapshot', [sFilesRun1, sFilesRun2], [], ... 393 'target', 2, ... % SSP projectors 394 'Comment', 'SSP projectors'); 395 396 397 398 %% ===== TUTORIAL #15: IMPORT EVENTS ================================================= 399 % =================================================================================== 400 disp([10 'DEMO> Tutorial #15: Import events' 10]); 401 % Process: Import MEG/EEG: Events (Run01) 402 sFilesEpochs1 = bst_process('CallProcess', 'process_import_data_event', sFilesRun1, [], ... 403 'subjectname', SubjectName, ... 404 'condition', '', ... 405 'eventname', 'standard, deviant', ... 406 'timewindow', [], ... 407 'epochtime', [-0.100, 0.500], ... 408 'createcond', 0, ... 409 'ignoreshort', 1, ... 410 'usectfcomp', 1, ... 411 'usessp', 1, ... 412 'freq', [], ... 413 'baseline', [-0.1, -0.0017]); 414 % Process: Import MEG/EEG: Events (Run02) 415 sFilesEpochs2 = bst_process('CallProcess', 'process_import_data_event', sFilesRun2, [], ... 416 'subjectname', SubjectName, ... 417 'condition', '', ... 418 'eventname', 'standard, deviant', ... 419 'timewindow', [], ... 420 'epochtime', [-0.100, 0.500], ... 421 'createcond', 0, ... 422 'ignoreshort', 1, ... 423 'usectfcomp', 1, ... 424 'usessp', 1, ... 425 'freq', [], ... 426 'baseline', [-0.1, -0.0017]); 427 % Display raster plot 428 hFigRaster = view_erpimage({sFilesEpochs1.FileName}, 'erpimage', 'MEG'); 429 panel_display(); 430 bst_report('Snapshot', hFigRaster, sFilesEpochs1(1).FileName, 'ERP image'); 431 close(hFigRaster); 432 433 434 %% ===== TUTORIAL #16: AVERAGE ======================================================= 435 % =================================================================================== 436 disp([10 'DEMO> Tutorial #16: Average' 10]); 437 % Process: Average: By trial group (folder average) 438 sFilesAvg = bst_process('CallProcess', 'process_average', [sFilesEpochs1, sFilesEpochs2], [], ... 439 'avgtype', 5, ... % By trial groups (folder average) 440 'avg_func', 1, ... % Arithmetic average: mean(x) 441 'weighted', 0, ... 442 'keepevents', 1); 443 % Process: Delete events 'cardiac' 444 bst_process('CallProcess', 'process_evt_delete', sFilesAvg, [], ... 445 'eventname', 'cardiac'); 446 % Process: Snapshot: Recordings time series 447 bst_process('CallProcess', 'process_snapshot', sFilesAvg, [], ... 448 'target', 5, ... % Recordings time series 449 'modality', 1, ... % MEG (All) 450 'Comment', 'Evoked response'); 451 % Set colormap: global color scale 452 bst_colormaps('SetMaxMode', 'meg', 'global'); 453 % Process: Snapshot: Recordings topography (contact sheet) 454 bst_process('CallProcess', 'process_snapshot', sFilesAvg, [], ... 455 'target', 7, ... % Recordings topography (contact sheet) 456 'modality', 1, ... % MEG 457 'contact_time', [0, 0.350], ... 458 'contact_nimage', 15, ... 459 'Comment', 'Evoked response'); 460 461 % Process: Average+Stderr: By trial group (subject average) 462 sFilesAvgAll = bst_process('CallProcess', 'process_average', [sFilesEpochs1, sFilesEpochs2], [], ... 463 'avgtype', 6, ... % By trial group (subject average) 464 'avg_func', 7, ... % Arithmetic average + Standard error 465 'weighted', 0, ... 466 'keepevents', 1); 467 % Process: Delete events 'cardiac' 468 bst_process('CallProcess', 'process_evt_delete', sFilesAvgAll, [], ... 469 'eventname', 'cardiac'); 470 % Process: Delete events 'saccade' 471 bst_process('CallProcess', 'process_evt_delete', sFilesAvgAll, [], ... 472 'eventname', 'saccade'); 473 % Process: Snapshot: Recordings time series 474 bst_process('CallProcess', 'process_snapshot', sFilesAvgAll, [], ... 475 'target', 5, ... % Recordings time series 476 'modality', 1, ... % MEG (All) 477 'Comment', 'Evoked response'); 478 479 480 %% ===== TUTORIAL #17: EXPLORATION =================================================== 481 % =================================================================================== 482 disp([10 'DEMO> Tutorial #17: Bad segments' 10]); 483 % View averages 484 hFigMeg1 = view_timeseries(sFilesAvg(1).FileName, 'MEG'); 485 hFigMeg2 = view_timeseries(sFilesAvg(2).FileName, 'MEG'); 486 hFigEeg1 = view_timeseries(sFilesAvg(1).FileName, 'Misc'); 487 hFigEeg2 = view_timeseries(sFilesAvg(2).FileName, 'Misc'); 488 hFigTopo1 = view_topography(sFilesAvg(1).FileName, 'MEG', '2DSensorCap'); 489 hFigTopo2 = view_topography(sFilesAvg(2).FileName, 'MEG', '2DSensorCap'); 490 hFigTp2 = view_topography(sFilesAvg(3).FileName, 'MEG', '3DSensorCap'); 491 hFigTp3 = view_topography(sFilesAvg(3).FileName, 'MEG', '2DDisc'); 492 hFigTp4 = view_topography(sFilesAvg(3).FileName, 'MEG', '2DLayout'); 493 % Set time: 90ms 494 panel_time('SetCurrentTime', 0.090); 495 % Set filters: 40Hz low-pass, no high-pass 496 panel_filter('SetFilters', 1, 40, 0, [], 0, [], 0, 0); 497 % View selected sensors 498 SelectedChannels = {'MLC31', 'MLC32'}; 499 bst_figures('SetSelectedRows', SelectedChannels); 500 view_timeseries(sFilesAvg(4).FileName, [], SelectedChannels); 501 % Select time window 502 figure_timeseries('SetTimeSelectionManual', hFigMeg1, [0.070, 0.130]); 503 % Show sensors on 2DSensorCap topography 504 isMarkers = 1; 505 isLabels = 0; 506 figure_3d('ViewSensors', hFigTopo1, isMarkers, isLabels); 507 % Display time contact sheet for a figure 508 pause(0.5); 509 hContactFig = view_contactsheet( hFigTopo2, 'time', 'fig', [], 12, [0 0.120] ); 510 pause(0.5); 511 close(hContactFig); 512 513 514 515 %% ===== TUTORIAL #18: COLORMAPS ===================================================== 516 % =================================================================================== 517 disp([10 'DEMO> Tutorial #18: Colormaps' 10]); 518 % Set 'Meg' colormap to 'jet' 519 bst_colormaps('SetColormapName', 'meg', 'jet'); 520 pause(0.5); 521 % Set 'Meg' colormap to 'rwb' 522 bst_colormaps('SetColormapName', 'meg', 'cmap_rbw'); 523 % Set colormap to display absolute values 524 bst_colormaps('SetColormapAbsolute', 'meg', 1); 525 % Normalize colormap for each time frame 526 bst_colormaps('SetMaxMode', 'meg', 'local'); 527 % Hide colorbar 528 bst_colormaps('SetDisplayColorbar', 'meg', 0); 529 pause(0.5); 530 % Restore colormap to default values 531 bst_colormaps('RestoreDefaults', 'meg'); 532 % Edit good/bad channel for current file 533 gui_edit_channelflag(sFilesAvg(1).FileName); 534 % Close figures 535 pause(0.5); 536 bst_memory('UnloadAll', 'Forced'); 537 538 539 540 %% ===== TUTORIAL #20: HEAD MODEL ==================================================== 541 % =================================================================================== 542 disp([10 'DEMO> Tutorial #20: Head model' 10]); 543 % Process: Compute head model 544 bst_process('CallProcess', 'process_headmodel', sFilesAvg, [], ... 545 'comment', '', ... 546 'sourcespace', 1, ... 547 'meg', 3); % Overlapping spheres 548 % Get study structure 549 sStudy = bst_get('Study', sFilesAvg(1).iStudy); 550 % Show spheres 551 hFig = view_spheres(sStudy.HeadModel(sStudy.iHeadModel).FileName, sStudy.Channel.FileName, sSubject); 552 pause(0.5); 553 close(hFig); 554 555 556 557 %% ===== TUTORIAL #21: NOISE COVARIANCE ============================================== 558 % =================================================================================== 559 disp([10 'DEMO> Tutorial #21: Noise covariance' 10]); 560 % Process: Compute covariance (noise or data) 561 bst_process('CallProcess', 'process_noisecov', sFilesNoise, [], ... 562 'baseline', [], ... 563 'sensortypes', 'MEG, EEG, SEEG, ECOG', ... 564 'target', 1, ... % Noise covariance (covariance over baseline time window) 565 'dcoffset', 1, ... % Block by block, to avoid effects of slow shifts in data 566 'identity', 0, ... 567 'copycond', 1, ... 568 'copysubj', 0, ... 569 'replacefile', 1); % Replace 570 % Process: Snapshot: Noise covariance 571 bst_process('CallProcess', 'process_snapshot', sFilesNoise, [], ... 572 'target', 3, ... % Noise covariance 573 'Comment', 'Noise covariance'); 574 575 % Process: Compute covariance (noise or data) [Run01] 576 bst_process('CallProcess', 'process_noisecov', sFilesEpochs1, [], ... 577 'baseline', [-0.1, -0.0017], ... 578 'datatimewindow', [0, 0.5], ... 579 'sensortypes', 'MEG, EEG, SEEG, ECOG', ... 580 'target', 2, ... % Data covariance (covariance over data time window) 581 'dcoffset', 1, ... % Block by block, to avoid effects of slow shifts in data 582 'identity', 0, ... 583 'copycond', 0, ... 584 'copysubj', 0, ... 585 'replacefile', 1); % Replace 586 % Process: Compute covariance (noise or data) [Run02] 587 bst_process('CallProcess', 'process_noisecov', sFilesEpochs2, [], ... 588 'baseline', [-0.1, -0.0017], ... 589 'datatimewindow', [0, 0.5], ... 590 'sensortypes', 'MEG, EEG, SEEG, ECOG', ... 591 'target', 2, ... % Data covariance (covariance over data time window) 592 'dcoffset', 1, ... % Block by block, to avoid effects of slow shifts in data 593 'identity', 0, ... 594 'copycond', 0, ... 595 'copysubj', 0, ... 596 'replacefile', 1); % Replace 597 % Process: Snapshot: Data covariance 598 bst_process('CallProcess', 'process_snapshot', [sFilesEpochs1(1), sFilesEpochs2(1)], [], ... 599 'target', 12, ... % Data covariance 600 'Comment', 'Data covariance'); 601 602 603 604 %% ===== TUTORIAL #22: SOURCE ESTIMATION ============================================= 605 % =================================================================================== 606 disp([10 'DEMO> Tutorial #22: Source estimation' 10]); 607 % === GET DEVIANT AVERAGE RUN01 === 608 % Process: Select recordings in: Subject01/S01_AEF_20131218_01_600Hz_notch 609 sFiles01 = bst_process('CallProcess', 'process_select_files_data', [], [], ... 610 'subjectname', SubjectName, ... 611 'condition', 'S01_AEF_20131218_01_600Hz_notch', ... 612 'includebad', 0); 613 % Process: Select file comments with tag: deviant 614 sFilesAvgDeviant01 = bst_process('CallProcess', 'process_select_tag', sFiles01, [], ... 615 'tag', 'Avg: deviant', ... 616 'search', 2, ... % Search the file comments 617 'select', 1); % Select only the files with the tag 618 619 % === CONSTRAINED EXAMPLE === 620 % Minimum norm options 621 InverseOptions = struct(... 622 'Comment', 'MN: MEG', ... 623 'InverseMethod', 'minnorm', ... 624 'InverseMeasure', 'amplitude', ... 625 'SourceOrient', {{'fixed'}}, ... 626 'Loose', 0.2, ... 627 'UseDepth', 1, ... 628 'WeightExp', 0.5, ... 629 'WeightLimit', 10, ... 630 'NoiseMethod', 'reg', ... 631 'NoiseReg', 0.1, ... 632 'SnrMethod', 'fixed', ... 633 'SnrRms', 0.001, ... 634 'SnrFixed', 3, ... 635 'ComputeKernel', 1, ... 636 'DataTypes', {{'MEG'}}); 637 % Process: Compute sources [2018] 638 sFilesSrcDeviant01 = bst_process('CallProcess', 'process_inverse_2018', sFilesAvgDeviant01, [], ... 639 'output', 2, ... % Kernel only: one per file 640 'inverse', InverseOptions); 641 642 % === DISPLAY SOURCES MANUALLY === 643 % View time series 644 hFigSrc1 = view_timeseries(sFilesAvgDeviant01(1).FileName, 'MEG'); 645 % View on the cortex surface 646 hFigSrc2 = script_view_sources(sFilesSrcDeviant01.FileName, 'cortex'); 647 % Set current time to 90ms 648 panel_time('SetCurrentTime', 0.090); 649 % Set orientation 650 figure_3d('SetStandardView', hFigSrc2, 'left'); 651 % Set surface threshold 652 iSurf = 1; 653 thresh = .30; 654 panel_surface('SetDataThreshold', hFigSrc2, iSurf, thresh); 655 % Set surface smoothing 656 panel_surface('SetSurfaceSmooth', hFigSrc2, iSurf, .6, 0); 657 % Show sulci 658 panel_surface('SetShowSulci', hFigSrc2, iSurf, 1); 659 660 % View sources on MRI (3D orthogonal slices) 661 hFigSrc3 = script_view_sources(sFilesSrcDeviant01.FileName, 'mri3d'); 662 panel_surface('SetDataThreshold', hFigSrc3, iSurf, thresh); 663 % Set the position of the cuts in the 3D figure 664 cutsPosVox = [74 93 159]; 665 panel_surface('PlotMri', hFigSrc3, cutsPosVox); 666 667 % View sources with MRI Viewer 668 hFigSrc4 = script_view_sources(sFilesSrcDeviant01.FileName, 'mriviewer'); 669 panel_surface('SetDataThreshold', hFigSrc4, iSurf, thresh); 670 % Set the position of the cuts in the MRI Viewer (values in millimeters) 671 figure_mri('SetLocation', 'voxel', hFigSrc4, [], cutsPosVox); 672 % Close figures 673 close([hFigSrc1 hFigSrc2 hFigSrc3 hFigSrc4]); 674 675 % === UNCONSTRAINED EXAMPLE === 676 % Unconstrained minnorm 677 InverseOptions.Comment = 'MN: MEG'; 678 InverseOptions.InverseMeasure = 'amplitude'; 679 InverseOptions.SourceOrient = {'free'}; 680 % Process: Compute sources [2018] 681 sFilesSrcUnconst = bst_process('CallProcess', 'process_inverse_2018', sFilesAvgDeviant01, [], ... 682 'output', 2, ... % Kernel only: one per file 683 'inverse', InverseOptions); 684 685 686 % === NORMALIZED SOURCES === 687 % dSPM 688 InverseOptions.Comment = 'dSPM: MEG'; 689 InverseOptions.InverseMeasure = 'dspm2018'; 690 InverseOptions.SourceOrient = {'fixed'}; 691 sFilesSrcDspm = bst_process('CallProcess', 'process_inverse_2018', sFilesAvgDeviant01, [], ... 692 'output', 2, ... % Kernel only: one per file 693 'inverse', InverseOptions); 694 % sLORETA (old function) 695 sFilesSrcSloreta = bst_process('CallProcess', 'process_inverse', sFilesAvgDeviant01, [], ... 696 'comment', '', ... 697 'method', 3, ... % sLORETA 698 'wmne', struct(... 699 'SourceOrient', {{'fixed'}}, ... 700 'loose', 0.2, ... 701 'SNR', 3, ... 702 'pca', 1, ... 703 'diagnoise', 0, ... 704 'regnoise', 1, ... 705 'magreg', 0.1, ... 706 'gradreg', 0.1, ... 707 'depth', 1, ... 708 'weightexp', 0.5, ... 709 'weightlimit', 10), ... 710 'sensortypes', 'MEG, MEG MAG, MEG GRAD, EEG', ... 711 'output', 2); % Kernel only: one per file 712 % Process: Z-score normalization: [-100ms,-2ms] 713 sFilesSrcZscore = bst_process('CallProcess', 'process_baseline_norm', sFilesSrcDeviant01, [], ... 714 'baseline', [-0.100, -0.002], ... 715 'source_abs', 0, ... 716 'method', 'zscore'); % Z-score transformation: x_std = (x - μ) / σ 717 718 % Process: Snapshot: Sources (one time) 719 bst_process('CallProcess', 'process_snapshot', sFilesSrcDeviant01, [], ... 720 'target', 8, ... % Sources (one time) 721 'orient', 1, ... % left 722 'time', 0.09, ... 723 'threshold', 30, ... 724 'Comment', 'Current density map (Constrained)'); 725 bst_process('CallProcess', 'process_snapshot', sFilesSrcDspm, [], ... 726 'target', 8, ... % Sources (one time) 727 'orient', 1, ... % left 728 'time', 0.09, ... 729 'threshold', 60, ... 730 'Comment', 'dSPM'); 731 bst_process('CallProcess', 'process_snapshot', sFilesSrcSloreta, [], ... 732 'target', 8, ... % Sources (one time) 733 'orient', 1, ... % left 734 'time', 0.09, ... 735 'threshold', 60, ... 736 'Comment', 'sLORETA'); 737 bst_process('CallProcess', 'process_snapshot', sFilesSrcZscore, [], ... 738 'target', 8, ... % Sources (one time) 739 'orient', 1, ... % left 740 'time', 0.09, ... 741 'threshold', 60, ... 742 'Comment', 'Z-score'); 743 bst_process('CallProcess', 'process_snapshot', sFilesSrcUnconst, [], ... 744 'target', 8, ... % Sources (one time) 745 'orient', 1, ... % left 746 'time', 0.0917, ... 747 'threshold', 0, ... 748 'Comment', 'Current density map (Unconstrained)'); 749 750 % === DELETE EXPERIMENTS === 751 % Process: Delete constrained example 752 bst_process('CallProcess', 'process_delete', [sFilesSrcDeviant01, sFilesSrcDspm, sFilesSrcSloreta, sFilesSrcZscore, sFilesSrcUnconst], [], ... 753 'target', 1); % Delete selected files 754 755 756 % === SHARED KERNEL === 757 % Constrained minnorm 758 InverseOptions.Comment = 'MN: MEG'; 759 InverseOptions.InverseMeasure = 'amplitude'; 760 InverseOptions.SourceOrient = {'fixed'}; 761 % Process: Compute sources [2018] 762 sFilesAvgSrc = bst_process('CallProcess', 'process_inverse_2018', sFilesAvg, [], ... 763 'output', 1, ... % Kernel only: shared 764 'inverse', InverseOptions); 765 766 767 % === AVERAGE SOURCES ACROSS RUNS === 768 % Process: Average: By trial group (subject average) 769 sFilesIntraSrc = bst_process('CallProcess', 'process_average', sFilesAvgSrc, [], ... 770 'avgtype', 6, ... % By trial group (subject average) 771 'avg_func', 1, ... % Arithmetic average: mean(x) 772 'weighted', 1, ... 773 'scalenormalized', 0); 774 % Process: Low-pass:40Hz 775 sFilesIntraSrcLow = bst_process('CallProcess', 'process_bandpass', sFilesIntraSrc, [], ... 776 'highpass', 0, ... 777 'lowpass', 40, ... 778 'attenuation', 'strict', ... % 60dB 779 'mirror', 0, ... 780 'overwrite', 0); 781 % Process: Z-score normalization: [-100ms,-2ms] 782 sFilesIntraZscore = bst_process('CallProcess', 'process_baseline_norm', sFilesIntraSrcLow, [], ... 783 'baseline', [-0.100, -0.0017], ... 784 'source_abs', 0, ... 785 'method', 'zscore'); % Z-score transformation: x_std = (x - μ) / σ 786 787 % Process: Delete intermediate results 788 bst_process('CallProcess', 'process_delete', sFilesIntraSrcLow, [], ... 789 'target', 1); % Delete selected files 790 % Screen captures 791 bst_process('CallProcess', 'process_snapshot', sFilesIntraZscore, [], ... 792 'target', 8, ... % Sources (one time) 793 'orient', 1, ... % left 794 'time', 0.09, ... 795 'threshold', 40, ... 796 'Comment', 'Average across runs (left)'); 797 bst_process('CallProcess', 'process_snapshot', sFilesIntraZscore, [], ... 798 'target', 8, ... % Sources (one time) 799 'orient', 2, ... % right 800 'time', 0.09, ... 801 'threshold', 40, ... 802 'Comment', 'Average across runs (right)'); 803 bst_process('CallProcess', 'process_snapshot', sFilesIntraZscore(1), [], ... 804 'target', 9, ... % Sources (contact sheet) 805 'orient', 1, ... % left 806 'contact_time', [0, 0.35], ... 807 'contact_nimage', 15, ... 808 'threshold', 20, ... 809 'Comment', 'Average across runs'); 810 811 812 %% ===== TUTORIAL #23: SCOUTS ======================================================== 813 % =================================================================================== 814 disp([10 'DEMO> Tutorial #23: Scouts' 10]); 815 % Load surface file 816 sCortex = in_tess_bst(CortexFile); 817 % Add new scouts in first atlas 818 sCortex.iAtlas = find(strcmpi({sCortex.Atlas.Name}, 'Destrieux')); 819 % Save file 820 bst_save(file_fullpath(CortexFile), sCortex, 'v7'); 821 % Unload everything 822 bst_memory('UnloadAll', 'Forced'); 823 % Find scouts indices to display: {'A1L', 'A1R', 'IFGL', 'M1L'} 824 [tmp,iScouts,tmp] = intersect({sCortex.Atlas(sCortex.iAtlas).Scouts.Label}, {'G_temp_sup-G_T_transv L', 'G_temp_sup-G_T_transv R', 'G_front_inf-Opercular L', 'G_precentral L'}); 825 826 % View cortex 827 hFigSurf1 = view_surface(CortexFile, [], [], 'NewFigure'); 828 hFigSurf2 = view_surface(CortexFile, [], [], 'NewFigure'); 829 figure_3d('SetStandardView', hFigSurf1, 'left'); 830 figure_3d('SetStandardView', hFigSurf2, 'right'); 831 panel_surface('SetSurfaceSmooth', hFigSurf1, 1, .6, 0); 832 panel_surface('SetSurfaceSmooth', hFigSurf2, 1, .6, 0); 833 panel_surface('SetShowSulci', hFigSurf1, 1, 1); 834 panel_surface('SetShowSulci', hFigSurf2, 1, 1); 835 % Configure scouts display 836 panel_scout('SetScoutsOptions', 0, 1, 1, 'all', 0.7, 1, 1, 0); 837 % View scouts 838 hFigScouts = view_scouts({sFilesIntraZscore.FileName}, iScouts); 839 hLegend = findobj(hFigScouts, 'Type', 'legend'); 840 if ~isempty(hLegend) && ishandle(hLegend(1)) 841 set(hLegend(1), 'Units', 'pixels'); 842 pos = get(hLegend(1), 'Position'); 843 set(hLegend(1), 'Position', [1, 1, pos(3), pos(4)]); 844 end 845 % Save figures 846 bst_report('Snapshot', hFigScouts, sFilesIntraZscore(1).FileName, 'Scouts', [100 100 670 250]); 847 % Close all 848 pause(1); 849 bst_memory('UnloadAll', 'Forced'); 850 851 852 853 %% ===== TUTORIAL #24: TIME-FREQUENCY ================================================ 854 % =================================================================================== 855 disp([10 'DEMO> Tutorial #24: Time-frequency' 10]); 856 % Process: Simulate generic signals 857 sSim = bst_process('CallProcess', 'process_simulate_matrix', [], [], ... 858 'subjectname', 'Test', ... 859 'condition', 'Simulation', ... 860 'samples', 6000, ... 861 'srate', 1000, ... 862 'matlab', ['f1 = 2; f2 = 20; f3 = 50;' 10 'i =2000:6000;' 10 'Data(1,i) = sin(f1*2*pi*t(i)) + 0.4 * cos(f2*2*pi*t(i));' 10 'Data = Data + 0.2 * sin(f3*2*pi*t) + 0.4 * rand(1,6000);']); 863 864 % Time-frequency options 865 TfOptions = struct(... 866 'Comment', 'Power,1-60Hz', ... 867 'TimeBands', [], ... 868 'Freqs', [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60], ... 869 'MorletFc', 1, ... 870 'MorletFwhmTc', 3, ... 871 'ClusterFuncTime', 'none', ... 872 'Measure', 'power', ... 873 'Output', 'average', ... 874 'SaveKernel', 0); 875 % Process: Time-frequency (Morlet wavelets) 876 sSimTf1 = bst_process('CallProcess', 'process_timefreq', sSim, [], ... 877 'edit', TfOptions, ... 878 'normalize', 'multiply'); % 1/f compensation: Multiply output values by frequency 879 880 % === NORMALIZATION === 881 % Process: Time-frequency (Morlet wavelets) 882 sSimTf2 = bst_process('CallProcess', 'process_timefreq', sSim, [], ... 883 'edit', TfOptions, ... 884 'normalize', 'none'); % None: Save non-standardized time-frequency maps 885 % Process: Spectral flattening 886 sSimTf2Multi = bst_process('CallProcess', 'process_tf_norm', sSimTf2, [], ... 887 'normalize', 'multiply', ... % 1/f compensation (multiple by frequency) 888 'overwrite', 0); 889 % Process: Event-related perturbation (ERS/ERD): [750ms,1250ms] 890 sSimTf2Ersd = bst_process('CallProcess', 'process_baseline_norm', sSimTf2, [], ... 891 'baseline', [0.75, 1.25], ... 892 'method', 'ersd', ... % Event-related perturbation (ERS/ERD): x_std = (x - μ) / μ * 100 893 'overwrite', 0); 894 895 % Process: Snapshot: Time-frequency maps 896 bst_process('CallProcess', 'process_snapshot', sSimTf2, [], ... 897 'target', 14, ... % Time-frequency maps 898 'Comment', 'Not normalized'); 899 % Process: Snapshot: Time-frequency maps 900 bst_process('CallProcess', 'process_snapshot', sSimTf2Multi, [], ... 901 'target', 14, ... % Time-frequency maps 902 'Comment', 'Spectral flattening: 1/f compensation'); 903 % Process: Snapshot: Time-frequency maps 904 bst_process('CallProcess', 'process_snapshot', sSimTf2Ersd, [], ... 905 'target', 14, ... % Time-frequency maps 906 'Comment', 'ERS/ERD'); 907 % Spectrum/Time series 908 hFigTf1 = view_spectrum(sSimTf2.FileName, 'Spectrum'); 909 hFigTf2 = view_spectrum(sSimTf2.FileName, 'TimeSeries'); 910 panel_time('SetCurrentTime', 0.5); 911 panel_freq('SetCurrentFreq', 2, 0); 912 bst_report('Snapshot', [hFigTf1 hFigTf2], sSimTf2.FileName, 'Not normalized: 2s/20Hz', [200, 200, 400, 250]); 913 panel_time('SetCurrentTime', 2.02); 914 panel_freq('SetCurrentFreq', 20, 0); 915 bst_report('Snapshot', [hFigTf1 hFigTf2], sSimTf2.FileName, 'Not normalized: 2s/20Hz', [200, 200, 400, 250]); 916 bst_memory('UnloadAll', 'Forced'); 917 918 919 % === HILBERT TRANSFORM === 920 % Process: Hilbert transform 921 sSimHilbert = bst_process('CallProcess', 'process_hilbert', sSim, [], ... 922 'edit', struct(... 923 'Comment', 'Power', ... 924 'TimeBands', [], ... 925 'Freqs', {{'delta', '2, 4', 'mean'; 'theta', '5, 7', 'mean'; 'alpha', '8, 12', 'mean'; 'beta', '15, 29', 'mean'; 'gamma1', '30, 59', 'mean'; 'gamma2', '60, 90', 'mean'}}, ... 926 'ClusterFuncTime', 'none', ... 927 'Measure', 'power', ... 928 'Output', 'all', ... 929 'SaveKernel', 0), ... 930 'normalize', 'none', ... % None: Save non-standardized time-frequency maps 931 'mirror', 0); 932 % Process: Spectral flattening 933 sSimHilbertMulti = bst_process('CallProcess', 'process_tf_norm', sSimHilbert, [], ... 934 'normalize', 'multiply', ... % 1/f compensation (multiple by frequency) 935 'overwrite', 0); 936 % Process: Event-related perturbation (ERS/ERD): [750ms,1250ms] 937 sSimHilbertErsd = bst_process('CallProcess', 'process_baseline_norm', sSimHilbert, [], ... 938 'baseline', [0.75, 1.25], ... 939 'method', 'ersd', ... % Event-related perturbation (ERS/ERD): x_std = (x - μ) / μ * 100 940 'overwrite', 0); 941 942 % Process: Snapshot: Time-frequency maps 943 bst_process('CallProcess', 'process_snapshot', sSimHilbert, [], ... 944 'target', 14, ... % Time-frequency maps 945 'Comment', 'Not normalized'); 946 % Process: Snapshot: Time-frequency maps 947 bst_process('CallProcess', 'process_snapshot', sSimHilbertMulti, [], ... 948 'target', 14, ... % Time-frequency maps 949 'Comment', 'Spectral flattening: 1/f compensation'); 950 % Process: Snapshot: Time-frequency maps 951 bst_process('CallProcess', 'process_snapshot', sSimHilbertErsd, [], ... 952 'target', 14, ... % Time-frequency maps 953 'Comment', 'ERS/ERD'); 954 955 % === SINGLE TRIALS === 956 TfOptions.Comment = 'Avg,Power,1-150Hz'; 957 TfOptions.Freqs = [1, 2, 3.1, 4.2, 5.4, 6.7, 8, 9.5, 11, 12.6, 14.3, 16.1, 18.1, 20.1, 22.3, 24.6, 27, 29.6, 32.4, 35.3, 38.4, 41.6, 45.1, 48.8, 52.7, 56.9, 61.3, 66, 70.9, 76.2, 81.8, 87.7, 94, 100.6, 107.7, 115.2, 123.1, 131.6, 140.5, 150]; 958 % Process: Time-frequency (Morlet wavelets) 959 sEpochsAvgTf = bst_process('CallProcess', 'process_timefreq', sFilesEpochs1, [], ... 960 'sensortypes', 'MEG, EEG', ... 961 'edit', TfOptions, ... 962 'normalize', 'none'); % None: Save non-standardized time-frequency maps 963 % Process: Event-related perturbation (ERS/ERD): [-75ms,0ms] 964 sEpochsAvgTfErsd = bst_process('CallProcess', 'process_baseline_norm', sEpochsAvgTf, [], ... 965 'baseline', [-0.075, 0], ... 966 'method', 'ersd', ... % Event-related perturbation (ERS/ERD): x_std = (x - μ) / μ * 100 967 'overwrite', 0); 968 969 % === DISPLAY === 970 % View time-frequency file 971 hFigTf1 = view_timefreq(sEpochsAvgTfErsd.FileName, 'SingleSensor'); 972 % Configure display 973 sOptions = panel_display('GetDisplayOptions'); 974 sOptions.HideEdgeEffects = 1; 975 sOptions.HighResolution = 1; 976 panel_display('SetDisplayOptions', sOptions); 977 % Other display modes 978 hFigTf2 = view_timefreq(sEpochsAvgTfErsd.FileName, 'AllSensors'); 979 hFigTf3 = view_timefreq(sEpochsAvgTfErsd.FileName, '2DLayout'); 980 hFigTf4 = view_timefreq(sEpochsAvgTfErsd.FileName, '2DLayoutOpt'); 981 bst_colormaps('SetColormapName', 'stat2', 'jet'); 982 bst_colormaps('SetColormapAbsolute', 'stat2', 1); 983 bst_report('Snapshot', hFigTf1, sEpochsAvgTfErsd.FileName, 'Time-frequency', [200, 200, 400, 250]); 984 bst_report('Snapshot', [hFigTf2 hFigTf3 hFigTf4], sEpochsAvgTfErsd.FileName, 'Time-frequency', [200, 200, 750, 400]); 985 close([hFigTf1 hFigTf2 hFigTf3 hFigTf4]); 986 % Image [channel x time] 987 hFigTf5 = view_erpimage(sEpochsAvgTfErsd.FileName, 'trialimage'); 988 % Topographies 989 hFigTf6 = view_topography(sEpochsAvgTfErsd.FileName, 'MEG', '3DSensorCap'); 990 hFigTf7 = view_topography(sEpochsAvgTfErsd.FileName, 'MEG', '2DDisc'); 991 hFigTf8 = view_topography(sEpochsAvgTfErsd.FileName, 'MEG', '2DSensorCap'); 992 hFigTf9 = view_topography(sEpochsAvgTfErsd.FileName, 'MEG', '2DLayout'); 993 panel_time('SetCurrentTime', 0.175); 994 panel_freq('SetCurrentFreq', 8, 0); 995 bst_report('Snapshot', [hFigTf5 hFigTf6 hFigTf7 hFigTf8 hFigTf9], sEpochsAvgTfErsd.FileName, 'Time-frequency: 8Hz', [200, 200, 400, 250]); 996 close([hFigTf5 hFigTf6 hFigTf7 hFigTf8 hFigTf9]); 997 998 999 % === AVERAGE FOR SCOUTS === 1000 % Select all sources for the single deviant epochs 1001 sFilesEpochDeviantSrc = bst_process('CallProcess', 'process_select_files_results', [], [], ... 1002 'subjectname', SubjectName, ... 1003 'condition', '', ... 1004 'includebad', 0); 1005 sFilesEpochDeviantSrc = bst_process('CallProcess', 'process_select_tag', sFilesEpochDeviantSrc, [], ... 1006 'tag', 'deviant', ... 1007 'search', 1, ... % Search the file names 1008 'select', 1); % Select only the files with the tag 1009 sFilesEpochDeviantSrc = bst_process('CallProcess', 'process_select_tag', sFilesEpochDeviantSrc, [], ... 1010 'tag', 'average', ... 1011 'search', 1, ... % Search the file names 1012 'select', 2); % Exclude the files with the tag 1013 1014 % Process: Time-frequency (Morlet wavelets) 1015 sFilesTfScout = bst_process('CallProcess', 'process_timefreq', sFilesEpochDeviantSrc, [], ... 1016 'clusters', {'Destrieux', {'G_temp_sup-G_T_transv L', 'G_temp_sup-G_T_transv R', 'G_front_inf-Opercular L', 'G_precentral L'}}, ... 1017 'scoutfunc', 1, ... % Mean 1018 'edit', struct(... 1019 'Comment', 'Deviant: Scouts,Avg,Power,1-150Hz', ... 1020 'TimeBands', [], ... 1021 'Freqs', [1, 2, 3.1, 4.2, 5.4, 6.7, 8, 9.5, 11, 12.6, 14.3, 16.1, 18.1, 20.1, 22.3, 24.6, 27, 29.6, 32.4, 35.3, 38.4, 41.6, 45.1, 48.8, 52.7, 56.9, 61.3, 66, 70.9, 76.2, 81.8, 87.7, 94, 100.6, 107.7, 115.2, 123.1, 131.6, 140.5, 150], ... 1022 'MorletFc', 1, ... 1023 'MorletFwhmTc', 3, ... 1024 'ClusterFuncTime', 'after', ... 1025 'Measure', 'power', ... 1026 'Output', 'average', ... 1027 'SaveKernel', 0), ... 1028 'normalize', 'none'); % None: Save non-standardized time-frequency maps) 1029 % Process: Event-related perturbation (ERS/ERD): [-75ms,0ms] 1030 sFilesTfScoutErsd = bst_process('CallProcess', 'process_baseline_norm', sFilesTfScout, [], ... 1031 'baseline', [-0.075, 0], ... 1032 'method', 'ersd', ... % Event-related perturbation (ERS/ERD): x_std = (x - μ) / μ * 100 1033 'overwrite', 0); 1034 % Process: Snapshot: Time-frequency maps 1035 bst_process('CallProcess', 'process_snapshot', sFilesTfScoutErsd, [], ... 1036 'target', 14, ... % Time-frequency maps 1037 'Comment', 'ERS/ERD'); 1038 1039 1040 % === FULL CORTEX / HILBERT === 1041 % Process: Hilbert transform 1042 sFilesHilbertCortex = bst_process('CallProcess', 'process_hilbert', sFilesEpochDeviantSrc, [], ... 1043 'clusters', {}, ... 1044 'scoutfunc', 1, ... % Mean 1045 'edit', struct(... 1046 'Comment', 'Deviant: Avg,Magnitude', ... 1047 'TimeBands', [], ... 1048 'Freqs', {{'delta', '2, 4', 'mean'; 'theta', '5, 7', 'mean'; 'alpha', '8, 12', 'mean'; 'beta', '15, 29', 'mean'; 'gamma1', '30, 59', 'mean'; 'gamma2', '60, 90', 'mean'}}, ... 1049 'ClusterFuncTime', 'none', ... 1050 'Measure', 'power', ... 1051 'Output', 'average', ... 1052 'RemoveEvoked', 0, ... 1053 'SaveKernel', 0), ... 1054 'normalize', 'none', ... % None: Save non-standardized time-frequency maps 1055 'mirror', 0); 1056 % Process: Event-related perturbation (ERS/ERD): [-75ms,0ms] 1057 sFilesHilbertCortexErsd = bst_process('CallProcess', 'process_baseline_norm', sFilesHilbertCortex, [], ... 1058 'baseline', [-0.075, 0], ... 1059 'method', 'ersd', ... % Event-related perturbation (ERS/ERD): x_std = (x - μ) / μ * 100 1060 'overwrite', 0); 1061 1062 % View results 1063 hFigTf1 = view_surface_data([], sFilesHilbertCortexErsd.FileName); 1064 hFigTf2 = view_timefreq(sFilesHilbertCortexErsd.FileName, 'SingleSensor', 362); 1065 figure_3d('SetStandardView', hFigTf1, 'left'); 1066 panel_surface('SetDataThreshold', hFigTf1, 1, 0.5); 1067 panel_time('SetCurrentTime', 0.175); 1068 panel_freq('SetCurrentFreq', 3); 1069 bst_colormaps('RestoreDefaults', 'stat2'); 1070 bst_report('Snapshot', [hFigTf1 hFigTf2], sFilesHilbertCortexErsd.FileName, 'Hilbert transform: Alpha band', [200, 200, 400, 250]); 1071 bst_memory('UnloadAll', 'Forced'); 1072 1073 1074 %% ===== TUTORIAL #25: DIFFRERENCE =================================================== 1075 % =================================================================================== 1076 disp([10 'DEMO> Tutorial #25: Difference' 10]); 1077 1078 % ===== SELECT TRIALS (DATA) ===== 1079 % Process: Select recordings in: Subject01/* 1080 sFilesAll = bst_process('CallProcess', 'process_select_files_data', [], [], ... 1081 'subjectname', SubjectName, ... 1082 'condition', '', ... 1083 'includebad', 0); 1084 % Process: Select file names with tag: deviant_trial 1085 sEpochDeviant = bst_process('CallProcess', 'process_select_tag', sFilesAll, [], ... 1086 'tag', 'deviant_trial', ... 1087 'search', 1, ... % Search the file names 1088 'select', 1); % Select only the files with the tag 1089 % Process: Select file names with tag: standard_trial 1090 sEpochStandard = bst_process('CallProcess', 'process_select_tag', sFilesAll, [], ... 1091 'tag', 'standard_trial', ... 1092 'search', 1, ... % Search the file names 1093 'select', 1); % Select only the files with the tag 1094 1095 % ===== SELECT TRIALS (SOURCE) ===== 1096 % Process: Select recordings in: Subject01/* 1097 sFilesAllSrc = bst_process('CallProcess', 'process_select_files_results', [], [], ... 1098 'subjectname', SubjectName, ... 1099 'condition', '', ... 1100 'includebad', 0); 1101 % Process: Select file names with tag: deviant_trial 1102 sEpochDeviantSrc = bst_process('CallProcess', 'process_select_tag', sFilesAllSrc, [], ... 1103 'tag', 'deviant_trial', ... 1104 'search', 1, ... % Search the file names 1105 'select', 1); % Select only the files with the tag 1106 % Process: Select file names with tag: standard_trial 1107 sEpochStandardSrc = bst_process('CallProcess', 'process_select_tag', sFilesAllSrc, [], ... 1108 'tag', 'standard_trial', ... 1109 'search', 1, ... % Search the file names 1110 'select', 1); % Select only the files with the tag 1111 1112 % ===== ABSOLUTE DIFFERENCE ====== 1113 % Process: Difference: A-B, abs 1114 sDiffSrc = bst_process('CallProcess', 'process_diff_ab', sFilesIntraSrc(1).FileName, sFilesIntraSrc(2).FileName, ... 1115 'source_abs', 1); 1116 % Process: Set comment: deviant|abs - standard|abs 1117 sDiffSrc = bst_process('CallProcess', 'process_set_comment', sDiffSrc, [], ... 1118 'tag', 'deviant|abs - standard|abs', ... 1119 'isindex', 1); 1120 % Process: Low-pass:40Hz 1121 sDiffSrcZscore = bst_process('CallProcess', 'process_bandpass', sDiffSrc, [], ... 1122 'highpass', 0, ... 1123 'lowpass', 40, ... 1124 'attenuation', 'strict', ... % 60dB 1125 'mirror', 0, ... 1126 'overwrite', 1); 1127 % Process: Z-score transformation: [-100ms,-2ms] 1128 sDiffSrcZscore = bst_process('CallProcess', 'process_baseline_norm', sDiffSrcZscore, [], ... 1129 'baseline', [-0.1, -0.002], ... 1130 'source_abs', 0, ... 1131 'method', 'zscore', ... % Z-score transformation: x_std = (x - μ) / σ 1132 'overwrite', 0); 1133 % Process: Snapshot: Sources (contact sheet) 1134 bst_process('CallProcess', 'process_snapshot', sDiffSrcZscore, [], ... 1135 'target', 9, ... % Sources (contact sheet) 1136 'modality', 1, ... % MEG (All) 1137 'orient', 1, ... % left 1138 'contact_time', [0, 0.35], ... 1139 'contact_nimage', 15, ... 1140 'threshold', 30, ... 1141 'Comment', 'Difference deviant - standard (absolute)'); 1142 1143 % ===== RELATIVE DIFFERENCE ===== 1144 % Process: Difference: A-B 1145 sDiffSrcRel = bst_process('CallProcess', 'process_diff_ab', sFilesIntraSrc(1).FileName, sFilesIntraSrc(2).FileName, ... 1146 'source_abs', 0); 1147 % Process: Set comment: deviant - standard 1148 sDiffSrcRel = bst_process('CallProcess', 'process_set_comment', sDiffSrcRel, [], ... 1149 'tag', 'deviant - standard', ... 1150 'isindex', 1); 1151 % Process: Low-pass:40Hz 1152 sDiffSrcRelZscore = bst_process('CallProcess', 'process_bandpass', sDiffSrcRel, [], ... 1153 'highpass', 0, ... 1154 'lowpass', 40, ... 1155 'attenuation', 'strict', ... % 60dB 1156 'mirror', 0, ... 1157 'overwrite', 1); 1158 % Process: Z-score transformation: [-100ms,-2ms] 1159 sDiffSrcRelZscore = bst_process('CallProcess', 'process_baseline_norm', sDiffSrcRelZscore, [], ... 1160 'baseline', [-0.1, -0.002], ... 1161 'source_abs', 0, ... 1162 'method', 'zscore', ... % Z-score transformation: x_std = (x - μ) / σ 1163 'overwrite', 0); 1164 % Configure colormap: hot/absolute 1165 bst_colormaps('SetColormapName', 'stat2', 'hot'); 1166 bst_colormaps('SetColormapAbsolute', 'stat2', 1); 1167 % Process: Snapshot: Sources (contact sheet) 1168 bst_process('CallProcess', 'process_snapshot', sDiffSrcRelZscore, [], ... 1169 'target', 9, ... % Sources (contact sheet) 1170 'modality', 1, ... % MEG (All) 1171 'orient', 1, ... % left 1172 'contact_time', [0, 0.35], ... 1173 'contact_nimage', 15, ... 1174 'threshold', 30, ... 1175 'Comment', 'Difference deviant - standard (relative)'); 1176 % Restore colormap: rwb/relative 1177 bst_colormaps('SetColormapName', 'stat2', 'cmap_rbw'); 1178 bst_colormaps('SetColormapAbsolute', 'stat2', 0); 1179 1180 % ===== DIFFERENCE OF MEANS ===== 1181 % Process: Select uniform number of files [uniform] 1182 [sEpochDeviantUni, sEpochStandardUni] = bst_process('CallProcess', 'process_select_uniform2', sEpochDeviant, sEpochStandard, ... 1183 'nfiles', 0, ... 1184 'method', 4); % Uniformly distributed 1185 % Process: Difference of means [mean] 1186 sDiffMean = bst_process('CallProcess', 'process_diff_mean', sEpochDeviantUni, sEpochStandardUni, ... 1187 'avg_func', 1, ... % Arithmetic average mean(A) - mean(B) 1188 'weighted', 0); 1189 % Process: Snapshot: Recordings time series 1190 bst_process('CallProcess', 'process_snapshot', sDiffMean, [], ... 1191 'target', 5, ... % Recordings time series 1192 'modality', 1, ... % MEG (All) 1193 'Comment', 'Difference of means'); 1194 1195 1196 1197 %% ===== TUTORIAL #26: STATISTICS ==================================================== 1198 % =================================================================================== 1199 disp([10 'DEMO> Tutorial #26: Statistics' 10]); 1200 % ===== HISTOGRAMS ===== 1201 % Process: Extract values: [160ms] MLP57 1202 sHistoDeviant = bst_process('CallProcess', 'process_extract_values', sEpochDeviant, [], ... 1203 'timewindow', [0.16, 0.16], ... 1204 'sensortypes', 'MLP57', ... 1205 'isabs', 0, ... 1206 'avgtime', 0, ... 1207 'avgrow', 0, ... 1208 'dim', 2, ... % Concatenate time (dimension 2) 1209 'Comment', ''); 1210 % Process: Extract values: [160ms] MLP57 1211 sHistoStandard = bst_process('CallProcess', 'process_extract_values', sEpochStandard, [], ... 1212 'timewindow', [0.16, 0.16], ... 1213 'sensortypes', 'MLP57', ... 1214 'isabs', 0, ... 1215 'avgtime', 0, ... 1216 'avgrow', 0, ... 1217 'dim', 2, ... % Concatenate time (dimension 2) 1218 'Comment', ''); 1219 % Display histograms 1220 hFigHisto = view_histogram({sHistoDeviant.FileName, sHistoStandard.FileName}); 1221 bst_report('Snapshot', hFigHisto, sHistoDeviant.FileName, 'Histograms for MLP57/160ms'); 1222 close(hFigHisto); 1223 1224 % ===== EXEMPLE #1: PARAMETRIC/DATA ===== 1225 % Process: t-test [equal] [-100ms,500ms] H0:(A-B = 0) 1226 sTestParamData = bst_process('CallProcess', 'process_test_parametric2', sEpochDeviant, sEpochStandard, ... 1227 'timewindow', [-0.1, 0.5], ... 1228 'sensortypes', '', ... 1229 'isabs', 0, ... 1230 'avgtime', 0, ... 1231 'avgrow', 0, ... 1232 'Comment', '', ... 1233 'test_type', 'ttest_equal', ... % Student's t-test (equal variance) t = (mean(A)-mean(B)) / (Sx * sqrt(1/nA + 1/nB))Sx = sqrt(((nA-1)*var(A) + (nB-1)*var(B)) / (nA+nB-2)) df = nA + nB - 2 1234 'tail', 'two'); % Two-tailed 1235 % Set display properties 1236 StatThreshOptions = bst_get('StatThreshOptions'); 1237 StatThreshOptions.pThreshold = 0.05; 1238 StatThreshOptions.Correction = 'fdr'; 1239 StatThreshOptions.Control = [1 2 3]; 1240 bst_set('StatThreshOptions', StatThreshOptions); 1241 % Process: Snapshot: Recordings time series 1242 bst_process('CallProcess', 'process_snapshot', sTestParamData, [], ... 1243 'target', 5, ... % Recordings time series 1244 'modality', 1, ... % MEG (All) 1245 'time', 0.16, ... 1246 'Comment', 'Parametric t-test (p<0.05, FDR)'); 1247 % Process: Snapshot: Recordings topography (one time) 1248 bst_process('CallProcess', 'process_snapshot', sTestParamData, [], ... 1249 'target', 6, ... % Recordings topography (one time) 1250 'modality', 1, ... % MEG (All) 1251 'time', 0.16, ... 1252 'Comment', 'Parametric t-test (p<0.05, FDR)'); 1253 1254 % ===== EXEMPLE #2: NON-PARAMETRIC/DATA ===== 1255 % Process: Perm t-test equal [-100ms,500ms MEG] H0:(A=B), H1:(A<>B) 1256 sTestPermData = bst_process('CallProcess', 'process_test_permutation2', sEpochDeviant, sEpochStandard, ... 1257 'timewindow', [-0.1, 0.5], ... 1258 'sensortypes', 'MEG', ... 1259 'isabs', 0, ... 1260 'avgtime', 0, ... 1261 'avgrow', 0, ... 1262 'iszerobad', 1, ... 1263 'Comment', '', ... 1264 'test_type', 'ttest_equal', ... % Student's t-test (equal variance) t = (mean(A)-mean(B)) / (Sx * sqrt(1/nA + 1/nB))Sx = sqrt(((nA-1)*var(A) + (nB-1)*var(B)) / (nA+nB-2)) 1265 'randomizations', 1000, ... 1266 'tail', 'two'); % Two-tailed 1267 % Process: Snapshot: Recordings time series 1268 bst_process('CallProcess', 'process_snapshot', sTestPermData, [], ... 1269 'target', 5, ... % Recordings time series 1270 'modality', 1, ... % MEG (All) 1271 'time', 0.16, ... 1272 'Comment', 'Non-parametric t-test (p<0.05, FDR)'); 1273 % Process: Snapshot: Recordings topography (one time) 1274 bst_process('CallProcess', 'process_snapshot', sTestPermData, [], ... 1275 'target', 6, ... % Recordings topography (one time) 1276 'modality', 1, ... % MEG (All) 1277 'time', 0.16, ... 1278 'Comment', 'Non-parametric t-test (p<0.05, FDR)'); 1279 1280 % ===== EXEMPLE #3: CLUSTER/DATA ===== 1281 % Process: FT t-test unequal cluster [-100ms,500ms MEG] H0:(A=B), H1:(A<>B) 1282 sTestClustData = bst_process('CallProcess', 'process_ft_timelockstatistics', sEpochDeviant, sEpochStandard, ... 1283 'sensortypes', 'MEG', ... 1284 'timewindow', [-0.1, 0.5], ... 1285 'isabs', 0, ... 1286 'avgtime', 0, ... 1287 'avgchan', 0, ... 1288 'randomizations', 1000, ... 1289 'statistictype', 1, ... % Independent t-test 1290 'tail', 'two', ... % Two-tailed 1291 'correctiontype', 2, ... % cluster 1292 'minnbchan', 0, ... 1293 'clusteralpha', 0.05); 1294 % Process: Snapshot: Recordings time series 1295 bst_process('CallProcess', 'process_snapshot', sTestClustData, [], ... 1296 'target', 5, ... % Recordings time series 1297 'time', 0.16, ... 1298 'modality', 1, ... % MEG (All) 1299 'Comment', 'Cluster-based permutation test'); 1300 % Process: Snapshot: Recordings topography (one time) 1301 bst_process('CallProcess', 'process_snapshot', sTestClustData, [], ... 1302 'target', 6, ... % Recordings topography (one time) 1303 'modality', 1, ... % MEG (All) 1304 'time', 0.16, ... 1305 'Comment', 'Cluster-based permutation test'); 1306 1307 % ===== EXAMPLE #4: PARAMETRIC/SOURCES ===== 1308 % Process: t-test [equal] [-100ms,500ms] H0:(A-B = 0) 1309 sTestParamSrc = bst_process('CallProcess', 'process_test_parametric2', sEpochDeviantSrc, sEpochStandardSrc, ... 1310 'timewindow', [-0.1, 0.5], ... 1311 'scoutsel', {}, ... 1312 'scoutfunc', 1, ... % Mean 1313 'isnorm', 0, ... 1314 'avgtime', 0, ... 1315 'Comment', '', ... 1316 'test_type', 'ttest_equal', ... % Student's t-test (equal variance) t = (mean(A)-mean(B)) / (Sx * sqrt(1/nA + 1/nB))Sx = sqrt(((nA-1)*var(A) + (nB-1)*var(B)) / (nA+nB-2)) df = nA + nB - 2 1317 'tail', 'two'); % Two-tailed 1318 % Process: Difference of means [abs(mean)] 1319 sDiffMeanSrc = bst_process('CallProcess', 'process_diff_mean', sEpochDeviantSrc, sEpochStandardSrc, ... 1320 'avg_func', 2, ... % Absolute value of average abs(mean(A)) - abs(mean(B)) 1321 'weighted', 0); 1322 % Process: Apply statistic threshold: p<0.05 (FDR:1,2,3) 1323 sDiffMeanSrcThresh = bst_process('CallProcess', 'process_extract_pthresh2', sTestParamSrc, sDiffMeanSrc, ... 1324 'pthresh', 0.05, ... 1325 'correction', 3, ... % False discovery rate (FDR) 1326 'control1', 1, ... 1327 'control2', 1, ... 1328 'control3', 1); 1329 % Process: Snapshot: Sources (one time) 1330 bst_process('CallProcess', 'process_snapshot', sTestParamSrc, [], ... 1331 'target', 8, ... % Sources (one time) 1332 'orient', 1, ... % left 1333 'time', 0.148, ... 1334 'threshold', 40, ... 1335 'Comment', 'Parametric t-test (p<0.05, FDR)'); 1336 % Process: Snapshot: Sources (one time) 1337 bst_process('CallProcess', 'process_snapshot', sDiffMeanSrc, [], ... 1338 'target', 8, ... % Sources (one time) 1339 'orient', 1, ... % left 1340 'time', 0.148, ... 1341 'threshold', 40, ... 1342 'Comment', 'abs(average(deviant)) - abs(average(standard))'); 1343 % Process: Snapshot: Sources (one time) 1344 bst_process('CallProcess', 'process_snapshot', sDiffMeanSrcThresh, [], ... 1345 'target', 8, ... % Sources (one time) 1346 'orient', 1, ... % left 1347 'time', 0.148, ... 1348 'threshold', 0, ... 1349 'Comment', 'Different of mean thresholded with t-test results'); 1350 1351 % ===== EXAMPLE #5: PARAMETRIC/SCOUTS ===== 1352 % Process: t-test equal [-100ms,500ms] H0:(A=B), H1:(A<>B) 1353 sTestParamScout = bst_process('CallProcess', 'process_test_parametric2', sEpochDeviantSrc, sEpochStandardSrc, ... 1354 'timewindow', [-0.1, 0.5], ... 1355 'scoutsel', {'Destrieux', {'G_front_inf-Opercular L', 'G_precentral L', 'G_temp_sup-G_T_transv L'}}, ... 1356 'scoutfunc', 1, ... % Mean 1357 'isnorm', 0, ... 1358 'avgtime', 0, ... 1359 'Comment', '', ... 1360 'test_type', 'ttest_equal', ... % Student's t-test (equal variance) t = (mean(A)-mean(B)) / (Sx * sqrt(1/nA + 1/nB))Sx = sqrt(((nA-1)*var(A) + (nB-1)*var(B)) / (nA+nB-2)) df = nA + nB - 2 1361 'tail', 'two'); % Two-tailed 1362 % Process: Apply statistic threshold: p<0.05 (FDR:1,2,3) 1363 sTestParamScoutThresh = bst_process('CallProcess', 'process_extract_pthresh', sTestParamScout, [], ... 1364 'pthresh', 0.05, ... 1365 'correction', 3, ... % False discovery rate (FDR) 1366 'control1', 1, ... 1367 'control2', 1, ... 1368 'control3', 1); 1369 % Process: Compute head model 1370 bst_process('CallProcess', 'process_headmodel', sTestParamScoutThresh, [], ... 1371 'sourcespace', 1, ... % Cortex surface 1372 'meg', 3); % Overlapping spheres 1373 % Process: Simulate recordings from scouts 1374 sSimulData = bst_process('CallProcess', 'process_simulate_recordings', sTestParamScoutThresh, [], ... 1375 'scouts', {'Destrieux', {'G_front_inf-Opercular L', 'G_precentral L', 'G_temp_sup-G_T_transv L'}}, ... 1376 'savesources', 1); 1377 1378 % Get corresponding source file 1379 [sStudy,iStudy,iRes] = bst_get('ResultsForDataFile', sSimulData.FileName); 1380 sSimulSrc = sStudy.Result(iRes).FileName; 1381 % Reset visualization filters 1382 panel_filter('SetFilters', 0, [], 0, [], 0, [], 0, 0); 1383 % Process: Snapshot: Recordings time series 1384 bst_process('CallProcess', 'process_snapshot', sTestParamScout, [], ... 1385 'target', 5, ... % Recordings time series 1386 'time', 0.148, ... 1387 'Comment', 'Parametric t-test (p<0.05, FDR)'); 1388 % Process: Snapshot: Sources (one time) 1389 bst_process('CallProcess', 'process_snapshot', sSimulSrc, [], ... 1390 'target', 8, ... % Sources (one time) 1391 'orient', 1, ... % left 1392 'time', 0.148, ... 1393 'threshold', 0, ... 1394 'Comment', 'Simulated sources'); 1395 % Process: Snapshot: Recordings time series 1396 bst_process('CallProcess', 'process_snapshot', sSimulData, [], ... 1397 'target', 5, ... % Recordings time series 1398 'modality', 1, ... % MEG (All) 1399 'time', 0.148, ... 1400 'Comment', 'Simulated MEG recordings'); 1401 1402 1403 % ===== EXAMPLE #6: NON-PARAMETRIC/TIMEFREQ ===== 1404 TfOptions.Output = 'all'; 1405 % Process: Time-frequency (Morlet wavelets) / DEVIANT 1406 sEpochDeviantTf = bst_process('CallProcess', 'process_timefreq', sEpochDeviant, [], ... 1407 'sensortypes', 'MLP57', ... 1408 'edit', TfOptions, ... 1409 'normalize', 'none'); % None: Save non-standardized time-frequency maps 1410 % Process: Time-frequency (Morlet wavelets) / STANDARD 1411 sEpochStandardTf = bst_process('CallProcess', 'process_timefreq', sEpochStandard, [], ... 1412 'sensortypes', 'MLP57', ... 1413 'edit', TfOptions, ... 1414 'normalize', 'none'); % None: Save non-standardized time-frequency maps 1415 % Process: Perm t-test equal [-100ms,500ms 1-150Hz] H0:(A=B), H1:(A<>B) 1416 sTestTf = bst_process('CallProcess', 'process_test_permutation2', sEpochDeviantTf, sEpochStandardTf, ... 1417 'timewindow', [-0.1, 0.5], ... 1418 'freqrange', [1, 150], ... 1419 'rows', '', ... 1420 'isabs', 0, ... 1421 'avgtime', 0, ... 1422 'avgrow', 0, ... 1423 'avgfreq', 0, ... 1424 'matchrows', 0, ... 1425 'iszerobad', 1, ... 1426 'Comment', '', ... 1427 'test_type', 'ttest_equal', ... % Student's t-test (equal variance) t = (mean(A)-mean(B)) / (Sx * sqrt(1/nA + 1/nB))Sx = sqrt(((nA-1)*var(A) + (nB-1)*var(B)) / (nA+nB-2)) 1428 'randomizations', 1000, ... 1429 'tail', 'two'); % Two-tailed 1430 % Set stat threshold 1431 StatThreshOptions.pThreshold = 0.05; 1432 StatThreshOptions.Correction = 'none'; 1433 bst_set('StatThreshOptions', StatThreshOptions); 1434 % Process: Snapshot: Time-frequency maps 1435 bst_process('CallProcess', 'process_snapshot', sTestTf, [], ... 1436 'target', 14, ... % Time-frequency maps 1437 'Comment', 'Non-parametric t-test (p<0.05, Uncorrected)'); 1438 % Process: Delete intermediate results 1439 bst_process('CallProcess', 'process_delete', [sEpochDeviantTf, sEpochStandardTf], [], ... 1440 'target', 1); % Delete selected files 1441 1442 1443 1444 %% ===== SAVE REPORT ===== 1445 % Save and display report 1446 ReportFile = bst_report('Save', []); 1447 if ~isempty(reports_dir) && ~isempty(ReportFile) 1448 bst_report('Export', ReportFile, reports_dir); 1449 else 1450 bst_report('Open', ReportFile); 1451 end 1452 1453 disp([10 'DEMO> Done.' 10]);

For an example of a script illustrating how to create loops, look at the tutorial MEG visual: group study: brainstorm3/toolbox/script/tutorial_visual_group.m

1 function tutorial_visual_group(ProtocolName, reports_dir) 2 % TUTORIAL_VISUAL_GROUP: Runs the Brainstorm/SPM group analysis pipeline (group analysis, BIDS VERSION) 3 % 4 % ONLINE TUTORIALS: 5 % - https://neuroimage.usc.edu/brainstorm/Tutorials/VisualSingle 6 % - https://neuroimage.usc.edu/brainstorm/Tutorials/VisualGroup 7 % 8 % INPUTS: 9 % - ProtocolName : Name of the protocol in which the recordings for the 19 subjects have been imported (TutorialVisual or TutorialGroup) 10 % This folder must include: the run-level averages (recordings and time-frequency) and the source kernels for each run 11 % - reports_dir : If defined, exports all the reports as HTML to this folder 12 13 % @============================================================================= 14 % This function is part of the Brainstorm software: 15 % https://neuroimage.usc.edu/brainstorm 16 % 17 % Copyright (c) University of Southern California & McGill University 18 % This software is distributed under the terms of the GNU General Public License 19 % as published by the Free Software Foundation. Further details on the GPLv3 20 % license can be found at http://www.gnu.org/copyleft/gpl.html. 21 % 22 % FOR RESEARCH PURPOSES ONLY. THE SOFTWARE IS PROVIDED "AS IS," AND THE 23 % UNIVERSITY OF SOUTHERN CALIFORNIA AND ITS COLLABORATORS DO NOT MAKE ANY 24 % WARRANTY, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO WARRANTIES OF 25 % MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, NOR DO THEY ASSUME ANY 26 % LIABILITY OR RESPONSIBILITY FOR THE USE OF THIS SOFTWARE. 27 % 28 % For more information type "brainstorm license" at command prompt. 29 % =============================================================================@ 30 % 31 % Author: Francois Tadel, Elizabeth Bock, 2016-2018 32 33 % ===== CHECK PROTOCOL ===== 34 % Start Brainstorm without the GUI 35 if ~brainstorm('status') 36 brainstorm nogui 37 end 38 % Output folder for reports 39 if (nargin < 2) || isempty(reports_dir) || ~isdir(reports_dir) 40 reports_dir = []; 41 end 42 % You have to specify the folder in which the tutorial dataset is unzipped 43 if (nargin < 1) || isempty(ProtocolName) 44 ProtocolName = 'TutorialGroup'; 45 end 46 % Select current protocol 47 iProtocol = bst_get('Protocol', ProtocolName); 48 if isempty(iProtocol) 49 error(['Unknown protocol: ' ProtocolName]); 50 end 51 if (iProtocol ~= bst_get('iProtocol')) 52 gui_brainstorm('SetCurrentProtocol', iProtocol); 53 end 54 % Process: Select results files in: sub-01/sub-01_ses-meg_task-facerecognition_run-01_proc-sss_meg_notch 55 sFiles = bst_process('CallProcess', 'process_select_files_results', [], [], ... 56 'subjectname', 'sub-01', ... 57 'condition', 'sub-01_ses-meg_task-facerecognition_run-01_proc-sss_meg_notch'); 58 if isempty(sFiles) 59 error(['No source files available in folder: sub-01/sub-01_ses-meg_task-facerecognition_run-01_proc-sss_meg_notch.' 10 ... 60 'You should run tutorial_visual_single first, or download the ' 10 ... 61 'protocol TutorialGroup.zip from the Brainstorm website.']); 62 end 63 % Process: Select file comments with tag: Unconstr 64 sFiles = bst_process('CallProcess', 'process_select_tag', sFiles, [], 'tag', 'Unconstr'); 65 if ~isempty(sFiles) 66 isUnconstrained = 1; 67 else 68 isUnconstrained = 0; 69 end 70 % Start a new report 71 bst_report('Start'); 72 73 74 % ===== DELETE ALL THE EXISTING FILES ===== 75 GroupSubjectName = bst_get('NormalizedSubjectName'); 76 DirIntra = bst_get('DirAnalysisIntra'); 77 % Delete group analysis 78 [sSubject, iSubject] = bst_get('Subject', GroupSubjectName); 79 if ~isempty(sSubject) 80 db_delete_subjects(iSubject); 81 end 82 % Get all the files in the intra subjects 83 sAvgSrc = bst_process('CallProcess', 'process_select_files_data', [], [], ... 84 'subjectname', 'All', ... 85 'condition', DirIntra, ... 86 'includeintra', 1); 87 % If there are existing files: delete the contents of the folders 88 if ~isempty(sAvgSrc) 89 % Get the base folder of the protocol 90 ProtocolInfo = bst_get('ProtocolInfo', iProtocol); 91 % Get all the folders 92 AllIntraDir = unique(cellfun(@bst_fileparts, {sAvgSrc.FileName}, 'UniformOutput', 0)); 93 % Loop over the intra folders 94 for iSubj = 1:length(AllIntraDir) 95 % Get all the files 96 subjIntraDir = dir(fullfile(ProtocolInfo.STUDIES, AllIntraDir{iSubj}, '*.mat')); 97 subjIntraFiles = setdiff({subjIntraDir.name}, 'brainstormstudy.mat'); 98 % Delete all the files 99 for iFile = 1:length(subjIntraFiles) 100 delete(fullfile(ProtocolInfo.STUDIES, AllIntraDir{iSubj}, subjIntraFiles{iFile})); 101 end 102 end 103 % Reload protocol 104 db_reload_database(iProtocol); 105 end 106 107 108 %% ===== SUBJECT AVERAGES: FAMOUS, UNFAMILIAR, SCRAMBLED ============================== 109 % ==================================================================================== 110 111 % ===== SUBJECT AVERAGE: MEG+EEG ===== 112 % Process: Select data files in: */*/Avg 113 sAvgRun = bst_process('CallProcess', 'process_select_files_data', [], [], ... 114 'subjectname', 'All', ... 115 'tag', 'Avg'); 116 % Process: Weighted Average: By trial group (subject average) 117 sAvgSubj = bst_process('CallProcess', 'process_average', sAvgRun, [], ... 118 'avgtype', 6, ... % By trial group (subject average) 119 'avg_func', 1, ... % Arithmetic average: mean(x) 120 'weighted', 1, ... 121 'keepevents', 1); 122 % Process: Average: By trial group (grand average) 123 sAvgGroup = bst_process('CallProcess', 'process_average', sAvgSubj, [], ... 124 'avgtype', 7, ... % By trial group (grand average) 125 'avg_func', 1, ... % Arithmetic average: mean(x) 126 'weighted', 0, ... 127 'keepevents', 0); 128 129 % ===== SUBJECT AVERAGES: SOURCES (EEG) ===== 130 % Process: Select source files in: */*/EEG 131 sAvgRunSrcEeg = bst_process('CallProcess', 'process_select_files_results', [], [], 'tag', 'EEG'); 132 % Process: Weighted Average: By trial group (subject average) - EEG 133 sAvgSubjSrcEeg = bst_process('CallProcess', 'process_average', sAvgRunSrcEeg, [], ... 134 'avgtype', 6, ... % By trial group (subject average) 135 'avg_func', 1, ... % Arithmetic average: mean(x) 136 'weighted', 1, ... 137 'scalenormalized', 0); 138 % Process: Add tag: EEG 139 sAvgSubjSrcEeg = bst_process('CallProcess', 'process_add_tag', sAvgSubjSrcEeg, [], ... 140 'tag', 'EEG', ... 141 'output', 1); % Add to comment 142 143 % ===== SUBJECT AVERAGES: SOURCES (MEG) ===== 144 % Process: Select source files in: */*/MEG 145 sAvgRunSrcMeg = bst_process('CallProcess', 'process_select_files_results', [], [], 'tag', 'MEG'); 146 % Process: Weighted Average: By trial group (subject average) - MEG 147 sAvgSubjSrcMeg = bst_process('CallProcess', 'process_average', sAvgRunSrcMeg, [], ... 148 'avgtype', 6, ... % By trial group (subject average) 149 'avg_func', 1, ... % Arithmetic average: mean(x) 150 'weighted', 1, ... 151 'scalenormalized', 1); 152 % Process: Add tag: MEG 153 sAvgSubjSrcMeg = bst_process('CallProcess', 'process_add_tag', sAvgSubjSrcMeg, [], ... 154 'tag', 'MEG', ... 155 'output', 1); % Add to comment 156 157 % ===== SUBJECT AVERAGES: TIMEFREQ ===== 158 % Process: Select time-frequency files in: */* 159 sAvgRunTf = bst_process('CallProcess', 'process_select_files_timefreq', [], []); 160 % Process: Weighted Average: By trial group (subject average) 161 sAvgSubjTf = bst_process('CallProcess', 'process_average', sAvgRunTf, [], ... 162 'avgtype', 6, ... % By trial group (subject average) 163 'avg_func', 1, ... % Arithmetic average: mean(x) 164 'weighted', 1, ... 165 'scalenormalized', 1); 166 167 168 169 %% ===== SUBJECT AVERAGES: FACES ====================================================== 170 % ==================================================================================== 171 172 % ===== FACES AVERAGE: MEG+EEG ====== 173 % Process: Select data files in: */* 174 sAllData = bst_process('CallProcess', 'process_select_files_data', [], [], ... 175 'subjectname', 'All', ... 176 'includeintra', 1); 177 % Process: Select file comments with tag: "WAvg: Avg: Famous (6 files)" 178 sAvgSubjFamous = bst_process('CallProcess', 'process_select_tag', sAllData, [], 'tag', 'WAvg: Avg: Famous (6 files)'); 179 % Process: Select file comments with tag: "WAvg: Avg: Unfamiliar (6 files)" 180 sAvgSubjUnfamiliar = bst_process('CallProcess', 'process_select_tag', sAllData, [], 'tag', 'WAvg: Avg: Unfamiliar (6 files)'); 181 % Process: Weighted Average A&B 182 sAvgSubjFaces = bst_process('CallProcess', 'process_average_ab', sAvgSubjFamous, sAvgSubjUnfamiliar, ... 183 'weighted', 1); 184 % Process: Add tag: Faces 185 sAvgSubjFaces = bst_process('CallProcess', 'process_add_tag', sAvgSubjFaces, [], ... 186 'tag', 'Faces', ... 187 'output', 2); % Add to file name 188 % Process: Set comment 189 sAvgSubjFaces = bst_process('CallProcess', 'process_set_comment', sAvgSubjFaces, [], ... 190 'tag', 'WAvg: Avg: Faces', ... 191 'isindex', 0); 192 % Process: Average: By trial group (grand average) 193 sAvgGroupFaces = bst_process('CallProcess', 'process_average', sAvgSubjFaces, [], ... 194 'avgtype', 7, ... % By trial group (grand average) 195 'avg_func', 1, ... % Arithmetic average: mean(x) 196 'weighted', 0, ... 197 'keepevents', 0); 198 199 % ===== FACES AVERAGE: SOURCES (EEG) ====== 200 % Process: Select source files in: */*/EEG 201 sSrcEeg = bst_process('CallProcess', 'process_select_files_results', [], [], ... 202 'subjectname', 'All', ... 203 'tag', 'EEG', ... 204 'includeintra', 1); 205 % Process: Select file comments with tag: "WAvg: Avg: Famous" 206 sAvgSubjFamousSrcEeg = bst_process('CallProcess', 'process_select_tag', sSrcEeg, [], 'tag', 'WAvg: Avg: Famous'); 207 % Process: Select file comments with tag: "WAvg: Avg: Unfamiliar" 208 sAvgSubjUnfamiliarSrcEeg = bst_process('CallProcess', 'process_select_tag', sSrcEeg, [], 'tag', 'WAvg: Avg: Unfamiliar'); 209 % Process: Weighted Average A&B 210 sAvgSubjFacesSrcEeg = bst_process('CallProcess', 'process_average_ab', sAvgSubjFamousSrcEeg, sAvgSubjUnfamiliarSrcEeg, ... 211 'weighted', 1, ... 212 'scalenormalized', 1); 213 % Process: Add tag: Faces 214 sAvgSubjFacesSrcEeg = bst_process('CallProcess', 'process_add_tag', sAvgSubjFacesSrcEeg, [], ... 215 'tag', 'Faces', ... 216 'output', 2); % Add to file name 217 % Process: Set comment 218 sAvgSubjFacesSrcEeg = bst_process('CallProcess', 'process_set_comment', sAvgSubjFacesSrcEeg, [], ... 219 'tag', 'WAvg: Avg: Faces | EEG', ... 220 'isindex', 0); 221 222 % ===== FACES AVERAGE: SOURCES (MEG) ====== 223 % Process: Select source files in: */* 224 sSrcMeg = bst_process('CallProcess', 'process_select_files_results', [], [], ... 225 'subjectname', 'All', ... 226 'tag', 'MEG', ... 227 'includeintra', 1); 228 % Process: Select file comments with tag: Famous 229 sAvgSubjFamousSrcMeg = bst_process('CallProcess', 'process_select_tag', sSrcMeg, [], 'tag', 'Famous'); 230 % Process: Select file comments with tag: Unfamiliar 231 sAvgSubjUnfamiliarSrcMeg = bst_process('CallProcess', 'process_select_tag', sSrcMeg, [], 'tag', 'Unfamiliar'); 232 % Process: Weighted Average A&B 233 sAvgSubjFacesSrcMeg = bst_process('CallProcess', 'process_average_ab', sAvgSubjFamousSrcMeg, sAvgSubjUnfamiliarSrcMeg, ... 234 'weighted', 1, ... 235 'scalenormalized', 1); 236 % Process: Add tag: Faces 237 sAvgSubjFacesSrcMeg = bst_process('CallProcess', 'process_add_tag', sAvgSubjFacesSrcMeg, [], ... 238 'tag', 'Faces', ... 239 'output', 2); % Add to file name 240 % Process: Set comment 241 sAvgSubjFacesSrcMeg = bst_process('CallProcess', 'process_set_comment', sAvgSubjFacesSrcMeg, [], ... 242 'tag', 'WAvg: Avg: Faces | MEG', ... 243 'isindex', 0); 244 245 % ===== FACES AVERAGE: TIMEFREQ ====== 246 % Process: Select time-frequency files in: */* 247 sAllTf = bst_process('CallProcess', 'process_select_files_timefreq', [], [], ... 248 'subjectname', 'All', ... 249 'includeintra', 1); 250 % Process: Select file comments with tag: "WAvg: Avg: Famous" 251 sAvgSubjFamousTf = bst_process('CallProcess', 'process_select_tag', sAllTf, [], 'tag', 'WAvg: Avg: Famous'); 252 % Process: Select file comments with tag: "WAvg: Avg: Unfamiliar" 253 sAvgSubjUnfamiliarTf = bst_process('CallProcess', 'process_select_tag', sAllTf, [], 'tag', 'WAvg: Avg: Unfamiliar'); 254 % Process: Weighted Average A&B 255 sAvgSubjFacesTf = bst_process('CallProcess', 'process_average_ab', sAvgSubjFamousTf, sAvgSubjUnfamiliarTf, ... 256 'weighted', 1); 257 % Process: Add tag: Faces 258 sAvgSubjFacesTf = bst_process('CallProcess', 'process_add_tag', sAvgSubjFacesTf, [], ... 259 'tag', 'Faces', ... 260 'output', 2); % Add to file name 261 % Process: Set comment 262 sAvgSubjFacesTf = bst_process('CallProcess', 'process_set_comment', sAvgSubjFacesTf, [], ... 263 'tag', 'WAvg: Avg: Faces', ... 264 'isindex', 0); 265 266 267 268 %% ===== SUBJECT AVERAGES: WITHIN-SUBJECT DIFFERENCES ================================= 269 % ==================================================================================== 270 271 % ===== FACES - SCRAMBLED: SOURCES EEG ===== 272 % Process: Select source files in: */*/EEG 273 sSrcEeg = bst_process('CallProcess', 'process_select_files_results', [], [], ... 274 'subjectname', 'All', ... 275 'tag', 'EEG', ... 276 'includeintra', 1); 277 % Process: Select file comments with tag: "WAvg: Avg: Faces" 278 sAvgSubjFacesSrcEeg = bst_process('CallProcess', 'process_select_tag', sSrcEeg, [], 'tag', 'WAvg: Avg: Faces'); 279 % Process: Select file comments with tag: "WAvg: Avg: Scrambled" 280 sAvgSubjScrambledSrcEeg = bst_process('CallProcess', 'process_select_tag', sSrcEeg, [], 'tag', 'WAvg: Avg: Scrambled'); 281 % Process: Difference: A-B 282 sDiffFacesSubjSrcEeg = bst_process('CallProcess', 'process_diff_ab', sAvgSubjFacesSrcEeg, sAvgSubjScrambledSrcEeg, ... 283 'source_abs', 0); 284 % Process: Set comment 285 sDiffFacesSubjSrcEeg = bst_process('CallProcess', 'process_set_comment', sDiffFacesSubjSrcEeg, [], ... 286 'tag', 'Faces - Scrambled | EEG', ... 287 'isindex', 0); 288 289 % ===== FACES - SCRAMBLED: SOURCES MEG ===== 290 % Process: Select source files in: */*/MEG 291 sSrcMeg = bst_process('CallProcess', 'process_select_files_results', [], [], ... 292 'subjectname', 'All', ... 293 'tag', 'MEG', ... 294 'includeintra', 1); 295 % Process: Select file comments with tag: "WAvg: Avg: Faces" 296 sAvgSubjFacesSrcMeg = bst_process('CallProcess', 'process_select_tag', sSrcMeg, [], 'tag', 'WAvg: Avg: Faces'); 297 % Process: Select file comments with tag: "WAvg: Avg: Scrambled" 298 sAvgSubjScrambledSrcMeg = bst_process('CallProcess', 'process_select_tag', sSrcMeg, [], 'tag', 'WAvg: Avg: Scrambled'); 299 % Process: Difference: A-B 300 sDiffFacesSubjSrcMeg = bst_process('CallProcess', 'process_diff_ab', sAvgSubjFacesSrcMeg, sAvgSubjScrambledSrcMeg, ... 301 'source_abs', 0); 302 % Process: Set comment 303 sDiffFacesSubjSrcMeg = bst_process('CallProcess', 'process_set_comment', sDiffFacesSubjSrcMeg, [], ... 304 'tag', 'Faces - Scrambled | MEG', ... 305 'isindex', 0); 306 307 % ===== FAMOUS - UNFAMILIAR: SOURCE EEG ===== 308 % Process: Select file comments with tag: "WAvg: Avg: Famous" 309 sAvgSubjFamousSrcEeg = bst_process('CallProcess', 'process_select_tag', sSrcEeg, [], 'tag', 'WAvg: Avg: Famous'); 310 % Process: Select file comments with tag: "WAvg: Avg: Unfamiliar" 311 sAvgSubjUnfamiliarSrcEeg = bst_process('CallProcess', 'process_select_tag', sSrcEeg, [], 'tag', 'WAvg: Avg: Unfamiliar'); 312 % Process: Difference: A-B 313 sDiffFamousSubjSrcEeg = bst_process('CallProcess', 'process_diff_ab', sAvgSubjFamousSrcEeg, sAvgSubjUnfamiliarSrcEeg, ... 314 'source_abs', 0); 315 % Process: Set comment 316 sDiffFamousSubjSrcEeg = bst_process('CallProcess', 'process_set_comment', sDiffFamousSubjSrcEeg, [], ... 317 'tag', 'Famous - Unfamiliar | EEG', ... 318 'isindex', 0); 319 320 % ===== FAMOUS - UNFAMILIAR: SOURCE MEG ===== 321 % Process: Select file comments with tag: "WAvg: Avg: Famous" 322 sAvgSubjFamousSrcMeg = bst_process('CallProcess', 'process_select_tag', sSrcMeg, [], 'tag', 'WAvg: Avg: Famous'); 323 % Process: Select file comments with tag: "WAvg: Avg: Unfamiliar" 324 sAvgSubjUnfamiliarSrcMeg = bst_process('CallProcess', 'process_select_tag', sSrcMeg, [], 'tag', 'WAvg: Avg: Unfamiliar'); 325 % Process: Difference: A-B 326 sDiffFamousSubjSrcMeg = bst_process('CallProcess', 'process_diff_ab', sAvgSubjFamousSrcMeg, sAvgSubjUnfamiliarSrcMeg, ... 327 'source_abs', 0); 328 % Process: Set comment 329 sDiffFamousSubjSrcMeg = bst_process('CallProcess', 'process_set_comment', sDiffFamousSubjSrcMeg, [], ... 330 'tag', 'Famous - Unfamiliar | MEG', ... 331 'isindex', 0); 332 333 334 335 %% ===== SUBJECT AVERAGES: FILTER AND NORMALIZE ======================================= 336 % ==================================================================================== 337 338 % ===== FILTER: MEG+EEG ===== 339 % Process: Select data files in: */Intra 340 sAvgData = bst_process('CallProcess', 'process_select_files_data', [], [], ... 341 'subjectname', 'All', ... 342 'condition', DirIntra, ... 343 'includeintra', 1); 344 % Process: Low-pass:32Hz 345 sAvgData = bst_process('CallProcess', 'process_bandpass', sAvgData, [], ... 346 'sensortypes', 'MEG, EEG', ... 347 'highpass', 0, ... 348 'lowpass', 32, ... 349 'attenuation', 'strict', ... % 60dB 350 'mirror', 0, ... 351 'overwrite', 1); 352 % Process: Extract time: [-202ms,900ms] 353 sAvgData = bst_process('CallProcess', 'process_extract_time', sAvgData, [], ... 354 'timewindow', [-0.2018, 0.9], ... 355 'overwrite', 1); 356 357 % ===== FILTER AND NORMALIZE: SOURCES ===== 358 % Process: Select source files in: */Intra 359 sAvgSrc = bst_process('CallProcess', 'process_select_files_results', [], [], ... 360 'subjectname', 'All', ... 361 'condition', DirIntra, ... 362 'includeintra', 1); 363 % Process: Low-pass:32Hz 364 sAvgSrc = bst_process('CallProcess', 'process_bandpass', sAvgSrc, [], ... 365 'highpass', 0, ... 366 'lowpass', 32, ... 367 'attenuation', 'strict', ... % 60dB 368 'mirror', 0, ... 369 'overwrite', 1); 370 % Process: Extract time: [-202ms,900ms] 371 sAvgSrc = bst_process('CallProcess', 'process_extract_time', sAvgSrc, [], ... 372 'timewindow', [-0.2018, 0.9], ... 373 'overwrite', 1); 374 % Process: Z-score transformation: [-202ms,-5ms] 375 sAvgSrc = bst_process('CallProcess', 'process_baseline_norm', sAvgSrc, [], ... 376 'baseline', [-0.2018, -0.005], ... 377 'source_abs', 0, ... 378 'method', 'zscore', ... % Z-score transformation: x_std = (x - μ) / σ 379 'overwrite', 1); 380 381 % ===== NORMALIZE: TIMEFREQ ===== 382 % Process: Select time-frequency files in: */Intra 383 sAvgTf = bst_process('CallProcess', 'process_select_files_timefreq', [], [], ... 384 'subjectname', 'All', ... 385 'condition', DirIntra, ... 386 'includeintra', 1); 387 % Process: Event-related perturbation (ERS/ERD): [-200ms,-4ms] 388 sAvgTf = bst_process('CallProcess', 'process_baseline_norm', sAvgTf, [], ... 389 'baseline', [-0.2, -0.004], ... 390 'method', 'ersd', ... % Event-related perturbation (ERS/ERD): x_std = (x - μ) / μ * 100 391 'overwrite', 1); 392 393 % Save report 394 ReportFile = bst_report('Save', []); 395 if ~isempty(reports_dir) && ~isempty(ReportFile) 396 bst_report('Export', ReportFile, bst_fullfile(reports_dir, ['report_' ProtocolName '_1prepare.html'])); 397 end 398 399 400 401 %% ===== SUBJECT AVERAGES: SCREEN CAPTURES ============================================ 402 % ==================================================================================== 403 % Start a new report 404 bst_report('Start'); 405 % Set colormap: local color scale 406 bst_colormaps('SetMaxMode', 'stat2', 'local'); 407 408 % ===== FACES: RECORDINGS ====== 409 % Process: Select data files in: */*/WAvg: Avg: Faces | 410 sAvgFaces = bst_process('CallProcess', 'process_select_files_data', [], [], ... 411 'subjectname', 'All', ... 412 'tag', 'WAvg: Avg: Faces |', ... 413 'includeintra', 1); 414 % Process: Snapshot: Recordings time series 415 bst_process('CallProcess', 'process_snapshot', sAvgFaces, [], ... 416 'target', 1, ... % Sensor/MRI registration 417 'modality', 1); % MEG 418 % Process: Snapshot: Recordings time series 419 bst_process('CallProcess', 'process_snapshot', sAvgFaces, [], ... 420 'target', 5, ... % Recordings time series 421 'modality', 4, ... % EEG 422 'time', 0.1055); 423 424 % ===== FACES: SOURCES MEG ====== 425 % Process: Select source files in: */*/WAvg: Avg: Faces | MEG 426 sAvgFacesSrcMeg = bst_process('CallProcess', 'process_select_files_results', [], [], ... 427 'subjectname', 'All', ... 428 'condition', DirIntra, ... 429 'tag', 'WAvg: Avg: Faces | MEG', ... 430 'includeintra', 1); 431 % Process: Snapshot: Sources (one time) 432 bst_process('CallProcess', 'process_snapshot', sAvgFacesSrcMeg, [], ... 433 'target', 8, ... % Sources (one time) 434 'orient', 6, ... % back 435 'time', 0.1055, ... 436 'threshold', 10); 437 438 % ===== FACES: TIME-FREQ EEG ====== 439 % Process: Select time-frequency files in: */*/WAvg: Avg: Faces 440 sAvgFacesTf = bst_process('CallProcess', 'process_select_files_timefreq', [], [], ... 441 'subjectname', 'All', ... 442 'tag', 'WAvg: Avg: Faces', ... 443 'includeintra', 1); 444 % Process: Snapshot: Time-frequency maps 445 bst_process('CallProcess', 'process_snapshot', sAvgFacesTf, [], ... 446 'target', 14, ... % Time-frequency maps 447 'time', 0.1055, ... 448 'rowname', 'EEG070'); 449 450 % Save report 451 ReportFile = bst_report('Save', []); 452 if ~isempty(reports_dir) && ~isempty(ReportFile) 453 bst_report('Export', ReportFile, bst_fullfile(reports_dir, ['report_' ProtocolName '_2snapshots.html'])); 454 end 455 456 457 458 %% ===== GROUP ANALYSIS: FACES-SCRAMBLED: MEG/EEG ===================================== 459 % ==================================================================================== 460 % Start a new report 461 bst_report('Start'); 462 % Set colormap: global color scale 463 bst_colormaps('SetMaxMode', 'meg', 'global'); 464 bst_colormaps('SetMaxMode', 'eeg', 'global'); 465 % Set display properties for statistics: p<0.05, FDR-corrected 466 StatThreshOptions = bst_get('StatThreshOptions'); 467 StatThreshOptions.pThreshold = 0.05; 468 StatThreshOptions.Correction = 'fdr'; 469 StatThreshOptions.Control = [1 2 3]; 470 bst_set('StatThreshOptions', StatThreshOptions); 471 472 % ===== GRAND AVERAGES ===== 473 % Process: Select data files in: Group_analysis/Intra 474 sAvgData = bst_process('CallProcess', 'process_select_files_data', [], [], ... 475 'subjectname', GroupSubjectName, ... 476 'condition', DirIntra, ... 477 'includeintra', 1); 478 % Take screen captures 479 DataScreenCapture(sAvgData(1), 1, 1, ''); 480 DataScreenCapture(sAvgData(2), 1, 1, ''); 481 DataScreenCapture(sAvgData(3), 1, 1, ''); 482 DataScreenCapture(sAvgData(4), 1, 1, ''); 483 % Define function for screen captures 484 function DataScreenCapture(sFiles, isMEG, isEEG, Comment) 485 if isMEG 486 % Process: Snapshot: Recordings time series 487 bst_process('CallProcess', 'process_snapshot', sFiles, [], ... 488 'target', 5, ... % Recordings time series 489 'modality', 3, ... % MEG (Magnetometers) 490 'time', 0.1055, ... 491 'Comment', ['MEG MAG: ' Comment]); 492 % Process: Snapshot: Recordings topography (contact sheet) 493 bst_process('CallProcess', 'process_snapshot', sFiles, [], ... 494 'target', 7, ... % Recordings topography (contact sheet) 495 'modality', 3, ... % MEG (Magnetometers) 496 'contact_time', [0.050, 0.300], ... 497 'contact_nimage', 6, ... 498 'Comment', ['MEG MAG: ' Comment]); 499 end 500 if isEEG 501 % Process: Snapshot: Recordings time series 502 bst_process('CallProcess', 'process_snapshot', sFiles, [], ... 503 'target', 5, ... % Recordings time series 504 'modality', 4, ... % EEG 505 'time', 0.1055, ... 506 'Comment', ['EEG: ' Comment]); 507 % Process: Snapshot: Recordings topography (contact sheet) 508 bst_process('CallProcess', 'process_snapshot', sFiles, [], ... 509 'target', 7, ... % Recordings topography (contact sheet) 510 'modality', 4, ... % EEG 511 'contact_time', [0.050, 0.300], ... 512 'contact_nimage', 6, ... 513 'Comment', ['EEG: ' Comment]); 514 end 515 end 516 517 % ===== FACES-SCRAMBLED: MEAN DIFF ===== 518 % Process: Select data files in: */DirIntra/"Avg: Faces | " 519 sSubjAvgFacesData = bst_process('CallProcess', 'process_select_files_data', [], [], ... 520 'subjectname', 'All', ... 521 'condition', DirIntra, ... 522 'tag', 'Avg: Faces | ', ... 523 'includeintra', 1); 524 % Process: Select data files in: */DirIntra/"WAvg: Avg: Scrambled (6 files)" 525 sSubjAvgScrambledData = bst_process('CallProcess', 'process_select_files_data', [], [], ... 526 'subjectname', 'All', ... 527 'condition', DirIntra, ... 528 'tag', 'WAvg: Avg: Scrambled (6 files)', ... 529 'includeintra', 1); 530 % Process: Difference of means [mean] 531 sDiffFacesData = bst_process('CallProcess', 'process_diff_mean', sSubjAvgFacesData, sSubjAvgScrambledData, ... 532 'avg_func', 1, ... % Arithmetic average mean(A) - mean(B) 533 'weighted', 0); 534 % Process: Set comment 535 sDiffFacesData = bst_process('CallProcess', 'process_set_comment', sDiffFacesData, [], ... 536 'tag', 'Faces - Scrambled', ... 537 'isindex', 0); 538 % Take screen captures 539 DataScreenCapture(sDiffFacesData, 1, 1, sDiffFacesData.Comment); 540 541 % ===== FACES-SCRAMBLED: PARAMETRIC T-TEST ===== 542 % Process: Select data files in: */DirIntra/"Avg: Faces | " 543 sSubjAvgFacesData = bst_process('CallProcess', 'process_select_files_data', [], [], ... 544 'subjectname', 'All', ... 545 'condition', DirIntra, ... 546 'tag', 'Avg: Faces | ', ... 547 'includeintra', 1); 548 % Process: Select data files in: */DirIntra/"WAvg: Avg: Scrambled (6 files)" 549 sSubjAvgScrambledData = bst_process('CallProcess', 'process_select_files_data', [], [], ... 550 'subjectname', 'All', ... 551 'condition', DirIntra, ... 552 'tag', 'WAvg: Avg: Scrambled (6 files)', ... 553 'includeintra', 1); 554 % Process: t-test paired [ALL] H0:(A=B), H1:(A<>B) 555 sStatParamFacesData = bst_process('CallProcess', 'process_test_parametric2p', sSubjAvgFacesData, sSubjAvgScrambledData, ... 556 'timewindow', [], ... 557 'sensortypes', '', ... 558 'isabs', 0, ... 559 'avgtime', 0, ... 560 'avgrow', 0, ... 561 'Comment', '', ... 562 'test_type', 'ttest_paired', ... % Paired Student's t-test (A-B)~N(m,v)t = mean(A-B) / std(A-B) * sqrt(n) df=n-1 563 'tail', 'two'); % Two-tailed 564 % Process: Set comment 565 sStatParamFacesData = bst_process('CallProcess', 'process_set_comment', sStatParamFacesData, [], ... 566 'tag', 'Faces - Scrambled: Parametric t-test', ... 567 'isindex', 0); 568 % Take screen captures 569 DataScreenCapture(sStatParamFacesData, 1, 1, sStatParamFacesData.Comment); 570 571 % ===== FACES-SCRAMBLED: NON-PARAMETRIC T-TEST ===== 572 % Process: Perm t-test paired [All] H0:(A=B), H1:(A<>B) 573 sStatPermFacesData = bst_process('CallProcess', 'process_test_permutation2p', sSubjAvgFacesData, sSubjAvgScrambledData, ... 574 'timewindow', [], ... 575 'sensortypes', '', ... 576 'isabs', 0, ... 577 'avgtime', 0, ... 578 'avgrow', 0, ... 579 'iszerobad', 0, ... 580 'Comment', '', ... 581 'test_type', 'ttest_paired', ... % Paired Student's t-test T = mean(A-B) / std(A-B) * sqrt(n) 582 'randomizations', 1000, ... 583 'tail', 'two'); % Two-tailed 584 % Process: Set comment 585 sStatPermFacesData = bst_process('CallProcess', 'process_set_comment', sStatPermFacesData, [], ... 586 'tag', 'Faces - Scrambled: Permutation t-test', ... 587 'isindex', 0); 588 % Take screen captures 589 DataScreenCapture(sStatPermFacesData, 1, 1, sStatPermFacesData.Comment); 590 591 % ===== FACES-SCRAMBLED: CLUSTER T-TEST ===== 592 % Process: FT t-test paired cluster [all EEG] H0:(A=B), H1:(A<>B) 593 sStatClustFacesData = bst_process('CallProcess', 'process_ft_timelockstatistics', sSubjAvgFacesData, sSubjAvgScrambledData, ... 594 'sensortypes', 'EEG', ... 595 'timewindow', [], ... 596 'isabs', 0, ... 597 'avgtime', 0, ... 598 'avgchan', 0, ... 599 'randomizations', 1000, ... 600 'statistictype', 2, ... % Paired t-test 601 'tail', 'two', ... % Two-tailed 602 'correctiontype', 2, ... % cluster 603 'minnbchan', 0, ... 604 'clusteralpha', 0.05); 605 % If FieldTrip is available and results were returned 606 if ~isempty(sStatClustFacesData) 607 % Process: Set comment 608 sStatClustFacesData = bst_process('CallProcess', 'process_set_comment', sStatClustFacesData, [], ... 609 'tag', 'Faces - Scrambled: Cluster t-test EEG', ... 610 'isindex', 0); 611 % Take screen captures 612 DataScreenCapture(sStatClustFacesData, 0, 1, sStatClustFacesData.Comment); 613 end 614 615 616 %% ===== GROUP ANALYSIS: FACES-SCRAMBLED: MEG/EEG ===================================== 617 % ==================================================================================== 618 % ===== FAMOUS-UNFAMILIAR: MEAN DIFF ===== 619 % Process: Select data files in: */DirIntra/"Avg: Famous | " 620 sSubjAvgFamousData = bst_process('CallProcess', 'process_select_files_data', [], [], ... 621 'subjectname', 'All', ... 622 'condition', DirIntra, ... 623 'tag', 'Avg: Famous (6 files)', ... 624 'includeintra', 1); 625 % Process: Select data files in: */DirIntra/"WAvg: Avg: Unfamiliar (6 files)" 626 sSubjAvgUnfamiliarData = bst_process('CallProcess', 'process_select_files_data', [], [], ... 627 'subjectname', 'All', ... 628 'condition', DirIntra, ... 629 'tag', 'WAvg: Avg: Unfamiliar (6 files)', ... 630 'includeintra', 1); 631 % Process: Difference of means [mean] 632 sDiffFamousData = bst_process('CallProcess', 'process_diff_mean', sSubjAvgFamousData, sSubjAvgUnfamiliarData, ... 633 'avg_func', 1, ... % Arithmetic average mean(A) - mean(B) 634 'weighted', 0); 635 % Process: Set comment 636 sDiffFamousData = bst_process('CallProcess', 'process_set_comment', sDiffFamousData, [], ... 637 'tag', 'Famous - Unfamiliar', ... 638 'isindex', 0); 639 % Take screen captures 640 DataScreenCapture(sDiffFamousData, 1, 1, sDiffFamousData.Comment); 641 642 % ===== FAMOUS-UNFAMILIAR: PARAMETRIC T-TEST ===== 643 % Process: t-test paired [All] H0:(A=B), H1:(A<>B) 644 sStatParamFamousData = bst_process('CallProcess', 'process_test_parametric2p', sSubjAvgFamousData, sSubjAvgUnfamiliarData, ... 645 'timewindow', [], ... 646 'sensortypes', '', ... 647 'isabs', 0, ... 648 'avgtime', 0, ... 649 'avgrow', 0, ... 650 'Comment', '', ... 651 'test_type', 'ttest_paired', ... % Paired Student's t-test (A-B)~N(m,v)t = mean(A-B) / std(A-B) * sqrt(n) df=n-1 652 'tail', 'two'); % Two-tailed 653 % Process: Set comment 654 sStatParamFamousData = bst_process('CallProcess', 'process_set_comment', sStatParamFamousData, [], ... 655 'tag', 'Famous - Unfamiliar: Parametric t-test', ... 656 'isindex', 0); 657 % Take screen captures 658 DataScreenCapture(sStatParamFamousData, 1, 1, sStatParamFamousData.Comment); 659 660 % ===== FAMOUS-UNFAMILIAR: CLUSTER T-TEST ===== 661 % Process: FT t-test paired cluster [all EEG] H0:(A=B), H1:(A<>B) 662 sStatClustFamousData = bst_process('CallProcess', 'process_ft_timelockstatistics', sSubjAvgFamousData, sSubjAvgUnfamiliarData, ... 663 'sensortypes', 'EEG', ... 664 'timewindow', [], ... 665 'isabs', 0, ... 666 'avgtime', 0, ... 667 'avgchan', 0, ... 668 'randomizations', 1000, ... 669 'statistictype', 2, ... % Paired t-test 670 'tail', 'two', ... % Two-tailed 671 'correctiontype', 2, ... % cluster 672 'minnbchan', 0, ... 673 'clusteralpha', 0.05); 674 % If FieldTrip is available and results were returned 675 if ~isempty(sStatClustFacesData) 676 % Process: Set comment 677 sStatClustFamousData = bst_process('CallProcess', 'process_set_comment', sStatClustFamousData, [], ... 678 'tag', 'Famous - Unfamiliar: Cluster t-test EEG', ... 679 'isindex', 0); 680 % Take screen captures 681 DataScreenCapture(sStatClustFamousData, 0, 1, sStatClustFamousData.Comment); 682 end 683 684 % Save report 685 ReportFile = bst_report('Save', []); 686 if ~isempty(reports_dir) && ~isempty(ReportFile) 687 bst_report('Export', ReportFile, bst_fullfile(reports_dir, ['report_' ProtocolName '_3meeg.html'])); 688 end 689 690 691 692 %% ===== FACES-SCRAMBLED: SOURCES (MEG) =============================================== 693 % ==================================================================================== 694 % Start a new report 695 bst_report('Start'); 696 697 % ===== PROJECT SOURCES ===== 698 % Process: Select source files in: */Intra 699 sAvgSrc = bst_process('CallProcess', 'process_select_files_results', [], [], ... 700 'subjectname', 'All', ... 701 'condition', DirIntra, ... 702 'includeintra', 1); 703 % Process: Absolute values / Norm for unconstrained 704 if isUnconstrained 705 % Process: Unconstrained to flat map 706 sAvgSrc = bst_process('CallProcess', 'process_source_flat', sAvgSrc, [], ... 707 'method', 1); % Norm: sqrt(x^2+y^2+z^2) 708 else 709 sAvgSrc = bst_process('CallProcess', 'process_absolute', sAvgSrc, [], ... 710 'overwrite', 1); 711 end 712 % Process: Project on default anatomy 713 sProjSrc = bst_process('CallProcess', 'process_project_sources', sAvgSrc, [], ... 714 'headmodeltype', 'surface'); % Cortex surface 715 % Process: Spatial smoothing (3.00,abs) 716 sProjSrc = bst_process('CallProcess', 'process_ssmooth_surfstat', sProjSrc, [], ... 717 'fwhm', 3, ... 718 'overwrite', 1, ... 719 'source_abs', 1); 720 721 % ===== MOVE FILES ===== 722 % File tags / destination folder 723 ListTags = {'Avg: Faces | EEG', 'Faces_EEG'; ... 724 'Avg: Faces | MEG', 'Faces_MEG'; ... 725 'Avg: Famous (6 files) | EEG', 'Famous_EEG'; ... 726 'Avg: Famous (6 files) | MEG', 'Famous_MEG'; ... 727 'Avg: Unfamiliar (6 files) | EEG', 'Unfamiliar_EEG'; ... 728 'Avg: Unfamiliar (6 files) | MEG', 'Unfamiliar_MEG'; ... 729 'Avg: Scrambled (6 files) | EEG', 'Scrambled_EEG'; ... 730 'Avg: Scrambled (6 files) | MEG', 'Scrambled_MEG'; ... 731 'Faces - Scrambled | EEG', 'Faces-Scrambled_EEG'; ... 732 'Faces - Scrambled | MEG', 'Faces-Scrambled_MEG'; ... 733 'Famous - Unfamiliar | EEG', 'Famous-Unfamiliar_EEG'; ... 734 'Famous - Unfamiliar | MEG', 'Famous-Unfamiliar_MEG'}; 735 % Move group by group 736 for i = 1:size(ListTags,1) 737 % Find names of files to move 738 iTag = find(~cellfun(@(c)isempty(strfind(c, ListTags{i,1})), {sProjSrc.Comment})); 739 if isempty(iTag) 740 continue; 741 end 742 % Process: Move files 743 bst_process('CallProcess', 'process_movefile', {sProjSrc(iTag).FileName}, [], ... 744 'subjectname', GroupSubjectName, ... 745 'folder', ListTags{i,2}); 746 end 747 748 749 % ===== LOOP ON MODALITY ===== 750 AllModalities = {'MEG','EEG'}; 751 for iMod = 1:length(AllModalities) 752 % Selected modality for this loop 753 Mod = AllModalities{iMod}; 754 755 % ===== |FACES-SCRAMBLED|: MEAN DIFF ====== 756 % Process: Select source files in: Group_analysis/Faces-Scrambled_MOD 757 sDiffFaces = bst_process('CallProcess', 'process_select_files_results', [], [], ... 758 'subjectname', GroupSubjectName, ... 759 'condition', ['Faces-Scrambled_', Mod]); 760 % Process: Weighted Average: By folder (grand average) 761 sAbsDiffFaces = bst_process('CallProcess', 'process_average', sDiffFaces, [], ... 762 'avgtype', 4, ... % By folder (grand average) 763 'avg_func', 1, ... % Arithmetic average: mean(x) 764 'weighted', 0, ... 765 'scalenormalized', 0); 766 % Process: Set comment 767 sAbsDiffFaces = bst_process('CallProcess', 'process_set_comment', sAbsDiffFaces, [], ... 768 'tag', ['mean(|Faces-Scrambled|) | ' Mod], ... 769 'isindex', 0); 770 % Process: Move files 771 sAbsDiffFaces = bst_process('CallProcess', 'process_movefile', sAbsDiffFaces, [], ... 772 'subjectname', GroupSubjectName, ... 773 'folder', DirIntra); 774 % Colormap: Custom max: [-10,+10] Z 775 bst_colormaps('SetMaxCustom', 'stat2', [], -10, 10); 776 % Process: Snapshot: Sources (contact sheet) 777 bst_process('CallProcess', 'process_snapshot', sAbsDiffFaces, [], ... 778 'target', 9, ... % Sources (contact sheet) 779 'orient', 4, ... % bottom 780 'contact_time', [0.05, 0.4], ... 781 'contact_nimage', 16, ... 782 'threshold', 30, ... 783 'Comment', [Mod ': mean(|Faces-Scrambled|)']); 784 % Process: Snapshot: Sources (contact sheet) 785 bst_process('CallProcess', 'process_snapshot', sAbsDiffFaces, [], ... 786 'target', 9, ... % Sources (contact sheet) 787 'orient', 2, ... % right 788 'contact_time', [0.05, 0.4], ... 789 'contact_nimage', 16, ... 790 'threshold', 30, ... 791 'Comment', [Mod ': mean(|Faces-Scrambled|)']); 792 793 % ===== |FACES-SCRAMBLED|: PARAMETRIC CHI2-TEST ====== 794 % Process: Select source files in: Group_analysis/Faces-Scrambled_MOD 795 sDiffFaces = bst_process('CallProcess', 'process_select_files_results', [], [], ... 796 'subjectname', GroupSubjectName, ... 797 'condition', ['Faces-Scrambled_' Mod]); 798 % Process: Chi2-test [all] H0:(|Zi| = 0) 799 sChiParamFaces = bst_process('CallProcess', 'process_test_parametric1', sDiffFaces, [], ... 800 'timewindow', [], ... 801 'scoutsel', {}, ... 802 'scoutfunc', 1, ... % Mean 803 'isnorm', 0, ... 804 'avgtime', 0, ... 805 'Comment', '', ... 806 'test_type', 'chi2_onesample', ... % One-sample Chi2 test Zi~N(0,1), i=1..nQ = sum(|Zi|^2) Q~Chi2(n) 807 'tail', 'two'); % Two-tailed 808 % Process: Set comment 809 sChiParamFaces = bst_process('CallProcess', 'process_set_comment', sChiParamFaces, [], ... 810 'tag', ['|Faces-Scrambled|=0: Parametric Chi2 test | ' Mod], ... 811 'isindex', 0); 812 % Process: Move files 813 sChiParamFaces = bst_process('CallProcess', 'process_movefile', sChiParamFaces, [], ... 814 'subjectname', GroupSubjectName, ... 815 'folder', DirIntra); 816 % Set colormap: global color scale 817 bst_colormaps('SetMaxMode', 'stat2', 'global'); 818 % Process: Snapshot: Sources (contact sheet) 819 bst_process('CallProcess', 'process_snapshot', sChiParamFaces, [], ... 820 'target', 9, ... % Sources (contact sheet) 821 'orient', 4, ... % bottom 822 'contact_time', [0.05, 0.4], ... 823 'contact_nimage', 16, ... 824 'Comment', [Mod ': |Faces-Scrambled|=0: Parametric Chi2-test']); 825 % Process: Snapshot: Sources (contact sheet) 826 bst_process('CallProcess', 'process_snapshot', sChiParamFaces, [], ... 827 'target', 9, ... % Sources (contact sheet) 828 'orient', 2, ... % right 829 'contact_time', [0.05, 0.4], ... 830 'contact_nimage', 16, ... 831 'Comment', [Mod ': |Faces-Scrambled|=0: Parametric Chi2-test']); 832 833 % ===== LOG(|FACES-SCRAMBLED|): PARAMETRIC CHI2-TEST ====== 834 % Process: Select source files in: Group_analysis/Faces-Scrambled_MOD 835 sDiffFaces = bst_process('CallProcess', 'process_select_files_results', [], [], ... 836 'subjectname', GroupSubjectName, ... 837 'condition', ['Faces-Scrambled_' Mod]); 838 % Process: Duplicate folders: Add tag "_log" 839 sDiffFacesLog = bst_process('CallProcess', 'process_duplicate', sDiffFaces, [], ... 840 'target', 2, ... % Duplicate folders 841 'tag', '_log'); 842 % Process: Run Matlab command 843 sDiffFacesLog = bst_process('CallProcess', 'process_matlab_eval', sDiffFacesLog, [], ... 844 'matlab', 'Data = log(Data);', ... 845 'overwrite', 1); 846 % Process: Chi2-test [all] H0:(|Zi| = 0) 847 sChiParamFacesLog = bst_process('CallProcess', 'process_test_parametric1', sDiffFacesLog, [], ... 848 'timewindow', [], ... 849 'scoutsel', {}, ... 850 'scoutfunc', 1, ... % Mean 851 'isnorm', 0, ... 852 'avgtime', 0, ... 853 'Comment', '', ... 854 'test_type', 'chi2_onesample', ... % One-sample Chi2 test Zi~N(0,1), i=1..nQ = sum(|Zi|^2) Q~Chi2(n) 855 'tail', 'two'); % Two-tailed 856 % Process: Set comment 857 sChiParamFacesLog = bst_process('CallProcess', 'process_set_comment', sChiParamFacesLog, [], ... 858 'tag', ['log(|Faces-Scrambled|)=0: Parametric Chi2 test | ' Mod], ... 859 'isindex', 0); 860 % Process: Move files 861 sChiParamFacesLog = bst_process('CallProcess', 'process_movefile', sChiParamFacesLog, [], ... 862 'subjectname', GroupSubjectName, ... 863 'folder', DirIntra); 864 % Set colormap: global color scale 865 bst_colormaps('SetMaxMode', 'stat2', 'global'); 866 % Process: Snapshot: Sources (contact sheet) 867 bst_process('CallProcess', 'process_snapshot', sChiParamFacesLog, [], ... 868 'target', 9, ... % Sources (contact sheet) 869 'orient', 4, ... % bottom 870 'contact_time', [0.05, 0.4], ... 871 'contact_nimage', 16, ... 872 'Comment', [Mod ': log(|Faces-Scrambled|)=0: Parametric Chi2-test']); 873 % Process: Snapshot: Sources (contact sheet) 874 bst_process('CallProcess', 'process_snapshot', sChiParamFacesLog, [], ... 875 'target', 9, ... % Sources (contact sheet) 876 'orient', 2, ... % right 877 'contact_time', [0.05, 0.4], ... 878 'contact_nimage', 16, ... 879 'Comment', [Mod ': log(|Faces-Scrambled|)=0: Parametric Chi2-test']); 880 881 882 % ===== |FACES|-|SCRAMBLED|: MEAN DIFF ====== 883 % Process: Select source files in: Group_analysis/Faces_MOD 884 sAvgFaces = bst_process('CallProcess', 'process_select_files_results', [], [], ... 885 'subjectname', GroupSubjectName, ... 886 'condition', ['Faces_' Mod]); 887 % Process: Select source files in: Group_analysis/Scrambled_MOD 888 sAvgScrambled = bst_process('CallProcess', 'process_select_files_results', [], [], ... 889 'subjectname', GroupSubjectName, ... 890 'condition', ['Scrambled_' Mod]); 891 % Process: Difference of means [abs(mean)] 892 sDiffAbsFaces = bst_process('CallProcess', 'process_diff_mean', sAvgFaces, sAvgScrambled, ... 893 'avg_func', 2, ... % Absolute value of average abs(mean(A)) - abs(mean(B)) 894 'weighted', 0); 895 % Process: Set comment 896 sDiffAbsFaces = bst_process('CallProcess', 'process_set_comment', sDiffAbsFaces, [], ... 897 'tag', ['mean(|Faces|)-mean(|Scrambled|) | ' Mod], ... 898 'isindex', 0); 899 % Colormap: Custom max: [-10,+10] Z 900 bst_colormaps('SetMaxCustom', 'stat2', [], -10, 10); 901 % Process: Snapshot: Sources (contact sheet) 902 bst_process('CallProcess', 'process_snapshot', sDiffAbsFaces, [], ... 903 'target', 9, ... % Sources (contact sheet) 904 'orient', 4, ... % bottom 905 'contact_time', [0.05, 0.4], ... 906 'contact_nimage', 16, ... 907 'threshold', 30, ... 908 'Comment', [Mod ': mean(|Faces|)-mean(|Scrambled|)']); 909 % Process: Snapshot: Sources (contact sheet) 910 bst_process('CallProcess', 'process_snapshot', sDiffAbsFaces, [], ... 911 'target', 9, ... % Sources (contact sheet) 912 'orient', 2, ... % right 913 'contact_time', [0.05, 0.4], ... 914 'contact_nimage', 16, ... 915 'threshold', 30, ... 916 'Comment', [Mod ': mean(|Faces|)-mean(|Scrambled|)']); 917 918 % ===== |FACES|-|SCRAMBLED|: PARAMETRIC T-TEST ====== 919 % Process: Select source files in: Group_analysis/Faces_MOD 920 sAvgFaces = bst_process('CallProcess', 'process_select_files_results', [], [], ... 921 'subjectname', GroupSubjectName, ... 922 'condition', ['Faces_' Mod]); 923 % Process: Select source files in: Group_analysis/Scrambled_MOD 924 sAvgScrambled = bst_process('CallProcess', 'process_select_files_results', [], [], ... 925 'subjectname', GroupSubjectName, ... 926 'condition', ['Scrambled_' Mod]); 927 % Process: t-test paired [-202ms,900ms] H0:(A=B), H1:(A<>B) 928 sTtestParamFaces = bst_process('CallProcess', 'process_test_parametric2p', sAvgFaces, sAvgScrambled, ... 929 'timewindow', [], ... 930 'scoutsel', {}, ... 931 'scoutfunc', 1, ... % Mean 932 'isnorm', 0, ... 933 'avgtime', 0, ... 934 'Comment', '', ... 935 'test_type', 'ttest_paired', ... % Paired Student's t-test (A-B)~N(m,v)t = mean(A-B) / std(A-B) * sqrt(n) df=n-1 936 'tail', 'two'); % Two-tailed 937 % Process: Set comment 938 sTtestParamFaces = bst_process('CallProcess', 'process_set_comment', sTtestParamFaces, [], ... 939 'tag', ['|Faces|=|Scrambled|: Parametric t-test | ' Mod], ... 940 'isindex', 0); 941 % Set colormap: global color scale 942 bst_colormaps('SetMaxMode', 'stat2', 'global'); 943 % Process: Snapshot: Sources (contact sheet) 944 bst_process('CallProcess', 'process_snapshot', sTtestParamFaces, [], ... 945 'target', 9, ... % Sources (contact sheet) 946 'orient', 4, ... % bottom 947 'contact_time', [0.05, 0.4], ... 948 'contact_nimage', 16, ... 949 'Comment', [Mod ': |Faces|-|Scrambled|: Parametric t-test']); 950 % Process: Snapshot: Sources (contact sheet) 951 bst_process('CallProcess', 'process_snapshot', sTtestParamFaces, [], ... 952 'target', 9, ... % Sources (contact sheet) 953 'orient', 2, ... % right 954 'contact_time', [0.05, 0.4], ... 955 'contact_nimage', 16, ... 956 'Comment', [Mod ': |Faces|-|Scrambled|: Parametric t-test']); 957 end 958 959 960 961 % Save report 962 ReportFile = bst_report('Save', []); 963 if ~isempty(reports_dir) && ~isempty(ReportFile) 964 bst_report('Export', ReportFile, bst_fullfile(reports_dir, ['report_' ProtocolName '_4sources.html'])); 965 end 966 967 968 969 970 end 971

Report viewer

Click on Run to start the script.

As this process is taking screen captures, do not use your computer for something else at the same time: if another window covers the Brainstorm figures, it will not capture the right images.

At the end, the report viewer is opened to show the status of all the processes, the information messages, the list of input and output files, and the screen captures. The report is saved in your home folder ($home/.brainstorm/reports). If you close this window, you can get it back with the menu File > Report viewer.








Feedback: Comments, bug reports, suggestions, questions
Email address (if you expect an answer):


Tutorials/Scripting (last edited 2016-07-20 22:35:14 by FrancoisTadel)