
%% start with data that was preprocessed in FieldTrip
subj = 15;
prefix = sprintf('Sub%02d', subj);
load([prefix '_raw']);  % this is called "data" rather than "raw"

% load the results from beamformer_part1
load vol
load sens
load mri_realigned


%% deal with maxfilter

% the data has been maxfiltered and subsequently contatenated
% this results in an ill-conditioned estimate of covariance or CSD

cfg             = [];
cfg.method      = 'pca';
cfg.updatesens  = 'no';
cfg.channel     = 'MEGMAG';
comp = ft_componentanalysis(cfg, data);

cfg             = [];
cfg.updatesens  = 'no';
cfg.component   = comp.label(51:end);
data_fix = ft_rejectcomponent(cfg, comp);


%%

pos1 = [21 -64 30];
pos2 = [0 35 83];

cfg = [];
cfg.location = pos1;
figure; ft_sourceplot(cfg, mri_realigned);
cfg.location = pos2;
figure; ft_sourceplot(cfg, mri_realigned);

%%

cfg             = [];
cfg.vol         = vol;
cfg.grad        = sens;
cfg.senstype    = 'meg';
cfg.method      = 'lcmv';
cfg.lcmv.keepfilter = 'yes';
cfg.lcmv.projectmom = 'yes';
cfg.grid.unit   = 'mm';
cfg.grid.pos    = pos1;
source1 = ft_sourceanalysis(cfg, timelock2);

cfg.grid.pos = pos2;
source2 = ft_sourceanalysis(cfg, timelock2);


%% construct single-trial virtual channel representation

virtualchannel_raw = [];
virtualchannel_raw.label = {'pos1'; 'pos2'};
virtualchannel_raw.trialinfo = data_fix.trialinfo;
for i=1:882
  % note that this is the non-filtered raw data
  virtualchannel_raw.time{i}       = data_fix.time{i};
  virtualchannel_raw.trial{i}(1,:) = source1.avg.filter{1} * data_fix.trial{i}(:,:);
  virtualchannel_raw.trial{i}(2,:) = source2.avg.filter{1} * data_fix.trial{i}(:,:);
end

%%

cfg                 = [];
cfg.keeptrials      = 'yes';
cfg.preproc.demean  = 'yes';
cfg.preproc.baselinewindow = [-inf 0];
virtualchannel_avg = ft_timelockanalysis(cfg, virtualchannel_raw);

figure
plot(virtualchannel_avg.time, virtualchannel_avg.avg)
legend(virtualchannel_avg.label);


%%

cfg         = [];
cfg.method  = 'wavelet';
cfg.output  = 'powandcsd';
cfg.foi     = 4:2:70;
cfg.toi     = -0.200:0.020:1.000;
virtualchannel_wavelet = ft_freqanalysis(cfg, virtualchannel_raw);

cfg                 = [];
cfg.baselinetype    = 'relative';
cfg.baseline        = [-inf 0];
cfg.channel         = {'pos1'};
figure; ft_singleplotTFR(cfg, virtualchannel_wavelet);

cfg.channel         = {'pos2'};
figure; ft_singleplotTFR(cfg, virtualchannel_wavelet);


%%

cfg = [];
cfg.method = 'coh';
coherence = ft_connectivityanalysis(cfg, virtualchannel_wavelet);


figure
imagesc(coherence.time, coherence.freq, squeeze(coherence.cohspctrm(1,:,:)));
axis xy
