addpath(genpath('./support_routines'));
addpath(genpath('./support_tools'));

clear all; close all; clc;

sSubjectName='monkey 2';
sVisualArea='V1V4';

subjectInfos=getSubjectInfos(sSubjectName);
subjectDataV1=getSubjectData(sSubjectName,'V1','ALL','LFPs','ALL');
subjectDataV4=getSubjectData(sSubjectName,'V4','ALL','LFPs','ALL');

SAMPLING_FREQ=1017.38;
NUM_TIME_SAMPLES=512;
NUM_FREQ_SAMPLES=NUM_TIME_SAMPLES/2+1;
timeAxis=linspace(0,NUM_TIME_SAMPLES/SAMPLING_FREQ,NUM_TIME_SAMPLES);
freqAxis=linspace(0,SAMPLING_FREQ/2,NUM_FREQ_SAMPLES);
freqIndices2Plot=find(freqAxis<=100);

lamLayerTag={'Supra','Granr','Infra'};
gratCondTag={'RF','OUT1','OUT2'};

plotPmtsPenAvgOutMeanPoolLayers=1;

% REMOVE DATA IN TIME WINDOWS OTHER THAN PREFIRSTDIM TO SAVE MEMORY
rmpp=[];
for pp=1:length(subjectDataV1.LfpStruct)
    if ~isempty(subjectDataV1.LfpStruct(pp).Sorted)
        rmfn=fieldnames(subjectDataV1.LfpStruct(pp).Sorted.Full.Supra.RF);
        rmfn(cellfun(@(fnm) strcmpi(fnm,'PreFirstDimDataBiZSc'),rmfn))=[];
        for ll=1:3
            subjectDataV1.LfpStruct(pp).Sorted.Full.(lamLayerTag{ll}).RF=rmfield(subjectDataV1.LfpStruct(pp).Sorted.Full.(lamLayerTag{ll}).RF,rmfn);
            subjectDataV1.LfpStruct(pp).Sorted.Full.(lamLayerTag{ll}).OUT1=rmfield(subjectDataV1.LfpStruct(pp).Sorted.Full.(lamLayerTag{ll}).OUT1,rmfn);
            subjectDataV1.LfpStruct(pp).Sorted.Full.(lamLayerTag{ll}).OUT2=rmfield(subjectDataV1.LfpStruct(pp).Sorted.Full.(lamLayerTag{ll}).OUT2,rmfn);
        end
    else
        rmpp=cat(1,rmpp,pp);
    end
end
subjectDataV1.LfpStruct(rmpp)=[]; subjectDataV1.penInfos(rmpp)=[]; subjectDataV1.penIDs(rmpp)=[];
rmpp=[];
for pp=1:length(subjectDataV4.LfpStruct)
    if ~isempty(subjectDataV4.LfpStruct(pp).Sorted)
        rmfn=fieldnames(subjectDataV4.LfpStruct(pp).Sorted.Full.Supra.RF);
        rmfn(cellfun(@(fnm) strcmpi(fnm,'PreFirstDimDataBiZSc'),rmfn))=[];
        for ll=1:3
            subjectDataV4.LfpStruct(pp).Sorted.Full.(lamLayerTag{ll}).RF=rmfield(subjectDataV4.LfpStruct(pp).Sorted.Full.(lamLayerTag{ll}).RF,rmfn);
            subjectDataV4.LfpStruct(pp).Sorted.Full.(lamLayerTag{ll}).OUT1=rmfield(subjectDataV4.LfpStruct(pp).Sorted.Full.(lamLayerTag{ll}).OUT1,rmfn);
            subjectDataV4.LfpStruct(pp).Sorted.Full.(lamLayerTag{ll}).OUT2=rmfield(subjectDataV4.LfpStruct(pp).Sorted.Full.(lamLayerTag{ll}).OUT2,rmfn);
        end
    else
        rmpp=cat(1,rmpp,pp);
    end
end
subjectDataV4.LfpStruct(rmpp)=[]; subjectDataV4.penInfos(rmpp)=[]; subjectDataV4.penIDs(rmpp)=[];


%% Chronux settings (support_tools) - http://chronux.org/
% "Observed Brain Dynamics", P. Mitra and H. Bokil, Oxford University Press, New York, 2008.
paramsMT=[];
paramsMT.Fs=SAMPLING_FREQ;
paramsMT.tapers = [2 3]; % TW K K=2*TW-1
paramsMT.trialave = 1;
paramsMT.pad=0;

