Differences between revisions 1 and 55 (spanning 54 versions)

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.

Requirements: You need a license for the Matlab environment in order to use these tools and execute custom scripts. If you are running the compiled version of Brainstorm with the MCR library, the only custom code you can run is through the menu File > Matlab console and the process "Run Matlab command".

Contents

Starting a new script
Line by line: Header
Line by line: Body
Line by line: Footer
Simplify the calls
Interaction with Matlab
Naming conventions
Running the script
Running the script again
Starting Brainstorm
Database requests [TODO]
File manipulation [TODO]
Additional quality control
Loop over subject and runs [TODO]
How to process an entire study
Final script
Report viewer

Starting a new script

The easiest way to get started with a new Brainstorm script is to use the script generator, already introduced in the tutorial Select files and run processes. Select some files in the Process1 or Process2 tabs, select a list of processes, and use the menu Generate .m script. The example below should work on the the protocol "TutorialIntroduction" created during the introduction tutorials.

In the Process1 tab, leave the selection box empty and click on [Run]. Instead of selecting the files from the Brainstorm interface, we will select them directly from the database using a script.
Select process File > Select files: Recordings (do not execute immediately)
Subject=Subject01, Condition=[Empty], File comment=Avg: deviant (the space is important).
This process selects all the recordings with a comment including the string "Avg: deviant", from all the folders in Subject01 (except "Intra-subject" and "Common files"). We expect to get two files: the averages of the deviant condition for both runs.
Add process Pre-process > Band-pass filter: Lower cutoff=0Hz, Upper cutoff=30Hz, Mirror.
Add process File > Save snapshot: Recordings time series, Sensors=MEG.
This will apply a low-pass filter at 30Hz and save a screen capture of the signals in the report.
Do not run the pipeline, select the menu Generate .m script instead. It saves a new .m file and opens it in the Matlab editor. Close the pipeline editor window and look at the script.
The script you just generated can be the starting point to your own custom script. The following sections explain line by line how they work and how to edit them.

Line by line: Header

% Script generated by Brainstorm (19-Jul-2016)

All the lines starting with a "%" are comments, they are never executed.

% Input files
sFiles = [];
SubjectNames = {...
    'Subject01'};

These lines define the script inputs:

sFiles: The list of files in input. Currently empty because we did not select anything in input in the Process1 list. If we had selected files, it would contain a cell array of strings with relative file paths.
SubjectNames: List of subject names that are used in the script. Most of the times, the generated script would contain only one entry, but it written as a cell array so that it is easier to extend it to multiple subjects with a loop (described further in this tutorial).

% Start a new report
bst_report('Start', sFiles);

Start a new report of activity: Clears all the previous logs and gets ready to record new messages. The report will collect all the messages that are generated during the execution of the script by the various processes. You can explicitely add screen captures and additional messages to the current report with the function bst_report. This report will remain open until the function bst_report('Start') is called again. To display the current report, use the menu File > Report viewer.

The syntax function_name('SubFunction', arguments) is used a lot in Brainstorm, to call a subfunction available inside a .m file. This line above calls the function Report() in the file brainstorm3/toolbox/process/bst_report.m. This is made possible with the use of the short script "macro_methodcall" you will see at the top of many files. Many of the Brainstorm .m files are actually libraries of functions, rather than simple "scripts" or "functions".

Line by line: Body

% Process: Select data files in: Subject01/*/Avg: deviant
sFiles = bst_process('CallProcess', 'process_select_files_data', sFiles, [], ...
    'subjectname',   SubjectNames{1}, ...
    'condition',     '', ...
    'tag',           'Avg: deviant', ...
    'includebad',    0, ...
    'includeintra',  0, ...
    'includecommon', 0);

% Process: Low-pass:30Hz
sFiles = bst_process('CallProcess', 'process_bandpass', sFiles, [], ...
    'highpass',    0, ...
    'lowpass',     30, ...
    'mirror',      1, ...
    'sensortypes', 'MEG', ...
    'overwrite',   0);

% Process: Snapshot: Recordings time series
sFiles = bst_process('CallProcess', 'process_snapshot', sFiles, [], ...
    'target',         5, ...  % Recordings time series
    'modality',       1, ...  % MEG (All)
    'orient',         4, ...  % bottom
    'time',           0.11, ...
    'contact_time',   [0, 0.1], ...
    'contact_nimage', 12, ...
    'threshold',      20, ...
    'Comment',        'Run');

You will find one block per process you selected in the pipeline editor. They all have the same syntax:
output_files = bst_process('CallProcess', process_name, input_files_A, input_files_B, options_list);

process_name: String indicating the function corresponding to the process to execute. To know from the pipeline editor what is the path to the process function: hover your mouse over the selected process, as illustrated in this tutorial.
input_files_A: List of input files in Process1, or FilesA in Process2. It can be a cell array of files names (full path, or relative path from the protocol folder), or an array of structures describing the files in the database (returned by a previous call to bst_process).
input_files_B: Empty for Process1, or FilesB in Process2. Cell array of strings or array of struct.
options_list: Pairs of (option_name, option_values), one for each option of the process.
output_files: Array of structures describing the files in output of the process. If the process created new files, this variable contains the new files. If the process didn't create new files or was modifying exiting files, most of the time this variable would contain the same files as the input list.

Line by line: Footer

% Save and display report
ReportFile = bst_report('Save', sFiles);

Closes the current report and saves it in the user Brainstorm report folder ($HOME/.brainstorm/reports). These reports are in .mat format and contain all the information necessary to re-run the execution exactly in the same way, but they not easy to read.

The parameter "sFiles" is optional, it indicates what are the files that are considered as the final results of the script. You can remove it without breaking your script: ReportFile = bst_report('Save');

bst_report('Open', ReportFile);

Opens the report viewer to display what happened during the execution. This is equivalent to using the menu File > Report viewer. You can comment this line (ie. add a "%" at the beginning of the line) if you don't want to show the report at the end of the execution.

% bst_report('Export', ReportFile, ExportDir);

This function exports the report in readable format, as an HTML that includes all the screen captures embedded in it. It is disabled by default. If you want to use this feature: remove the "%" at the beginning of the line, and define the variable ExportDir.

ExportDir must be a string that defines where to save the HTML report. It can be wither the absolute path to a HTML file (eg. 'C:\Users\myuser\Documents\report_example.html') or just a folder (eg. 'C:\Users\myuser\Documents'). If you enter only a path to a folder, a default file name including the protocol name and a date tag is generated (report_ProtocolName_YYMMDD_HHMMSS.html).

Simplify the calls

This script you generated is like any Matlab script: you can edit it, rename the variables, add tests and loops, etc. The first important thing to understand is how to edit the options or change the inputs/outputs of a process. The script generator uses only one variable for all the file lists (sFiles) and the output process is always the input of the following process. This is usually too restrictive to write a full analysis script: we commonly need to have multiple lists of files or to run two different operations on the same file.

Let's consider the first process call, which selects the averages for the Deviant condition in both runs.

sFiles = bst_process('CallProcess', 'process_select_files_data', sFiles, [], ...
    'subjectname',   SubjectNames{1}, ...
    'condition',     '', ...
    'tag',           'Avg: deviant', ...
    'includebad',    0, ...
    'includeintra',  0, ...
    'includecommon', 0);

There is no need to set the parameter sFiles, because there is no input, you can replace it with an empty matrix []. You can therefore remove the line "sFiles = [];". We can also rename the output variable "sAvgData", to be a bit more specific.

sAvgData = bst_process('CallProcess', 'process_select_files_data', [], [], ...

You can omit all the options that are not defined, not used, or kept to their default values:

sAvgData = bst_process('CallProcess', 'process_select_files_data', [], [], ...
    'subjectname',   SubjectNames{1}, ...
    'tag',           'Avg: deviant');

Edit low-pass filter call: Change the input to sAvgData and the output to sAvgDataLow, this way you will be able to keep track of the two files if you need to use them independently later in the script.

sAvgDataLow = bst_process('CallProcess', 'process_bandpass', sAvgData, [], ...
    'highpass',    0, ...
    'lowpass',     30, ...
    'sensortypes', 'MEG');

Edit the call to the snapshot process: Change the input to sAvgDataLow, and remove the output parameter (we are not expecting output from it).

bst_process('CallProcess', 'process_snapshot', sAvgDataLow, [], ...
    'target',   5, ...  % Recordings time series
    'modality', 1);     % MEG (All)

Replace the last lines with the following, to export the report instead of opening in the report viewer (edit the file path to point at your own user folder instead).

ReportFile = bst_report('Save');
bst_report('Export', ReportFile, 'C:\Users\franc\Documents\report_test.html');

Interaction with Matlab

Select the code for the first process in the Matlab editor, right-click > Evaluate selection (or press F9).

If you haven't executed your script yet, you will get the following error in the Matlab command window:

Undefined variable "SubjectNames" or class "SubjectNames".

The variable SubjectNames is not defined yet: Execute the first two lines "SubjectNames = {'Subject01'}", then try again. You should now have a new variable in your Matlab workspace, which points at the two average files. Type "sAvgData(1)" in your command window to display the first element:

>> sAvgData(1)
ans =
          iStudy: 6
           iItem: 1
        FileName: 'Subject01/S01_AEF_20131218_01_600Hz_notch/data_deviant_average_160513_1329.mat'
        FileType: 'data'
         Comment: 'Avg: deviant (39 files)'
       Condition: 'S01_AEF_20131218_01_600Hz_notch'
     SubjectFile: 'Subject01/brainstormsubject.mat'
     SubjectName: 'Subject01'
        DataFile: ''
     ChannelFile: 'Subject01/S01_AEF_20131218_01_600Hz_notch/channel_ctf_acc1.mat'
    ChannelTypes: {'ADC A'  'ADC V'  'DAC'  'ECG'  'EOG'  'FitErr'  'HLU'  'MEG'  'MEG REF' ...}

The field "sAvgData(1).FileName" contains the relative path to the to the Deviant average in the first run. This structure sAvgData contains also a lot of information that can be helpful in your script:

iStudy / iItem: Reference of the file in the database (described later in this tutorial).
FileType: 'raw' (continuous files), 'data' (recordings), 'results' (sources), 'timefreq' (time-frequency, spectrum and connectivity), or 'matrix' (any time series extracted from other files).
Comment: Comment field of the file (what is displayed in the database explorer)
Condition: Name of the condition/folder in which the file is located
SubjectFile: Relative path to the subject file (brainstormsubject.mat)
SubjectName: Name of the subject
DataFile: For types 'results' or 'timefreq', path of the parent file in the database explorer
ChannelFile: Relative path to the channel file
ChannelTypes: Cell array of channel types available for the input file

Naming conventions

To help you navigate in the Brainstorm code, here are some naming conventions:

Structures: Name starting with a "s" followed by a capital letter (eg. sFiles).
Indices: Either loop variables or array indices, name starting with a "i" (eg. iSubject, iStudy, iTime).
Counts: Number of elements in a group, name starting with a "n" (eg. nAvg, nTrials, nSubjects).
File names: Scripts and functions, only lower case, separation with "_" (eg. process_average)
Sub-functions: Inside a .m file, name starting with a capital, camel-case (eg. CallProcess, Start).

Running the script

The simplified script looks like this:

% Input files
SubjectNames = {'Subject01'};
% Start a new report
bst_report('Start');

% Process: Select data files in: Subject01/*/Avg: deviant
sAvgData = bst_process('CallProcess', 'process_select_files_data', [], [], ...
    'subjectname',   SubjectNames{1}, ...
    'tag',           'Avg: deviant');
% Process: Low-pass:30Hz
sAvgDataLow = bst_process('CallProcess', 'process_bandpass', sAvgData, [], ...
    'highpass',    0, ...
    'lowpass',     30, ...
    'sensortypes', 'MEG');
% Process: Snapshot: Recordings time series
bst_process('CallProcess', 'process_snapshot', sAvgDataLow, [], ...
    'target',   5, ...  % Recordings time series
    'modality', 1);     % MEG (All)

% Save and display report
ReportFile = bst_report('Save');
bst_report('Export', ReportFile, 'C:\Users\franc\Documents\report_test.html');

You have three ways to execute it:

Select all the lines (Ctrl+A) and evaluate it in Matlab (F9).
In the Editor toolbar of the Matlab environment, click on the button [Run].
Save the file, go to this folder with Matlab (or add it to your path) and type the name of the script in the command window (without the ".m" at the end).

At the end of the execution, nothing happens, because we indicated we wanted to just export the report instead of opening it. To check out the report of execution: use the menu File > Report viewer from the Brainstorm window, or open the report_test.html that was saved somewhere on your computer.

On this page, you can get review everything that happened in the script: when it was executed, how long it took, what are the processes that were executed, some produced messages (two files were selected with the first process), and at the end you have the screen captures taken by process_snapshot.

Running the script again

If you execute the script again, it will not behave as expected anymore. The selection process we used assumes that there is only one file per folder with a comment that includes "Avg: deviant". This is not the case anymore after the execution, because the low-pass filtered files also contain the same tags. The execution of the first process of the script now returns 4 files.

>> sAvgData = bst_process('CallProcess', 'process_select_files_data', [], [], ...
    'subjectname',   SubjectNames{1}, ...
    'tag',           'Avg: deviant')

sAvgData =
1x4 struct array with fields:
    iStudy
    iItem
    ...

In order to exclude the low-pass filtered files from this selection, you can add another process that will refine the selection. Use the script generator again to create a template call for another process, then copy-paste it in your script.

In Process1: Select any recordings file (we will not run anything, just generate a script).
Select process File > Select files: By tag: Search="low", Search the comment, Ignore the files.
Select the menu Generate .m script (make sure you do not overwrite the script you are currently working on), then close the pipeline editor.
Copy-paste and edit the call to process_select_tag to your main script.

Now the file selection part of your script should look like this, and should return only two files:

