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

 

%% Localisation of max beta suppression

[minval, minpowindx] = min(source_diff.avg.pow(:));

pos_cm = source_diff.pos(minpowindx,:) ;

 

%% You may choose yourself a location of interest

pos_cm = [0.0 5.0 9.0 ] ; %% CS

%% pos_cm = [0.0 -4.0 5.0 ] ; %% Others loc --- checked where !

 

 

%% Used to store the EP computed at each virtual sensor location

evoked_virtual_zscore = {} ;

evoked_virtual = {} ;

index = 1 ;

 

for i_scan = -60:5:60

 

 

    %% Compute filter for corresponding virtual sensor

    cfg                    = [];

    cfg.method             = 'lcmv';

    cfg.headmodel          = vol ;

    cfg.grid.pos           = [pos_cm(1)+i_scan/10 pos_cm(2) pos_cm(3)] ;

    %% cfg.grid.pos           = [pos_cm(1) pos_cm(2)+i_scan/10 pos_cm(3)] ;

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

    pow_data = [];

    pow_data.label = {'pow_x', 'pow_y', 'pow_z'};

    pow_data.time = trials.time;

    for i=1:length(trials.trial)

        pow_data.trial{i} = beamformer * trials.trial{i}(chansel,:);

    end

    

    %     %% construct a single virtual channel in the maximum power orientation (not usefull if fixed orig above)

    %     timeseries = cat(2, pow_data.trial{:});

    %     [u, s, v] = svd(timeseries, 'econ');

    %

    %     % this is equal to the first column of matrix V, apart from the scaling with s(1,1)

    %     timeseriesmaxproj = u(:,1)' * timeseries;

    %     virtualchanneldata = [];

    %     virtualchanneldata.label = {'cortex'};

    %     virtualchanneldata.time = trials.time;

    %     for i=1:length(trials.trial)

    %         virtualchanneldata.trial{i} = u(:,1)' * beamformer * trials.trial{i}(chansel,:);

    %     end

    

    %% Together with option fixedori (fixed oreintation) this is suffisent (we have only one vector not three)

    virtualchanneldata = [];

    virtualchanneldata.label = {'cortex'};

    virtualchanneldata.time = trials.time;

    for i=1:length(trials.trial)

        virtualchanneldata.trial{i} = pow_data.trial{i} ;

    end

    

    

%     %% Time-frequency analysis

%     cfg              = [];

%     cfg.output       = 'pow';

%     cfg.channel      = 'cortex';

%     cfg.method       = 'mtmconvol';

%     cfg.taper        = 'hanning';

%     cfg.foi          = 7:2:100;                         % 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 -0.5 to 1.5 sec in steps of 0.05 sec (50 ms)

%     TFRvirtualchanneldata = ft_freqanalysis(cfg, virtualchanneldata);

%     title(['X = ' num2str(pos_cm(1)+i_scan/10)]) ;

%     print(['TF-virtualChannel_alpha_' num2str(pos_cm(1)+i_scan/10)],'-dpng')

% 

%     %% Display the results as a TF map

%     cfg = [];

%     cfg.baseline     = reference_time_window ;

%     cfg.baselinetype = 'db';

%     cfg.channelname   = 'cortex'; % top figure

%     cfg.zlim = [-5 5] ;

%     figure;ft_singleplotTFR(cfg, TFRvirtualchanneldata);

    

    

    %% Evoked potential on the virtual sensor

    cfg.demean             = 'yes';

    cfg.baselinewindow     = reference_time_window;

    evoked_virtual{index}  = ft_timelockanalysis(cfg, virtualchanneldata);

    

    evoked_virtual_zscore{index} = evoked_virtual{index} ;

    mean_bl = mean(evoked_virtual{index}.avg(1,nearest(evoked_virtual{1}.time, -2.0):nearest(evoked_virtual{1}.time, -1.8)) ,2) ;

    std_bl = std(evoked_virtual{index}.avg(1,nearest(evoked_virtual{1}.time, -2.0):nearest(evoked_virtual{1}.time, -1.8)) , 0, 2) ;

    evoked_virtual_zscore{index}.avg = (evoked_virtual{index}.avg - mean_bl)/std_bl ;

    

    %% Display the evoked potential

    cfg.ylim = [-0.2 0.2] ;

    figure ; ft_singleplotER(cfg, evoked_virtual{index}) ; 

    title(['X = ' num2str(pos_cm(1)+i_scan/10)]) ;

    print(['EV-mag-virtualChannel_beta_' num2str(pos_cm(2)+i_scan/10)],'-dpng')

    

    cfg.ylim = [-10 10] ;

    figure ; ft_singleplotER(cfg, evoked_virtual_zscore{index}) ; 

    title(['X = ' num2str(pos_cm(1)+i_scan/10)]) ;

    print(['EV-mag-zscore-virtualChannel_beta_' num2str(pos_cm(2)+i_scan/10)],'-dpng')

    

    index = index + 1 ;

end

 

figure

for i = 1:index-1

   plot(evoked_virtual_zscore{1}.time, evoked_virtual_zscore{i}.avg);

   hold on

end

hold off

 

 

%% Coherence with EMG

 

%% %% First using EMG

channelsofinterest = {'MEG*2', 'MEG*3'}; %% Gradiometers on an Elekta system

 

 

%% Compute cross spectrum matrix

cfg            = [];

cfg.method     = 'mtmfft';

cfg.output     = 'powandcsd';

cfg.taper           = 'dpss';

cfg.tapsmofrq       = 5;

cfg.foi             = 10;

cfg.keeptrials = 'yes';

cfg.channel    = {'MEG*2' 'MEG*3' 'BIO005'};

cfg.channelcmb = {'MEG*2' 'MEG*2'; 'MEG*3' 'MEG*3';, 'MEG*2' 'MEG*3' ; 'MEG*2' 'BIO005' ; 'MEG*3' 'BIO005'};

 

powcsd           = ft_freqanalysis(cfg, data_timewindow_active);

 

 

% Create leadfield grid

cfg                 = [];

cfg.channel         = channelsofinterest ;

cfg.grad            = powcsd.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);

 

cfg                 = [];

cfg.method          = 'dics';

cfg.refchan         = 'BIO005';

cfg.frequency       = freqofinterest ;

cfg.grid            = grid;

cfg.vol             = vol;

cfg.grid.unit       = 'cm';

cfg.dics.lambda  = beamformer_lambda_normalization;

%% cfg.dics.reducerank = 2; % default for MEG is 2, for EEG is 3

 

source              = ft_sourceanalysis(cfg, powcsd);

 

cfg              = [];

cfg.method       = 'ortho';

cfg.funparameter = 'coh';

ft_sourceplot(cfg, source);