%% FIND PENs W/SAME TRIALS IN V1 AND V4
penIDsV1=arrayfun(@(jj) subjectDataV1.LfpStruct(jj).penNum,1:length(subjectDataV1.LfpStruct));
penIDsV4=arrayfun(@(jj) subjectDataV4.LfpStruct(jj).penNum,1:length(subjectDataV4.LfpStruct));
penIDsV1V4=intersect(penIDsV1,penIDsV4);% penIDsV1(arrayfun(@(jj) any(penIDsV4==penIDsV1(jj)), 1:length(penIDsV1)));
penIndicesV1=arrayfun(@(jj) find(penIDsV1==penIDsV1V4(jj)), 1:length(penIDsV1V4));
penIndicesV4=arrayfun(@(jj) find(penIDsV4==penIDsV1V4(jj)), 1:length(penIDsV1V4));

counter=0;

for pp=1:length(penIDsV1V4)
    for cnd=1:3 % LOOP OVER ATTENTIONAL CONDITIONS
        % FIND TRIALS WITH SAME NLX EVENT IDs FOR V1 AND V4
        currPenCndTrialIDsV1=subjectDataV1.LfpStruct(penIndicesV1(pp)).Sorted.Labels.selectedNlxEvents(remod(ceil(subjectDataV1.LfpStruct(penIndicesV1(pp)).Sorted.Labels.GratCondition/6),3)==cnd);
        currPenCndTrialIDsV4=subjectDataV4.LfpStruct(penIndicesV4(pp)).Sorted.Labels.selectedNlxEvents(remod(ceil(subjectDataV4.LfpStruct(penIndicesV4(pp)).Sorted.Labels.GratCondition/6),3)==cnd);
        currPenCndTrialIDsV1V4=intersect(currPenCndTrialIDsV1,currPenCndTrialIDsV4);
        currPenCndTrialIndicesV1=arrayfun(@(jj) find(currPenCndTrialIDsV1==currPenCndTrialIDsV1V4(jj)), 1:length(currPenCndTrialIDsV1V4));
        currPenCndTrialIndicesV4=arrayfun(@(jj) find(currPenCndTrialIDsV4==currPenCndTrialIDsV1V4(jj)), 1:length(currPenCndTrialIDsV1V4));
        % check that trial IDs actually match
        if ~all(currPenCndTrialIDsV1(currPenCndTrialIndicesV1) == currPenCndTrialIDsV4(currPenCndTrialIndicesV4))
            error('Some of the trials IDs for V1 do not match trial IDs for V4.');
        end
        % FIND TRIALS WITH SAME FIRST DIM DELAY FOR V1 AND V4
        currPenCndFirstDimDelayV1=subjectDataV1.LfpStruct(penIndicesV1(pp)).Sorted.Delays.FirstDim(remod(ceil(subjectDataV1.LfpStruct(penIndicesV1(pp)).Sorted.Labels.GratCondition/6),3)==cnd);
        currPenCndFirstDimDelayV4=subjectDataV4.LfpStruct(penIndicesV4(pp)).Sorted.Delays.FirstDim(remod(ceil(subjectDataV4.LfpStruct(penIndicesV4(pp)).Sorted.Labels.GratCondition/6),3)==cnd);
        if ~all(currPenCndFirstDimDelayV1(currPenCndTrialIndicesV1)==currPenCndFirstDimDelayV4(currPenCndTrialIndicesV4))
            error('FirstDim Delays not matching in V1 and V4 trials');
        end
        % Computing Sp Coh V1
        for ll1=1:3 % LOOP OVER LAMINAR LAYER 1
            for ll2=ll1:3 % LOOP OVER LAMINAR LAYER 2
                currNChll1=size(subjectDataV1.LfpStruct(penIndicesV1(pp)).Sorted.Full.(lamLayerTag{ll1}).(gratCondTag{cnd}).PreFirstDimDataBiZSc,1);
                currNChll2=size(subjectDataV1.LfpStruct(penIndicesV1(pp)).Sorted.Full.(lamLayerTag{ll2}).(gratCondTag{cnd}).PreFirstDimDataBiZSc,1);
                currSpCoh=nan(currNChll1*currNChll2,NUM_FREQ_SAMPLES);
                currSpPhi=nan(currNChll1*currNChll2,NUM_FREQ_SAMPLES);
                if currNChll1>0
                    for ii=1:currNChll1
                        currData1ch=squeeze(subjectDataV1.LfpStruct(penIndicesV1(pp)).Sorted.Full.(lamLayerTag{ll1}).(gratCondTag{cnd}).PreFirstDimDataBiZSc(ii,currPenCndTrialIndicesV1,end-NUM_TIME_SAMPLES+1:end))';
                        for jj=((ii+1)*(ll1==ll2)+1*(ll1~=ll2)):currNChll2
                            currData2ch=squeeze(subjectDataV1.LfpStruct(penIndicesV1(pp)).Sorted.Full.(lamLayerTag{ll2}).(gratCondTag{cnd}).PreFirstDimDataBiZSc(jj,currPenCndTrialIndicesV1,end-NUM_TIME_SAMPLES+1:end))';
                            [currSpCoh((ii-1)*currNChll2+jj,:),currSpPhi((ii-1)*currNChll2+jj,:),~,~,~]=coherencyc(currData1ch,currData2ch,paramsMT);
                        end
                    end
                end
                sCohV1Struct(pp).([lamLayerTag{ll1} 'to' lamLayerTag{ll2}]).(gratCondTag{cnd}).PreFirstDimSpCoh=currSpCoh;
                sCohV1Struct(pp).([lamLayerTag{ll1} 'to' lamLayerTag{ll2}]).(gratCondTag{cnd}).PreFirstDimSpPhi=currSpPhi;
                eval(clearoutmem('currN* currD* currS*')); % saves memory at runtime
            end
        end
        % Computing Sp Coh V4
        for ll1=1:3 % LOOP OVER LAMINAR LAYER 1
            for ll2=ll1:3 % LOOP OVER LAMINAR LAYER 2
                currNChll1=size(subjectDataV4.LfpStruct(penIndicesV4(pp)).Sorted.Full.(lamLayerTag{ll1}).(gratCondTag{cnd}).PreFirstDimDataBiZSc,1);
                currNChll2=size(subjectDataV4.LfpStruct(penIndicesV4(pp)).Sorted.Full.(lamLayerTag{ll2}).(gratCondTag{cnd}).PreFirstDimDataBiZSc,1);
                currSpCoh=nan(currNChll1*currNChll2,NUM_FREQ_SAMPLES);
                currSpPhi=nan(currNChll1*currNChll2,NUM_FREQ_SAMPLES);
                if currNChll1>0
                    for ii=1:currNChll1
                        currData1ch=squeeze(subjectDataV4.LfpStruct(penIndicesV4(pp)).Sorted.Full.(lamLayerTag{ll1}).(gratCondTag{cnd}).PreFirstDimDataBiZSc(ii,currPenCndTrialIndicesV4,end-NUM_TIME_SAMPLES+1:end))';
                        for jj=((ii+1)*(ll1==ll2)+1*(ll1~=ll2)):currNChll2
                            currData2ch=squeeze(subjectDataV4.LfpStruct(penIndicesV4(pp)).Sorted.Full.(lamLayerTag{ll2}).(gratCondTag{cnd}).PreFirstDimDataBiZSc(jj,currPenCndTrialIndicesV4,end-NUM_TIME_SAMPLES+1:end))';
                            [currSpCoh((ii-1)*currNChll2+jj,:),currSpPhi((ii-1)*currNChll2+jj,:),~,~,~]=coherencyc(currData1ch,currData2ch,paramsMT);
                        end
                    end
                end
                sCohV4Struct(pp).([lamLayerTag{ll1} 'to' lamLayerTag{ll2}]).(gratCondTag{cnd}).PreFirstDimSpCoh=currSpCoh;
                sCohV4Struct(pp).([lamLayerTag{ll1} 'to' lamLayerTag{ll2}]).(gratCondTag{cnd}).PreFirstDimSpPhi=currSpPhi;
                eval(clearoutmem('currN* currD* currS*')); % saves memory at runtime
            end
        end
        % Computing Sp Coh V1V4
        for ll1=1:3 % LOOP OVER LAMINAR LAYER 1
            for ll2=ll1:3 % LOOP OVER LAMINAR LAYER 2
                currNChll1=size(subjectDataV1.LfpStruct(penIndicesV1(pp)).Sorted.Full.(lamLayerTag{ll1}).(gratCondTag{cnd}).PreFirstDimDataBiZSc,1);
                currNChll2=size(subjectDataV4.LfpStruct(penIndicesV4(pp)).Sorted.Full.(lamLayerTag{ll2}).(gratCondTag{cnd}).PreFirstDimDataBiZSc,1);
                currSpCoh=nan(currNChll1*currNChll2,NUM_FREQ_SAMPLES);
                currSpPhi=nan(currNChll1*currNChll2,NUM_FREQ_SAMPLES);
                if currNChll1>0 && currNChll2>0
                    for ii=1:currNChll1
                        currData1ch=squeeze(subjectDataV1.LfpStruct(penIndicesV1(pp)).Sorted.Full.(lamLayerTag{ll1}).(gratCondTag{cnd}).PreFirstDimDataBiZSc(ii,currPenCndTrialIndicesV1,end-NUM_TIME_SAMPLES+1:end))';
                        for jj=1:currNChll2
                            currData2ch=squeeze(subjectDataV4.LfpStruct(penIndicesV4(pp)).Sorted.Full.(lamLayerTag{ll2}).(gratCondTag{cnd}).PreFirstDimDataBiZSc(jj,currPenCndTrialIndicesV4,end-NUM_TIME_SAMPLES+1:end))';
                            if size(currData1ch,2)~=size(currData2ch,2)
                                warning(['#trials for V1 (' num2str(size(currData1ch,2)) ') do not match #trials for V4 (' num2str(size(currData2ch,2)) ')']);
                            end
                            [currSpCoh((ii-1)*currNChll2+jj,:),currSpPhi((ii-1)*currNChll2+jj,:),~,~,~]=coherencyc(currData1ch,currData2ch,paramsMT);
                        end
                    end
                end
                sCohV1V4Struct(pp).([lamLayerTag{ll1} 'V1to' lamLayerTag{ll2} 'V4']).(gratCondTag{cnd}).PreFirstDimSpCoh=currSpCoh;
                sCohV1V4Struct(pp).([lamLayerTag{ll1} 'V1to' lamLayerTag{ll2} 'V4']).(gratCondTag{cnd}).PreFirstDimSpPhi=currSpPhi;
                eval(clearoutmem('currN* currD* currS*')); % saves memory at runtime
            end
        end
    end
