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

%%

cfg = [];
cfg.channel = 'MEGMAG';
cfg.method = 'wavelet';
cfg.output = 'fourier';
cfg.keeptrials = 'yes';
cfg.foi = 16;
cfg.toi = 0.150;
freq = ft_freqanalysis(cfg, data_fix);


%%

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

% save grid grid


%%

cfg           = [];
cfg.vol       = vol;  % from FT
cfg.grad      = sens; % from FT
cfg.senstype  = 'meg';
cfg.grid      = grid;
cfg.method    = 'pcc';
cfg.pcc.fixedori = 'yes';
cfg.latency   = [0.140 0.160];
cfg.frequency = [14 18];

source = ft_sourceanalysis(cfg, freq);

figure
plot(source.avg.mom{source.inside(1)}, '.')
xlabel('real');
ylabel('imag');


%%

pos = [21 -64 30];

% compute the nearest grid location
dif = grid.pos;
dif(:,1) = dif(:,1)-pos(1);
dif(:,2) = dif(:,2)-pos(2);
dif(:,3) = dif(:,3)-pos(3);
dif = sqrt(sum(dif.^2,2));
[distance, refindx] = min(dif);

cfg = [];
cfg.method    = 'coh';
% cfg.complex   = 'abs';
cfg.complex   = 'absimag';
cfg.refindx   = refindx;
conn = ft_connectivityanalysis(cfg, source);

% the output contains both the actual source position, as well as the position of the reference
% this is ugly and will probably change in future FieldTrip versions
orgpos = conn.pos(:,1:3);
refpos = conn.pos(:,4:6);
conn.pos = orgpos;

%% visualise the seed-based connectivity results

cfg               = [];
cfg.funparameter  = 'cohspctrm';
figure; ft_sourceplot(cfg, conn);

cfg             = [];
cfg.parameter   = 'cohspctrm';
sourceI = ft_sourceinterpolate(cfg, conn, mri_realigned);

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

%% look at connectivity difference

cfg         = [];
cfg.trials  = find(freq.trialinfo==1 | freq.trialinfo==2); % 1=Famous, 2=Unfamiliar
source1 = ft_selectdata(cfg, source);

cfg.trials  = find(freq.trialinfo==3); % 3=Scrambled
source2 = ft_selectdata(cfg, source);

cfg         = [];
cfg.method  = 'coh';
cfg.complex = 'absimag';
cfg.refindx = refindx;
conn1 = ft_connectivityanalysis(cfg, source1);
conn2 = ft_connectivityanalysis(cfg, source2);

conn1.pos = conn1.pos(:,1:3);
conn2.pos = conn2.pos(:,1:3);

cfg           = [];
cfg.parameter = 'cohspctrm';
cfg.operation = 'x1-x2';
conn_dif = ft_math(cfg, conn1, conn2);

cfg           = [];
cfg.parameter = 'cohspctrm';
source1int    = ft_sourceinterpolate(cfg, conn1, mri_realigned);
source2int    = ft_sourceinterpolate(cfg, conn2, mri_realigned);
source_difint = ft_sourceinterpolate(cfg, conn_dif, mri_realigned);

cfg               = [];
cfg.funparameter  = 'cohspctrm';
cfg.funcolorlim   = [-0.1 0.1];
cfg.maskparameter = 'cohspctrm';
cfg.opacitylim    = [-0.15 0.15];
figure; ft_sourceplot(cfg, source_difint);

%% look at the analysis history

% for saving to disk
prefix = sprintf('Sub%02d', subj);

cfg           = [];
cfg.filename  = [prefix '_source_difint.html'];
ft_analysispipeline(cfg, source_difint);

