|
Size: 6406
Comment:
|
Size: 8853
Comment:
|
| Deletions are marked like this. | Additions are marked like this. |
| Line 1: | Line 1: |
| = Tutorial 20: Advanced scripting = | = Tutorial 28: Scripting = '''[TUTORIAL UNDER DEVELOPMENT: NOT READY FOR PUBLIC USE] ''' |
| Line 3: | Line 5: |
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. |
|
| Line 6: | Line 10: |
| = From CTF = The main window includes a graphical batching interface that directly benefits from the database explorer: files are organized as a tree of subjects and conditions, and simple drag-and-drop operations readily select files for subsequent batch processing. Most of the Brainstorm features are available through this interface, including pre-processing of the recordings, averaging, time-frequency decompositions, and computing statistics. A full analysis pipeline can be created in a few minutes, saved in the user’s preferences and reloaded in one click, executed directly or exported as a Matlab script. |
== 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). |
| Line 9: | Line 13: |
| The available processes are organized in a plug-in structure. Any Matlab script that is added to the plug-in folder (brainstorm3/toolbox/process/functions/) and has the right format will be automatically detected and made available in the GUI. This mechanism makes the contribution from other developers to Brainstorm very easy. | * '''Prototype''': Start by processing one or two subjects completely '''interactively''' (exactly like in the introduction tutorials). Use the few pilot subjects that you have for your study to prototype the analysis pipeline and check manually all the intermediate stages. Take notes of what you're doing along the way, so that you can later write a script that reproduces the same operations. * '''Anatomical fiducials''': Set NAS/LPA/RPA and compute the MNI transformation for each subject. * '''Segmentation''': Run FreeSurfer/BrainSuite to get surfaces and atlases for all the subjects. * '''File > Batch MRI fiducials''': This menu prompts for the selection of the fiducials for all the subjects and saves a file __fiducials.m__ in each segmentation folder. You will not have to redo this even if you have to start over your analysis from the beginning. * '''Script''': Write a loop that calls the process "Import anatomy folder" for all the subjects. * '''Alternatives''': Create and import the subjects one by one and set the fiducials at the import time. Or use the default anatomy for all the subjects (or use [[Tutorials/TutWarping|warped templates]]). |
| Line 11: | Line 20: |
| == Creating a pipeline == | * '''Script #1''': Pre-processing: Loop on the subjects and the acquisition runs. * '''Create link to raw files''': Link the subject and noise recordings to the database. * '''Event markers''': Read and group triggers from digital and analog channel, fix stimulation delays * '''Evaluation''': Power spectrum density of the recordings to evaluate their quality. * '''Pre-processing''': Notch filter, sinusoid removal, band-pass filter. * '''Evaluation''': Power spectrum density of the recordings to make sure the filters worked well. * '''Cleanup''': Delete the links to the original files (the filtered ones are copied in the database). * '''Detect artifacts''': Detect heartbeats, Detect eye blinks, Remove simultaneous. * '''Compute SSP''': Heartbeats, Blinks (this selects the first component of each decomposition) * '''Compute ICA''': If you have some artifacts you'd like to remove with ICA (no default selection). * '''Screenshots''': Check the MRI/sensors registration, PSD before and after corrections, SSP. * '''Export the report''': One report per subject, or one report for all the subjects, saved in HTML. * '''Manual inspection #1''': * '''Check the reports''': Information messages (number of events, errors and warnings) and screen captures (registration problems, obvious noisy channels, incorrect SSP topographies). * '''Mark bad channels''': Open the recordings, select the channels and mark them as bad. Or use the process "Set bad channels" to mark the same bad channels in multiple files. * '''Fix the SSP/ICA''': For the suspicious runs: Open the file, adjust the list of blink and cardiac events, remove and recompute the SSP decompositions, manually select the components. * '''Detect other artifacts''': Run the process on all the runs of all the subjects at once (select all the files in Process1 and run the process, or generate the equivalent script). * '''Mark bad segments''': Review the artifacts detected in 1-7Hz and 40-240Hz, keep only the ones you really want to remove, then mark the event categories as bad. Review quickly the rest of the file and check that there are no other important artifacts. * '''Additional SSP''': If you find one type of artifact that repeats (typically saccades and SQUID jumps), you can create additional SSP projectors, either with the process "SSP: Generic" or directly from a topography figure (right-click on the figure > Snapshot> Use as SSP projector). * '''Script #2''': Subject-level analysis: Epoching, averaging, sources, time-frequency. * '''Importing''': Process "Import MEG/EEG: Events" and "Pre-process > Remove DC offset". * '''Averaging''': Average trials by run, average runs by subject (registration problem in MEG). * '''Noise covariance''': Compute from empty room or resting recordings, copy to other folders. * '''Head model''': Compute for each run, or compute once and copy if the runs are co-registered. * '''Sources''': Compute for each run, average across runs and subjects in source space for MEG. * '''Time-frequency''': Computation with Hilbert transform or Morlet wavelets, then normalize. * '''Screenshots''': Check the quality of all the averages (time series, topographies, sources). * '''Export the report''': One report per subject, or one report for all the subjects, saved in HTML. * '''Manual inspection #2''': * '''Check the reports''': Check the number of epochs imported and averaged in each condition, check the screen capture of the averages (all the primary responses should be clearly visible). * '''Regions of interest''': If not using predefined regions from an atlas, define the scouts on the anatomy of each subject (or on the template and then project them to the subjects). * '''Script #3''': Group analysis, ROI-based analysis, etc. * '''Averaging''': Group averages for the sensor data, the sources and the time-frequency maps. * '''Statistics''': Contrast between conditions or groups of subjects. * '''Regions of interest''': Any operation that involve scouts. |
| Line 13: | Line 56: |
| === List of processes === 1. Click on Run. The Process selection window appears, with which you can create an analysis pipeline (ie. a list of process that are applied on the selected files one after the other). The first button in the toolbar shows the list of processed that are currently available. If you click on a menu, it's added to the list. |
== Script generation == http://neuroimage.usc.edu/brainstorm/Tutorials/PipelineEditor |
| Line 16: | Line 59: |
| {{http://neuroimage.usc.edu/brainstorm/Tutorials/TutProcesses?action=AttachFile&do=get&target=pipeline1.gif|pipeline1.gif|class="attachment"}} 1. Some menus appear in grey (example: Sources > Spatial smoothing). This means that they are not meant to be applied to the type of data that you have in input, or at the end of the current pipeline. The "spatial smoothing" process may only be run on source files. 1. When you select a process, a list of options specific to this process is shown in the window. * To delete a process: Select it and press the ''Delete'' key, or the big cross in the toolbar. * With the "up arrow" and "down arrow" buttons in the toolbar, you can move up/down a process in the pipeline. 1. Now add the following processes, and set their options: * '''Pre-process > Band-pass filter''': 2Hz - 30Hz * In some processes, you can specify the type(s) of sensors on which you want to apply the process. This way you can for instance apply different filters on the EEG and the MEG, if you have both in the same files. * '''Extract > Extract time''': 40.00ms - 49.60ms, overwrite initial file * This will extract from each file a small time window around the main response peak. * Selecting the overwrite option will replace the previous file (bandpass), with the output of this process (bandpass+extract). This option is usually unselected for the first process in the list, then selected automatically. * '''Average > Average over time''': Overwrite initial file * Compute the average over this small time window. |
== Script edition == Add loops, load files, ... |
| Line 30: | Line 62: |
| {{http://neuroimage.usc.edu/brainstorm/Tutorials/TutProcesses?action=AttachFile&do=get&target=pipeline2.gif|pipeline2.gif|class="attachment"}} 1. Save your pipeline: Click on the last button in the toolbar > Save > New... > Type "process_avg45". |
Loops: http://neuroimage.usc.edu/forums/showthread.php?2429-Problem-using-tags |
| Line 33: | Line 64: |
| === Saving/exporting a pipeline === The last button in the the toolbar offers a list of menus to save, load and export the pipelines. |
== File manipulation == * Modify a structure manually: Export to Matlab/Import from Matlab * File manipulation: file_short, file_fullpath, in_bst_*... * Documentation of all file structures: point at the appropriate tutorials * Select files from the database (with bst_get and processes) |
| Line 36: | Line 70: |
| . {{http://neuroimage.usc.edu/brainstorm/Tutorials/TutProcesses?action=AttachFile&do=get&target=pipeline3.gif|pipeline3.gif|class="attachment"}} | == Final script == The following script from the Brainstorm distribution reproduces the introduction tutorials ("Get started"): '''brainstorm3/toolbox/script/tutorial_introduction.m''' |
| Line 38: | Line 73: |
| * '''Load''': List of processes that are saved in the user preferences * '''Load from file''': Import a pipeline from a pipeline_...mat file (previously saved with the menu "Save as Matlab matrix") * '''Save''': Save the pipeline in the user preferences, to be able to access it really fast after * '''Save as Matlab matrix''': Exports the pipeline for a Matlab structure in a .mat file. Allows different users to exchange their analysis pipelines (or a single user between different computers). * '''Generate .m script''': This option generates a human-readable Matlab script that can be re-used for other purposes or modified. * '''Delete''': Remove a pipeline that is saved in the user preferences. * '''Reset options''': Brainstorm saves automatically for each user the options of all the processes. This menu removes all the saved options, and set them back to the default values. |
<<HTML(<div style="border:1px solid black; background-color:#EEEEFF; width:720px; height:500px; overflow:scroll; padding:10px; font-family: Consolas,Menlo,Monaco,Lucida Console,Liberation Mono,DejaVu Sans Mono,Bitstream Vera Sans Mono,Courier New,monospace,sans-serif; font-size: 13px; white-space: pre;">)>><<EmbedContent("http://neuroimage.usc.edu/bst/viewcode.php?f=tutorial_introduction.m")>><<HTML(</div >)>> |
| Line 46: | Line 75: |
| Here is the Matlab script that is generated automatically for this pipeline. | <<BR>>For an example of a script illustrating how to create loops, look at the tutorial [[Tutorials/VisualGroup|MEG visual: group study]]. '''brainstorm3/toolbox/script/tutorial_visual_''''''group.m''' |
| Line 48: | Line 77: |
| . {{http://neuroimage.usc.edu/brainstorm/Tutorials/TutProcesses?action=AttachFile&do=get&target=script.gif|script.gif|class="attachment"}} . |
<<HTML(<div style="border:1px solid black; background-color:#EEEEFF; width:720px; height:500px; overflow:scroll; padding:10px; font-family: Consolas,Menlo,Monaco,Lucida Console,Liberation Mono,DejaVu Sans Mono,Bitstream Vera Sans Mono,Courier New,monospace,sans-serif; font-size: 13px; white-space: pre;">)>><<EmbedContent("http://neuroimage.usc.edu/bst/viewcode.php?f=tutorial_visual_group.m")>><<HTML(</div >)>> |
| Line 51: | Line 79: |
| Click on Ok, in the pipeline window. After a few seconds, you will see two new files in the database, and the "Report viewer" window. | == Report viewer == Click on Run to start the script. |
| Line 53: | Line 82: |
| === Report viewer === Each time the pipeline editor is used to executed to run a list of processes, a report is generated and saved in the user home folder (/home/username/reports/). The report viewer shows as an HTML page some of the information saved in this report structure: the date and duration of execution, the list of processes, the input and output files. It reports all the warning and errors that happen during the execution. |
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. |
| Line 56: | Line 84: |
| The report viewer does not necessarily appear automatically at the end of the last process: it is shown only when more than one processes were executed, or when any of the processes returned an error or a warning. | 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. |
| Line 58: | Line 86: |
| When running processes manually from a script, the calls bst_report(Start, Save, Open) explicitly indicate when the logging of the events should start and stop. You can add images to the reports for quality control using the process "File > Save snapshot". {{http://neuroimage.usc.edu/brainstorm/Tutorials/TutProcesses?action=AttachFile&do=get&target=report.gif|output.gif|class="attachment"}} After you close the report window, you can re-open the last report with the main menu of the Brainstorm window: '''File > Report viewer'''. With the buttons in the toolbar, you can go back to the previous reports saved from the same protocol. <<EmbedContent("http://neuroimage.usc.edu/bst/get_prevnext.php?prev=Tutorials/Connectivity")>> |
<<EmbedContent("http://neuroimage.usc.edu/bst/get_prevnext.php?prev=Tutorials/Workflows")>> |
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.
Contents
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).
Prototype: Start by processing one or two subjects completely interactively (exactly like in the introduction tutorials). Use the few pilot subjects that you have for your study to prototype the analysis pipeline and check manually all the intermediate stages. Take notes of what you're doing along the way, so that you can later write a script that reproduces the same operations.
Anatomical fiducials: Set NAS/LPA/RPA and compute the MNI transformation for each subject.
Segmentation: Run FreeSurfer/BrainSuite to get surfaces and atlases for all the subjects.
File > Batch MRI fiducials: This menu prompts for the selection of the fiducials for all the subjects and saves a file fiducials.m in each segmentation folder. You will not have to redo this even if you have to start over your analysis from the beginning.
Script: Write a loop that calls the process "Import anatomy folder" for all the subjects.
Alternatives: Create and import the subjects one by one and set the fiducials at the import time. Or use the default anatomy for all the subjects (or use warped templates).
Script #1: Pre-processing: Loop on the subjects and the acquisition runs.
Create link to raw files: Link the subject and noise recordings to the database.
Event markers: Read and group triggers from digital and analog channel, fix stimulation delays
Evaluation: Power spectrum density of the recordings to evaluate their quality.
Pre-processing: Notch filter, sinusoid removal, band-pass filter.
Evaluation: Power spectrum density of the recordings to make sure the filters worked well.
Cleanup: Delete the links to the original files (the filtered ones are copied in the database).
Detect artifacts: Detect heartbeats, Detect eye blinks, Remove simultaneous.
Compute SSP: Heartbeats, Blinks (this selects the first component of each decomposition)
Compute ICA: If you have some artifacts you'd like to remove with ICA (no default selection).
Screenshots: Check the MRI/sensors registration, PSD before and after corrections, SSP.
Export the report: One report per subject, or one report for all the subjects, saved in HTML.
Manual inspection #1:
Check the reports: Information messages (number of events, errors and warnings) and screen captures (registration problems, obvious noisy channels, incorrect SSP topographies).
Mark bad channels: Open the recordings, select the channels and mark them as bad. Or use the process "Set bad channels" to mark the same bad channels in multiple files.
Fix the SSP/ICA: For the suspicious runs: Open the file, adjust the list of blink and cardiac events, remove and recompute the SSP decompositions, manually select the components.
Detect other artifacts: Run the process on all the runs of all the subjects at once (select all the files in Process1 and run the process, or generate the equivalent script).
Mark bad segments: Review the artifacts detected in 1-7Hz and 40-240Hz, keep only the ones you really want to remove, then mark the event categories as bad. Review quickly the rest of the file and check that there are no other important artifacts.
Additional SSP: If you find one type of artifact that repeats (typically saccades and SQUID jumps), you can create additional SSP projectors, either with the process "SSP: Generic" or directly from a topography figure (right-click on the figure > Snapshot> Use as SSP projector).
Script #2: Subject-level analysis: Epoching, averaging, sources, time-frequency.
Importing: Process "Import MEG/EEG: Events" and "Pre-process > Remove DC offset".
Averaging: Average trials by run, average runs by subject (registration problem in MEG).
Noise covariance: Compute from empty room or resting recordings, copy to other folders.
Head model: Compute for each run, or compute once and copy if the runs are co-registered.
Sources: Compute for each run, average across runs and subjects in source space for MEG.
Time-frequency: Computation with Hilbert transform or Morlet wavelets, then normalize.
Screenshots: Check the quality of all the averages (time series, topographies, sources).
Export the report: One report per subject, or one report for all the subjects, saved in HTML.
Manual inspection #2:
Check the reports: Check the number of epochs imported and averaged in each condition, check the screen capture of the averages (all the primary responses should be clearly visible).
Regions of interest: If not using predefined regions from an atlas, define the scouts on the anatomy of each subject (or on the template and then project them to the subjects).
Script #3: Group analysis, ROI-based analysis, etc.
Averaging: Group averages for the sensor data, the sources and the time-frequency maps.
Statistics: Contrast between conditions or groups of subjects.
Regions of interest: Any operation that involve scouts.
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
- Modify a structure manually: Export to Matlab/Import from Matlab
- File manipulation: file_short, file_fullpath, in_bst_*...
- Documentation of all file structures: point at the appropriate tutorials
- Select files from the database (with bst_get and processes)
Final script
The following script from the Brainstorm distribution reproduces the introduction tutorials ("Get started"): 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.