end


%% Pool Sp. Coh. Magnitude / Phase across sessions
PoolSpCohStructV1.AllLayers.RF.PreFirstDimSpCoh=[];
PoolSpCohStructV1.AllLayers.OUT1.PreFirstDimSpCoh=[];
PoolSpCohStructV1.AllLayers.OUT2.PreFirstDimSpCoh=[];
PoolSpCohStructV4.AllLayers.RF.PreFirstDimSpCoh=[];
PoolSpCohStructV4.AllLayers.OUT1.PreFirstDimSpCoh=[];
PoolSpCohStructV4.AllLayers.OUT2.PreFirstDimSpCoh=[];
PoolSpCohStructV1V4.AllLayers.RF.PreFirstDimSpCoh=[];
PoolSpCohStructV1V4.AllLayers.OUT1.PreFirstDimSpCoh=[];
PoolSpCohStructV1V4.AllLayers.OUT2.PreFirstDimSpCoh=[];
for ll1=1:3
    for ll2=ll1:3
        for cnd=1:3
            currPoolPreFirstDimSpCoh1=[];
            currPoolPreFirstDimSpPhi1=[];
            currPoolPreFirstDimSpCoh2=[];
            currPoolPreFirstDimSpPhi2=[];
            currPoolPreFirstDimSpCoh3=[];
            currPoolPreFirstDimSpPhi3=[];
            for pp=1:length(sCohV1V4Struct)
                currPoolPreFirstDimSpCoh1=cat(1, currPoolPreFirstDimSpCoh1, sCohV1Struct(pp).([lamLayerTag{ll1} 'to' lamLayerTag{ll2}]).(gratCondTag{cnd}).PreFirstDimSpCoh);
                currPoolPreFirstDimSpPhi1=cat(1, currPoolPreFirstDimSpPhi1, sCohV1Struct(pp).([lamLayerTag{ll1} 'to' lamLayerTag{ll2}]).(gratCondTag{cnd}).PreFirstDimSpPhi);
                currPoolPreFirstDimSpCoh2=cat(1, currPoolPreFirstDimSpCoh2, sCohV4Struct(pp).([lamLayerTag{ll1} 'to' lamLayerTag{ll2}]).(gratCondTag{cnd}).PreFirstDimSpCoh);
                currPoolPreFirstDimSpPhi2=cat(1, currPoolPreFirstDimSpPhi2, sCohV4Struct(pp).([lamLayerTag{ll1} 'to' lamLayerTag{ll2}]).(gratCondTag{cnd}).PreFirstDimSpPhi);
                currPoolPreFirstDimSpCoh3=cat(1, currPoolPreFirstDimSpCoh3, sCohV1V4Struct(pp).([lamLayerTag{ll1} 'V1to' lamLayerTag{ll2} 'V4']).(gratCondTag{cnd}).PreFirstDimSpCoh);
                currPoolPreFirstDimSpPhi3=cat(1, currPoolPreFirstDimSpPhi3, sCohV1V4Struct(pp).([lamLayerTag{ll1} 'V1to' lamLayerTag{ll2} 'V4']).(gratCondTag{cnd}).PreFirstDimSpPhi);
            end
            currPoolPreFirstDimSpCoh1(all(isnan(currPoolPreFirstDimSpCoh1')),:)=[];
            currPoolPreFirstDimSpPhi1(all(isnan(currPoolPreFirstDimSpPhi1')),:)=[];
            currPoolPreFirstDimSpCoh2(all(isnan(currPoolPreFirstDimSpCoh2')),:)=[];
            currPoolPreFirstDimSpPhi2(all(isnan(currPoolPreFirstDimSpPhi2')),:)=[];
            currPoolPreFirstDimSpCoh3(all(isnan(currPoolPreFirstDimSpCoh3')),:)=[];
            currPoolPreFirstDimSpPhi3(all(isnan(currPoolPreFirstDimSpPhi3')),:)=[];
            PoolSpCohStructV1.AllLayers.(gratCondTag{cnd}).PreFirstDimSpCoh=cat(1,PoolSpCohStructV1.AllLayers.(gratCondTag{cnd}).PreFirstDimSpCoh,currPoolPreFirstDimSpCoh1);
            PoolSpCohStructV4.AllLayers.(gratCondTag{cnd}).PreFirstDimSpCoh=cat(1,PoolSpCohStructV4.AllLayers.(gratCondTag{cnd}).PreFirstDimSpCoh,currPoolPreFirstDimSpCoh2);
            PoolSpCohStructV1V4.AllLayers.(gratCondTag{cnd}).PreFirstDimSpCoh=cat(1,PoolSpCohStructV1V4.AllLayers.(gratCondTag{cnd}).PreFirstDimSpCoh,currPoolPreFirstDimSpCoh3);
            PoolSpCohStructV1.([lamLayerTag{ll1} 'to' lamLayerTag{ll2}]).(gratCondTag{cnd}).PreFirstDimSpCoh=currPoolPreFirstDimSpCoh1;
            PoolSpCohStructV1.([lamLayerTag{ll1} 'to' lamLayerTag{ll2}]).(gratCondTag{cnd}).PreFirstDimSpPhi=currPoolPreFirstDimSpPhi1;
            PoolSpCohStructV4.([lamLayerTag{ll1} 'to' lamLayerTag{ll2}]).(gratCondTag{cnd}).PreFirstDimSpCoh=currPoolPreFirstDimSpCoh2;
            PoolSpCohStructV4.([lamLayerTag{ll1} 'to' lamLayerTag{ll2}]).(gratCondTag{cnd}).PreFirstDimSpPhi=currPoolPreFirstDimSpPhi2;
            PoolSpCohStructV1V4.([lamLayerTag{ll1} 'V1to' lamLayerTag{ll2} 'V4']).(gratCondTag{cnd}).PreFirstDimSpCoh=currPoolPreFirstDimSpCoh3;
            PoolSpCohStructV1V4.([lamLayerTag{ll1} 'V1to' lamLayerTag{ll2} 'V4']).(gratCondTag{cnd}).PreFirstDimSpPhi=currPoolPreFirstDimSpPhi3;
            eval(clearoutmem('curr*'));  % saves memory at runtime
        end
    end
end

%% PLOT MEAN ACROSS PENS PRE-FIRST-DIM OUT=OUT1/2+OUT2/2, LAYERS POOLED
if plotPmtsPenAvgOutMeanPoolLayers
    figure('units','normalized','outerposition',[0 0 1 .6]);
    subplot(1,3,1);
    currSpectCohRF=PoolSpCohStructV1.AllLayers.(gratCondTag{1}).PreFirstDimSpCoh;
    currSpectCohOUT=[nanmean(cat(4,PoolSpCohStructV1.AllLayers.(gratCondTag{2}).PreFirstDimSpCoh,...
        PoolSpCohStructV1.AllLayers.(gratCondTag{3}).PreFirstDimSpCoh),4)];
    plotmsem(freqAxis,currSpectCohRF,'r');
    plotmsem(freqAxis,currSpectCohOUT,'b');
    xlim([0 100]);
    xlabel('Frequency [Hz]'); 
    pValuesVec=nan(1,sum(freqAxis<100));
    for ff=1:sum(freqAxis<100)
        pValuesVec(ff)=signrank(currSpectCohRF(:,ff),currSpectCohOUT(:,ff));
    end
    [~,~,~,pValuesVecFDR]=fdr_bh(pValuesVec);
    pValuesBar05FDR=nan(1,sum(freqAxis<=100));
    pValuesBar05FDR(pValuesVecFDR<=.05)=0;
    plot(freqAxis(freqAxis<100),pValuesBar05FDR,'color',[0 0 0],'linewidth',2.5);
    set(gca,'XTick',0:20:100); 
    title('field-field Sp. Coh. within V1 depths');
    
    subplot(1,3,2);
    currSpectCohRF=PoolSpCohStructV4.AllLayers.(gratCondTag{1}).PreFirstDimSpCoh;
    currSpectCohOUT=[nanmean(cat(4,PoolSpCohStructV4.AllLayers.(gratCondTag{2}).PreFirstDimSpCoh,...
        PoolSpCohStructV4.AllLayers.(gratCondTag{3}).PreFirstDimSpCoh),4)];
    plotmsem(freqAxis,currSpectCohRF,'r');
    plotmsem(freqAxis,currSpectCohOUT,'b');
    xlim([0 100]);
    xlabel('Frequency [Hz]'); 
    pValuesVec=nan(1,sum(freqAxis<=100));
    for ff=1:sum(freqAxis<100)
        pValuesVec(ff)=signrank(currSpectCohRF(:,ff),currSpectCohOUT(:,ff));
    end
    [~,~,~,pValuesVecFDR]=fdr_bh(pValuesVec);
    pValuesBar05FDR=nan(1,sum(freqAxis<100));
    pValuesBar05FDR(pValuesVecFDR<=.05)=0;
    plot(freqAxis(freqAxis<100),pValuesBar05FDR,'color',[0 0 0],'linewidth',2.5);
    set(gca,'XTick',0:20:100); 
    title('field-field Sp. Coh. within V4 depths');
    
    subplot(1,3,3);
    currSpectCohRF=PoolSpCohStructV1V4.AllLayers.(gratCondTag{1}).PreFirstDimSpCoh;
    currSpectCohOUT=[nanmean(cat(4,PoolSpCohStructV1V4.AllLayers.(gratCondTag{2}).PreFirstDimSpCoh,...
        PoolSpCohStructV1V4.AllLayers.(gratCondTag{3}).PreFirstDimSpCoh),4)];
    plotmsem(freqAxis,currSpectCohRF,'r');
    plotmsem(freqAxis,currSpectCohOUT,'b');
    xlim([0 100]);
    xlabel('Frequency [Hz]'); 
    pValuesVec=nan(1,sum(freqAxis<=100));
    for ff=1:sum(freqAxis<100)
        pValuesVec(ff)=signrank(currSpectCohRF(:,ff),currSpectCohOUT(:,ff));
    end
    [~,~,~,pValuesVecFDR]=fdr_bh(pValuesVec);
    pValuesBar05FDR=nan(1,sum(freqAxis<100));
    pValuesBar05FDR(pValuesVecFDR<=.05)=0;
    plot(freqAxis(freqAxis<100),pValuesBar05FDR,'color',[0 0 0],'linewidth',2.5);
    set(gca,'XTick',0:20:100); 
    title('field-field Sp. Coh. V1-V4');
    supertitle([sSubjectName ' - Avg across PENs - LFPs Pre First Dim Spectral Coherence']);
end