% Process: Select data files in: Subject01/*/Avg: Deviant
sAvgData = bst_process('CallProcess', 'process_select_files_data', [], [], ...
    'subjectname',   SubjectNames{1}, ...
    'tag',           'Avg: deviant');
% Process: Ignore file comments with tag: low
sAvgData = bst_process('CallProcess', 'process_select_tag', sAvgData, [], ...
    'tag',    'low', ...
    'search', 2, ...  % Search the file comments
    'select', 2);  % Ignore the files with the tag

With this last modification, your script is more robust. It can be executed multiple times without completely changing its behavior. When you are fetching files from the database using comment tags or file names, always pay attention to this aspect: the database grows and the further you go, the more specific your requests have to be.

A good practice can be to tag explicitely the output files your script generates if you need to fetch them later. You can use the process File > Add tag and File > Set comment.

Starting Brainstorm

Brainstorm must be running in the background for these scripts to run properly. The interface doesn't have to be visible on the screen, but the database engine must be running for processing requests. At the beginning of your script, you can explicitely start or restart Brainstorm.

brainstorm: Start Brainstorm with the regular GUI.

brainstorm nogui: Start in silent mode, without the GUI, without progress bar, without user interactions.

If you want to start Brainstorm only if Brainstorm is not already running, you can use the following code:

if ~brainstorm('status')
    brainstorm nogui
end

To select a specific protocol at the beginning of your script:

ProtocolName = 'TutorialIntroduction';
% Get the protocol index
iProtocol = bst_get('Protocol', ProtocolName);
if isempty(iProtocol)
    error(['Unknown protocol: ' ProtocolName]);
end
% Select the current procotol
gui_brainstorm('SetCurrentProtocol', iProtocol);

To delete the protocol and start over:

% Delete existing protocol
gui_brainstorm('DeleteProtocol', ProtocolName);
% Create new protocol
gui_brainstorm('CreateProtocol', ProtocolName, 0, 0);

Database requests [TODO]

bst_get / bst_set

iStudy: Index of the study, get the structure with sStudy=bst_get('Study', sInput.iStudy)
iItem: Index of file in its category (example for recordings: sStudy.Data(sInput.iItem))

sStudy = bst_get('Study', sAvgData(1).iStudy)

sStudy.Data(sAvgData(1).iItem)

File manipulation [TODO]

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)
Reading / writing / exporting functions

Additional quality control

Create figures and snapshot them

Generate warnings and errors

Loop over subject and runs [TODO]

Creating loops is not supported yet by the script generator, but relatively easy to do from a script without having to know too much about Matlab programming.

1) Fill the cell array SubjectNames with all your subjects names, with the same dimensions as the list of input raw files (sFiles)

2) Add a "for" loop that includes all the bst_process() calls (leave the bst_report() calls and input definition outside)
3) Inside the loop, replace SubjectNames with SubjectNames{i} and sFiles with sFiles(i)

How to process an entire study

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.

Final script

The following script from the Brainstorm distribution reproduces the introduction tutorials ("Get started"): brainstorm3/toolbox/script/tutorial_introduction.m

function tutorial_introduction(tutorial_dir, reports_dir)
% TUTORIAL_INTRODUCTION: Script that runs all the Brainstorm introduction tutorials.
%
% INPUTS: 
%    - tutorial_dir : Directory where the sample_introduction.zip file has been unzipped
%    - reports_dir  : Directory where to save the execution report (instead of displaying it)

% @=============================================================================
% This function is part of the Brainstorm software:
% https://neuroimage.usc.edu/brainstorm
% 
% Copyright (c) University of Southern California & McGill University
% This software is distributed under the terms of the GNU General Public License
% as published by the Free Software Foundation. Further details on the GPLv3
% license can be found at http://www.gnu.org/copyleft/gpl.html.
% 
% FOR RESEARCH PURPOSES ONLY. THE SOFTWARE IS PROVIDED "AS IS," AND THE
% UNIVERSITY OF SOUTHERN CALIFORNIA AND ITS COLLABORATORS DO NOT MAKE ANY
% WARRANTY, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO WARRANTIES OF
% MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, NOR DO THEY ASSUME ANY
% LIABILITY OR RESPONSIBILITY FOR THE USE OF THIS SOFTWARE.
%
% For more information type "brainstorm license" at command prompt.
% =============================================================================@
%
% Author: Francois Tadel, 2016-2017


% ===== FILES TO IMPORT =====
% Output folder for reports
if (nargin < 2) || isempty(reports_dir) || ~isdir(reports_dir)
    reports_dir = [];
end
% You have to specify the folder in which the tutorial dataset is unzipped
if (nargin == 0) || isempty(tutorial_dir) || ~file_exist(tutorial_dir)
    error('The first argument must be the full path to the dataset folder.');
end
% Subject name
SubjectName = 'Subject01';
% Build the path of the files to import
AnatDir    = fullfile(tutorial_dir, 'sample_introduction', 'anatomy');
Run1File   = fullfile(tutorial_dir, 'sample_introduction', 'data', 'S01_AEF_20131218_01_600Hz.ds');
Run2File   = fullfile(tutorial_dir, 'sample_introduction', 'data', 'S01_AEF_20131218_02_600Hz.ds');
NoiseFile  = fullfile(tutorial_dir, 'sample_introduction', 'data', 'S01_Noise_20131218_02_600Hz.ds');
% Check if the folder contains the required files
if ~file_exist(Run1File)
    error(['The folder ' tutorial_dir ' does not contain the folder from the file sample_introduction.zip.']);
end
% Re-inialize random number generator
if (bst_get('MatlabVersion') >= 712)
    rng('default');
end


%% ===== TUTORIAL #1: CREATE PROTOCOL ================================================
%  ===================================================================================
disp([10 'DEMO> Tutorial #1: Create protocol' 10]);
% The protocol name has to be a valid folder name (no spaces, no weird characters...)
ProtocolName = 'TutorialIntroduction';
% Start brainstorm without the GUI
if ~brainstorm('status')
    brainstorm nogui
end
% Delete existing protocol
gui_brainstorm('DeleteProtocol', ProtocolName);
% Create new protocol
gui_brainstorm('CreateProtocol', ProtocolName, 0, 0);
% Start a new report
bst_report('Start');
% Reset colormaps
bst_colormaps('RestoreDefaults', 'meg');
% Set the current display mode to 'butterfly'
bst_set('TSDisplayMode', 'butterfly');


%% ===== TUTORIAL #2: IMPORT ANATOMY =================================================
%  ===================================================================================
disp([10 'DEMO> Tutorial #2: Import anatomy' 10]);
% Process: Import FreeSurfer folder
bst_process('CallProcess', 'process_import_anatomy', [], [], ...
    'subjectname', SubjectName, ...
    'mrifile',     {AnatDir, 'FreeSurfer'}, ...
    'nvertices',   15000, ...
    'nas', [127, 213, 139], ...
    'lpa', [ 52, 113,  96], ...
    'rpa', [202, 113,  91]);
% This automatically calls the SPM registration procedure because the AC/PC/IH points are not defined



%% ===== TUTORIAL #3: EXPLORE ANATOMY ================================================
%  ===================================================================================
disp([10 'DEMO> Tutorial #3: Explore anatomy' 10]);
% Get subject definition
sSubject = bst_get('Subject', SubjectName);
% Get MRI file and surface files
MriFile    = sSubject.Anatomy(sSubject.iAnatomy).FileName;
CortexFile = sSubject.Surface(sSubject.iCortex).FileName;
HeadFile   = sSubject.Surface(sSubject.iScalp).FileName;
% Display MRI
hFigMri1 = view_mri(MriFile);
hFigMri3 = view_mri_3d(MriFile, [], [], 'NewFigure');
hFigMri2 = view_mri_slices(MriFile, 'x', 20); 
pause(0.5);
% Close figures
close([hFigMri1 hFigMri2 hFigMri3]);
% Display scalp and cortex
hFigSurf = view_surface(HeadFile);
hFigSurf = view_surface(CortexFile, [], [], hFigSurf);
hFigMriSurf = view_mri(MriFile, CortexFile);
% Figure configuration
iTess = 2;
panel_surface('SetShowSulci',     hFigSurf, iTess, 1);
panel_surface('SetSurfaceColor',  hFigSurf, iTess, [1 0 0]);
panel_surface('SetSurfaceSmooth', hFigSurf, iTess, 0.5, 0);
panel_surface('SetSurfaceTransparency', hFigSurf, iTess, 0.8);
figure_3d('SetStandardView', hFigSurf, 'left');
pause(0.5);
% Close figures
close([hFigSurf hFigMriSurf]);



%% ===== TUTORIAL #4: CHANNEL FILE ===================================================
%  ===================================================================================
disp([10 'DEMO> Tutorial #4: Channel file' 10]);
% Process: Create link to raw files
sFilesRun1 = bst_process('CallProcess', 'process_import_data_raw', [], [], ...
    'subjectname',  SubjectName, ...
    'datafile',     {Run1File, 'CTF'}, ...
    'channelalign', 1);
sFilesRun2 = bst_process('CallProcess', 'process_import_data_raw', [], [], ...
    'subjectname',  SubjectName, ...
    'datafile',     {Run2File, 'CTF'}, ...
    'channelalign', 1);
sFilesNoise = bst_process('CallProcess', 'process_import_data_raw', [], [], ...
    'subjectname',  SubjectName, ...
    'datafile',     {NoiseFile, 'CTF'}, ...
    'channelalign', 0);
sFilesRaw = [sFilesRun1, sFilesRun2, sFilesNoise];
% Process: Snapshot: Sensors/MRI registration
bst_process('CallProcess', 'process_snapshot', [sFilesRun1, sFilesRun2], [], ...
    'target',   1, ...  % Sensors/MRI registration
    'modality', 1, ...  % MEG (All)
    'orient',   1, ...  % left
    'Comment',  'MEG/MRI Registration');

% View sensors
hFig = view_surface(HeadFile);
hFig = view_channels(sFilesRun1.ChannelFile, 'MEG', 1, 1, hFig);
% Hide sensors
pause(0.5);
hFig = view_channels(sFilesRun1.ChannelFile, 'MEG', 0, 0, hFig);
% View coils
hFig = view_channels(sFilesRun1.ChannelFile, 'CTF', 1, 1, hFig);
% View helmet
pause(0.5);
hFig = view_helmet(sFilesRun1.ChannelFile, hFig);
pause(0.5);
close(hFig);
% Edit good/bad channel for current file
gui_edit_channel(sFilesRun1.ChannelFile);
pause(0.5);
% Unload everything
bst_memory('UnloadAll', 'Forced');



%% ===== TUTORIAL #5: REVIEW RAW =====================================================
%  ===================================================================================
disp([10 'DEMO> Tutorial #5: Review raw' 10]);
% Process: Convert to continuous (CTF): Continuous
bst_process('CallProcess', 'process_ctf_convert', sFilesRaw, [], ...
    'rectype', 2);  % Continuous

% View recordings
hFigMeg = view_timeseries(sFilesRun1.FileName, 'MEG');
hFigEeg = view_timeseries(sFilesRun1.FileName, 'Misc');
hFigSel = view_timeseries(sFilesRun1.FileName, 'MEG', {'MLT11','MLT12','MLT13'});
% Figure configuration
pause(0.5);
panel_record('SetTimeLength', 3);
panel_record('SetStartTime', 100);
panel_record('SetDisplayMode', hFigMeg, 'column');
panel_montage('SetCurrentMontage', hFigMeg, 'CTF LT');
% Set filters: panel_filter('SetFilters', LowPassEnabled, LowPassValue, HighPassEnabled, HighPassValue, SinRemovalEnabled, SinRemovalValue, MirrorEnabled, FullSourcesEnabled)
panel_filter('SetFilters', 1, 100, 1, 1, 0, [], 0, 0);
pause(0.5);
panel_record('SetDisplayMode', hFigMeg, 'butterfly');
panel_montage('SetCurrentMontage', hFigMeg, '');
% Close figures
close([hFigMeg hFigEeg hFigSel]);



%% ===== TUTORIAL #8: STIM DELAYS ====================================================
%  ===================================================================================
disp([10 'DEMO> Tutorial #8: Stim delays' 10]);
% Process: Detect: standard_fix
bst_process('CallProcess', 'process_evt_detect_analog', [sFilesRun1, sFilesRun2], [], ...
    'eventname',   'standard_fix', ...
    'channelname', 'UADC001', ...
    'timewindow',  [], ...
    'threshold',   1.2, ...
    'blanking',    0.2, ...
    'highpass',    0, ...
    'lowpass',     0, ...
    'refevent',    'standard', ...
    'isfalling',   0, ...
    'ispullup',    0, ...
    'isclassify',  0);
% Process: Detect: deviant_fix
bst_process('CallProcess', 'process_evt_detect_analog', [sFilesRun1, sFilesRun2], [], ...
    'eventname',   'deviant_fix', ...
    'channelname', 'UADC001', ...
    'timewindow',  [], ...
    'threshold',   1.2, ...
    'blanking',    0.2, ...
    'highpass',    0, ...
    'lowpass',     0, ...
    'refevent',    'deviant', ...
    'isfalling',   0, ...
    'ispullup',    0, ...
    'isclassify',  0);
% Process: Read from channel
bst_process('CallProcess', 'process_evt_read', [sFilesRun1, sFilesRun2], [], ...
    'stimchan',  'UDIO001', ...
    'trackmode', 1, ...  % Value: detect the changes of channel value
    'zero',      0);

% Process: Delete events
bst_process('CallProcess', 'process_evt_delete', [sFilesRun1, sFilesRun2], [], ...
    'eventname', 'standard, deviant, button');
% Process: Rename event (standard_fix>standard)
bst_process('CallProcess', 'process_evt_rename', [sFilesRun1, sFilesRun2], [], ...
    'src',  'standard_fix', ...
    'dest', 'standard');
% Process: Rename event (deviant_fix>deviant)
bst_process('CallProcess', 'process_evt_rename', [sFilesRun1, sFilesRun2], [], ...
    'src',  'deviant_fix', ...
    'dest', 'deviant');
% Process: Rename event (64>button)
bst_process('CallProcess', 'process_evt_rename', [sFilesRun1, sFilesRun2], [], ...
    'src',  '64', ...
    'dest', 'button');



%% ===== TUTORIAL #10: FREQUENCY FILTERS =============================================
%  ===================================================================================
disp([10 'DEMO> Tutorial #10: Frequency filters' 10]);
% Process: Sinusoid removal: 60Hz 120Hz 180Hz 300Hz
sFilesNotch = bst_process('CallProcess', 'process_notch', sFilesRaw, [], ...
    'freqlist',    [60, 120, 180], ...
    'sensortypes', 'MEG', ...
    'read_all',    0);
% Process: Power spectrum density (Welch)
sFilesPsd = bst_process('CallProcess', 'process_psd', [sFilesRaw, sFilesNotch], [], ...
    'timewindow',  [], ...
    'win_length',  4, ...
    'win_overlap', 50, ...
    'clusters',    {}, ...
    'sensortypes', 'MEG', ...
    'edit', struct(...
         'Comment',    'Power', ...
         'TimeBands',  [], ...
         'Freqs',      [], ...
         'ClusterFuncTime', 'none', ...
         'Measure',    'power', ...
         'Output',     'all', ...
         'SaveKernel', 0));
% Process: Snapshot: Frequency spectrum
bst_process('CallProcess', 'process_snapshot', sFilesPsd, [], ...
    'target',   10, ...  % Frequency spectrum
    'modality', 1, ...   % MEG (All)
    'Comment',  'Power spectrum density');
% Process: Delete folders
bst_process('CallProcess', 'process_delete', sFilesRaw, [], ...
    'target', 2);  % Delete folders
% Separate the three outputs
sFilesRun1  = {sFilesNotch(1).FileName};
sFilesRun2  = {sFilesNotch(2).FileName};
sFilesNoise = {sFilesNotch(3).FileName};



%% ===== TUTORIAL #11: BAD CHANNELS ==================================================
%  ===================================================================================
% % Process: Set bad channels
% sFiles = bst_process('CallProcess', 'process_channel_setbad', sFilesRun2, [], ...
%     'sensortypes', 'MRT51, MLO52, MLO42, MLO43');



%% ===== TUTORIAL #12: ARTIFACTS DETECTION ===========================================
%  ===================================================================================
disp([10 'DEMO> Tutorial #12: Artifacts detection' 10]);
% Process: Detect heartbeats
bst_process('CallProcess', 'process_evt_detect_ecg', [sFilesRun1, sFilesRun2], [], ...
    'channelname', 'ECG', ...
    'timewindow',  [], ...
    'eventname',   'cardiac');
% Process: Detect eye blinks
bst_process('CallProcess', 'process_evt_detect_eog',  [sFilesRun1, sFilesRun2], [], ...
    'channelname', 'VEOG', ...
    'timewindow',  [], ...
    'eventname',   'blink');
% Process: Remove simultaneous
bst_process('CallProcess', 'process_evt_remove_simult', [sFilesRun1, sFilesRun2], [], ...
    'remove', 'cardiac', ...
    'target', 'blink', ...
    'dt',     0.25, ...
    'rename', 0);



%% ===== TUTORIAL #13: SSP ===========================================================
%  ===================================================================================
disp([10 'DEMO> Tutorial #13: SSP' 10]);
% Process: SSP ECG: cardiac
bst_process('CallProcess', 'process_ssp_ecg', [sFilesRun1, sFilesRun2], [], ...
    'eventname',   'cardiac', ...
    'sensortypes', 'MEG', ...
    'usessp',      0, ...
    'select',      1);
% Process: SSP EOG: blink
bst_process('CallProcess', 'process_ssp_eog', sFilesRun1, [], ...
    'eventname',   'blink', ...
    'sensortypes', 'MEG', ...
    'usessp',      0, ...
    'select',      [1 2]);
bst_process('CallProcess', 'process_ssp_eog', sFilesRun2, [], ...
    'eventname',   'blink', ...
    'sensortypes', 'MEG', ...
    'usessp',      0, ...
    'select',      1);


%% ===== TUTORIAL #14: BAD SEGMENTS ==================================================
%  ===================================================================================
disp([10 'DEMO> Tutorial #14: Bad segments' 10]);
% Process: Detect other artifacts
bst_process('CallProcess', 'process_evt_detect_badsegment', [sFilesRun1, sFilesRun2], [], ...
    'timewindow',  [], ...
    'sensortypes', 'MEG', ...
    'threshold',   3, ...  % 3
    'isLowFreq',   1, ...
    'isHighFreq',  1);

% Process: Rename event (1-7Hz > saccade)  (Run02 only)
bst_process('CallProcess', 'process_evt_rename', sFilesRun2, [], ...
    'src',  '1-7Hz', ...
    'dest', 'saccade');

% Manual selection of saccades (cannot be done from the pipeline editor: manual edition of the structures)
sMatRun2 = in_bst_data(sFilesRun2{1}, 'F');
iEvtSaccade = find(strcmpi({sMatRun2.F.events.label}, 'saccade'));
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];
sMatRun2.F.events(iEvtSaccade).epochs  = ones(1, size(sMatRun2.F.events(iEvtSaccade).times, 2));
sMatRun2.F.events(iEvtSaccade).channels = [];
sMatRun2.F.events(iEvtSaccade).notes    = [];
bst_save(file_fullpath(sFilesRun2{1}), sMatRun2, 'v6', 1);

% Process: SSP: saccade  (Run02 only)
bst_process('CallProcess', 'process_ssp', sFilesRun2, [], ...
    'timewindow',  [], ...
    'eventname',   'saccade', ...
    'eventtime',   [-0.2, 0.2], ...
    'bandpass',    [1.5, 7], ...
    'sensortypes', 'MEG', ...
    'usessp',      1, ...
    'saveerp',     0, ...
    'method',      1, ...
    'select',      1);
% Process: Detect other artifacts  (Run02 only)
bst_process('CallProcess', 'process_evt_detect_badsegment', sFilesRun2, [], ...
    'timewindow',  [], ...
    'sensortypes', 'MEG', ...
    'threshold',   3, ...  % 3
    'isLowFreq',   1, ...
    'isHighFreq',  1);

% Process: Rename event (1-7Hz > bad_1-7Hz)
bst_process('CallProcess', 'process_evt_rename', [sFilesRun1, sFilesRun2], [], ...
    'src',  '1-7Hz', ...
    'dest', 'bad_1-7Hz');
% Process: Rename event (40-240Hz > bad_40-240Hz)
bst_process('CallProcess', 'process_evt_rename', [sFilesRun1, sFilesRun2], [], ...
    'src',  '40-240Hz', ...
    'dest', 'bad_40-240Hz');
% Process: Snapshot: SSP projectors
bst_process('CallProcess', 'process_snapshot', [sFilesRun1, sFilesRun2], [], ...
    'target',  2, ...  % SSP projectors
    'Comment', 'SSP projectors');



%% ===== TUTORIAL #15: IMPORT EVENTS =================================================
%  ===================================================================================
disp([10 'DEMO> Tutorial #15: Import events' 10]);
% Process: Import MEG/EEG: Events (Run01)
sFilesEpochs1 = bst_process('CallProcess', 'process_import_data_event', sFilesRun1, [], ...
    'subjectname', SubjectName, ...
    'condition',   '', ...
    'eventname',   'standard, deviant', ...
    'timewindow',  [], ...
    'epochtime',   [-0.100, 0.500], ...
    'createcond',  0, ...
    'ignoreshort', 1, ...
    'usectfcomp',  1, ...
    'usessp',      1, ...
    'freq',        [], ...
    'baseline',    [-0.1, -0.0017]);
% Process: Import MEG/EEG: Events (Run02)
sFilesEpochs2 = bst_process('CallProcess', 'process_import_data_event', sFilesRun2, [], ...
    'subjectname', SubjectName, ...
    'condition',   '', ...
    'eventname',   'standard, deviant', ...
    'timewindow',  [], ...
    'epochtime',   [-0.100, 0.500], ...
    'createcond',  0, ...
    'ignoreshort', 1, ...
    'usectfcomp',  1, ...
    'usessp',      1, ...
    'freq',        [], ...
    'baseline',    [-0.1, -0.0017]);
% Display raster plot
hFigRaster = view_erpimage({sFilesEpochs1.FileName}, 'erpimage', 'MEG');
panel_display();
bst_report('Snapshot', hFigRaster, sFilesEpochs1(1).FileName, 'ERP image');
close(hFigRaster);


%% ===== TUTORIAL #16: AVERAGE =======================================================
%  ===================================================================================
disp([10 'DEMO> Tutorial #16: Average' 10]);
% Process: Average: By trial group (folder average)
sFilesAvg = bst_process('CallProcess', 'process_average', [sFilesEpochs1, sFilesEpochs2], [], ...
    'avgtype',    5, ...  % By trial groups (folder average)
    'avg_func',   1, ...  % Arithmetic average: mean(x)
    'weighted',   0, ...
    'keepevents', 1);
% Process: Delete events 'cardiac'
bst_process('CallProcess', 'process_evt_delete', sFilesAvg, [], ...
    'eventname', 'cardiac');
% Process: Snapshot: Recordings time series
bst_process('CallProcess', 'process_snapshot', sFilesAvg, [], ...
    'target',     5, ...  % Recordings time series
    'modality',   1, ...  % MEG (All)
    'Comment',    'Evoked response');
% Set colormap: global color scale
bst_colormaps('SetMaxMode', 'meg', 'global');
% Process: Snapshot: Recordings topography (contact sheet)
bst_process('CallProcess', 'process_snapshot', sFilesAvg, [], ...
    'target',         7, ...  % Recordings topography (contact sheet)
    'modality',       1, ...  % MEG
    'contact_time',   [0, 0.350], ...
    'contact_nimage', 15, ...
    'Comment',        'Evoked response');

% Process: Average+Stderr: By trial group (subject average)
sFilesAvgAll = bst_process('CallProcess', 'process_average', [sFilesEpochs1, sFilesEpochs2], [], ...
    'avgtype',    6, ...  % By trial group (subject average)
    'avg_func',   7, ...  % Arithmetic average + Standard error
    'weighted',   0, ...
    'keepevents', 1);
% Process: Delete events 'cardiac'
bst_process('CallProcess', 'process_evt_delete', sFilesAvgAll, [], ...
    'eventname', 'cardiac');
% Process: Delete events 'saccade'
bst_process('CallProcess', 'process_evt_delete', sFilesAvgAll, [], ...
    'eventname', 'saccade');
% Process: Snapshot: Recordings time series
bst_process('CallProcess', 'process_snapshot', sFilesAvgAll, [], ...
    'target',     5, ...  % Recordings time series
    'modality',   1, ...  % MEG (All)
    'Comment',    'Evoked response');


%% ===== TUTORIAL #17: EXPLORATION ===================================================
%  ===================================================================================
disp([10 'DEMO> Tutorial #17: Bad segments' 10]);
% View averages
hFigMeg1  = view_timeseries(sFilesAvg(1).FileName, 'MEG');
hFigMeg2  = view_timeseries(sFilesAvg(2).FileName, 'MEG');
hFigEeg1  = view_timeseries(sFilesAvg(1).FileName, 'Misc');
hFigEeg2  = view_timeseries(sFilesAvg(2).FileName, 'Misc');
hFigTopo1 = view_topography(sFilesAvg(1).FileName, 'MEG', '2DSensorCap');
hFigTopo2 = view_topography(sFilesAvg(2).FileName, 'MEG', '2DSensorCap');
hFigTp2   = view_topography(sFilesAvg(3).FileName, 'MEG', '3DSensorCap');
hFigTp3   = view_topography(sFilesAvg(3).FileName, 'MEG', '2DDisc');
hFigTp4   = view_topography(sFilesAvg(3).FileName, 'MEG', '2DLayout');
% Set time: 90ms
panel_time('SetCurrentTime', 0.090);
% Set filters: 40Hz low-pass, no high-pass
panel_filter('SetFilters', 1, 40, 0, [], 0, [], 0, 0);
% View selected sensors
SelectedChannels = {'MLC31', 'MLC32'};
bst_figures('SetSelectedRows', SelectedChannels);
view_timeseries(sFilesAvg(4).FileName, [], SelectedChannels);
% Select time window
figure_timeseries('SetTimeSelectionManual', hFigMeg1, [0.070, 0.130]);
% Show sensors on 2DSensorCap topography
isMarkers = 1;
isLabels = 0;
figure_3d('ViewSensors', hFigTopo1, isMarkers, isLabels);
% Display time contact sheet for a figure
pause(0.5);
hContactFig = view_contactsheet( hFigTopo2, 'time', 'fig', [], 12, [0 0.120] );
pause(0.5);
close(hContactFig);



%% ===== TUTORIAL #18: COLORMAPS =====================================================
%  ===================================================================================
disp([10 'DEMO> Tutorial #18: Colormaps' 10]);
% Set 'Meg' colormap to 'jet'
bst_colormaps('SetColormapName', 'meg', 'jet');
pause(0.5);
% Set 'Meg' colormap to 'rwb'
bst_colormaps('SetColormapName', 'meg', 'cmap_rbw');
% Set colormap to display absolute values
bst_colormaps('SetColormapAbsolute', 'meg', 1);
% Normalize colormap for each time frame
bst_colormaps('SetMaxMode', 'meg', 'local');
% Hide colorbar
bst_colormaps('SetDisplayColorbar', 'meg', 0);
pause(0.5);
% Restore colormap to default values
bst_colormaps('RestoreDefaults', 'meg');
% Edit good/bad channel for current file
gui_edit_channelflag(sFilesAvg(1).FileName);
% Close figures
pause(0.5);
bst_memory('UnloadAll', 'Forced');



%% ===== TUTORIAL #20: HEAD MODEL ====================================================
%  ===================================================================================
disp([10 'DEMO> Tutorial #20: Head model' 10]);
% Process: Compute head model
bst_process('CallProcess', 'process_headmodel', sFilesAvg, [], ...
    'comment',      '', ...
    'sourcespace',  1, ...
    'meg',          3);  % Overlapping spheres
% Get study structure
sStudy = bst_get('Study', sFilesAvg(1).iStudy);
% Show spheres
hFig = view_spheres(sStudy.HeadModel(sStudy.iHeadModel).FileName, sStudy.Channel.FileName, sSubject);
pause(0.5);
close(hFig);



%% ===== TUTORIAL #21: NOISE COVARIANCE ==============================================
%  ===================================================================================
disp([10 'DEMO> Tutorial #21: Noise covariance' 10]);
% Process: Compute covariance (noise or data)
bst_process('CallProcess', 'process_noisecov', sFilesNoise, [], ...
    'baseline',       [], ...
    'sensortypes',    'MEG, EEG, SEEG, ECOG', ...
    'target',         1, ...  % Noise covariance     (covariance over baseline time window)
    'dcoffset',       1, ...  % Block by block, to avoid effects of slow shifts in data
    'identity',       0, ...
    'copycond',       1, ...
    'copysubj',       0, ...
    'replacefile',    1);  % Replace
% Process: Snapshot: Noise covariance
bst_process('CallProcess', 'process_snapshot', sFilesNoise, [], ...
    'target',  3, ...  % Noise covariance
    'Comment', 'Noise covariance');

% Process: Compute covariance (noise or data)  [Run01]
bst_process('CallProcess', 'process_noisecov', sFilesEpochs1, [], ...
    'baseline',       [-0.1, -0.0017], ...
    'datatimewindow', [0, 0.5], ...
    'sensortypes',    'MEG, EEG, SEEG, ECOG', ...
    'target',         2, ...  % Data covariance      (covariance over data time window)
    'dcoffset',       1, ...  % Block by block, to avoid effects of slow shifts in data
    'identity',       0, ...
    'copycond',       0, ...
    'copysubj',       0, ...
    'replacefile',    1);  % Replace
% Process: Compute covariance (noise or data)  [Run02]
bst_process('CallProcess', 'process_noisecov', sFilesEpochs2, [], ...
    'baseline',       [-0.1, -0.0017], ...
    'datatimewindow', [0, 0.5], ...
    'sensortypes',    'MEG, EEG, SEEG, ECOG', ...
    'target',         2, ...  % Data covariance      (covariance over data time window)
    'dcoffset',       1, ...  % Block by block, to avoid effects of slow shifts in data
    'identity',       0, ...
    'copycond',       0, ...
    'copysubj',       0, ...
    'replacefile',    1);  % Replace
% Process: Snapshot: Data covariance
bst_process('CallProcess', 'process_snapshot', [sFilesEpochs1(1), sFilesEpochs2(1)], [], ...
    'target',  12, ...  % Data covariance
    'Comment', 'Data covariance');



%% ===== TUTORIAL #22: SOURCE ESTIMATION =============================================
%  ===================================================================================
disp([10 'DEMO> Tutorial #22: Source estimation' 10]);
% === GET DEVIANT AVERAGE RUN01 ===
% Process: Select recordings in: Subject01/S01_AEF_20131218_01_600Hz_notch
sFiles01 = bst_process('CallProcess', 'process_select_files_data', [], [], ...
    'subjectname', SubjectName, ...
    'condition',   'S01_AEF_20131218_01_600Hz_notch', ...
    'includebad',  0);
% Process: Select file comments with tag: deviant
sFilesAvgDeviant01 = bst_process('CallProcess', 'process_select_tag', sFiles01, [], ...
    'tag',    'Avg: deviant', ...
    'search', 2, ...  % Search the file comments
    'select', 1);     % Select only the files with the tag

% === CONSTRAINED EXAMPLE ===
% Minimum norm options
InverseOptions = struct(...
    'Comment',        'MN: MEG', ...
    'InverseMethod',  'minnorm', ...
    'InverseMeasure', 'amplitude', ...
    'SourceOrient',   {{'fixed'}}, ...
    'Loose',          0.2, ...
    'UseDepth',       1, ...
    'WeightExp',      0.5, ...
    'WeightLimit',    10, ...
    'NoiseMethod',    'reg', ...
    'NoiseReg',       0.1, ...
    'SnrMethod',      'fixed', ...
    'SnrRms',         0.001, ...
    'SnrFixed',       3, ...
    'ComputeKernel',  1, ...
    'DataTypes',      {{'MEG'}});
% Process: Compute sources [2018]
sFilesSrcDeviant01 = bst_process('CallProcess', 'process_inverse_2018', sFilesAvgDeviant01, [], ...
    'output',  2, ...  % Kernel only: one per file
    'inverse', InverseOptions);

% === DISPLAY SOURCES MANUALLY ===
% View time series
hFigSrc1 = view_timeseries(sFilesAvgDeviant01(1).FileName, 'MEG');
% View on the cortex surface
hFigSrc2 = script_view_sources(sFilesSrcDeviant01.FileName, 'cortex');
% Set current time to 90ms
panel_time('SetCurrentTime', 0.090);
% Set orientation
figure_3d('SetStandardView', hFigSrc2, 'left');
% Set surface threshold
iSurf = 1;
thresh = .30; 
panel_surface('SetDataThreshold', hFigSrc2, iSurf, thresh);
% Set surface smoothing
panel_surface('SetSurfaceSmooth', hFigSrc2, iSurf, .6, 0);
% Show sulci
panel_surface('SetShowSulci', hFigSrc2, iSurf, 1);

% View sources on MRI (3D orthogonal slices)
hFigSrc3 = script_view_sources(sFilesSrcDeviant01.FileName, 'mri3d');
panel_surface('SetDataThreshold', hFigSrc3, iSurf, thresh);
% Set the position of the cuts in the 3D figure
cutsPosVox = [74 93 159];
panel_surface('PlotMri', hFigSrc3, cutsPosVox);

% View sources with MRI Viewer
hFigSrc4 = script_view_sources(sFilesSrcDeviant01.FileName, 'mriviewer');
panel_surface('SetDataThreshold', hFigSrc4, iSurf, thresh);
% Set the position of the cuts in the MRI Viewer (values in millimeters)
figure_mri('SetLocation', 'voxel', hFigSrc4, [], cutsPosVox);
% Close figures
close([hFigSrc1 hFigSrc2 hFigSrc3 hFigSrc4]);

% === UNCONSTRAINED EXAMPLE ===
% Unconstrained minnorm
InverseOptions.Comment        = 'MN: MEG';
InverseOptions.InverseMeasure = 'amplitude';
InverseOptions.SourceOrient   = {'free'};
% Process: Compute sources [2018]
sFilesSrcUnconst = bst_process('CallProcess', 'process_inverse_2018', sFilesAvgDeviant01, [], ...
    'output',  2, ...  % Kernel only: one per file
    'inverse', InverseOptions);


% === NORMALIZED SOURCES ===
% dSPM
InverseOptions.Comment        = 'dSPM: MEG';
InverseOptions.InverseMeasure = 'dspm2018';
InverseOptions.SourceOrient   = {'fixed'};
sFilesSrcDspm = bst_process('CallProcess', 'process_inverse_2018', sFilesAvgDeviant01, [], ...
    'output',  2, ...  % Kernel only: one per file
    'inverse', InverseOptions);
% sLORETA (old function)
sFilesSrcSloreta = bst_process('CallProcess', 'process_inverse', sFilesAvgDeviant01, [], ...
    'comment', '', ...
    'method',  3, ...  % sLORETA
    'wmne', struct(...
         'SourceOrient', {{'fixed'}}, ...
         'loose',        0.2, ...
         'SNR',          3, ...
         'pca',          1, ...
         'diagnoise',    0, ...
         'regnoise',     1, ...
         'magreg',       0.1, ...
         'gradreg',      0.1, ...
         'depth',        1, ...
         'weightexp',    0.5, ...
         'weightlimit',  10), ...
    'sensortypes', 'MEG, MEG MAG, MEG GRAD, EEG', ...
    'output', 2);  % Kernel only: one per file
% Process: Z-score normalization: [-100ms,-2ms]
sFilesSrcZscore = bst_process('CallProcess', 'process_baseline_norm', sFilesSrcDeviant01, [], ...
    'baseline',   [-0.100, -0.002], ...
    'source_abs', 0, ...
    'method',     'zscore');   % Z-score transformation:    x_std = (x - μ) / σ

% Process: Snapshot: Sources (one time)
bst_process('CallProcess', 'process_snapshot', sFilesSrcDeviant01, [], ...
    'target',    8, ...  % Sources (one time)
    'orient',    1, ...  % left
    'time',      0.09, ...
    'threshold', 30, ...
    'Comment',   'Current density map (Constrained)');
bst_process('CallProcess', 'process_snapshot', sFilesSrcDspm, [], ...
    'target',    8, ...  % Sources (one time)
    'orient',    1, ...  % left
    'time',      0.09, ...
    'threshold', 60, ...
    'Comment',   'dSPM');
bst_process('CallProcess', 'process_snapshot', sFilesSrcSloreta, [], ...
    'target',    8, ...  % Sources (one time)
    'orient',    1, ...  % left
    'time',      0.09, ...
    'threshold', 60, ...
    'Comment',   'sLORETA');
bst_process('CallProcess', 'process_snapshot', sFilesSrcZscore, [], ...
    'target',    8, ...  % Sources (one time)
    'orient',    1, ...  % left
    'time',      0.09, ...
    'threshold', 60, ...
    'Comment',   'Z-score');
bst_process('CallProcess', 'process_snapshot', sFilesSrcUnconst, [], ...
    'target',         8, ...  % Sources (one time)
    'orient',         1, ...  % left
    'time',           0.0917, ...
    'threshold',      0, ...
    'Comment',        'Current density map (Unconstrained)');

% === DELETE EXPERIMENTS ===
% Process: Delete constrained example
bst_process('CallProcess', 'process_delete', [sFilesSrcDeviant01, sFilesSrcDspm, sFilesSrcSloreta, sFilesSrcZscore, sFilesSrcUnconst], [], ...
    'target', 1);  % Delete selected files


% === SHARED KERNEL ===
% Constrained minnorm
InverseOptions.Comment        = 'MN: MEG';
InverseOptions.InverseMeasure = 'amplitude';
InverseOptions.SourceOrient   = {'fixed'};
% Process: Compute sources [2018]
sFilesAvgSrc = bst_process('CallProcess', 'process_inverse_2018', sFilesAvg, [], ...
    'output',  1, ...  % Kernel only: shared
    'inverse', InverseOptions);


% === AVERAGE SOURCES ACROSS RUNS ===
% Process: Average: By trial group (subject average)
sFilesIntraSrc = bst_process('CallProcess', 'process_average', sFilesAvgSrc, [], ...
    'avgtype',         6, ...  % By trial group (subject average)
    'avg_func',        1, ...  % Arithmetic average:  mean(x)
    'weighted',        1, ...
    'scalenormalized', 0);
% Process: Low-pass:40Hz
sFilesIntraSrcLow = bst_process('CallProcess', 'process_bandpass', sFilesIntraSrc, [], ...
    'highpass',    0, ...
    'lowpass',     40, ...
    'attenuation', 'strict', ...  % 60dB
    'mirror',      0, ...
    'overwrite',   0);
% Process: Z-score normalization: [-100ms,-2ms]
sFilesIntraZscore = bst_process('CallProcess', 'process_baseline_norm', sFilesIntraSrcLow, [], ...
    'baseline',   [-0.100, -0.0017], ...
    'source_abs', 0, ...
    'method',     'zscore');  % Z-score transformation:    x_std = (x - μ) / σ

% Process: Delete intermediate results
bst_process('CallProcess', 'process_delete', sFilesIntraSrcLow, [], ...
    'target', 1);  % Delete selected files
% Screen captures
bst_process('CallProcess', 'process_snapshot', sFilesIntraZscore, [], ...
    'target',    8, ...  % Sources (one time)
    'orient',    1, ...  % left
    'time',      0.09, ...
    'threshold', 40, ...
    'Comment',   'Average across runs (left)');
bst_process('CallProcess', 'process_snapshot', sFilesIntraZscore, [], ...
    'target',    8, ...  % Sources (one time)
    'orient',    2, ...  % right
    'time',      0.09, ...
    'threshold', 40, ...
    'Comment',   'Average across runs (right)');
bst_process('CallProcess', 'process_snapshot', sFilesIntraZscore(1), [], ...
    'target',         9, ...  % Sources (contact sheet)
    'orient',         1, ...  % left
    'contact_time',   [0, 0.35], ...
    'contact_nimage', 15, ...
    'threshold',      20, ...
    'Comment',        'Average across runs');


%% ===== TUTORIAL #23: SCOUTS ========================================================
%  ===================================================================================
disp([10 'DEMO> Tutorial #23: Scouts' 10]);
% Load surface file
sCortex = in_tess_bst(CortexFile);
% Add new scouts in first atlas
sCortex.iAtlas = find(strcmpi({sCortex.Atlas.Name}, 'Destrieux'));
% Save file
bst_save(file_fullpath(CortexFile), sCortex, 'v7');
% Unload everything
bst_memory('UnloadAll', 'Forced');
% Find scouts indices to display: {'A1L', 'A1R', 'IFGL', 'M1L'}
[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'});

% View cortex
hFigSurf1 = view_surface(CortexFile, [], [], 'NewFigure');
hFigSurf2 = view_surface(CortexFile, [], [], 'NewFigure');
figure_3d('SetStandardView', hFigSurf1, 'left');
figure_3d('SetStandardView', hFigSurf2, 'right');
panel_surface('SetSurfaceSmooth', hFigSurf1, 1, .6, 0);
panel_surface('SetSurfaceSmooth', hFigSurf2, 1, .6, 0);
panel_surface('SetShowSulci', hFigSurf1, 1, 1);
panel_surface('SetShowSulci', hFigSurf2, 1, 1);
% Configure scouts display
panel_scout('SetScoutsOptions', 0, 1, 1, 'all', 0.7, 1, 1, 0);
% View scouts
hFigScouts = view_scouts({sFilesIntraZscore.FileName}, iScouts);
hLegend = findobj(hFigScouts, 'Type', 'legend');
if ~isempty(hLegend) && ishandle(hLegend(1))
    set(hLegend(1), 'Units', 'pixels');
    pos = get(hLegend(1), 'Position');
    set(hLegend(1), 'Position', [1, 1, pos(3), pos(4)]);
end
% Save figures
bst_report('Snapshot', hFigScouts, sFilesIntraZscore(1).FileName, 'Scouts', [100 100 670 250]);
% Close all
pause(1);
bst_memory('UnloadAll', 'Forced');



%% ===== TUTORIAL #24: TIME-FREQUENCY ================================================
%  ===================================================================================
disp([10 'DEMO> Tutorial #24: Time-frequency' 10]);
% Process: Simulate generic signals
sSim = bst_process('CallProcess', 'process_simulate_matrix', [], [], ...
    'subjectname', 'Test', ...
    'condition',   'Simulation', ...
    'samples',     6000, ...
    'srate',       1000, ...
    '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);']);

% Time-frequency options
TfOptions = struct(...
    'Comment',         'Power,1-60Hz', ...
    'TimeBands',       [], ...
    '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], ...
    'MorletFc',        1, ...
    'MorletFwhmTc',    3, ...
    'ClusterFuncTime', 'none', ...
    'Measure',         'power', ...
    'Output',          'average', ...
    'SaveKernel',       0);
% Process: Time-frequency (Morlet wavelets)
sSimTf1 = bst_process('CallProcess', 'process_timefreq', sSim, [], ...
    'edit',      TfOptions, ...
    'normalize', 'multiply');  % 1/f compensation: Multiply output values by frequency

% === NORMALIZATION ===
% Process: Time-frequency (Morlet wavelets)
sSimTf2 = bst_process('CallProcess', 'process_timefreq', sSim, [], ...
    'edit',      TfOptions, ...
    'normalize', 'none');  % None: Save non-standardized time-frequency maps
% Process: Spectral flattening
sSimTf2Multi = bst_process('CallProcess', 'process_tf_norm', sSimTf2, [], ...
    'normalize', 'multiply', ...  % 1/f compensation (multiple by frequency)
    'overwrite', 0);
% Process: Event-related perturbation (ERS/ERD): [750ms,1250ms]
sSimTf2Ersd = bst_process('CallProcess', 'process_baseline_norm', sSimTf2, [], ...
    'baseline',  [0.75, 1.25], ...
    'method',    'ersd', ...  % Event-related perturbation (ERS/ERD):    x_std = (x - μ) / μ * 100
    'overwrite', 0);

% Process: Snapshot: Time-frequency maps
bst_process('CallProcess', 'process_snapshot', sSimTf2, [], ...
    'target',   14, ...  % Time-frequency maps
    'Comment',  'Not normalized');
% Process: Snapshot: Time-frequency maps
bst_process('CallProcess', 'process_snapshot', sSimTf2Multi, [], ...
    'target',   14, ...  % Time-frequency maps
    'Comment',  'Spectral flattening: 1/f compensation');
% Process: Snapshot: Time-frequency maps
bst_process('CallProcess', 'process_snapshot', sSimTf2Ersd, [], ...
    'target',   14, ...  % Time-frequency maps
    'Comment',  'ERS/ERD');
% Spectrum/Time series
hFigTf1 = view_spectrum(sSimTf2.FileName, 'Spectrum');
hFigTf2 = view_spectrum(sSimTf2.FileName, 'TimeSeries');
panel_time('SetCurrentTime', 0.5);
panel_freq('SetCurrentFreq', 2, 0);
bst_report('Snapshot', [hFigTf1 hFigTf2], sSimTf2.FileName, 'Not normalized: 2s/20Hz', [200, 200, 400, 250]);
panel_time('SetCurrentTime', 2.02);
panel_freq('SetCurrentFreq', 20, 0);
bst_report('Snapshot', [hFigTf1 hFigTf2], sSimTf2.FileName, 'Not normalized: 2s/20Hz', [200, 200, 400, 250]);
bst_memory('UnloadAll', 'Forced');


% === HILBERT TRANSFORM ===
% Process: Hilbert transform
sSimHilbert = bst_process('CallProcess', 'process_hilbert', sSim, [], ...
    'edit', struct(...
         'Comment',         'Power', ...
         'TimeBands',       [], ...
         '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'}}, ...
         'ClusterFuncTime', 'none', ...
         'Measure',         'power', ...
         'Output',          'all', ...
         'SaveKernel',      0), ...
    'normalize', 'none', ...  % None: Save non-standardized time-frequency maps
    'mirror',    0);
% Process: Spectral flattening
sSimHilbertMulti = bst_process('CallProcess', 'process_tf_norm', sSimHilbert, [], ...
    'normalize', 'multiply', ...  % 1/f compensation (multiple by frequency)
    'overwrite', 0);
% Process: Event-related perturbation (ERS/ERD): [750ms,1250ms]
sSimHilbertErsd = bst_process('CallProcess', 'process_baseline_norm', sSimHilbert, [], ...
    'baseline',  [0.75, 1.25], ...
    'method',    'ersd', ...  % Event-related perturbation (ERS/ERD):    x_std = (x - μ) / μ * 100
    'overwrite', 0);

% Process: Snapshot: Time-frequency maps
bst_process('CallProcess', 'process_snapshot', sSimHilbert, [], ...
    'target',   14, ...  % Time-frequency maps
    'Comment',  'Not normalized');
% Process: Snapshot: Time-frequency maps
bst_process('CallProcess', 'process_snapshot', sSimHilbertMulti, [], ...
    'target',   14, ...  % Time-frequency maps
    'Comment',  'Spectral flattening: 1/f compensation');
% Process: Snapshot: Time-frequency maps
bst_process('CallProcess', 'process_snapshot', sSimHilbertErsd, [], ...
    'target',   14, ...  % Time-frequency maps
    'Comment',  'ERS/ERD');

% === SINGLE TRIALS ===
TfOptions.Comment = 'Avg,Power,1-150Hz';
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];
% Process: Time-frequency (Morlet wavelets)
sEpochsAvgTf = bst_process('CallProcess', 'process_timefreq', sFilesEpochs1, [], ...
    'sensortypes', 'MEG, EEG', ...
    'edit',        TfOptions, ...
    'normalize',   'none');  % None: Save non-standardized time-frequency maps
% Process: Event-related perturbation (ERS/ERD): [-75ms,0ms]
sEpochsAvgTfErsd = bst_process('CallProcess', 'process_baseline_norm', sEpochsAvgTf, [], ...
    'baseline',  [-0.075, 0], ...
    'method',    'ersd', ...  % Event-related perturbation (ERS/ERD):    x_std = (x - μ) / μ * 100
    'overwrite', 0);

% === DISPLAY ===
% View time-frequency file
hFigTf1 = view_timefreq(sEpochsAvgTfErsd.FileName, 'SingleSensor');
% Configure display
sOptions = panel_display('GetDisplayOptions');
sOptions.HideEdgeEffects = 1;
sOptions.HighResolution = 1;
panel_display('SetDisplayOptions', sOptions);
% Other display modes
hFigTf2 = view_timefreq(sEpochsAvgTfErsd.FileName, 'AllSensors');
hFigTf3 = view_timefreq(sEpochsAvgTfErsd.FileName, '2DLayout');
hFigTf4 = view_timefreq(sEpochsAvgTfErsd.FileName, '2DLayoutOpt');
bst_colormaps('SetColormapName', 'stat2', 'jet');
bst_colormaps('SetColormapAbsolute', 'stat2', 1);
bst_report('Snapshot', hFigTf1, sEpochsAvgTfErsd.FileName, 'Time-frequency', [200, 200, 400, 250]);
bst_report('Snapshot', [hFigTf2 hFigTf3 hFigTf4], sEpochsAvgTfErsd.FileName, 'Time-frequency', [200, 200, 750, 400]);
close([hFigTf1 hFigTf2 hFigTf3 hFigTf4]);
% Image [channel x time]
hFigTf5 = view_erpimage(sEpochsAvgTfErsd.FileName, 'trialimage');
% Topographies
hFigTf6 = view_topography(sEpochsAvgTfErsd.FileName, 'MEG', '3DSensorCap');
hFigTf7 = view_topography(sEpochsAvgTfErsd.FileName, 'MEG', '2DDisc');
hFigTf8 = view_topography(sEpochsAvgTfErsd.FileName, 'MEG', '2DSensorCap');
hFigTf9 = view_topography(sEpochsAvgTfErsd.FileName, 'MEG', '2DLayout');
panel_time('SetCurrentTime', 0.175);
panel_freq('SetCurrentFreq', 8, 0);
bst_report('Snapshot', [hFigTf5 hFigTf6 hFigTf7 hFigTf8 hFigTf9], sEpochsAvgTfErsd.FileName, 'Time-frequency: 8Hz', [200, 200, 400, 250]);
close([hFigTf5 hFigTf6 hFigTf7 hFigTf8 hFigTf9]);


% === AVERAGE FOR SCOUTS ===
% Select all sources for the single deviant epochs
sFilesEpochDeviantSrc = bst_process('CallProcess', 'process_select_files_results', [], [], ...
    'subjectname', SubjectName, ...
    'condition',   '', ...
    'includebad',  0);
sFilesEpochDeviantSrc = bst_process('CallProcess', 'process_select_tag', sFilesEpochDeviantSrc, [], ...
    'tag',    'deviant', ...
    'search', 1, ...  % Search the file names
    'select', 1);     % Select only the files with the tag
sFilesEpochDeviantSrc = bst_process('CallProcess', 'process_select_tag', sFilesEpochDeviantSrc, [], ...
    'tag',    'average', ...
    'search', 1, ...  % Search the file names
    'select', 2);     % Exclude the files with the tag

% Process: Time-frequency (Morlet wavelets)
sFilesTfScout = bst_process('CallProcess', 'process_timefreq', sFilesEpochDeviantSrc, [], ...
    'clusters',   {'Destrieux', {'G_temp_sup-G_T_transv L', 'G_temp_sup-G_T_transv R', 'G_front_inf-Opercular L', 'G_precentral L'}}, ...
    'scoutfunc',  1, ...  % Mean
    'edit', struct(...
         'Comment',         'Deviant: Scouts,Avg,Power,1-150Hz', ...
         'TimeBands',       [], ...
         '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], ...
         'MorletFc',        1, ...
         'MorletFwhmTc',    3, ...
         'ClusterFuncTime', 'after', ...
         'Measure',         'power', ...
         'Output',          'average', ...
         'SaveKernel',      0), ...
    'normalize',  'none');   % None: Save non-standardized time-frequency maps)
% Process: Event-related perturbation (ERS/ERD): [-75ms,0ms]
sFilesTfScoutErsd = bst_process('CallProcess', 'process_baseline_norm', sFilesTfScout, [], ...
    'baseline',  [-0.075, 0], ...
    'method',    'ersd', ...  % Event-related perturbation (ERS/ERD):    x_std = (x - μ) / μ * 100
    'overwrite', 0);
% Process: Snapshot: Time-frequency maps
bst_process('CallProcess', 'process_snapshot', sFilesTfScoutErsd, [], ...
    'target',   14, ...  % Time-frequency maps
    'Comment',  'ERS/ERD');


% === FULL CORTEX / HILBERT ===
% Process: Hilbert transform
sFilesHilbertCortex = bst_process('CallProcess', 'process_hilbert', sFilesEpochDeviantSrc, [], ...
    'clusters',  {}, ...
    'scoutfunc', 1, ...  % Mean
    'edit',      struct(...
         'Comment',         'Deviant: Avg,Magnitude', ...
         'TimeBands',       [], ...
         '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'}}, ...
         'ClusterFuncTime', 'none', ...
         'Measure',         'power', ...
         'Output',          'average', ...
         'RemoveEvoked',    0, ...
         'SaveKernel',      0), ...
    'normalize', 'none', ...  % None: Save non-standardized time-frequency maps
    'mirror',    0);
% Process: Event-related perturbation (ERS/ERD): [-75ms,0ms]
sFilesHilbertCortexErsd = bst_process('CallProcess', 'process_baseline_norm', sFilesHilbertCortex, [], ...
    'baseline',  [-0.075, 0], ...
    'method',    'ersd', ...  % Event-related perturbation (ERS/ERD):    x_std = (x - μ) / μ * 100
    'overwrite', 0);

% View results
hFigTf1 = view_surface_data([], sFilesHilbertCortexErsd.FileName);
hFigTf2 = view_timefreq(sFilesHilbertCortexErsd.FileName, 'SingleSensor', 362);
figure_3d('SetStandardView', hFigTf1, 'left');
panel_surface('SetDataThreshold', hFigTf1, 1, 0.5);
panel_time('SetCurrentTime', 0.175);
panel_freq('SetCurrentFreq', 3);
bst_colormaps('RestoreDefaults', 'stat2');
bst_report('Snapshot', [hFigTf1 hFigTf2], sFilesHilbertCortexErsd.FileName, 'Hilbert transform: Alpha band', [200, 200, 400, 250]);
bst_memory('UnloadAll', 'Forced');


%% ===== TUTORIAL #25: DIFFRERENCE ===================================================
%  ===================================================================================
disp([10 'DEMO> Tutorial #25: Difference' 10]);

% ===== SELECT TRIALS (DATA) =====
% Process: Select recordings in: Subject01/*
sFilesAll = bst_process('CallProcess', 'process_select_files_data', [], [], ...
    'subjectname', SubjectName, ...
    'condition',   '', ...
    'includebad',  0);
% Process: Select file names with tag: deviant_trial
sEpochDeviant = bst_process('CallProcess', 'process_select_tag', sFilesAll, [], ...
    'tag',    'deviant_trial', ...
    'search', 1, ...  % Search the file names
    'select', 1);     % Select only the files with the tag
% Process: Select file names with tag: standard_trial
sEpochStandard = bst_process('CallProcess', 'process_select_tag', sFilesAll, [], ...
    'tag',    'standard_trial', ...
    'search', 1, ...  % Search the file names
    'select', 1);     % Select only the files with the tag

% ===== SELECT TRIALS (SOURCE) =====
% Process: Select recordings in: Subject01/*
sFilesAllSrc = bst_process('CallProcess', 'process_select_files_results', [], [], ...
    'subjectname', SubjectName, ...
    'condition',   '', ...
    'includebad',  0);
% Process: Select file names with tag: deviant_trial
sEpochDeviantSrc = bst_process('CallProcess', 'process_select_tag', sFilesAllSrc, [], ...
    'tag',    'deviant_trial', ...
    'search', 1, ...  % Search the file names
    'select', 1);     % Select only the files with the tag
% Process: Select file names with tag: standard_trial
sEpochStandardSrc = bst_process('CallProcess', 'process_select_tag', sFilesAllSrc, [], ...
    'tag',    'standard_trial', ...
    'search', 1, ...  % Search the file names
    'select', 1);     % Select only the files with the tag

% ===== ABSOLUTE DIFFERENCE ======
% Process: Difference: A-B, abs
sDiffSrc = bst_process('CallProcess', 'process_diff_ab', sFilesIntraSrc(1).FileName, sFilesIntraSrc(2).FileName, ...
    'source_abs', 1);
% Process: Set comment: deviant|abs - standard|abs
sDiffSrc = bst_process('CallProcess', 'process_set_comment', sDiffSrc, [], ...
    'tag',     'deviant|abs - standard|abs', ...
    'isindex', 1);
% Process: Low-pass:40Hz
sDiffSrcZscore = bst_process('CallProcess', 'process_bandpass', sDiffSrc, [], ...
    'highpass',    0, ...
    'lowpass',     40, ...
    'attenuation', 'strict', ...  % 60dB
    'mirror',      0, ...
    'overwrite',   1);
% Process: Z-score transformation: [-100ms,-2ms]
sDiffSrcZscore = bst_process('CallProcess', 'process_baseline_norm', sDiffSrcZscore, [], ...
    'baseline',   [-0.1, -0.002], ...
    'source_abs', 0, ...
    'method',     'zscore', ...  % Z-score transformation:    x_std = (x - μ) / σ
    'overwrite',  0);
% Process: Snapshot: Sources (contact sheet)
bst_process('CallProcess', 'process_snapshot', sDiffSrcZscore, [], ...
    'target',         9, ...  % Sources (contact sheet)
    'modality',       1, ...  % MEG (All)
    'orient',         1, ...  % left
    'contact_time',   [0, 0.35], ...
    'contact_nimage', 15, ...
    'threshold',      30, ...
    'Comment',        'Difference deviant - standard (absolute)');

% ===== RELATIVE DIFFERENCE =====
% Process: Difference: A-B
sDiffSrcRel = bst_process('CallProcess', 'process_diff_ab', sFilesIntraSrc(1).FileName, sFilesIntraSrc(2).FileName, ...
    'source_abs', 0);
% Process: Set comment: deviant - standard
sDiffSrcRel = bst_process('CallProcess', 'process_set_comment', sDiffSrcRel, [], ...
    'tag',     'deviant - standard', ...
    'isindex', 1);
% Process: Low-pass:40Hz
sDiffSrcRelZscore = bst_process('CallProcess', 'process_bandpass', sDiffSrcRel, [], ...
    'highpass',    0, ...
    'lowpass',     40, ...
    'attenuation', 'strict', ...  % 60dB
    'mirror',      0, ...
    'overwrite',   1);
% Process: Z-score transformation: [-100ms,-2ms]
sDiffSrcRelZscore = bst_process('CallProcess', 'process_baseline_norm', sDiffSrcRelZscore, [], ...
    'baseline',   [-0.1, -0.002], ...
    'source_abs', 0, ...
    'method',     'zscore', ...  % Z-score transformation:    x_std = (x - μ) / σ
    'overwrite',  0);
% Configure colormap: hot/absolute
bst_colormaps('SetColormapName', 'stat2', 'hot');
bst_colormaps('SetColormapAbsolute', 'stat2', 1);
% Process: Snapshot: Sources (contact sheet)
bst_process('CallProcess', 'process_snapshot', sDiffSrcRelZscore, [], ...
    'target',         9, ...  % Sources (contact sheet)
    'modality',       1, ...  % MEG (All)
    'orient',         1, ...  % left
    'contact_time',   [0, 0.35], ...
    'contact_nimage', 15, ...
    'threshold',      30, ...
    'Comment',        'Difference deviant - standard (relative)');
% Restore colormap: rwb/relative
bst_colormaps('SetColormapName', 'stat2', 'cmap_rbw');
bst_colormaps('SetColormapAbsolute', 'stat2', 0);

% ===== DIFFERENCE OF MEANS =====
% Process: Select uniform number of files  [uniform]
[sEpochDeviantUni, sEpochStandardUni] = bst_process('CallProcess', 'process_select_uniform2', sEpochDeviant, sEpochStandard, ...
    'nfiles', 0, ...
    'method', 4);  % Uniformly distributed
% Process: Difference of means [mean]
sDiffMean = bst_process('CallProcess', 'process_diff_mean', sEpochDeviantUni, sEpochStandardUni, ...
    'avg_func', 1, ...  % Arithmetic average mean(A) - mean(B)
    'weighted', 0);
% Process: Snapshot: Recordings time series
bst_process('CallProcess', 'process_snapshot', sDiffMean, [], ...
    'target',     5, ...  % Recordings time series
    'modality',   1, ...  % MEG (All)
    'Comment',    'Difference of means');



%% ===== TUTORIAL #26: STATISTICS ====================================================
%  ===================================================================================
disp([10 'DEMO> Tutorial #26: Statistics' 10]);
% ===== HISTOGRAMS =====
% Process: Extract values: [160ms] MLP57
sHistoDeviant = bst_process('CallProcess', 'process_extract_values', sEpochDeviant, [], ...
    'timewindow',  [0.16, 0.16], ...
    'sensortypes', 'MLP57', ...
    'isabs',       0, ...
    'avgtime',     0, ...
    'avgrow',      0, ...
    'dim',         2, ...  % Concatenate time (dimension 2)
    'Comment',     '');
% Process: Extract values: [160ms] MLP57
sHistoStandard = bst_process('CallProcess', 'process_extract_values', sEpochStandard, [], ...
    'timewindow',  [0.16, 0.16], ...
    'sensortypes', 'MLP57', ...
    'isabs',       0, ...
    'avgtime',     0, ...
    'avgrow',      0, ...
    'dim',         2, ...  % Concatenate time (dimension 2)
    'Comment',     '');
% Display histograms
hFigHisto = view_histogram({sHistoDeviant.FileName, sHistoStandard.FileName});
bst_report('Snapshot', hFigHisto, sHistoDeviant.FileName, 'Histograms for MLP57/160ms');
close(hFigHisto);

% ===== EXEMPLE #1: PARAMETRIC/DATA =====
% Process: t-test [equal] [-100ms,500ms]          H0:(A-B = 0)
sTestParamData = bst_process('CallProcess', 'process_test_parametric2', sEpochDeviant, sEpochStandard, ...
    'timewindow',    [-0.1, 0.5], ...
    'sensortypes',   '', ...
    'isabs',         0, ...
    'avgtime',       0, ...
    'avgrow',        0, ...
    'Comment',       '', ...
    '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
    'tail',          'two');  % Two-tailed
% Set display properties
StatThreshOptions = bst_get('StatThreshOptions');
StatThreshOptions.pThreshold = 0.05;
StatThreshOptions.Correction = 'fdr';
StatThreshOptions.Control    = [1 2 3];
bst_set('StatThreshOptions', StatThreshOptions);
% Process: Snapshot: Recordings time series
bst_process('CallProcess', 'process_snapshot', sTestParamData, [], ...
    'target',     5, ...  % Recordings time series
    'modality',   1, ...  % MEG (All)
    'time',       0.16, ...
    'Comment',    'Parametric t-test (p<0.05, FDR)');
% Process: Snapshot: Recordings topography (one time)
bst_process('CallProcess', 'process_snapshot', sTestParamData, [], ...
    'target',         6, ...  % Recordings topography (one time)
    'modality',       1, ...  % MEG (All)
    'time',           0.16, ...
    'Comment',        'Parametric t-test (p<0.05, FDR)');

% ===== EXEMPLE #2: NON-PARAMETRIC/DATA =====
% Process: Perm t-test equal [-100ms,500ms MEG]          H0:(A=B), H1:(A<>B)
sTestPermData = bst_process('CallProcess', 'process_test_permutation2', sEpochDeviant, sEpochStandard, ...
    'timewindow',     [-0.1, 0.5], ...
    'sensortypes',    'MEG', ...
    'isabs',          0, ...
    'avgtime',        0, ...
    'avgrow',         0, ...
    'iszerobad',      1, ...
    'Comment',        '', ...
    '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))
    'randomizations', 1000, ...
    'tail',           'two');  % Two-tailed
% Process: Snapshot: Recordings time series
bst_process('CallProcess', 'process_snapshot', sTestPermData, [], ...
    'target',     5, ...  % Recordings time series
    'modality',   1, ...  % MEG (All)
    'time',       0.16, ...
    'Comment',    'Non-parametric t-test (p<0.05, FDR)');
% Process: Snapshot: Recordings topography (one time)
bst_process('CallProcess', 'process_snapshot', sTestPermData, [], ...
    'target',         6, ...  % Recordings topography (one time)
    'modality',       1, ...  % MEG (All)
    'time',           0.16, ...
    'Comment',        'Non-parametric t-test (p<0.05, FDR)');

% ===== EXEMPLE #3: CLUSTER/DATA =====
% Process: FT t-test unequal cluster [-100ms,500ms MEG]          H0:(A=B), H1:(A<>B)
sTestClustData = bst_process('CallProcess', 'process_ft_timelockstatistics', sEpochDeviant, sEpochStandard, ...
    'sensortypes',    'MEG', ...
    'timewindow',     [-0.1, 0.5], ...
    'isabs',          0, ...
    'avgtime',        0, ...
    'avgchan',        0, ...
    'randomizations', 1000, ...
    'statistictype',  1, ...      % Independent t-test
    'tail',           'two', ...  % Two-tailed
    'correctiontype', 2, ...      % cluster
    'minnbchan',      0, ...
    'clusteralpha',   0.05);
% Process: Snapshot: Recordings time series
bst_process('CallProcess', 'process_snapshot', sTestClustData, [], ...
    'target',    5, ...  % Recordings time series
    'time',      0.16, ...
    'modality',  1, ...  % MEG (All)
    'Comment',   'Cluster-based permutation test');
% Process: Snapshot: Recordings topography (one time)
bst_process('CallProcess', 'process_snapshot', sTestClustData, [], ...
    'target',    6, ...  % Recordings topography (one time)
    'modality',  1, ...  % MEG (All)
    'time',      0.16, ...
    'Comment',   'Cluster-based permutation test');

% ===== EXAMPLE #4: PARAMETRIC/SOURCES =====
% Process: t-test [equal] [-100ms,500ms]          H0:(A-B = 0)
sTestParamSrc = bst_process('CallProcess', 'process_test_parametric2', sEpochDeviantSrc, sEpochStandardSrc, ...
    'timewindow',    [-0.1, 0.5], ...
    'scoutsel',      {}, ...
    'scoutfunc',     1, ...  % Mean
    'isnorm',        0, ...
    'avgtime',       0, ...
    'Comment',       '', ...
    '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
    'tail',          'two');  % Two-tailed
% Process: Difference of means [abs(mean)]
sDiffMeanSrc = bst_process('CallProcess', 'process_diff_mean', sEpochDeviantSrc, sEpochStandardSrc, ...
    'avg_func',   2, ...  % Absolute value of average abs(mean(A)) - abs(mean(B))
    'weighted',   0);
% Process: Apply statistic threshold: p<0.05 (FDR:1,2,3)
sDiffMeanSrcThresh = bst_process('CallProcess', 'process_extract_pthresh2', sTestParamSrc, sDiffMeanSrc, ...
    'pthresh',    0.05, ...
    'correction', 3, ...  % False discovery rate (FDR)
    'control1',   1, ...
    'control2',   1, ...
    'control3',   1);
% Process: Snapshot: Sources (one time)
bst_process('CallProcess', 'process_snapshot', sTestParamSrc, [], ...
    'target',    8, ...  % Sources (one time)
    'orient',    1, ...  % left
    'time',      0.148, ...
    'threshold', 40, ...
    'Comment',   'Parametric t-test (p<0.05, FDR)');
% Process: Snapshot: Sources (one time)
bst_process('CallProcess', 'process_snapshot', sDiffMeanSrc, [], ...
    'target',    8, ...  % Sources (one time)
    'orient',    1, ...  % left
    'time',      0.148, ...
    'threshold', 40, ...
    'Comment',   'abs(average(deviant)) - abs(average(standard))');
% Process: Snapshot: Sources (one time)
bst_process('CallProcess', 'process_snapshot', sDiffMeanSrcThresh, [], ...
    'target',    8, ...  % Sources (one time)
    'orient',    1, ...  % left
    'time',      0.148, ...
    'threshold', 0, ...
    'Comment',   'Different of mean thresholded with t-test results');

% ===== EXAMPLE #5: PARAMETRIC/SCOUTS =====
% Process: t-test equal [-100ms,500ms]          H0:(A=B), H1:(A<>B)
sTestParamScout = bst_process('CallProcess', 'process_test_parametric2', sEpochDeviantSrc, sEpochStandardSrc, ...
    'timewindow',    [-0.1, 0.5], ...
    'scoutsel',      {'Destrieux', {'G_front_inf-Opercular L', 'G_precentral L', 'G_temp_sup-G_T_transv L'}}, ...
    'scoutfunc',     1, ...  % Mean
    'isnorm',        0, ...
    'avgtime',       0, ...
    'Comment',       '', ...
    '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
    'tail',          'two');  % Two-tailed
% Process: Apply statistic threshold: p<0.05 (FDR:1,2,3)
sTestParamScoutThresh = bst_process('CallProcess', 'process_extract_pthresh', sTestParamScout, [], ...
    'pthresh',    0.05, ...
    'correction', 3, ...  % False discovery rate (FDR)
    'control1',   1, ...
    'control2',   1, ...
    'control3',   1);
% Process: Compute head model
bst_process('CallProcess', 'process_headmodel', sTestParamScoutThresh, [], ...
    'sourcespace', 1, ...  % Cortex surface
    'meg',         3);     % Overlapping spheres
% Process: Simulate recordings from scouts
sSimulData = bst_process('CallProcess', 'process_simulate_recordings', sTestParamScoutThresh, [], ...
    'scouts',      {'Destrieux', {'G_front_inf-Opercular L', 'G_precentral L', 'G_temp_sup-G_T_transv L'}}, ...
    'savesources', 1);

% Get corresponding source file
[sStudy,iStudy,iRes] = bst_get('ResultsForDataFile', sSimulData.FileName);
sSimulSrc = sStudy.Result(iRes).FileName;
% Reset visualization filters
panel_filter('SetFilters', 0, [], 0, [], 0, [], 0, 0);
% Process: Snapshot: Recordings time series
bst_process('CallProcess', 'process_snapshot', sTestParamScout, [], ...
    'target',    5, ...  % Recordings time series
    'time',      0.148, ...
    'Comment',   'Parametric t-test (p<0.05, FDR)');
% Process: Snapshot: Sources (one time)
bst_process('CallProcess', 'process_snapshot', sSimulSrc, [], ...
    'target',    8, ...  % Sources (one time)
    'orient',    1, ...  % left
    'time',      0.148, ...
    'threshold', 0, ...
    'Comment',   'Simulated sources');
% Process: Snapshot: Recordings time series
bst_process('CallProcess', 'process_snapshot', sSimulData, [], ...
    'target',    5, ...  % Recordings time series
    'modality',  1, ...  % MEG (All)
    'time',      0.148, ...
    'Comment',   'Simulated MEG recordings');


% ===== EXAMPLE #6: NON-PARAMETRIC/TIMEFREQ =====
TfOptions.Output = 'all';
% Process: Time-frequency (Morlet wavelets) / DEVIANT
sEpochDeviantTf = bst_process('CallProcess', 'process_timefreq', sEpochDeviant, [], ...
    'sensortypes', 'MLP57', ...
    'edit',        TfOptions, ...
    'normalize',   'none');  % None: Save non-standardized time-frequency maps
% Process: Time-frequency (Morlet wavelets) / STANDARD
sEpochStandardTf = bst_process('CallProcess', 'process_timefreq', sEpochStandard, [], ...
    'sensortypes', 'MLP57', ...
    'edit',        TfOptions, ...
    'normalize',   'none');  % None: Save non-standardized time-frequency maps
% Process: Perm t-test equal [-100ms,500ms 1-150Hz]          H0:(A=B), H1:(A<>B)
sTestTf = bst_process('CallProcess', 'process_test_permutation2', sEpochDeviantTf, sEpochStandardTf, ...
    'timewindow',     [-0.1, 0.5], ...
    'freqrange',      [1, 150], ...
    'rows',           '', ...
    'isabs',          0, ...
    'avgtime',        0, ...
    'avgrow',         0, ...
    'avgfreq',        0, ...
    'matchrows',      0, ...
    'iszerobad',      1, ...
    'Comment',        '', ...
    '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))
    'randomizations', 1000, ...
    'tail',           'two');  % Two-tailed
% Set stat threshold
StatThreshOptions.pThreshold = 0.05;
StatThreshOptions.Correction = 'none';
bst_set('StatThreshOptions', StatThreshOptions);
% Process: Snapshot: Time-frequency maps
bst_process('CallProcess', 'process_snapshot', sTestTf, [], ...
    'target',   14, ...  % Time-frequency maps
    'Comment',  'Non-parametric t-test (p<0.05, Uncorrected)');
% Process: Delete intermediate results
bst_process('CallProcess', 'process_delete', [sEpochDeviantTf, sEpochStandardTf], [], ...
    'target', 1);  % Delete selected files



%% ===== SAVE REPORT =====
% Save and display report
ReportFile = bst_report('Save', []);
if ~isempty(reports_dir) && ~isempty(ReportFile)
    bst_report('Export', ReportFile, reports_dir);
else
    bst_report('Open', ReportFile);
end

disp([10 'DEMO> Done.' 10]);

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

function tutorial_visual_single(bids_dir, reports_dir)
% TUTORIAL_VISUAL_SINGLE: Runs the Brainstorm/SPM group analysis pipeline (single subject, BIDS version).
%
% ONLINE TUTORIALS: https://neuroimage.usc.edu/brainstorm/Tutorials/VisualSingle
%
% INPUTS:
%    - bids_dir: Path to folder ds000117 (https://openneuro.org/datasets/ds000117)
%       |- derivatives/freesurfer/sub-XX                               : Segmentation folders generated with FreeSurfer
%       |- derivatives/meg_derivatives/sub-XX/ses-meg/meg/*.fif        : MEG+EEG recordings (processed with MaxFilter's SSS)
%       |- derivatives/meg_derivatives/sub-emptyroom/ses-meg/meg/*.fif : Empty room measurements
%    - reports_dir: If defined, exports all the reports as HTML to this folder

% @=============================================================================
% This function is part of the Brainstorm software:
% https://neuroimage.usc.edu/brainstorm
% 
% Copyright (c) University of Southern California & McGill University
% This software is distributed under the terms of the GNU General Public License
% as published by the Free Software Foundation. Further details on the GPLv3
% license can be found at http://www.gnu.org/copyleft/gpl.html.
% 
% FOR RESEARCH PURPOSES ONLY. THE SOFTWARE IS PROVIDED "AS IS," AND THE
% UNIVERSITY OF SOUTHERN CALIFORNIA AND ITS COLLABORATORS DO NOT MAKE ANY
% WARRANTY, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO WARRANTIES OF
% MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, NOR DO THEY ASSUME ANY
% LIABILITY OR RESPONSIBILITY FOR THE USE OF THIS SOFTWARE.
%
% For more information type "brainstorm license" at command prompt.
% =============================================================================@
%
% Author: Francois Tadel, Elizabeth Bock, 2016-2018


%% ===== SCRIPT VARIABLES =====
% Full list of subjects to process
SubjectNames = {'sub-01', 'sub-02', 'sub-03', 'sub-04', 'sub-05', 'sub-06', 'sub-07', 'sub-08', ...
                'sub-09', 'sub-10', 'sub-11', 'sub-12', 'sub-13', 'sub-14', 'sub-15', 'sub-16'};
% Empty-room dates for each subject (so that we can match automatically recordings with empty-room)
EmptyRoomSubj = 'sub-emptyroom';
AcquisitionDates = {'09-Apr-2009', '06-May-2009', '11-May-2009', '18-May-2009', '15-May-2009', '15-May-2009', '15-May-2009', '15-May-2009', ...
                    '15-May-2009', '15-May-2009', '01-Jun-2009', '01-Jun-2009', '01-Jun-2009', '26-Nov-2009', '08-Dec-2009', '08-Dec-2009'};
% Bad channels {iSubj} = {Run01, Run02, Run03, Run04, Run05, Run06}
BadChannels{1}  = {'EEG016', 'EEG070', 'EEG050',{'EEG008','EEG050'}, [], []};
BadChannels{2}  = {{'EEG027', 'EEG030', 'EEG038'}, 'EEG010', 'EEG010', 'EEG010', 'EEG010', 'EEG010'};
BadChannels{3}  = {{'EEG008','EEG017'}, {'EEG008','EEG017'}, {'EEG008','EEG017'}, {'EEG008','EEG017'}, {'EEG008','EEG017','EEG001'}, {'EEG008','EEG017','EEG020'}};
BadChannels{4}  = {{'EEG038'}, {'EEG038','EEG001','EEG016'}, {'EEG038','EEG001','EEG016'}, {'EEG038','EEG001'}, {'EEG038','EEG001','EEG016'}, {'EEG038','EEG001','EEG016'}};
BadChannels{5}  = {'EEG001', 'EEG001', [], [], [], []};
BadChannels{6}  = {'EEG068', [], 'EEG004', [], [], []};
BadChannels{7}  = {[], [], {'EEG004','EEG008'}, {'EEG004','EEG008'},{'EEG004','EEG008','EEG043','EEG045','EEG047'}, {'EEG004','EEG008'}};
BadChannels{8}  = {[], [], [], [], [], []};
BadChannels{9}  = {[], 'EEG004', 'EEG004', [], 'EEG004', 'EEG004'};
BadChannels{10} = {[], [], [], [], [], []};
BadChannels{11} = {{'EEG010','EEG050'}, 'EEG050', 'EEG050', 'EEG050', 'EEG050', 'EEG050'};
BadChannels{12} = {{'EEG024','EEG057'}, {'EEG024','EEG057'}, {'EEG024','EEG057'}, {'EEG024','EEG057','EEG070'}, {'EEG024','EEG057'}, {'EEG024','EEG057','EEG070'}};
BadChannels{13} = {'EEG009', 'EEG009', {'EEG009','EEG057','EEG69'}, 'EEG009', {'EEG009','EEG044'}, {'EEG009','EEG044'}};
BadChannels{14} = {'EEG029', 'EEG029', 'EEG029', {'EEG004','EEG008','EEG016','EEG029'}, {'EEG004','EEG008','EEG016','EEG029'}, {'EEG004','EEG008','EEG016','EEG029'}};
BadChannels{15} = {'EEG038', 'EEG038', 'EEG038', 'EEG038', {'EEG054','EEG038'}, 'EEG038'};
BadChannels{16} = {'EEG008', 'EEG008', 'EEG008', 'EEG008', 'EEG008', 'EEG008'};
% SSP components to remove {iSubj} = {sRun01, sRun02, sRun03, sRun03, sRun04, sRun05, sRun06},   sRun0X={ECG_GRAD,ECG_MAG}
SspSelect{1}  = {{1,1},         {1,1},   {1,1},   {1,1},   {1,1},   {1,1}}; 
SspSelect{2}  = {{1,1},         {1,1},   {1,1},   {1,1},   {3,1},   {1,1}};
SspSelect{3}  = {{[],1},        {1,1},   {1,1},   {1,1},   {1,1},   {1,1}};
SspSelect{4}  = {{[],1},        {[],1},  {[],1},  {[],1},  {[],1},  {[],1}};
SspSelect{5}  = {{2,1},         {1,1},   {1,1},   {[],1},  {1,1},   {1,1}};
SspSelect{6}  = {{2,1},         {2,1},   {1,1},   {2,1},   {1,1},   {2,1}};
SspSelect{7}  = {{1,1},         {1,1},   {1,1},   {1,1},   {1,1},   {1,1}};
SspSelect{8}  = {{1,1},         {1,1},   {[],1},  {2,1},   {1,1},   {2,1}};
SspSelect{9}  = {{1,1},         {1,1},   {[],1},  {[],1},  {1,1},   {1,1}};
SspSelect{10} = {{1,1},         {1,1},   {1,1},   {1,1},   {1,1},   {1,1}};
SspSelect{11} = {{[],1},        {[],1},  {[],1},  {[],1},  {[],[]}, {[],[]}};
SspSelect{12} = {{[1,2],[1,2]}, {1,1},   {1,1},   {1,1},   {1,1},   {1,1}};
SspSelect{13} = {{[],[]},       {[],[]}, {[],[]}, {[],[]}, {[],[]}, {[],[]}};
SspSelect{14} = {{1,1},         {1,1},   {1,1},   {1,1},   {1,1},   {1,1}};
SspSelect{15} = {{1,1},         {1,1},   {1,1},   {1,1},   {1,1},   {1,1}};
SspSelect{16} = {{1,1},         {1,1},   {1,1},   {2,1},   {1,1},   {1,1}};


%% ===== CREATE PROTOCOL =====
% Start brainstorm without the GUI
if ~brainstorm('status')
    brainstorm nogui
end
% Output folder for reports
if (nargin < 2) || isempty(reports_dir) || ~isdir(reports_dir)
    reports_dir = [];
end
% You have to specify the folder in which the tutorial dataset is unzipped
if (nargin < 1) || isempty(bids_dir) || ~file_exist(bids_dir) || ~file_exist(bst_fullfile(bids_dir, 'derivatives')) || ~file_exist(bst_fullfile(bids_dir, 'dataset_description.json'))
    error('The first argument must be the full path to the tutorial folder.');
end
% The protocol name has to be a valid folder name (no spaces, no weird characters...)
ProtocolName = 'TutorialVisual';
% Delete existing protocol
gui_brainstorm('DeleteProtocol', ProtocolName);
% Create new protocol
gui_brainstorm('CreateProtocol', ProtocolName, 0, 0);
% Set visualization filters: 40Hz low-pass, no high-pass
panel_filter('SetFilters', 1, 40, 0, [], 0, [], 0, 0);
% Set colormap: local color scale
bst_colormaps('SetMaxMode', 'meg', 'local');
bst_colormaps('SetMaxMode', 'eeg', 'local');


%% ===== PRE-PROCESS AND IMPORT =====
for iSubj = 1:16
    % Start a new report (one report per subject)
    bst_report('Start');
    disp(sprintf('\n===== IMPORT: SUBJECT #%d =====\n', iSubj));
    
    % If subject already exists: delete it
    [sSubject, iSubject] = bst_get('Subject', SubjectNames{iSubj});
    if ~isempty(sSubject)
        db_delete_subjects(iSubject);
    end
    
    % ===== FILES TO IMPORT =====
    % Build the path of the files to import
    AnatDir    = fullfile(bids_dir, 'derivatives', 'freesurfer', SubjectNames{iSubj}, 'ses-mri', 'anat');
    DataDir    = fullfile(bids_dir, 'derivatives',  'meg_derivatives', SubjectNames{iSubj}, 'ses-meg', 'meg');
    % Check if the folder contains the required files
    if ~file_exist(AnatDir)
        error(['The folder "' AnatDir '" does not exist.']);
    end
    if ~file_exist(DataDir)
        error(['The folder "' DataDir '" does not exist.']);
    end

    % ===== ANATOMY =====
    % Process: Import anatomy folder
    bst_process('CallProcess', 'process_import_anatomy', [], [], ...
        'subjectname', SubjectNames{iSubj}, ...
        'mrifile',     {AnatDir, 'FreeSurfer'}, ...
        'nvertices',   15000);
    
    % ===== PROCESS EACH RUN =====
    for iRun = 1:6
        % Files to import
        FifFile = bst_fullfile(DataDir, sprintf('%s_ses-meg_task-facerecognition_run-%02d_proc-sss_meg.fif', SubjectNames{iSubj}, iRun));

        % ===== LINK CONTINUOUS FILE =====
        % Process: Create link to raw file
        sFileRaw = bst_process('CallProcess', 'process_import_data_raw', [], [], ...
            'subjectname',    SubjectNames{iSubj}, ...
            'datafile',       {FifFile, 'FIF'}, ...
            'channelreplace', 1, ...
            'channelalign',   0);
        % Set acquisition date
        panel_record('SetAcquisitionDate', sFileRaw.iStudy, AcquisitionDates{iSubj});

        % ===== PREPARE CHANNEL FILE =====
        % Process: Set channels type
        bst_process('CallProcess', 'process_channel_settype', sFileRaw, [], ...
            'sensortypes', 'EEG061, EEG064', ...
            'newtype',     'NOSIG');
        bst_process('CallProcess', 'process_channel_settype', sFileRaw, [], ...
            'sensortypes', 'EEG062', ...
            'newtype',     'EOG');
        bst_process('CallProcess', 'process_channel_settype', sFileRaw, [], ...
            'sensortypes', 'EEG063', ...
            'newtype',     'ECG');

        % Process: Remove head points
        sFileRaw = bst_process('CallProcess', 'process_headpoints_remove', sFileRaw, [], ...
            'zlimit', 0);
        % Process: Refine registration
        sFileRaw = bst_process('CallProcess', 'process_headpoints_refine', sFileRaw, []);
        % Process: Project electrodes on scalp
        sFileRaw = bst_process('CallProcess', 'process_channel_project', sFileRaw, []);

        % Process: Snapshot: Sensors/MRI registration
        bst_process('CallProcess', 'process_snapshot', sFileRaw, [], ...
            'target',   1, ...  % Sensors/MRI registration
            'modality', 1, ...  % MEG (All)
            'orient',   1, ...  % left
            'Comment',  sprintf('MEG/MRI Registration: Subject #%d, Run #%d', iSubj, iRun));
        bst_process('CallProcess', 'process_snapshot', sFileRaw, [], ...
            'target',   1, ...  % Sensors/MRI registration
            'modality', 4, ...  % EEG
            'orient',   1, ...  % left
            'Comment',  sprintf('EEG/MRI Registration: Subject #%d, Run #%d', iSubj, iRun));

        % ===== IMPORT TRIGGERS =====
        % Process: Read from channel
        bst_process('CallProcess', 'process_evt_read', sFileRaw, [], ...
            'stimchan',  'STI101', ...
            'trackmode', 2, ...  % Bit: detect the changes for each bit independently
            'zero',      0);
        % Process: Group by name
        bst_process('CallProcess', 'process_evt_groupname', sFileRaw, [], ...
            'combine', 'Unfamiliar=3,4', ...
            'dt',      0, ...
            'delete',  1);
        % Process: Rename event
        bst_process('CallProcess', 'process_evt_rename', sFileRaw, [], ...
            'src',  '3', ...
            'dest', 'Famous');
        % Process: Rename event
        bst_process('CallProcess', 'process_evt_rename', sFileRaw, [], ...
            'src',  '5', ...
            'dest', 'Scrambled');
        % Process: Add time offset
        bst_process('CallProcess', 'process_evt_timeoffset', sFileRaw, [], ...
            'info',      [], ...
            'eventname', 'Famous, Unfamiliar, Scrambled', ...
            'offset',    0.0345);
        % Process: Delete events
        bst_process('CallProcess', 'process_evt_delete', sFileRaw, [], ...
            'eventname', '1,2,6,7,8,9,10,11,12,13,14,15,16');
        % Process: Detect cHPI activity (Elekta):STI201
        bst_process('CallProcess', 'process_evt_detect_chpi', sFileRaw, [], ...
            'eventname',   'chpi_bad', ...
            'channelname', 'STI201', ...
            'method',      'off');  % Mark as bad when the HPI coils are OFF
        
        % ===== FREQUENCY FILTERS =====
        % Process: Notch filter: 50Hz 100Hz 150Hz 200Hz
        sFileClean = bst_process('CallProcess', 'process_notch', sFileRaw, [], ...
            'freqlist',    [50, 100, 150, 200], ...
            'sensortypes', 'MEG, EEG', ...
            'read_all',    0);
        % Process: Power spectrum density (Welch)
        sFilesPsd = bst_process('CallProcess', 'process_psd', [sFileRaw, sFileClean], [], ...
            'timewindow',  [], ...
            'win_length',  4, ...
            'win_overlap', 50, ...
            'sensortypes', 'MEG, EEG', ...
            'edit', struct(...
                 'Comment',         'Power', ...
                 'TimeBands',       [], ...
                 'Freqs',           [], ...
                 'ClusterFuncTime', 'none', ...
                 'Measure',         'power', ...
                 'Output',          'all', ...
                 'SaveKernel',      0));
        % Process: Snapshot: Frequency spectrum
        bst_process('CallProcess', 'process_snapshot', sFilesPsd, [], ...
            'target',   10, ...  % Frequency spectrum
            'Comment',  sprintf('Power spctrum: Subject #%d, Run #%d', iSubj, iRun));

        % ===== BAD CHANNELS =====
        if ~isempty(BadChannels{iSubj}{iRun})
            % Process: Set bad channels
            bst_process('CallProcess', 'process_channel_setbad', sFileClean, [], ...
                'sensortypes', BadChannels{iSubj}{iRun});
        end

        % ===== EEG REFERENCE =====
        % Process: Re-reference EEG
        bst_process('CallProcess', 'process_eegref', sFileClean, [], ...
            'eegref',      'AVERAGE', ...
            'sensortypes', 'EEG');
    
        % ===== DETECT ARTIFACTS ======
        % Process: Detect heartbeats
        bst_process('CallProcess', 'process_evt_detect_ecg', sFileClean, [], ...
            'channelname', 'EEG063', ...
            'timewindow',  [], ...
            'eventname',   'cardiac');
        % Different amplitude thresholds for different subjects
        if strcmpi(SubjectNames{iSubj}, 'sub-05')
            thresholdMAX = 50;
        else
            thresholdMAX = 100;
        end
        % Process: Detect: blink_BAD    -   Detects all events where the amplitude exceeds 100uV
        bst_process('CallProcess', 'process_evt_detect_threshold', sFileClean, [], ...
            'eventname',    'blink_BAD', ...
            'channelname',  'EEG062', ...
            'timewindow',   [], ...
            'thresholdMAX', thresholdMAX, ...
            'units',        3, ...  % uV (10^-6)
            'bandpass',     [0.3, 20], ...
            'isAbsolute',   1, ...
            'isDCremove',   0);

        % ===== SSP COMPUTATION =====
        % Process: SSP ECG: cardiac
        bst_process('CallProcess', 'process_ssp_ecg', sFileClean, [], ...
            'eventname',   'cardiac', ...
            'sensortypes', 'MEG GRAD', ...
            'usessp',      1, ...
            'select',      SspSelect{iSubj}{iRun}{1}); 
        bst_process('CallProcess', 'process_ssp_ecg', sFileClean, [], ...
            'eventname',   'cardiac', ...
            'sensortypes', 'MEG MAG', ...
            'usessp',      1, ...
            'select',      SspSelect{iSubj}{iRun}{2});
        % Process: Snapshot: SSP projectors
        bst_process('CallProcess', 'process_snapshot', sFileClean, [], ...
            'target',   2, ...
            'Comment',  sprintf('Subject #%d, Run #%d', iSubj, iRun));   % SSP projectors
        
        % ===== IMPORT BAD EVENTS =====
        % Get bad segments: this is typically done manually, not from a script
        BadSegments = GetBadSegments(iSubj, iRun);
        % Process: Import from file
        bst_process('CallProcess', 'process_evt_import', sFileClean, [], ...
            'evtfile', {BadSegments, 'ARRAY-TIMES'}, ...
            'evtname', 'BAD');
        
        % ===== IMPORT TRIALS =====
        % Process: Import MEG/EEG: Events
        sFilesEpochs = bst_process('CallProcess', 'process_import_data_event', sFileClean, [], ...
            'subjectname', SubjectNames{iSubj}, ...
            'condition',   '', ...
            'eventname',   'Famous, Scrambled, Unfamiliar', ...
            'timewindow',  [], ...
            'epochtime',   [-0.5, 1.2], ...
            'createcond',  0, ...
            'ignoreshort', 1, ...
            'usectfcomp',  1, ...
            'usessp',      1, ...
            'freq',        [], ...
            'baseline',    [-0.5, -0.0009]);
        
        % ===== AVERAGE: RUN =====
        % Process: Average: By trial group (folder average)
        sFilesAvg = bst_process('CallProcess', 'process_average', sFilesEpochs, [], ...
            'avgtype',    5, ...  % By trial group (folder average)
            'avg_func',   1, ...  % Arithmetic average:  mean(x)
            'weighted',   0, ...
            'keepevents', 0);
        % Process: Snapshot: Recordings time series
        bst_process('CallProcess', 'process_snapshot', sFilesAvg, [], ...
            'target',   5, ...  % Recordings time series
            'modality', 4, ...  % EEG
            'time',     0.11, ...
            'Comment',  sprintf('Subject #%d, Run #%d', iSubj, iRun));
        % Process: Snapshot: Recordings topography
        bst_process('CallProcess', 'process_snapshot', sFilesAvg, [], ...
            'target',   6, ...  % Recordings topography (one time)
            'modality', 4, ...  % EEG
            'time',     0.11, ...
            'Comment',  sprintf('Subject #%d, Run #%d', iSubj, iRun));

        % ===== COMPUTE NOISECOV: EEG =====
        % Process: Compute covariance (noise or data)
        bst_process('CallProcess', 'process_noisecov', sFilesEpochs, [], ...
            'baseline',       [-0.5, -0.0009], ...
            'sensortypes',    'EEG', ...
            'target',         1, ...  % Noise covariance     (covariance over baseline time window)
            'dcoffset',       1, ...  % Block by block, to avoid effects of slow shifts in data
            'identity',       0, ...
            'copycond',       0, ...
            'copysubj',       0, ...
            'replacefile',    1);  % Replace
    end
    
    % Save report
    ReportFile = bst_report('Save', []);
    if ~isempty(reports_dir) && ~isempty(ReportFile)
        bst_report('Export', ReportFile, bst_fullfile(reports_dir, ['report_' ProtocolName '_' SubjectNames{iSubj} '.html']));
    end
end


%% ===== EMPTY ROOM RECORDINGS =====
disp(sprintf('\n===== IMPORT: EMPTY-ROOM =====\n'));
% Loop on all the noise sessions
NoiseFiles = {};
for ses = {'20090409', '20090506', '20090511', '20090515', '20090518', '20090601', '20091126', '20091208'}
    NoiseFiles{end+1} = fullfile(bids_dir, 'derivatives', 'meg_derivatives', EmptyRoomSubj, ['ses-' ses{1}], 'meg', ['sub-emptyroom_ses-' ses{1} '_task-noise_proc-sss_meg.fif']);
end
% Process: Create link to raw file
sFilesNoise = bst_process('CallProcess', 'process_import_data_raw', [], [], ...
    'subjectname',    EmptyRoomSubj, ...
    'datafile',       {NoiseFiles, 'FIF'}, ...
    'channelreplace', 1, ...
    'channelalign',   0);
% Process: Notch filter: 50Hz 100Hz 150Hz 200Hz
sFileNoiseClean = bst_process('CallProcess', 'process_notch', sFilesNoise, [], ...
        'freqlist',    [50, 100, 150, 200], ...
        'sensortypes', 'MEG, EEG', ...
        'read_all',    0);
% Process: Compute noise covariance
bst_process('CallProcess', 'process_noisecov', sFileNoiseClean, [], ...
    'baseline',    [], ...
    'sensortypes', 'MEG', ...
    'target',      1, ...  % Noise covariance     (covariance over baseline time window)
    'dcoffset',    1, ...  % Block by block, to avoid effects of slow shifts in data
    'identity',    0, ...
    'copycond',    1, ...
    'copysubj',    1, ...
    'copymatch',   1, ...
    'replacefile', 2);  % Merge
       

%% ===== SOURCE ESTIMATION =====
% Start a new report (one report for the source estimation of all the subjects)
bst_report('Start');
% Loop on the subjects: This loop is separated from the previous one, because we should 
% compute the BEM surfaces after importing all the runs, so that the registration is done 
% using the high resolution head surface, instead of the smooth scalp BEM layer.
for iSubj = 1:length(SubjectNames)
    disp(sprintf('\n===== SOURCES: SUBJECT #%d =====\n', iSubj));
    
    % ===== BEM SURFACES =====
    % Process: Generate BEM surfaces
    bst_process('CallProcess', 'process_generate_bem', [], [], ...
        'subjectname', SubjectNames{iSubj}, ...
        'nscalp',      1082, ...
        'nouter',      642, ...
        'ninner',      642, ...
        'thickness',   4, ...
        'method',      'brainstorm');

    % ===== SELECT ALL AVERAGES =====
    % Process: Select data files in: */*
    sFilesAvg = bst_process('CallProcess', 'process_select_files_data', [], [], ...
        'subjectname',   SubjectNames{iSubj});
    % Process: Select file comments with tag: Avg
    sFilesAvg = bst_process('CallProcess', 'process_select_tag', sFilesAvg, [], ...
        'tag',    'Avg');  % Select only the files with the tag

    % ===== COMPUTE HEAD MODELS =====
    % Process: Compute head model (only for the first run of the subject)
    bst_process('CallProcess', 'process_headmodel', sFilesAvg(1), [], ...
        'sourcespace', 1, ...  % Cortex surface
        'meg',         3, ...  % Overlapping spheres
        'eeg',         3, ...  % OpenMEEG BEM
        'ecog',        1, ...  % 
        'seeg',        1, ...  % 
        'openmeeg',    struct(...
             'BemSelect',    [1, 1, 1], ...
             'BemCond',      [1, 0.0125, 1], ...
             'BemNames',     {{'Scalp', 'Skull', 'Brain'}}, ...
             'BemFiles',     {{}}, ...
             'isAdjoint',    0, ...
             'isAdaptative', 1, ...
             'isSplit',      0, ...
             'SplitLength',  4000));
    % Get all the runs for this subject (ie the list of the study indices)
    iStudyOther = setdiff(unique([sFilesAvg.iStudy]), sFilesAvg(1).iStudy);
    % Copy the forward model file to the other runs
    sHeadmodel = bst_get('HeadModelForStudy', sFilesAvg(1).iStudy);
    for iStudy = iStudyOther
        db_add(iStudy, sHeadmodel.FileName);
    end
    
    % ===== COMPUTE SOURCES: MEG =====
    % Process: Compute sources [2018]
    sAvgSrcMeg = bst_process('CallProcess', 'process_inverse_2018', sFilesAvg, [], ...
        'output',  1, ...  % Kernel only: shared
        'inverse', struct(...
             'Comment',        'MN: MEG ALL', ...
             'InverseMethod',  'minnorm', ...
             'InverseMeasure', 'amplitude', ...
             'SourceOrient',   {{'fixed'}}, ...
             'Loose',          0.2, ...
             'UseDepth',       1, ...
             'WeightExp',      0.5, ...
             'WeightLimit',    10, ...
             'NoiseMethod',    'reg', ...
             'NoiseReg',       0.1, ...
             'SnrMethod',      'fixed', ...
             'SnrRms',         1e-06, ...
             'SnrFixed',       3, ...
             'ComputeKernel',  1, ...
             'DataTypes',      {{'MEG GRAD', 'MEG MAG'}}));
    % Process: Snapshot: Sources (one time)   -   Loop only to get a correct comment for the report
    for i = 1:length(sAvgSrcMeg)
        bst_process('CallProcess', 'process_snapshot', sAvgSrcMeg(i), [], ...
            'target',         8, ...  % Sources (one time)
            'orient',         4, ...  % bottom
            'time',           0.11, ...
            'threshold',      20, ...
            'Comment',        ['MEG sources: ' sFilesAvg(i).FileName]);
    end
    
    % ===== COMPUTE SOURCES: EEG =====
    % Process: Compute sources [2018]
    sAvgSrcEeg = bst_process('CallProcess', 'process_inverse_2018', sFilesAvg, [], ...
        'output',  1, ...  % Kernel only: shared
        'inverse', struct(...
             'Comment',        'MN: EEG', ...
             'InverseMethod',  'minnorm', ...
             'InverseMeasure', 'amplitude', ...
             'SourceOrient',   {{'fixed'}}, ...
             'Loose',          0.2, ...
             'UseDepth',       1, ...
             'WeightExp',      0.5, ...
             'WeightLimit',    10, ...
             'NoiseMethod',    'reg', ...
             'NoiseReg',       0.1, ...
             'SnrMethod',      'fixed', ...
             'SnrRms',         1e-06, ...
             'SnrFixed',       3, ...
             'ComputeKernel',  1, ...
             'DataTypes',      {{'EEG'}}));
    % Process: Snapshot: Sources (one time)   -   Loop only to get a correct comment for the report
    for i = 1:length(sAvgSrcEeg)
        bst_process('CallProcess', 'process_snapshot', sAvgSrcEeg(i), [], ...
            'target',         8, ...  % Sources (one time)
            'orient',         4, ...  % bottom
            'time',           0.11, ...
            'threshold',      10, ...
            'Comment',        ['EEG sources: ' sFilesAvg(i).FileName]);
    end
