%% 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 seedLet's choose a location of interest, here the central culcsus 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}) ; |