ft_defaults ; %% Load fieldtrip configuration

%% Read two runs, define trials around LEFT_MVT_CLEAN marker (-2s to 2s around the marker), read also BIO005 which is an EMG.
[ trials ] = trial_definition_ft( {'cardiocor_blinkcor_run03_tsss.fif' 'cardiocor_blinkcor_run04_tsss.fif'}, 'LEFT_MVT_CLEAN', 2.0, 2.0, 'BIO005' ) ;

%% Load the head model created during the MRI pre processing
load vol ; %% load head model located in the subject directory and obtained from the T1 MRI of the subject
 
%% Load the MRI data for visualisation purpose
load('mri_aligned.mat') ; %% MRI data needed for vizualisation
load('segmentedmri.mat') ;

overall_time_window = [-2 2] ; %% if you modify this modify the trials def above ! Is this used below ???

%% Define time window for active state and baseline
active_time_window = [0.2 0.6] ;
reference_time_window = [-1.4 -1.0] ;

%% Basic configuration of the beamformer
beamformer_lambda_normalization = '10%' ;
lead_field_depth_normalization = 'column' ;

%% Frequency of interest and associated width (we will examine here the 16 Hz to 28 Hz band)
freqofinterest = 22 ;
freqhalfwin = 6 ;

%% Select the channel of interest
channelsofinterest = {'MEG*2', 'MEG*3'}; %% Gradiometers on an Elekta system
%% channelsofinterest = {'MEG*1'} ; %% Magnetometer on an Elekta system

% Select time window of interest
cfg = [];
cfg.toilim = active_time_window ;
data_timewindow_active = ft_redefinetrial(cfg,trials);

% Select time window of control
cfg = [];
cfg.toilim = reference_time_window ;
data_timewindow_ref = ft_redefinetrial(cfg,trials);

% Concatenate for common spatial filter computation
data_all = ft_appenddata([], data_timewindow_active, data_timewindow_ref);
 
% Trials will full window
cfg = [];
cfg.toilim = overall_time_window ;
full_data_all = ft_redefinetrial(cfg,trials);

%% Compute covariance of the data (as a whole - both conditions)
cfg                   = [];
cfg.covariance        = 'yes';
cfg.channel           = channelsofinterest ;
cfg.vartrllength      = 0;
cfg.covariancewindow  = 'all';
cfg.trials            = 'all';
tlock                 = ft_timelockanalysis(cfg, trials);

% Create leadfield grid

cfg                 = [];
cfg.channel         = channelsofinterest ;
cfg.grad            = data_all.grad;
cfg.vol             = vol ;
cfg.dics.reducerank = 2; % default for MEG is 2, for EEG is 3
cfg.grid.resolution = 1.0;   % use a 3-D grid with a 0.5 cm resolution
cfg.grid.unit       = 'cm';
cfg.grid.tight      = 'yes';
cfg.normalize       =  lead_field_depth_normalization ; %% Depth normalization
[grid] = ft_prepare_leadfield(cfg);

%% LCMV ANALYSIS TO GET TEMPORAL INFORMATION with virtual sensors

%% Retrieve correct labelling of the trials
data_all.trialinfo = [zeros(length(data_timewindow_active.trial), 1); ones(length(data_timewindow_ref.trial), 1)];

%% Compute covariance of the data (as a whole - both conditions)
cfg                   = [];
cfg.covariance        = 'yes';
cfg.channel           = channelsofinterest ;
cfg.vartrllength      = 0;
cfg.covariancewindow  = 'all';
cfg.trials            = 'all';
tlock                 = ft_timelockanalysis(cfg, trials);

%% Let's choose a location of interest, here the central cursus we will define our virtual channels based on this seed WARNING: Unit = centimetre
pos_cm = [4.0 2.0 10.0 ] ;

%% Used to store the EP computed at each virtual sensor location
evoked_virtual_zscore = {} ;
evoked_virtual = {} ;

%% Compute filter for corresponding virtual sensor
cfg                    = [];
cfg.method             = 'lcmv';
cfg.headmodel          = vol ;
cfg.grid.pos           = pos_cm ;
cfg.grid.inside        = 1:size(cfg.grid.pos, 1);
cfg.grid.outside       = [];
cfg.grid.unit          = 'cm'; %% DO NOT FORGET THAT !
cfg.keepfilter         = 'yes';
cfg.lcmv.fixedori      = 'yes'; % project on axis of most variance using SVD
cfg.lcmv.projectnoise  = 'yes' ;
cfg.lcmv.lambda        = beamformer_lambda_normalization ;
cfg.lcmv.weightnorm    = 'nai' ;
cfg.lcmv.normalize     = lead_field_depth_normalization ; %% Depth normalization
source_idx       = ft_sourceanalysis(cfg, tlock);
beamformer = source_idx.avg.filter{1};

%% Apply filter to the raw data
chansel = ft_channelselection(channelsofinterest, data_all.label); % find MEG sensor names
chansel = match_str(data_all.label, chansel);
virtualchanneldata = [];
virtualchanneldata.label = {'virtual_sensor'};
virtualchanneldata.time = trials.time;
for i=1:length(trials.trial)
    virtualchanneldata.trial{i} = beamformer * trials.trial{i}(chansel,:);
end

%% Time-frequency analysis
cfg              = [];
cfg.output       = 'pow';
cfg.channel      = 'virtual_sensor';
cfg.method       = 'mtmconvol';
cfg.taper        = 'hanning';
cfg.foi          = 7:2:40;                         % analysis 2 to 40 Hz in steps of 2 Hz
cfg.t_ftimwin    = ones(length(cfg.foi),1).*0.5;   % length of time window = 0.5 sec
cfg.toi          = -2.0:0.05:2.0;                  % time window "slides" from -0.5 to 1.5 sec in steps of 0.05 sec (50 ms)
TFRvirtualchanneldata = ft_freqanalysis(cfg, virtualchanneldata);

%% Display the results as a TF map
cfg = [];
cfg.baseline     = reference_time_window ;
cfg.baselinetype = 'db';
cfg.channelname   = 'virtual_sensor'; % top figure
cfg.zlim = [-5 5] ;
figure;ft_singleplotTFR(cfg, TFRvirtualchanneldata);
title(['X = ' num2str(pos_cm(1)) ' Y = ' num2str(pos_cm(2)) ' Y = ' num2str(pos_cm(3))]) ;
 
%% Evoked potential on the virtual sensor
cfg.demean             = 'yes';
cfg.baselinewindow     = reference_time_window;
evoked_virtual{1}  = ft_timelockanalysis(cfg, virtualchanneldata);

%% Display the evoked potential
cfg.ylim = [-0.25 0.25] ;
figure ; ft_singleplotER(cfg, evoked_virtual{1}) ;