end
% Save report
ReportFile = bst_report('Save', []);
if ~isempty(reports_dir) && ~isempty(ReportFile)
    bst_report('Export', ReportFile, bst_fullfile(reports_dir, ['report_' ProtocolName '_sources.html']));
end


%% ===== TIME-FREQUENCY =====
% Start a new report (one report for the time-frequency of all the subjects)
bst_report('Start');
% List of conditions to process separately
AllConditions = {'Famous', 'Scrambled', 'Unfamiliar'};
% Channels to display in the screen capture, by order of preference (if the first channel is bad, use the following)
SelChannel = {'EEG070','EEG060','EEG065','EEG050','EEG003'};
% Compute one separate time-frequency average for each subject/run/condition
for iSubj = 1:length(SubjectNames)
    disp(sprintf('\n===== TIME-FREQUENCY: SUBJECT #%d =====\n', iSubj));
    for iRun = 1:6
        % Process: Select data files in: Subject/Run
        sTrialsAll = bst_process('CallProcess', 'process_select_files_data', [], [], ...
            'subjectname',   SubjectNames{iSubj}, ...
            'condition',     sprintf('sub-%02d_ses-meg_task-facerecognition_run-%02d_proc-sss_meg_notch', iSubj, iRun));
        % Loop on the conditions
        for iCond = 1:length(AllConditions)
            % Comment describing this average
            strComment = [SubjectNames{iSubj}, ' / ', sprintf('run_%02d', iRun), ' / ', AllConditions{iCond}];
            disp(['BST> ' strComment]);
            % Find the first good channel in the display list
            if isempty(BadChannels{iSubj}{iRun})
                iSel = 1;
            else 
                iSel = find(~ismember(SelChannel,BadChannels{iSubj}{iRun}), 1);
            end
            % Process: Select file comments with tag: Avg
            sTrialsCond = bst_process('CallProcess', 'process_select_tag', sTrialsAll, [], ...
                'tag',    [AllConditions{iCond}, '_trial'], ...
                'search', 1, ...  % Search the file names
                'select', 1);  % Select only the files with the tag
            % Process: Time-frequency (Morlet wavelets), averaged across trials
            sTimefreq = bst_process('CallProcess', 'process_timefreq', sTrialsCond, [], ...
                'sensortypes', 'MEG MAG, EEG', ...
                'edit',        struct(...
                     'Comment',         ['Avg: ' AllConditions{iCond} ', Power, 6-60Hz'], ...
                     'TimeBands',       [], ...
                     'Freqs',           [6, 6.8, 7.6, 8.6, 9.7, 11, 12.4, 14, 15.8, 17.9, 20.2, 22.8, 25.7, 29, 32.7, 37, 41.7, 47.1, 53.2, 60], ...
                     'MorletFc',        1, ...
                     'MorletFwhmTc',    3, ...
                     'ClusterFuncTime', 'none', ...
                     'Measure',         'power', ...
                     'Output',          'average', ...
                     'RemoveEvoked',    0, ...
                     'SaveKernel',      0), ...
                'normalize',   'none');  % None: Save non-standardized time-frequency maps
            % Process: Extract time: [-200ms,900ms]
            sTimefreq = bst_process('CallProcess', 'process_extract_time', sTimefreq, [], ...
                'timewindow', [-0.2, 0.9], ...
                'overwrite',  1);
            % Screen capture of one sensor
            hFigTf = view_timefreq(sTimefreq.FileName, 'SingleSensor', SelChannel{iSel});
            bst_report('Snapshot', hFigTf, strComment, 'Time-frequency', [200, 200, 400, 250]);
            close(hFigTf);
        end
    end
