function script_practicalMEEG_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
%         Takfarinas Medani, adapted for practicalMEEG, Oct 2025


%% ===== SCRIPT VARIABLES =====
% This script is adapted for the Zenodo data with the single subject sub-01
% Full list of subjects to process
% The protocol name has to be a valid folder name (no spaces, no weird characters...)
ProtocolName = 'PracticalMEEG2025_sub-01_Precomputed';
SubjectNames = {'sub-01'};
nbOfSubject = length(SubjectNames);
% Empty-room dates for each subject (so that we can match automatically recordings with empty-room)
EmptyRoomSubj = 'sub-emptyroom';
AcquisitionDates = {'09-Apr-2009'};
% Bad channels {iSubj} = {Run01, Run02, Run03, Run04, Run05, Run06}
BadChannels{1}  = {'EEG016'};
% SSP components to remove {iSubj} = {sRun01}, 
nbOfRun = 1; 
SspSelect{1}  = {{1,1}}; 

%% ===== INPUTS =====
% Parse command line
% 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

%% ===== START BRAINSTORM =====
% Start brainstorm without the GUI
if ~brainstorm('status')
    brainstorm nogui
end
% Disable visualization filters
panel_filter('SetFilters', 0, [], 0, [], 0, [], 0, 0);
% Restore all colormaps
bst_colormaps('RestoreDefaults', 'anatomy');
bst_colormaps('RestoreDefaults', 'meg');
bst_colormaps('RestoreDefaults', 'eeg');
bst_colormaps('RestoreDefaults', 'source');
bst_colormaps('RestoreDefaults', 'stat1');
bst_colormaps('RestoreDefaults', 'stat2');
bst_colormaps('RestoreDefaults', 'time');

% 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');
%% ===== CREATE PROTOCOL =====
% Delete existing protocol
gui_brainstorm('DeleteProtocol', ProtocolName);
% Create new protocol
gui_brainstorm('CreateProtocol', ProtocolName, 0, 0);


%% ===== PRE-PROCESS AND IMPORT =====
for iSubj = 1:nbOfSubject
    % 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: nbOfRun
        % 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:nbOfSubject
    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: 1 %length(SubjectNames)
    disp(sprintf('\n===== TIME-FREQUENCY: SUBJECT #%d =====\n', iSubj));
    for iRun = 1: nbOfRun
        % 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], ...
         };
        BadSeg = BadSegments{iSubj}{iRun}';
end