% Compute global coherence and primary eigenvecotrs
%
%   This code implements the global coherence and primary eigenvectors from cross-spectral matrix 
%   described in Cornelissen et al., 2015 eLife. 
%
%   Usage:
%   Main_GlobalCoherence: Main code to compute global coherence and primary eigenvectors using simulated data  
%   GlobalCoherence: Compute multitaper cross-spectral matrix and its singular value decompostion to get the global coherence and primary eigenvectors
%
%   From the paper:
%  "Age-dependent electroencephalogram (EEG) patterns during sevoflurane general anesthesia in infants"
%   Cornelissen, L, Kim, S-E, Purdon, PL, Brown, EN & Berde, CB
%   eLife, 2015
%
%   Copyright 2015 The General Hospital Coporation, authored by Seong-Eun Kim, Ph.D.
%   This work is licensed under a Creative Commons Attribution-NonCommercial-ShareAlike 4.0 International License.
%   (http://creativecommons.org/licenses/by-nc-sa/4.0/)
%
%   Last modified 5/26/2015
%
%**************************************************************************

close all; clear all;
%load simulated data
load('coherence_example.mat');
% adult_data: have a high coherence at frontal areas
% child_data: no coherence anywhere

%smapling rate: 250 Hz
Fs = 250;
%Compute primary eigenvetors for the cross-spectral matrix
[PU_adult,freq, gcoh_adult] = GlobalCoherence(Fs,adult_data); 
[PU_child,freq, gcoh_child] = GlobalCoherence(Fs,child_data); 

%Choose specific frequencies
freqbands = [0,0.5; 3.5,4; 11.5,12; 19.5,20];
freqbands_legend = {'0.5 Hz','4 Hz','12 Hz','20 Hz'};
%% plot
% Compute topoplots using eeglab funtions (Make sure to include eeglab funtions)
figure('units','normalized','position',[0 0 0.7 0.6], 'color','w');
set(gcf,'renderer', 'painters')
for jepoch = 1:2
    for jf = 1: 4,
        if jepoch == 1
            PU = PU_adult;
        else
            PU = PU_child;
        end

        fcurr = (freq>=freqbands(jf,1))&(freq<freqbands(jf,2));
        Ucurr = squeeze(nanmean(abs(PU(:,fcurr)),2));

        [U_headplot, xmesh, ymesh] = surfaceplot_2D_interp(Ucurr);
        
        U_headplot_mean = nanmean(U_headplot,3);
        U_headplot_mean( sqrt(xmesh.^2+ymesh.^2)>0.5) = nan;
        
        if jepoch==1
            subplot(2,4,jf)
        else
            subplot(2,4,jf+4)
        end
        contourf(xmesh,ymesh,U_headplot_mean,[0:.05:1])
        hold on
        PlotHeadEyesEars
        hold off
        axis image off

        caxis([0,0.5])
        if jepoch==1
            if  jf == 3;
                ht = text(-0.8,0.7,'Adult');
                set(ht, 'FontSize', 12, 'FontWeight', 'bold')
            end
            title([freqbands_legend(jf)])
        else
            if  jf == 3;
                ht = text(-0.8,0.7,'Child');
                set(ht, 'FontSize', 12, 'FontWeight', 'bold')
            end
            title([freqbands_legend(jf)])
        end
        shading flat
        
        drawnow
        
    end
    hc = colorbar;
    set(hc,'position',[0.93    0.39    0.013    0.23])
end