end
% Save report
ReportFile = bst_report('Save', []);
if ~isempty(reports_dir) && ~isempty(ReportFile)
    bst_report('Export', ReportFile, bst_fullfile(reports_dir, ['report_' ProtocolName '_timefreq.html']));
end
end




%% ===== SUPPORT FUNCTIONS =====
function BadSeg = GetBadSegments(iSubj, iRun)
     BadSegments{1} = {...
         [247.867 248.185; 598.999 598.999; 598.999 598.999; 611.999 611.999; 612.999 612.999; 613.999 613.999; 616.999 616.999; 617.999 617.999; 623.999 623.999; 715.209 715.467], ...
         [84.791 85.166], ...
         [79.183 80.167], ...
         [64.309 65.185], ...
         [90.958 91.167; 178.005 178.355; 293.282 295.919; 312.298 316.479; 353.835 357.716], ...
         [60.292 66.802; 69.975 71.210; 105.233 107.586; 108.822 109.506; 376.225 376.325]};
     BadSegments{2} = {...
         [223.806 224.199; 279.772 279.895; 453.241 455.108; 692.423 692.593], ...
         [65.298 66.194; 304.727 306.178; 399.165 400.732], ...
         [203.141 205.085; 281.579 287.883; 420.395 421.128], ...
         [387.118 388.229; 440.318 441.900; 554.825 558.744], ...
         [71.000 80.999; 82.750 87.367; 149.528 149.667; 264.747 267.995; 368.415 371.973; 376.263 378.763; 398.334 401.551; 537.410 541.645], ...
         [38.000 47.999; 47.825 50.046; 61.298 61.384; 249.653 253.379; 282.917 283.820; 286.135 287.616; 298.167 300.196; 328.254 329.511; 335.957 337.817; 478.277 480.707]};
     BadSegments{3} = {...
         [406.312 407.207; 727.055 728.714], ...
         [84.894 85.156; 152.028 152.946; 297.835 298.915; 418.272 421.845; 554.084 554.794], ...
         [73.758 74.159; 378.212 378.536; 406.065 407.099; 470.541 471.698; 488.900 491.168; 529.596 530.453], ...
         [94.874 95.152; 317.385 321.374; 325.696 327.055; 439.220 439.829; 454.473 455.175; 486.196 486.829; 518.660 522.015; 524.400 525.249; 562.417 570.325], ...
         [96.208 97.181; 98.942 99.096; 135.005 135.754; 143.990 144.599; 250.139 250.247; 300.459 300.559; 338.265 339.322; 545.913 546.067], ...
         [91.415 92.156; 284.843 286.525; 297.886 298.404; 317.046 317.163; 332.698 332.791; 358.946 359.402; 428.405 428.775; 478.374 478.690; 549.866 550.128]};
     BadSegments{4} = {...
         [22.967 22.967; 50.036 50.098; 52.058 52.058; 156.653 156.653; 171.565 173.386; 239.544 242.105; 268.162 270.175; 268.992 268.992; 316.032 316.032; 338.283 339.000; 357.959 361.909; 370.871 370.871; 381.579 383.677; 437.731 437.731; 463.482 468.505; 476.135 479.838; 486.652 488.272; 504.860 508.999], ...
         [309.493 311.707; 342.681 344.525; 354.019 357.321; 390.023 391.225; 393.926 395.855; 404.221 405.069; 432.522 435.932; 459.048 460.715; 471.763 478.529; 549.387 551.999; 591.087 594.143; 608.541 611.079; 624.847 626.615; 649.648 651.570], ...
         [57.411 58.198; 88.346 88.955; 200.761 202.335; 227.016 227.688; 257.726 258.054; 356.798 359.005; 404.260 411.003], ...
         [46.000 54.823; 61.000 70.332; 203.005 207.125; 275.875 278.121; 313.500 314.824; 337.973 338.636; 422.505 426.239], ...
         [58.000 62.479; 78.250 85.166; 89.955 91.360; 116.322 117.888; 130.013 131.987; 149.509 150.489; 174.650 175.823; 182.030 183.334; 196.758 197.384; 204.458 204.697; 205.236 208.663; 311.028 316.383; 320.700 327.181; 332.437 335.354; 344.205 346.133; 374.208 374.865; 385.519 386.214; 441.942 444.241; 453.957 456.997; 486.039 487.004; 501.238 504.185; 512.962 514.675; 553.398 556.215], ...
         [41.406 45.743; 58.681 59.144; 108.086 108.896; 140.633 143.750; 196.110 199.474; 210.778 210.971; 234.649 235.143; 258.081 259.632; 339.101 340.805; 390.277 390.609; 438.935 442.122; 528.221 534.031]};
     BadSegments{5} = {...
         [265.539 265.778; 266.334 266.495; 268.479 268.965; 367.428 367.636; 439.655 442.779; 453.497 453.853; 504.997 505.329; 519.513 519.683; 595.674 595.982; 602.000 602.463], ...
         [121.113 121.499; 124.971 126.213; 253.735 254.075; 272.232 272.464; 272.895 273.104; 346.368 346.645; 368.812 369.052; 406.382 406.605; 452.920 453.113; 454.903 455.112; 507.655 507.840; 508.766 509.013; 584.853 585.030; 594.831 595.656; 597.261 602.249], ...
         [37.251 37.497; 38.825 39.056; 40.615 41.849; 43.624 44.758; 53.333 53.641; 54.698 55.076; 57.668 59.196; 79.129 79.360; 81.475 81.714; 122.658 123.375; 284.787 285.296; 288.754 288.993; 345.790 346.022; 421.212 421.459; 481.428 482.207; 503.408 503.685; 504.272 504.449; 524.451 524.714; 526.913 531.499], ...
         [87.322 88.178; 91.085 91.325; 95.121 95.491; 114.174 114.397; 129.874 130.113; 151.220 151.544; 281.689 281.959; 532.966 533.345], ...
         [59.176 60.218; 74.854 75.317; 308.180 309.877; 380.705 381.059], ...
         [182.382 183.245; 196.220 196.736; 276.018 276.327; 292.490 294.086; 370.755 370.847; 435.644 436.624; 467.535 468.460; 522.838 525.847]};
     BadSegments{6} = {...
         [141.690 142.424; 157.070 157.417; 355.138 356.025; 423.999 423.999; 424.999 424.999; 426.999 426.999; 427.999 427.999; 486.430 488.151; 493.999 493.999; 501.511 501.619; 501.549 501.549; 501.585 501.585; 502.999 502.999; 503.999 503.999; 540.999 540.999; 541.999 541.999; 555.999 555.999; 556.999 556.999; 561.999 561.999; 563.999 563.999; 564.999 564.999; 565.999 565.999; 567.999 567.999], ...
         [64.898 65.161; 71.700 72.718; 226.185 226.740; 324.124 324.425; 329.062 329.301; 486.143 486.975], ...
         [62.300 63.148; 266.254 266.639; 409.920 410.221], ...
         [54.048 55.214; 330.893 331.255], ...
         [185.616 186.495; 331.411 331.796; 386.843 387.028; 387.999 387.999; 389.999 389.999; 434.575 434.875; 519.802 519.995], ...
         [44.720 45.167; 211.446 211.964; 368.955 369.172]};
     BadSegments{7} = {...
         [154.966 155.144; 551.639 551.855], ...
         [88.774 89.167; 107.999 107.999; 109.999 109.999; 110.999 110.999; 112.999 112.999; 113.999 113.999; 114.999 114.999; 119.999 119.999; 121.999 121.999; 124.223 125.465; 137.011 138.871], ...
         [81.627 82.136; 377.953 381.046], ...
         [241.136 242.171; 543.849 544.196; 596.639 598.553; 600.227 601.075; 603.999 603.999; 605.999 605.999; 609.999 609.999; 611.305 612.809; 614.999 614.999; 615.803 616.937; 623.694 625.877; 653.999 653.999; 655.055 655.756; 663.999 663.999], ...
         [68.852 69.481; 74.034 75.200; 78.426 78.600; 104.497 105.963; 253.951 254.034; 256.964 257.038; 257.915 258.048; 323.254 324.156; 365.880 368.131; 369.952 370.060; 371.728 372.999; 430.931 431.965; 535.521 542.999], ...
         [70.271 71.205; 94.441 96.445; 98.613 99.126; 112.318 112.749; 131.686 132.265; 148.935 150.615; 161.120 161.600; 205.325 208.214; 215.035 215.863; 217.403 218.818; 286.178 287.171; 399.075 404.853]};
     BadSegments{8} = {...
         [238.505 238.546; 256.354 257.224; 316.116 316.167; 341.519 341.558; 356.493 356.566; 380.095 380.170; 391.906 392.048; 448.457 448.562; 469.931 470.028; 530.488 530.575; 555.180 558.136; 562.152 562.245; 588.403 588.502; 625.205 625.662; 638.019 638.438; 649.982 650.008; 650.691 651.357; 651.925 652.061; 665.445 665.472; 695.923 696.015; 706.528 706.720; 729.706 732.534], ...
         [99.546 106.114; 113.245 114.199; 150.885 154.989; 277.486 278.075; 333.645 335.574; 339.358 340.646; 371.860 375.394; 497.459 499.156], ...
         [49.008 50.205; 231.446 233.406; 329.590 329.659; 355.101 356.019; 360.733 360.973; 372.955 374.891; 389.283 392.068; 453.610 455.431; 464.632 465.265; 489.209 489.996; 514.777 515.295], ...
         [178.368 179.232; 293.865 294.521; 418.252 418.314; 450.124 450.209; 480.236 480.445; 492.725 493.975; 495.170 497.624; 500.382 500.459; 504.247 504.402; 628.017 628.079; 628.827 628.905; 630.209 630.332], ...
         [100.616 101.172; 107.691 108.539; 188.814 188.875; 193.119 193.281; 207.964 208.033; 432.254 432.385; 489.834 489.911; 517.960 518.037; 518.955 519.040; 520.938 521.062; 521.945 522.037; 523.959 524.052; 525.942 526.057; 526.945 527.015; 528.958 529.059; 531.138 531.192; 531.979 532.048; 532.951 533.028; 533.962 534.024], ...
         [135.997 137.232; 383.565 383.657; 418.763 418.955]};
     BadSegments{9} = {...
         [215.107 216.187; 262.388 262.388; 287.519 287.635; 289.895 290.135; 311.999 311.999; 350.161 351.179; 526.154 527.033; 564.999 564.999; 584.935 585.059; 587.999 587.999; 601.999 601.999; 603.999 603.999; 608.999 608.999; 612.999 612.999], ...
         [47.667 47.775; 49.172 54.681; 67.195 70.805; 78.988 80.477; 138.701 138.948; 138.727 138.727; 138.728 138.728; 138.728 138.728; 138.734 138.734; 138.734 138.734; 138.881 138.881; 138.881 138.881; 138.907 138.907; 140.993 141.093; 141.037 141.037; 141.045 141.045; 155.795 155.864; 155.822 155.822; 155.849 155.849; 168.999 168.999; 206.999 206.999; 219.999 219.999; 225.999 225.999; 226.999 226.999; 228.999 228.999; 236.999 236.999; 242.463 242.463; 242.463 242.463; 242.487 242.487; 247.999 247.999; 251.999 251.999; 252.483 252.483; 253.315 254.163; 265.999 265.999; 267.999 267.999; 272.358 272.358; 272.410 272.410; 298.999 298.999; 300.999 300.999; 314.037 314.037; 314.047 314.047; 321.813 321.813; 321.813 321.813; 321.833 321.833; 329.117 329.117; 329.117 329.117; 329.156 329.156; 346.999 346.999; 347.999 347.999; 349.320 349.320; 349.329 349.329; 352.528 355.869; 364.040 364.040; 396.278 397.420; 404.865 404.865; 407.905 407.905; 407.905 407.905; 407.954 407.954; 418.454 418.454; 418.454 418.454; 418.486 418.486; 441.999 441.999; 444.999 444.999; 447.999 447.999; 453.550 453.650; 454.931 455.055; 457.999 457.999; 479.964 480.079; 481.999 481.999; 482.949 483.073; 488.948 489.064; 511.999 511.999; 520.999 520.999; 523.999 523.999; 526.954 527.069; 533.999 533.999; 537.999 537.999], ...
         [82.889 83.198; 206.143 207.424; 403.873 404.096; 406.790 407.485; 413.936 414.075; 415.912 416.112; 536.621 537.532], ...
         [182.999 182.999; 182.999 182.999; 183.999 183.999; 195.999 195.999; 208.999 208.999; 209.999 209.999; 278.601 278.955; 413.913 414.067; 415.250 417.792; 419.940 420.087; 420.999 420.999; 512.999 512.999; 514.363 516.254; 521.917 522.134; 522.999 522.999; 523.915 524.101; 533.999 533.999; 538.999 538.999; 539.892 540.062; 543.999 543.999; 548.874 549.089; 549.900 550.062; 552.755 554.075; 585.999 585.999; 587.957 588.073; 588.999 588.999; 594.999 594.999; 603.971 604.056; 604.928 605.097; 615.986 618.495; 623.999 623.999], ...
         [53.215 54.203; 164.227 168.965; 201.433 202.775; 292.973 295.889; 303.961 304.817; 309.354 311.236; 313.123 313.639; 322.547 323.017; 331.141 334.852; 356.637 356.985; 367.855 368.079; 369.995 373.459; 377.849 378.142; 406.072 407.306; 436.164 436.665; 458.033 458.905; 516.889 518.656; 517.812 518.684], ...
         [113.945 115.225; 198.570 198.863; 264.162 264.795; 383.705 385.641; 396.477 397.064; 399.706 406.457; 452.261 453.179; 486.338 487.102; 498.869 499.579; 507.968 508.925; 546.266 547.285; 558.825 560.128]};
     BadSegments{10} = {...
         [235.104 236.184; 324.272 325.028; 330.799 331.401; 541.062 541.826; 564.747 565.072], ...
         [100.581 101.229; 285.644 285.644; 285.644 285.644; 297.979 298.033; 298.999 298.999; 300.949 301.026; 301.999 301.999; 303.999 303.999; 304.999 304.999; 306.999 306.999; 307.999 307.999; 310.999 310.999; 311.999 311.999; 313.971 314.025; 314.999 314.999; 316.995 317.035; 317.999 317.999; 319.999 319.999; 320.923 321.046; 326.971 327.048; 327.999 327.999; 329.999 329.999; 330.999 330.999; 332.974 333.028; 333.999 333.999; 335.971 336.025; 336.999 336.999; 338.973 339.035; 339.999 339.999; 341.951 342.044; 343.999 343.999; 344.964 345.033; 346.999 346.999; 347.999 347.999; 452.015 453.411; 458.735 459.776; 467.974 468.089], ...
         [41.266 41.675; 53.905 55.240; 159.159 159.345; 220.255 220.394], ...
         [66.751 67.190; 201.674 201.812; 294.119 295.168; 303.188 303.528; 316.992 317.037; 317.999 317.999; 317.999 317.999; 319.999 319.999; 320.980 321.042; 322.999 322.999; 325.999 325.999; 326.955 327.048; 328.999 328.999; 333.999 333.999; 334.999 334.999; 342.999 342.999; 343.999 343.999; 351.915 352.038; 352.999 352.999; 407.979 408.056; 409.999 409.999; 411.999 411.999; 415.405 416.154; 435.999 435.999; 436.962 437.046; 468.885 469.100; 469.999 469.999; 470.999 470.999; 516.210 516.357], ...
         [105.369 106.172; 141.468 142.178; 199.008 199.818; 201.269 201.716; 352.851 353.083; 449.865 452.041; 459.508 459.802], ...
         [229.948 230.033; 230.999 230.999; 230.999 230.999; 259.588 259.990; 343.970 345.798; 361.999 361.999; 362.970 363.046; 364.999 364.999; 366.999 366.999; 367.962 368.031; 372.521 374.179; 405.999 405.999; 409.999 409.999; 411.999 411.999; 412.999 412.999; 434.999 434.999; 435.970 436.031; 516.999 516.999; 519.999 519.999]};
     BadSegments{11} = {...
         [179.045 180.225; 323.999 323.999; 323.999 323.999; 324.999 324.999; 327.978 328.025; 328.966 329.028; 330.999 330.999; 365.999 365.999; 367.999 367.999; 370.999 370.999; 371.999 371.999; 373.999 373.999; 375.965 376.042; 377.999 377.999], ...
         [52.107 53.156; 521.077 521.155], ...
         [65.331 66.156], ...
         [81.579 82.143; 108.565 108.565; 108.566 108.566; 112.176 112.176; 122.377 122.377; 122.565 122.565; 213.882 213.882; 215.305 215.305; 224.851 224.851; 224.919 224.919; 255.815 255.815; 257.893 257.893; 359.952 359.952; 361.630 361.754; 370.285 370.285; 376.853 376.853; 511.600 511.600; 513.631 513.631; 513.988 513.988; 518.215 518.215], ...
         [63.205 64.154; 170.951 171.036; 176.841 176.841; 176.842 176.842; 177.282 177.282; 223.886 223.971; 259.963 259.963; 261.547 261.547; 341.185 341.302; 368.999 368.999; 370.999 370.999; 374.999 374.999; 382.971 383.033; 383.999 383.999; 386.999 386.999; 388.958 389.035; 390.999 390.999; 391.955 392.031; 394.955 395.040; 396.999 396.999; 398.999 398.999; 400.999 400.999; 402.999 402.999; 404.963 405.033; 406.999 406.999; 429.829 429.829; 430.346 430.346; 465.184 465.184; 470.290 470.290], ...
         [50.195 51.145; 158.850 159.012]};
     BadSegments{12} = {...
         [193.435 194.175; 500.000 501.000; 08.988 509.868; 528.999 528.999; 528.999 528.999; 529.999 529.999; 531.999 531.999; 547.999 547.999], ...
         [133.484 134.185; 134.911 135.065; 553.999 553.999; 553.999 553.999; 554.999 554.999; 557.999 557.999; 558.999 558.999; 564.999 564.999; 565.977 566.046; 568.999 568.999; 569.951 570.028; 571.999 571.999; 579.959 580.028; 580.977 581.055; 583.971 584.033; 583.999 583.999; 585.999 585.999; 586.999 586.999; 588.980 589.026; 590.999 590.999; 591.958 592.044; 595.999 595.999], ...
         [46.028 47.216; 129.557 130.236; 200.999 200.999; 200.999 200.999; 201.975 202.075; 203.999 203.999; 204.967 205.053; 213.627 214.515; 218.958 219.044; 222.959 228.252; 309.999 309.999; 311.897 317.537; 311.999 311.999; 351.968 353.234; 399.999 399.999; 446.508 451.454; 487.999 487.999; 512.278 512.926], ...
         [82.535 84.062; 102.247 102.355; 134.999 134.999; 134.999 134.999; 136.999 136.999; 138.973 139.042; 147.999 147.999; 148.999 148.999; 149.999 149.999; 193.999 193.999; 194.999 194.999; 224.984 226.095; 241.954 242.077; 243.948 244.087; 280.999 280.999; 281.999 281.999; 305.000 305.451; 310.945 311.068; 328.999 328.999; 352.999 352.999; 353.999 353.999; 397.025 397.411; 415.944 418.096; 415.999 415.999; 470.999 470.999; 477.971 478.041; 483.997 484.737; 492.958 493.089; 493.999 493.999; 494.999 494.999; 522.999 522.999; 523.999 523.999; 524.955 525.039], ...
         [175.183 176.155; 246.112 246.991; 412.005 413.340; 483.915 484.061; 485.959 486.028; 497.360 498.880; 594.457 594.944; 596.641 598.006; 615.277 618.710], ...
         [66.900 72.232; 75.510 78.688; 543.495 545.999]};
     BadSegments{13} = {...
         [307.246 308.179; 627.881 628.089; 629.941 630.049], ...
         [171.506 172.301; 172.999 172.999; 430.999 430.999; 432.999 432.999; 434.999 434.999; 448.999 448.999; 479.999 479.999; 489.999 489.999; 490.999 490.999; 491.999 491.999], ...
         [92.059 93.209], ...
         [52.791 53.231; 54.999 54.999; 65.985 68.138; 91.240 91.386; 92.137 92.207; 105.346 105.493; 121.120 121.398; 146.546 146.955; 194.025 197.652; 242.985 243.571; 247.101 247.556; 269.889 270.083; 270.946 271.070; 274.408 275.866; 294.565 295.267; 319.370 319.879; 365.999 365.999; 375.913 376.044; 376.869 377.055; 377.999 377.999; 390.130 393.225; 403.965 404.035; 404.953 405.075; 419.265 419.565; 427.015 427.154; 427.879 428.118; 427.999 427.999; 428.913 429.067; 480.945 482.095; 484.285 484.571; 516.929 517.099; 517.963 518.079], ...
         [115.921 116.060; 117.999 117.999; 120.999 120.999; 126.896 127.058; 130.955 131.039; 131.957 132.043; 132.961 133.015; 134.999 134.999; 135.950 136.043; 139.999 139.999; 143.999 143.999; 144.850 145.066; 145.999 145.999; 146.999 146.999; 160.929 161.045; 161.999 161.999; 189.086 189.310; 193.939 194.063; 195.975 196.045; 199.961 200.045; 201.952 202.059; 202.939 205.084; 210.999 210.999; 211.738 213.629; 217.954 218.061; 218.987 219.049; 219.975 220.028; 258.953 259.037; 259.952 260.045; 260.963 261.032; 261.958 262.043; 262.954 263.054; 262.999 262.999; 288.999 288.999; 295.477 295.655; 298.976 299.054; 299.979 300.018; 300.975 301.005; 336.952 337.044; 337.977 338.032; 338.988 339.027; 338.999 338.999; 339.991 340.045; 340.975 341.029; 341.978 342.016; 343.976 344.054; 361.961 362.045; 365.942 366.050; 370.999 370.999; 384.961 385.045; 409.977 410.047; 415.984 416.030; 426.944 427.059; 430.999 430.999; 436.999 436.999; 448.999 448.999; 449.500 452.131; 481.962 482.055; 482.973 483.027; 484.975 485.060; 503.999 503.999; 504.999 504.999; 507.967 508.028; 508.955 509.032; 510.999 510.999; 512.999 512.999; 518.865 520.292; 522.958 523.043; 523.999 523.999; 524.948 525.034; 525.959 526.044; 526.946 527.055; 527.942 528.027; 528.961 529.053; 529.945 530.022; 535.999 535.999; 536.957 537.058; 538.967 539.036; 539.963 540.040; 541.976 542.045; 543.999 543.999; 551.956 552.034; 552.990 553.028; 553.977 554.024; 555.968 556.022; 556.204 557.454; 557.963 558.032; 559.969 560.023; 563.973 564.035; 563.999 563.999], ...
         [87.308 87.639; 107.039 108.189; 427.943 428.035; 437.965 438.027; 439.973 440.026; 491.275 492.571]};
     BadSegments{14} = {...
         [365.982 367.194; 368.999 368.999; 368.999 368.999; 369.999 369.999; 371.999 371.999; 373.999 373.999; 375.942 376.050; 376.999 376.999; 378.999 378.999; 380.999 380.999; 382.999 382.999; 383.999 383.999; 386.999 386.999; 387.999 387.999; 418.999 418.999; 444.999 444.999; 690.178 693.218], ...
         [101.796 102.183; 218.999 218.999; 219.999 219.999; 221.999 221.999; 222.999 222.999; 227.999 227.999; 229.999 229.999; 230.999 230.999; 232.999 232.999; 233.999 233.999; 235.999 235.999; 237.999 237.999; 238.985 239.031; 240.999 240.999; 242.999 242.999; 243.999 243.999; 246.999 246.999; 247.999 247.999; 253.999 253.999; 258.999 258.999; 259.999 259.999; 260.227 261.123; 264.999 264.999; 265.999 265.999; 268.999 268.999; 269.999 269.999; 274.999 274.999; 276.999 276.999; 277.999 277.999; 280.999 280.999; 282.966 283.021; 284.999 284.999; 285.999 285.999; 360.999 360.999; 362.973 363.026; 364.999 364.999; 365.999 365.999; 367.999 367.999; 368.985 369.024; 370.999 370.999; 371.999 371.999; 378.999 378.999; 379.999 379.999; 381.999 381.999; 383.999 383.999; 385.999 385.999; 387.999 387.999; 388.999 388.999; 398.999 398.999; 415.999 415.999; 417.999 417.999; 418.999 418.999; 421.999 421.999; 422.999 422.999], ...
         [89.285 90.187], ...
         [85.991 87.179; 88.999 88.999; 89.999 89.999; 92.966 93.027; 93.999 93.999; 107.999 107.999; 110.999 110.999; 111.971 112.025; 118.999 118.999; 121.999 121.999; 122.999 122.999; 166.999 166.999; 168.957 169.035; 173.999 173.999; 174.999 174.999; 175.957 176.035; 198.999 198.999; 199.999 199.999; 202.999 202.999; 203.999 203.999; 205.999 205.999; 207.999 207.999; 208.999 208.999; 210.974 211.044; 211.999 211.999; 212.999 212.999; 232.984 233.023; 234.999 234.999; 235.999 235.999; 236.999 236.999; 239.958 240.044; 240.999 240.999], ...
         [71.132 72.227], ...
         [98.134 98.827; 110.497 110.814; 114.000 117.272; 279.999 279.999; 279.999 279.999; 280.999 280.999; 282.999 282.999; 283.999 283.999; 284.999 284.999; 286.999 286.999; 287.999 287.999; 288.999 288.999; 289.999 289.999; 291.966 292.044; 292.999 292.999; 346.999 346.999; 437.608 437.646; 518.988 519.104; 522.969 523.015; 534.999 534.999; 536.967 537.029; 562.725 563.426]};
     BadSegments{15} = {...
         [66.755 67.172; 257.988 258.042; 258.714 259.385; 260.958 261.044; 265.931 266.062; 267.968 268.022; 268.978 269.033; 270.992 271.031; 272.964 273.017; 273.968 274.030; 275.975 276.059; 276.970 277.046; 288.976 289.038; 289.984 290.015; 295.994 296.025; 296.982 297.013; 300.966 301.152; 364.730 364.800], ...
         [46.531 46.555; 65.156 66.148; 87.535 87.642; 183.896 183.934; 294.342 294.375; 330.004 330.528; 333.785 333.815; 345.177 345.203; 399.544 399.573; 480.964 480.980], ...
         [53.557 54.175; 148.851 149.136], ...
         [43.681 44.214; 361.255 361.587; 409.462 411.021], ...
         [90.868 91.130; 410.935 411.044], ...
         [64.786 65.141; 132.385 132.786; 177.735 178.252; 278.747 278.902; 313.959 314.028; 314.970 315.031]};
     BadSegments{16} = {...
         [38.958 44.335; 305.280 312.826; 324.778 326.067; 393.183 398.367; 403.410 422.854], ...
         [49.726 50.251; 58.675 65.465; 263.918 264.088; 321.795 322.195; 522.881 523.081], ...
         [61.245 62.325; 65.628 65.905; 315.686 317.275; 335.184 336.449; 387.870 388.117; 389.892 390.062], ...
         [56.020 57.255; 60.711 61.096; 64.353 64.924; 73.202 73.757; 76.134 76.844; 95.625 96.859; 171.436 172.625; 178.359 178.715; 212.828 213.415; 352.425 352.617; 367.755 369.097; 488.400 489.426], ...
         [71.665 72.513; 76.347 77.181; 127.134 128.862], ...
         [46.190 47.209; 53.461 55.752; 194.317 194.510; 213.101 213.240; 216.356 216.442; 225.962 226.286; 245.163 245.464; 252.041 253.422; 254.648 254.856; 292.966 294.340; 340.928 341.067; 346.298 346.452]};
    BadSeg = BadSegments{iSubj}{iRun}';
