%% get data from SPM

D = spm_eeg_load('../spm-source-demo/PapMcbdspmeeg_run_01_sss.mat');

disp(D.condlist)

% convert data from SPM
raw      = D.ftraw(D.indchantype('MEGMAG'), D.indsample(-0.1):D.indsample(0.3), D.indtrial(D.condlist{1}, 'GOOD'));
timelock = D.fttimelock(D.indchantype('MEGMAG'), D.indsample(-0.1):D.indsample(0.3), D.indtrial(D.condlist{1}, 'GOOD'));

% alternative method
raw      = spm2fieldtrip(D);
timelock = ft_timelockanalysis([], raw);


%% 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"
load([prefix '_avg_Faces_vs_Scrambled']);
load([prefix '_avg_Famous']);
load([prefix '_avg_Unfamiliar']);
load([prefix '_avg_Scrambled']);

% load the results from part 1
load vol
load sens


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


%%

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

% save wavelet wavelet

cfg = [];
cfg.layout = 'neuromag306mag.lay';
cfg.baseline = [-inf 0];
cfg.baselinetype = 'relative';
figure; ft_multiplotTFR(cfg, wavelet)

%%

cfg = [];
cfg.grid.resolution = 7;
% cfg.inwardshift = -7; % allow dipoles 10mm outside the brain, this improves interpolation at the edges
cfg.grid.unit = 'mm';
cfg.vol  = vol;  % from FT
cfg.grad = sens; % from FT
cfg.senstype = 'meg';
cfg.normalize = 'yes';
grid = ft_prepare_leadfield(cfg, wavelet);

% save grid grid


%% perform whole-brain source reconstruction

cfg = [];
cfg.vol       = vol;  % from FT
cfg.grad      = sens; % from FT
cfg.senstype  = 'meg';
cfg.grid      = grid;
cfg.method    = 'dics';

cfg.frequency = [14 18];
cfg.latency   = [0.140 0.160];
sourceA = ft_sourceanalysis(cfg, wavelet);
cfg.latency   = [-0.100 -0.080];
sourceB = ft_sourceanalysis(cfg, wavelet);

% cfg.frequency = [40 65];
% cfg.latency   = [0.090 0.140];
% sourceA = ft_sourceanalysis(cfg, wavelet);
% cfg.latency   = [-0.100 -0.050];
% sourceB = ft_sourceanalysis(cfg, wavelet);

% 
% cfg.frequency = [12 20];
% cfg.latency   = [0.090 0.140];
% sourceA = ft_sourceanalysis(cfg, wavelet);
% cfg.latency   = [-0.050 0.000];
% sourceB = ft_sourceanalysis(cfg, wavelet);

% FT_MATH requires the time axis needs to be the same
sourceA.time = 0; 
sourceB.time = 0;

cfg = [];
cfg.parameter = 'pow';
cfg.operation = 'log10(x1/x2)'; % sourceA divided by sourceB
sourceR = ft_math(cfg, sourceA, sourceB);

cfg = [];
cfg.funparameter = 'pow';
figure; ft_sourceplot(cfg, sourceR);


%% interpolate and plot on individual anatomical MRI

cfg = [];
cfg.parameter = 'pow';
sourceI = ft_sourceinterpolate(cfg, sourceR, mri_realigned);

cfg = [];
cfg.funparameter = 'pow';
figure; ft_sourceplot(cfg, sourceI);



