%% read the continuous data and segment into 2 seconds epochs
%eeglab

%EEG_data = eeglab2fieldtrip('pre1.set', 'preprocessing', 'none');

%hdr    = ft_read_header('EEG.mat'); 
%data   = ft_read_data('EEG.mat', 'header', hdr ); 
%events = ft_read_event('EEG.mat', 'header', hdr );

cfg            = [];
cfg.dataset    = './muse_eeg_first_treatment/post2.edf'; % note that you may need to add the full path to the .ds directory
cfg.continuous = 'yes';
cfg.channel    = {'EEG'};
data           = ft_preprocessing(cfg);

cfg         = [];
cfg.length  = 2;
cfg.overlap = 0.5;
data        = ft_redefinetrial(cfg, data);

% this step is needed to 1) remove the DC-component, and to 2) get rid of a few segments of data at
% the end of the recording, which contains only 0's.
cfg        = [];
cfg.demean = 'yes';
cfg.trials = 1:(numel(data.trial)-6);
data       = ft_preprocessing(cfg, data);




%% use ICA in order to identify cardiac and blink components
%cfg                 = [];
%cfg.method          = 'runica';
%cfg.runica.maxsteps = 50;
%cfg.randomseed      = 0; % this can be uncommented to match the data that has been stored on disk
%comp                = ft_componentanalysis(cfg, datads);

%% visualize components

% these were the indices of the bad components that were identified
% they may be different if you re-run the ICA decomposition with a random randomseed.
%badcomp = [2 3 7 16];

%cfg            = [];
%cfg.channel    = badcomp;
%cfg.layout     = 'EEG1005.mat';
%cfg.compscale  = 'local';
%cfg.continuous = 'yes';
%ft_databrowser(cfg, comp);

%cfg           = [];
%cfg.component = badcomp;
%dataica       = ft_rejectcomponent(cfg, comp);


%cfg            = [];
%cfg.component  = badcomp;
%dataica        = ft_rejectcomponent(cfg, comp);


%% compute the power spectrum
cfg              = [];
cfg.output       = 'pow';
cfg.method       = 'mtmfft';
cfg.taper        = 'dpss';
cfg.tapsmofrq    = 1;
cfg.keeptrials   = 'no';
datapow          = ft_freqanalysis(cfg, data);


%% load the required geometrical information
load standard_bem.mat
load standard_sourcemodel3d8mm.mat

headmodel=vol


%% compute the leadfield
cfg             = [];
cfg.sourcemodel = sourcemodel;
cfg.headmodel   = headmodel;
cfg.channel     = {'EEG'};
cfg.elec        = 'standard_1005.elc'
lf              = ft_prepare_leadfield(cfg, data);



%% compute sensor level Fourier spectra, to be used for cross-spectral density computation.
cfg            = [];
cfg.method     = 'mtmfft';
cfg.output     = 'fourier';
cfg.keeptrials = 'yes';
cfg.tapsmofrq  = 1;
cfg.foi        = 10;
freq           = ft_freqanalysis(cfg, data);


%% do the source reconstruction
cfg                   = [];
cfg.frequency         = freq.freq;
cfg.method            = 'pcc';
cfg.grid              = lf;
cfg.headmodel         = headmodel;
cfg.keeptrials        = 'yes';
cfg.pcc.lambda        = '10%';
cfg.pcc.projectnoise  = 'yes';
cfg.pcc.fixedori      = 'yes';
source = ft_sourceanalysis(cfg, freq);
source = ft_sourcedescriptives([], source); % to get the neural-activity-index



load('dkatlas.mat')

%% plot the neural activity index (power/noise)
cfg               = [];
cfg.method        = 'surface';
cfg.funparameter  = 'nai';
cfg.maskparameter = cfg.funparameter;
cfg.funcolorlim   = [0.0 8];
cfg.opacitylim    = [3 8];
cfg.opacitymap    = 'rampup';
cfg.funcolormap   = 'jet';
cfg.colorbar      = 'no';
ft_sourceplot(cfg, source, dkatlas);
view([-90 30]);
light;


%% compute sensor level single trial power spectra
cfg              = [];
cfg.output       = 'pow';
cfg.method       = 'mtmfft';
cfg.taper        = 'dpss';
cfg.foilim       = [9 11];
cfg.tapsmofrq    = 1;
cfg.keeptrials   = 'yes';
datapow           = ft_freqanalysis(cfg, data);