end


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.

-  ⇤ ← Revision 1 as of 2015-02-03 16:36:35 → 
  Size: 296
  Editor: FrancoisTadel
  Comment:
+   ← Revision 55 as of 2016-07-21 22:33:46 → ⇥
  Size: 29227
  Editor: FrancoisTadel
  Comment:
-Deletions are marked like this.
+Additions are marked like this.
 Line 1:
-= Tutorial 20: Advanced scripting =
''Authors: Francois Tadel, Elizabeth Bock''
+= 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.

Requirements: You need a license for the Matlab environment in order to use these tools and execute custom scripts. If you are running the compiled version of Brainstorm with the MCR library, the only custom code you can run is through the menu File > Matlab console and the process "Run Matlab command".
-Line 6:
+Line 12:
-<<EmbedContent("http://neuroimage.usc.edu/bst/get_prevnext.php?prev=Tutorials/Connectivity")>>
+== Starting a new script ==
The easiest way to get started with a new Brainstorm script is to use the script generator, already introduced in the tutorial [[http://neuroimage.usc.edu/brainstorm/Tutorials/PipelineEditor#Saving_a_pipeline|Select files and run processes]]. Select some files in the Process1 or Process2 tabs, select a list of processes, and use the menu '''Generate .m script'''. The example below should work on the the protocol "TutorialIntroduction" created during the introduction tutorials.

 * In the Process1 tab, leave the selection box empty and click on [Run]. Instead of selecting the files from the Brainstorm interface, we will select them directly from the database using a script.
 * Select process '''File > Select files: Recordings''' (do not execute immediately)<<BR>>Subject='''Subject01''', Condition='''[Empty]''', File comment='''Avg: deviant''' (the space is important). <<BR>><<BR>> {{attachment:start1.gif||height="334",width="400"}}
 * This process selects all the recordings with a comment including the string "Avg: deviant", from all the folders in Subject01 (except "Intra-subject" and "Common files"). We expect to get two files: the averages of the deviant condition for both runs.
 * Add process '''Pre-process > Band-pass filter''': Lower cutoff='''0Hz''', Upper cutoff='''30Hz''', Mirror. <<BR>> '''Add process File > Save snapshot''': Recordings time series, Sensors=MEG. <<BR>> <<BR>> {{attachment:start2.gif||height="294",width="518"}}
 * This will apply a low-pass filter at 30Hz and save a screen capture of the signals in the report.
 * Do not run the pipeline, select the menu '''Generate .m script '''instead. It saves a new .m file and opens it in the Matlab editor. Close the pipeline editor window and look at the script.<<BR>><<BR>> {{attachment:start3.gif||height="272",width="673"}}
 * The script you just generated can be the starting point to your own  custom script. The following sections explain line by line how they work  and how to edit them.

== Line by line: Header ==
{{{
% Script generated by Brainstorm (19-Jul-2016)
}}}
All the lines starting with a "%" are comments, they are never executed.

{{{
% Input files
sFiles = [];
SubjectNames = {...
    'Subject01'};
}}}
These lines define the script inputs:

 * '''sFiles''':  The list of files in input. Currently empty because we did not select  anything in input in the Process1 list. If we had selected files, it  would contain a cell array of strings with relative file paths.
 * '''SubjectNames''': List of subject names that are used in the script. Most of the times, the generated script would contain only one entry, but it written as a cell array so that it is easier to extend it to multiple subjects with a loop (described further in this tutorial).

{{{
% Start a new report
bst_report('Start', sFiles);
}}}
Start a new report of activity: Clears all the previous logs and gets ready to record new messages. The report will collect all the messages that are generated during the execution of the script by the various processes. You can explicitely add screen captures and additional messages to the current report with the function bst_report. This report will remain open until the function bst_report('Start') is called again. To display the current report, use the menu [[http://neuroimage.usc.edu/brainstorm/Tutorials/PipelineEditor#Report_viewer|File > Report viewer]].

The syntax '''function_name'''(''''SubFunction'''', arguments) is used a lot in Brainstorm, to call a subfunction available inside a .m file. This line above calls the function Report() in the file brainstorm3/toolbox/process/bst_report.m. This is made possible with the use of the short script "macro_methodcall" you will see at the top of many files. Many of the Brainstorm .m files are actually libraries of functions, rather than simple "scripts" or "functions".

== Line by line: Body ==
{{{
% Process: Select data files in: Subject01/*/Avg: deviant
sFiles = bst_process('CallProcess', 'process_select_files_data', sFiles, [], ...
    'subjectname',   SubjectNames{1}, ...
    'condition',     '', ...
    'tag',           'Avg: deviant', ...
    'includebad',    0, ...
    'includeintra',  0, ...
    'includecommon', 0);

% Process: Low-pass:30Hz
sFiles = bst_process('CallProcess', 'process_bandpass', sFiles, [], ...
    'highpass',    0, ...
    'lowpass',     30, ...
    'mirror',      1, ...
    'sensortypes', 'MEG', ...
    'overwrite',   0);

% Process: Snapshot: Recordings time series
sFiles = bst_process('CallProcess', 'process_snapshot', sFiles, [], ...
    'target',         5, ...  % Recordings time series
    'modality',       1, ...  % MEG (All)
    'orient',         4, ...  % bottom
    'time',           0.11, ...
    'contact_time',   [0, 0.1], ...
    'contact_nimage', 12, ...
    'threshold',      20, ...
    'Comment',        'Run');
}}}
You will find one block per process you selected in the pipeline editor. They all have the same syntax: <<BR>> output_files = '''bst_process'''('CallProcess', process_name, input_files_A, input_files_B, options_list);

 * '''process_name''':  String indicating the function corresponding to the process to execute.  To know from the pipeline editor what is the path to the process  function: hover your mouse over the selected process, as illustrated in [[http://neuroimage.usc.edu/brainstorm/Tutorials/PipelineEditor#Saving_a_pipeline|this tutorial]].
 * '''input_files_A''':  List of input files in Process1, or FilesA in Process2. It can be a  cell array of files names (full path, or relative path from the protocol  folder), or an array of structures describing the files in the database  (returned by a previous call to bst_process).
 * '''input_files_B''': Empty for Process1, or FilesB in Process2. Cell array of strings or array of struct.
 * '''options_list''': Pairs of (option_name, option_values), one for each option of the process.
 * '''output_files''':  Array of structures describing the files in output of the process. If  the process created new files, this variable contains the new files. If  the process didn't create new files or was modifying exiting files, most  of the time this variable would contain the same files as the input  list.

== Line by line: Footer ==
{{{
% Save and display report
ReportFile = bst_report('Save', sFiles);
}}}
Closes the current report and saves it in the user Brainstorm report folder ($HOME/.brainstorm/reports). These reports are in .mat format and contain all the information necessary to re-run the execution exactly in the same way, but they not easy to read.

The parameter "sFiles" is optional, it indicates what are the files that are considered as the final results of the script. You can remove it without breaking your script: ReportFile = bst_report('Save');

{{{
bst_report('Open', ReportFile);
}}}
Opens the report viewer to display what happened during the execution. This is equivalent to using the menu [[http://neuroimage.usc.edu/brainstorm/Tutorials/PipelineEditor#Report_viewer|File > Report viewer]]. You can comment this line (ie. add a "%" at the beginning of the line) if you don't want to show the report at the end of the execution.

{{{
% bst_report('Export', ReportFile, ExportDir);
}}}
This function exports the report in readable format, as an HTML that includes all the screen captures embedded in it. It is disabled by default. If you want to use this feature: remove the "%" at the beginning of the line, and define the variable ExportDir.

ExportDir must be a string that defines where to save the HTML report. It can be wither the absolute path to a HTML file (eg. 'C:\Users\myuser\Documents\report_example.html') or just a folder (eg. 'C:\Users\myuser\Documents'). If you enter only a path to a folder, a default file name including the protocol name and a date tag is generated (report_ProtocolName_YYMMDD_HHMMSS.html).

== Simplify the calls ==
This script you generated is like any Matlab script: you can edit it, rename the variables, add tests and loops, etc. The first important thing to understand is how to edit the options or change the inputs/outputs of a process. The script generator uses only one variable for all the file lists ('''sFiles''') and the output process is always the input of the following process. This is usually too restrictive to write a full analysis script: we commonly need to have multiple lists of files or to run two different operations on the same file.

Let's consider the first process call, which  selects the averages for the Deviant condition in both runs.

{{{
sFiles = bst_process('CallProcess', 'process_select_files_data', sFiles, [], ...
    'subjectname',   SubjectNames{1}, ...
    'condition',     '', ...
    'tag',           'Avg: deviant', ...
    'includebad',    0, ...
    'includeintra',  0, ...
    'includecommon', 0);
}}}
There is no need to set the parameter sFiles, because there is no input, you can replace it with an empty matrix []. You can therefore remove the line "sFiles = [];". We can also rename the output variable "sAvgData", to be a bit more specific.

{{{
sAvgData = bst_process('CallProcess', 'process_select_files_data', [], [], ...
}}}
You can omit all the options that are not defined, not used, or kept to their default values:

{{{
sAvgData = bst_process('CallProcess', 'process_select_files_data', [], [], ...
    'subjectname',   SubjectNames{1}, ...
    'tag',           'Avg: deviant');
}}}
Edit low-pass filter call: Change the input to sAvgData and the output to sAvgDataLow, this way you will be able to keep track of the two files if you need to use them independently later in the script.

{{{
sAvgDataLow = bst_process('CallProcess', 'process_bandpass', sAvgData, [], ...
    'highpass',    0, ...
    'lowpass',     30, ...
    'sensortypes', 'MEG');
}}}
Edit the call to the snapshot process: Change the input to sAvgDataLow, and remove the output parameter (we are not expecting output from it).

{{{
bst_process('CallProcess', 'process_snapshot', sAvgDataLow, [], ...
    'target',   5, ...  % Recordings time series
    'modality', 1);     % MEG (All)
}}}
Replace the last lines with the following, to export  the report instead of opening in the report viewer (edit the file path  to point at your own user folder instead).

{{{
ReportFile = bst_report('Save');
bst_report('Export', ReportFile, 'C:\Users\franc\Documents\report_test.html');
}}}
== Interaction with Matlab ==
Select the code for the first process in the Matlab editor, right-click > Evaluate selection (or press F9).

{{attachment:edit1.gif}}

If you haven't executed your script yet, you will get the following error in the Matlab command window:

{{{
Undefined variable "SubjectNames" or class "SubjectNames".
}}}
The variable SubjectNames is not defined yet: Execute the first two lines "SubjectNames = {'Subject01'}", then try again. You should now have a new variable in your Matlab workspace, which points at the two average files. Type "sAvgData(1)" in your command window to display the first element:

{{{
>> sAvgData(1)
ans =
          iStudy: 6
           iItem: 1
        FileName: 'Subject01/S01_AEF_20131218_01_600Hz_notch/data_deviant_average_160513_1329.mat'
        FileType: 'data'
         Comment: 'Avg: deviant (39 files)'
       Condition: 'S01_AEF_20131218_01_600Hz_notch'
     SubjectFile: 'Subject01/brainstormsubject.mat'
     SubjectName: 'Subject01'
        DataFile: ''
     ChannelFile: 'Subject01/S01_AEF_20131218_01_600Hz_notch/channel_ctf_acc1.mat'
    ChannelTypes: {'ADC A'  'ADC V'  'DAC'  'ECG'  'EOG'  'FitErr'  'HLU'  'MEG'  'MEG REF' ...}
}}}
The field "sAvgData(1).FileName" contains the relative path to the to the Deviant average in the first run. This structure sAvgData contains also a lot of information that can be helpful in your script:

 * '''iStudy / iItem''': Reference of the file in the database (described later in this tutorial).
 * '''FileType''': 'raw' (continuous files), 'data' (recordings), 'results' (sources), 'timefreq' (time-frequency, spectrum and connectivity), or 'matrix' (any time series extracted from other files).
 * '''Comment''': Comment field of the file (what is displayed in the database explorer)
 * '''Condition''': Name of the condition/folder in which the file is located
 * '''SubjectFile''': Relative path to the subject file (brainstormsubject.mat)
 * '''SubjectName''': Name of the subject
 * '''DataFile''': For types 'results' or 'timefreq', path of the parent file in the database explorer
 * '''ChannelFile''': Relative path to the channel file
 * '''ChannelTypes''': Cell array of channel types available for the input file

== Naming conventions ==
To help you navigate in the Brainstorm code, here are some naming conventions:

 * '''Structures''': Name starting with a "s" followed by a capital letter (eg. sFiles).
 * '''Indices''': Either loop variables or array indices, name starting with a "i" (eg. iSubject, iStudy, iTime).
 * '''Counts''': Number of elements in a group, name starting with a "n" (eg. nAvg, nTrials, nSubjects).
 * '''File names''': Scripts and functions, only lower case, separation with "_" (eg. process_average)
 * '''Sub-functions''': Inside a .m file, name starting with a capital, camel-case (eg. CallProcess, Start).

== Running the script ==
The simplified script looks like this:

{{{
% Input files
SubjectNames = {'Subject01'};
% Start a new report
bst_report('Start');

% Process: Select data files in: Subject01/*/Avg: deviant
sAvgData = bst_process('CallProcess', 'process_select_files_data', [], [], ...
    'subjectname',   SubjectNames{1}, ...
    'tag',           'Avg: deviant');
% Process: Low-pass:30Hz
sAvgDataLow = bst_process('CallProcess', 'process_bandpass', sAvgData, [], ...
    'highpass',    0, ...
    'lowpass',     30, ...
    'sensortypes', 'MEG');
% Process: Snapshot: Recordings time series
bst_process('CallProcess', 'process_snapshot', sAvgDataLow, [], ...
    'target',   5, ...  % Recordings time series
    'modality', 1);     % MEG (All)

% Save and display report
ReportFile = bst_report('Save');
bst_report('Export', ReportFile, 'C:\Users\franc\Documents\report_test.html');
}}}
You have three ways to execute it:

 * Select all the lines (Ctrl+A) and evaluate it in Matlab (F9).
 * In the Editor toolbar of the Matlab environment, click on the button [Run].
 * Save the file, go to this folder with Matlab (or add it to your path) and type the name of the script in the command window (without the ".m" at the end).

At the end of the execution, nothing happens, because we indicated we wanted to just export the report instead of opening it. To check out the report of execution: use the menu '''File > Report viewer''' from the Brainstorm window, or open the '''report_test.html''' that was saved somewhere on your computer.

On this page, you can get review everything that happened in the script: when it was executed, how long it took, what are the processes that were executed, some produced messages (two files were selected with the first process), and at the end you have the screen captures taken by process_snapshot.

 . {{attachment:report1.gif||height="455",width="593"}}

== Running the script again ==
If you execute the script again, it will not behave as expected anymore. The selection process we used assumes that there is only one file per folder with a comment that includes "Avg: deviant". This is not the case anymore after the execution, because the low-pass filtered files also contain the same tags. The execution of the first process of the script now '''returns 4 files'''.

{{{
>> sAvgData = bst_process('CallProcess', 'process_select_files_data', [], [], ...
    'subjectname',   SubjectNames{1}, ...
    'tag',           'Avg: deviant')

sAvgData =
1x4 struct array with fields:
    iStudy
    iItem
    ...
}}}
In order to exclude the low-pass filtered files from this selection, you can add another process that will '''refine the selection'''. Use the script generator again to create a template call for another process, then copy-paste it in your script.

 * In Process1: Select any recordings file (we will not run anything, just generate a script).
 * Select process '''File > Select files: By tag''': Search="low", Search the comment, '''Ignore '''the files. <<BR>><<BR>> {{attachment:select_tag.gif||height="252",width="261"}}
 * Select the menu '''Generate .m script''' (make sure you do not overwrite the script you are currently working on), then close the pipeline editor.
 * Copy-paste and edit the call to process_select_tag to your main script.

Now the file selection part of your script should look like this, and should return only two files:

{{{
% Process: Select data files in: Subject01/*/Avg: Deviant
sAvgData = bst_process('CallProcess', 'process_select_files_data', [], [], ...
    'subjectname',   SubjectNames{1}, ...
    'tag',           'Avg: deviant');
% Process: Ignore file comments with tag: low
sAvgData = bst_process('CallProcess', 'process_select_tag', sAvgData, [], ...
    'tag',    'low', ...
    'search', 2, ...  % Search the file comments
    'select', 2);  % Ignore the files with the tag
}}}
With this last modification, your script is more robust. It can be executed multiple times without completely changing its behavior. When you are fetching files from the database using comment tags or file names, always pay attention to this aspect: the database grows and the further you go, the more specific your requests have to be.

A good practice can be to tag explicitely the output files your script generates if you need to fetch them later. You can use the process '''File > Add tag''' and '''File > Set comment'''.

 . {{attachment:add_tag.gif||height="203",width="456"}}

== Starting Brainstorm ==
Brainstorm must be running in the background for these scripts to run properly. The interface doesn't have to be visible on the screen, but the database engine must be running for processing requests. At the beginning of your script, you can explicitely start or restart Brainstorm.

'''brainstorm''': Start Brainstorm with the regular GUI.

'''brainstorm nogui''': Start in silent mode, without the GUI, without progress bar, without user interactions.

If you want to start Brainstorm only if Brainstorm is not already running, you can use the following code:

{{{
if ~brainstorm('status')
    brainstorm nogui
end
}}}
To select a specific protocol at the beginning of your script:

{{{
ProtocolName = 'TutorialIntroduction';
% Get the protocol index
iProtocol = bst_get('Protocol', ProtocolName);
if isempty(iProtocol)
    error(['Unknown protocol: ' ProtocolName]);
end
% Select the current procotol
gui_brainstorm('SetCurrentProtocol', iProtocol);
}}}
To delete the protocol and start over:

{{{
% Delete existing protocol
gui_brainstorm('DeleteProtocol', ProtocolName);
% Create new protocol
gui_brainstorm('CreateProtocol', ProtocolName, 0, 0);
}}}
== Database requests [TODO] ==
bst_get / bst_set

 * '''iStudy''': Index of the study, get the structure with sStudy=bst_get('Study', sInput.iStudy)
 * '''iItem''': Index of file in its category (example for recordings: sStudy.Data(sInput.iItem))

sStudy = bst_get('Study', sAvgData(1).iStudy)

sStudy.Data(sAvgData(1).iItem)

== File manipulation [TODO] ==
 * 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)
 * Reading / writing / exporting functions

== Additional quality control ==
Create figures and snapshot them

Generate warnings and errors

== Loop over subject and runs [TODO] ==
Creating loops is not supported yet by the script generator, but relatively easy to do from a script without having to know too much about Matlab programming.

1) Fill the cell array SubjectNames with all your subjects names, with  the same dimensions as the list of input raw files (sFiles)

 . 2) Add a "for" loop that includes all the bst_process() calls (leave the bst_report() calls and input definition outside)
 3) Inside the loop, replace SubjectNames with SubjectNames{i} and sFiles with sFiles(i)

== How to process an entire study ==
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  [[Tutorials/TutWarping|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.

== Final script ==
The following script from the Brainstorm distribution reproduces the introduction tutorials ("Get started"): '''brainstorm3/toolbox/script/tutorial_introduction.m'''

<<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 >)>>

<<BR>>For an example of a script illustrating how to create loops, look at the tutorial [[Tutorials/VisualSingle|MEG visual: single subject]]. '''brainstorm3/toolbox/script/tutorial_visual_single.m'''

<<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_single.m")>><<HTML(</div >)>>

== 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.

<<EmbedContent("http://neuroimage.usc.edu/bst/get_prevnext.php?prev=Tutorials/Workflows")>>

Feedback: Comments, bug reports, suggestions, questions

Email address (if you expect an answer):