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 - &mu;) / &sigma;

% 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 - &mu;) / &sigma;

% 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 - &mu;) / &mu; * 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 - &mu;) / &mu; * 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 - &mu;) / &mu; * 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 - &mu;) / &mu; * 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 - &mu;) / &mu; * 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 - &mu;) / &sigma;
    '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 - &mu;) / &sigma;
    '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]);
