%% 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 warped MNI grid
load grid ;
%% 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);

cfg              = [];
cfg.output       = 'pow';
cfg.channel      = 'MEG';
cfg.method       = 'mtmconvol';
cfg.taper        = 'hanning';
cfg.foi          = 7:2:40;                         % analysis 2 to 30 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 -2.0 to 2.0 sec in steps of 0.05 sec (50 ms)
cfg.trials       = 'all';

TFRhann = ft_freqanalysis(cfg, full_data_all);
 
%% Visualization
cfg = [];
cfg.xlim = [0.1 0.2];                
cfg.ylim = [20 20];                  
cfg.zlim = [-1e-28 1e-28];           
cfg.baseline     = [-0.2 -0.0]; 
cfg.baselinetype = 'absolute';
cfg.layout       = 'neuromag306mag.lay';

figure; ft_topoplotTFR(cfg,TFRhann);

% for the multiple plots also:
cfg = [];
cfg.ylim = [20 20];                  
cfg.baseline     = [-0.2 -0.0]; 
cfg.baselinetype = 'absolute';
cfg.layout       = 'neuromag306mag.lay';
cfg.xlim = [-0.4:0.1:0.4];
cfg.comment = 'xlim';
cfg.commentpos = 'title';
figure; ft_topoplotTFR(cfg,TFRhann);

cfg = [];
cfg.baseline     = [-1.0 -0.8]; 
cfg.baselinetype = 'absolute'; 
%% cfg.zlim         = [-1e-27 1e-27];           
cfg.showlabels   = 'yes';   
cfg.layout       = 'neuromag306mag.lay';
figure; ft_multiplotTFR(cfg, TFRhann);

%% Compute cross-sepctral density matrices
powcsd_all = compute_crosprect(channelsofinterest, freqofinterest, freqhalfwin, data_all) ;
powcsd_active = compute_crosprect(channelsofinterest, freqofinterest, freqhalfwin, data_timewindow_active) ;

powcsd_ref = compute_crosprect(channelsofinterest, freqofinterest, freqhalfwin, data_timewindow_ref) ;

% Create leadfield grid
cfg                 = [];
cfg.channel         = channelsofinterest ;
cfg.grad            = powcsd_active.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);

%% Compute DICS filters for the concatenated data
cfg              = [];
cfg.channel      = channelsofinterest ;
cfg.method       = 'dics';
cfg.frequency    = freqofinterest ;
cfg.grid         = grid;
cfg.headmodel    = vol;
cfg.senstype     = 'MEG';
cfg.dics.keepfilter   = 'yes'; % We wish to use the calculated filter later on
cfg.dics.projectnoise = 'yes';
cfg.dics.lambda  = beamformer_lambda_normalization;
source_all = ft_sourceanalysis(cfg, powcsd_all);

%% Apply computed filters to each active and reference windows
cfg              = [];
cfg.channel      = channelsofinterest ;
cfg.method       = 'dics';
cfg.frequency    = freqofinterest ;
cfg.grid         = grid;
cfg.grid.filter  = source_all.avg.filter;
cfg.dics.projectnoise = 'yes';
cfg.headmodel          = vol;
cfg.dics.lambda  = beamformer_lambda_normalization ;
cfg.senstype     ='MEG' ;

%% Source analysis active
source_active = ft_sourceanalysis(cfg, powcsd_active);

%% Source analysis reference
source_reference = ft_sourceanalysis(cfg, powcsd_ref);

%% Compute differences (power) between active and reference
source_diff = source_active ;
source_diff.avg.pow  = (source_active.avg.pow - source_reference.avg.pow) ./ (source_active.avg.pow + source_reference.avg.pow);

%% Display the results
visuBF( source_diff, 'title' )
%% MRI reslicing

mri_resliced = ft_volumereslice([], mri_aligned);

%% Display results superimposed on MRI

cfg            = [];
cfg.parameter = 'avg.pow';
source_active_int  = ft_sourceinterpolate(cfg, source_active, mri_resliced);
source_reference_int  = ft_sourceinterpolate(cfg, source_reference, mri_resliced);
source_diff_int  = source_active_int;
source_diff_int.pow  = (source_active_int.pow - source_reference_int.pow) ./ (source_active_int.pow + source_reference_int.pow);

%% Display the results
visuBF( source_diff_int, 'title' )