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

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

%% compute covariance

cfg = [];
cfg.channel = 'MEGMAG';
cfg.covariance = 'yes'; % compute the covariance for single trials, then average
% cfg.preproc.bpfilter = 'yes';
% cfg.preproc.bpfreq = [5 70];
% cfg.preproc.hpfilter = 'yes';
% cfg.preproc.hpfreq = 1;
% cfg.preproc.derivative = 'yes';
cfg.preproc.demean = 'yes';             % the PCA cleanup shifted the baseline
cfg.preproc.baselinewindow = [-inf 0];  % reapply the baseline correction
cfg.keeptrials = 'yes';
timelock1 = ft_timelockanalysis(cfg, data_fix);

cfg = [];
cfg.covariance = 'yes'; % compute the covariance of the averaged ERF
timelock2 = ft_timelockanalysis(cfg, timelock1);

figure
plot(timelock2.time, timelock2.avg)

cfg = [];
cfg.layout = 'neuromag306mag.lay';
figure; ft_multiplotER(cfg, timelock2);


%%

pos = [21 -64 30];

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

figure
plot(source.time, source.avg.mom{1})

%% construct single-trial virtual channel data

virtualchannel_raw = [];
virtualchannel_raw.label = {'source'};
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,:) = source.avg.filter{1} * data_fix.trial{i}(:,:);
end

%% average the virtual channel ERP

cfg = [];
cfg.keeptrials = 'yes';
cfg.preproc.demean = 'yes';
cfg.preproc.baselinewindow = [-inf 0];
virtualchannel_avg = ft_timelockanalysis(cfg, virtualchannel_raw);
cfg.trials = virtualchannel_raw.trialinfo==1;
virtualchannel_avg1 = ft_timelockanalysis(cfg, virtualchannel_raw);
cfg.trials = virtualchannel_raw.trialinfo==2;
virtualchannel_avg2 = ft_timelockanalysis(cfg, virtualchannel_raw);
cfg.trials = virtualchannel_raw.trialinfo==3;
virtualchannel_avg3 = ft_timelockanalysis(cfg, virtualchannel_raw);

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

figure
plot(virtualchannel_avg.time, [virtualchannel_avg1.avg; virtualchannel_avg2.avg; virtualchannel_avg3.avg]);
legend({'1-Famous', '2-Unfamiliar', '3-Scrambled'})

figure
imagesc(squeeze(virtualchannel_avg.trial))


%% investigate the virtual channel spectrally

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

cfg = [];
cfg.baseline = [-inf 0];
cfg.baselinetype = 'relative';
cfg.interactive = 'no';
figure; ft_singleplotTFR(cfg, virtualchannel_wavelet);


