% This script is intended to calculate coherence between LFP's and EEGs bipolar channels 
% for conditions ntrl and emo 
% For more info about parameters look at:
% http://www.fieldtriptoolbox.org/tutorial/coherence

load Path_directory\Filename.mat;  % File saved after running Script 1.


refchans={'LFPL23','LFPR23'};
eegchans={'EEGT4Cz','EEGT3Cz','EEGPzCz','EEGOzPz','EEGFzCz'};

neigh = struct;
neigh(1).label = 'EEGT4Cz';
neigh(1).neighblabel = {'EEGT3Cz','EEGPzCz','EEGFzCz'};
neigh(2).label = 'EEGT3Cz';
neigh(2).neighblabel = {'EEGT4Cz','EEGPzCz','EEGFzCz'};
neigh(3).label = 'EEGPzCz';
neigh(3).neighblabel = {'EEGT4Cz','EEGT3Cz','EEGFzCz'};
neigh(4).label = 'EEGOzPz';
neigh(4).neighblabel = {'EEGPzCz'};
neigh(5).label = 'EEGFzCz';
neigh(5).neighblabel = {'EEGT4Cz','EEGT3Cz','EEGPzCz'};

cfg=[];
cfg.latency=[0 4];
ntrl=ft_selectdata(cfg,ntrl);
emo=ft_selectdata(cfg,emo);

cfg=[];
cfg.length=1;                    % Segments length
cfg.overlap=.5;                  % Overlap percentage among segments
ntrl_cohdata=ft_redefinetrial(cfg,ntrl);
emo_cohdata=ft_redefinetrial(cfg,emo);

%frequency analysis
cfg=[];
cfg.method='mtmfft';
cfg.output='fourier';
cfg.taper='hanning';            
cfg.foi=1:25;
cfg.keeptrials='yes';
freq_ntrl=ft_freqanalysis(cfg,ntrl_cohdata);
freq_emo=ft_freqanalysis(cfg,emo_cohdata);

cfg=[];
cfg.method='coh';            % Indicates that coherence is calculated
cfg.complex='real';          % Indicates whether real of imaginary part of coherence is taken       
coh_ntrl=ft_connectivityanalysis(cfg,freq_ntrl);
coh_emo=ft_connectivityanalysis(cfg,freq_emo);


stats=cell(1,2);
for r=1:2
    
%%%%%%Graphing Desired Data%%%%%%%%%%%%%%%%%%%%%%%%%
figure;
cfg=[];
cfg.refchannel=refchans{r};
cfg.parameter='cohspctrm';
cfg.xlim=[4 25];
%cfg.ylim=[0 0.21];
cfg.channel='EEGFzCz';
ft_singleplotER(cfg,coh_ntrl,coh_emo);
cfg.channel='EEGT3Cz';
ft_singleplotER(cfg,coh_ntrl,coh_emo);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

cfg = [];
cfg.refchan_zcoh_jan=refchans{r};
cfg.parameter        = 'fourierspctrm';
cfg.neighbours       = neigh;
%cfg.avgoverfreq     = 'yes';
cfg.frequency        = [4 25];
cfg.method           = 'montecarlo';
cfg.statistic        = 'indepsamplesT';
cfg.correctm         = 'cluster';   % Statistical analysis
cfg.clusteralpha     = 0.05;        % p-value
cfg.clusterstatistic = 'maxsum';
cfg.tail             = 0;
cfg.clustertail      = 0;
cfg.alpha            = 0.05;
cfg.correcttail    = 'alpha';
cfg.numrandomization = 1000;
cfg.design           = [zeros(1,size(freq_emo.fourierspctrm,1))+1,zeros(1,size(freq_ntrl.fourierspctrm,1))+2];
cfg.ivar             = 1;
stat_left = ft_freqstatistics(cfg, freq_emo,freq_ntrl);
stats{r}=stat_left;

if isfield(stats{r},'posclusters') && ~isempty(stats{r}.posclusters) && stats{r}.posclusters(1).prob<=0.025
i=find(stats{r}.posclusterslabelmat==1);
[c,f]=ind2sub(size(stats{r}.stat),i);
sig_freq=unique(stats{r}.freq(f));
sig_chan=stats{r}.label{c};
display(['sig. pos. cluster in channel ',sig_chan, ' at freq ',num2str(sig_freq)])
end

if isfield(stats{r},'negclusters') && ~isempty(stats{r}.negclusters) && stats{r}.negclusters(1).prob<=0.025
i=find(stats{r}.negclustersclusterslabelmat==1);
[c,f]=ind2sub(size(stats{r}.stat),i);
sig_freq=unique(stats{r}.freq(f));
sig_chan=stats{r}.label{c};
display(['sig. neg. cluster in channel ',sig_chan, ' at freq ',num2str(sig_freq)])
end

end