%% identify the indices of trials with high and low alpha power
freqind = nearest(datapow.freq, 10);
tmp     = datapow.powspctrm(:,:,freqind);
chanind = find(mean(tmp,1)==max(mean(tmp,1)));  % find the sensor where power is max
indlow  = find(tmp(:,chanind)<=median(tmp(:,chanind)));
indhigh = find(tmp(:,chanind)>=median(tmp(:,chanind)));



%% compute the power spectrum for the median splitted data
cfg              = [];
cfg.trials       = indlow;
datapow_low      = ft_freqdescriptives(cfg, datapow);

cfg.trials       = indhigh;
datapow_high     = ft_freqdescriptives(cfg, datapow);

%% compute the difference between high and low
cfg = [];
cfg.parameter = 'powspctrm';
cfg.operation = 'divide';
powratio      = ft_math(cfg, datapow_high, datapow_low);

%% compute fourier spectra for frequency of interest according to the trial split
cfg            = [];
cfg.method     = 'mtmfft';
cfg.output     = 'fourier';
cfg.keeptrials = 'yes';
cfg.tapsmofrq  = 1;
cfg.foi        = 10;

cfg.trials = indlow;
freq_low   = ft_freqanalysis(cfg, data);

cfg.trials = indhigh;
freq_high  = ft_freqanalysis(cfg, data);

%% compute the beamformer filters based on the entire data
cfg                   = [];
cfg.frequency         = freq.freq;
cfg.method            = 'pcc';
cfg.grid              = lf;
cfg.headmodel         = headmodel;
cfg.keeptrials        = 'yes';
cfg.pcc.lambda        = '10%';
cfg.pcc.projectnoise  = 'yes';
cfg.pcc.keepfilter    = 'yes';
cfg.pcc.fixedori      = 'yes';
source = ft_sourceanalysis(cfg, freq);

% use the precomputed filters
cfg                   = [];
cfg.frequency         = freq.freq;
cfg.method            = 'pcc';
cfg.grid              = lf;
cfg.sourcemodel.filter       = source.avg.filter;
cfg.headmodel         = headmodel;
cfg.keeptrials        = 'yes';
cfg.pcc.lambda        = '10%';
cfg.pcc.projectnoise  = 'yes';
source_low  = ft_sourcedescriptives([], ft_sourceanalysis(cfg, freq_low));
source_high = ft_sourcedescriptives([], ft_sourceanalysis(cfg, freq_high));

cfg           = [];
cfg.operation = 'log10(x1)-log10(x2)';
cfg.parameter = 'pow';
source_ratio  = ft_math(cfg, source_high, source_low);

%% sourceinterpolate
cfg = [];
cfg.parameter    = 'pow';
sourceint = ft_sourceinterpolate(cfg,source_ratio,dkatlas);
cfg=[];
sourceint = ft_sourceparcellate(cfg,sourceint,dkatlas);

% create a fancy mask
source_ratio.mask = (1+tanh(2.*(source_ratio.pow./max(source_ratio.pow(:))-0.5)))./2;

cfg = [];
cfg.method        = 'surface';
cfg.funparameter  = 'pow';
cfg.maskparameter = 'mask';
cfg.funcolorlim   = [-.3 .3];
cfg.funcolormap   = 'jet';
cfg.colorbar      = 'no';
ft_sourceplot(cfg, source_ratio);
view([-90 30]);
light('style','infinite','position',[0 -200 200]);


%% compute connectivity
cfg         = [];
cfg.method  ='coh';
cfg.complex = 'absimag';
source_conn = ft_connectivityanalysis(cfg, source);


figure;imagesc(source_conn.cohspctrm);


cfg = [];
cfg.parcellation = 'parcellation';
cfg.parameter    = 'cohspctrm';
parc_conn = ft_sourceparcellate(cfg, source_conn, dkatlas);

figure;imagesc(parc_conn.cohspctrm);


cfg           = [];
cfg.method    = 'degrees';
cfg.parameter = 'cohspctrm';
cfg.threshold = .1;
network_full = ft_networkanalysis(cfg,source_conn);

%% sourceinterpolate
cfg = [];
cfg.parameter    = 'degrees';
network_int = ft_sourceinterpolate(cfg,network_full,dkatlas);
cfg=[];
network_int = ft_sourceparcellate(cfg, network_int, dkatlas);
%%
% create a fancy mask

cfg = [];
cfg.method        = 'surface';
cfg.funparameter  = 'degrees';
cfg.colorbar      = 'no';
figure(8);ft_sourceplot(cfg, network_int);
view([-90 30]);
light('style','infinite','position',[0 -200 200]);
colorbar off
material dull
set(gcf,'color','w');