V1-V4-LFPs-and-Visual-Attention/Code/process_Spectral_Power_V1V4_AllTimeWins.m

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

clear all; close all; clc;

sSubjectName='monkey 2';
sVisualArea='V1';
subjectInfos=getSubjectInfos(sSubjectName);
subjectData=getSubjectData(sSubjectName,sVisualArea,'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);

lamLayerTag={'Supra','Granr','Infra'};
gratCondTag={'RF','OUT1','OUT2'};
timeWinTag={'PostStim','Stationary','PostCue','PreFirstDim','PreSecondDim','PreThirdDim'};

plotPmtsPenAvg=0;
plotPmtsPenAvgOutMean=1;
plotPmtsPenAvgOutMeanPreFirstDimOverlayed=1;
plotPmtsPenAvgOutMeanPoolLayersPreFirstDimOverlayed=1;

for reportPmtMinusMean=1%0:1
    % reportPmtMinusMean allows to choose how to apply baseline normalization
    % reportPmtMinusMean = 0 reports (Pmt) ./ (MeanPmtbaseline)
    % reportPmtMinusMean = 1 reports (Pmt - MeanPmtbaseline)./(stdPmtbaseline)
    
    %% 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];
    paramsMt.pad = 0;
    
    pr=1;
    for pp=1:length(subjectData.penIDs) % LOOP OVER PENS
        if ~isempty(subjectData.LfpStruct(pp).Sorted)
            for ll=1:3 % LOOP OVER LAMINAR LAYERS
                for cnd=1:3 % LOOP OVER GRAT CONDITIONS
                    currLFPsPreStim=subjectData.LfpStruct(pp).Sorted.Full.(lamLayerTag{ll}).(gratCondTag{cnd}).PreStimDataBi(:,:,end-NUM_TIME_SAMPLES+1:end);
                    currLFPsPostStim=subjectData.LfpStruct(pp).Sorted.Full.(lamLayerTag{ll}).(gratCondTag{cnd}).PostStimDataBi(:,:,1:NUM_TIME_SAMPLES);
                    currLFPsStationary=subjectData.LfpStruct(pp).Sorted.Full.(lamLayerTag{ll}).(gratCondTag{cnd}).StationaryDataBi(:,:,1:NUM_TIME_SAMPLES);
                    currLFPsPostCue=subjectData.LfpStruct(pp).Sorted.Full.(lamLayerTag{ll}).(gratCondTag{cnd}).PostCueDataBi(:,:,1:NUM_TIME_SAMPLES);
                    currLFPsPreFirstDim=subjectData.LfpStruct(pp).Sorted.Full.(lamLayerTag{ll}).(gratCondTag{cnd}).PreFirstDimDataBi(:,:,end-NUM_TIME_SAMPLES+1:end);
                    
                    % SELECT GRATCONDS 2, 3, 4, 5, 9, 11, 15, 17, 20, 21, 22, 23, 27, 29, 33, 35
                    % This stage is necessary to ensure that time windows before dimmings do not include
                    % trials with changes in power due to dimming of stimuli / contrast changes at RF location
                    if ~isempty(subjectData.LfpStruct(pp).Sorted.Full.(lamLayerTag{ll}).(gratCondTag{cnd}).PreSecondDimDataBi)
                        currNumTrials2=length(subjectData.LfpStruct(pp).Sorted.Full.(lamLayerTag{ll}).(gratCondTag{cnd}).gratCondSecondDim);
                        currCorrectSecondDimTrials=find(arrayfun(@(jj) any([2 3 4 5 20 21 22 23 9 11 15 17 27 29 33 35]==subjectData.LfpStruct(pp).Sorted.Full.(lamLayerTag{ll}).(gratCondTag{cnd}).gratCondSecondDim(jj)), 1:currNumTrials2));
                        currLFPsPreSecondDim=subjectData.LfpStruct(pp).Sorted.Full.(lamLayerTag{ll}).(gratCondTag{cnd}).PreSecondDimDataBi(:,currCorrectSecondDimTrials,end-NUM_TIME_SAMPLES+1:end);
                    else
                        currLFPsPreSecondDim=[];
                    end
                    
                    % Note by the time of THIRD DIM the contrast is fixed in the RF iff the trial is ATTEND RF
                    if ~isempty(subjectData.LfpStruct(pp).Sorted.Full.(lamLayerTag{ll}).(gratCondTag{cnd}).PreThirdDimDataBi)
                        currNumTrials3=length(subjectData.LfpStruct(pp).Sorted.Full.(lamLayerTag{ll}).(gratCondTag{cnd}).gratCondThirdDim);
                        currCorrectThirdDimTrials=find(arrayfun(@(jj) any([2 3 4 5 20 21 22 23 9 11 15 17 27 29 33 35]==subjectData.LfpStruct(pp).Sorted.Full.(lamLayerTag{ll}).(gratCondTag{cnd}).gratCondThirdDim(jj)), 1:currNumTrials3));
                        currLFPsPreThirdDim=subjectData.LfpStruct(pp).Sorted.Full.(lamLayerTag{ll}).(gratCondTag{cnd}).PreThirdDimDataBi(:,currCorrectThirdDimTrials,end-NUM_TIME_SAMPLES+1:end);
                    else
                        currLFPsPreThirdDim=[];
                    end
                    
                    currNumChs=size(currLFPsPreStim,1);  currNumTrials=size(currLFPsPreFirstDim,2);
                    currNumTrials2=size(currLFPsPreSecondDim,2); currNumTrials3=size(currLFPsPreThirdDim,2);
                    
                    currPmtPrSt=nan(currNumChs,currNumTrials,NUM_FREQ_SAMPLES);
                    currPmtPsSt=nan(currNumChs,currNumTrials,NUM_FREQ_SAMPLES);
                    currPmtStat=nan(currNumChs,currNumTrials,NUM_FREQ_SAMPLES);
                    currPmtPsCu=nan(currNumChs,currNumTrials,NUM_FREQ_SAMPLES);
                    currPmtPFDm=nan(currNumChs,currNumTrials,NUM_FREQ_SAMPLES);
                    
                    currPmtPostStimZSc=nan(currNumChs,currNumTrials,NUM_FREQ_SAMPLES);
                    currPmtStationaryZSc=nan(currNumChs,currNumTrials,NUM_FREQ_SAMPLES);
                    currPmtPostCueZSc=nan(currNumChs,currNumTrials,NUM_FREQ_SAMPLES);
                    currPmtPreFirstDimZSc=nan(currNumChs,currNumTrials,NUM_FREQ_SAMPLES);
                    currPmtPreSecondDimZSc=nan(currNumChs,currNumTrials2,NUM_FREQ_SAMPLES);
                    currPmtPreThirdDimZSc=nan(currNumChs,currNumTrials3,NUM_FREQ_SAMPLES);
                    
                    for ch=1:currNumChs
                        % FULL BAND CASE: Pmt of LFPs ZScored in frequency domain
                        currPmtPreStim=mtspectrumc(squeeze(currLFPsPreStim(ch,:,:))',paramsMT)';
                        currPmtPostStim=mtspectrumc(squeeze(currLFPsPostStim(ch,:,:))',paramsMT)';
                        currPmtStationary=mtspectrumc(squeeze(currLFPsStationary(ch,:,:))',paramsMT)';
                        currPmtPostCue=mtspectrumc(squeeze(currLFPsPostCue(ch,:,:))',paramsMT)';
                        currPmtPreFirstDim=mtspectrumc(squeeze(currLFPsPreFirstDim(ch,:,:))',paramsMT)';
                        currPmtPreSecondDim=mtspectrumc(squeeze(currLFPsPreSecondDim(ch,:,:))',paramsMT)';
                        currPmtPreThirdDim=mtspectrumc(squeeze(currLFPsPreThirdDim(ch,:,:))',paramsMT)';
                        currPmtPreStimRepMean=repmat(nanmean(currPmtPreStim,1),currNumTrials,1);
                        currPmtPreStimRepStd=repmat(nanstd(currPmtPreStim,[],1),currNumTrials,1);
                        
                        if reportPmtMinusMean
                            currPmtPostStimZSc(ch,:,:)=(currPmtPostStim-currPmtPreStimRepMean)./(currPmtPreStimRepStd+eps);
                            currPmtStationaryZSc(ch,:,:)=(currPmtStationary-currPmtPreStimRepMean)./(currPmtPreStimRepStd+eps);
                            currPmtPostCueZSc(ch,:,:)=(currPmtPostCue-currPmtPreStimRepMean)./(currPmtPreStimRepStd+eps);
                            currPmtPreFirstDimZSc(ch,:,:)=(currPmtPreFirstDim-currPmtPreStimRepMean)./(currPmtPreStimRepStd+eps);
                            currPmtPreSecondDimZSc(ch,:,:)=(currPmtPreSecondDim-currPmtPreStimRepMean(1:currNumTrials2,:))./(currPmtPreStimRepStd(1:currNumTrials2,:)+eps);
                            currPmtPreThirdDimZSc(ch,:,:)=(currPmtPreThirdDim-currPmtPreStimRepMean(1:currNumTrials3,:))./(currPmtPreStimRepStd(1:currNumTrials3,:)+eps);
                        else
                            currPmtPostStimZSc(ch,:,:)=(currPmtPostStim)./(currPmtPreStimRepMean+eps);
                            currPmtStationaryZSc(ch,:,:)=(currPmtStationary)./(currPmtPreStimRepMean+eps);
                            currPmtPostCueZSc(ch,:,:)=(currPmtPostCue)./(currPmtPreStimRepMean+eps);
                            currPmtPreFirstDimZSc(ch,:,:)=(currPmtPreFirstDim)./(currPmtPreStimRepMean+eps);
                            currPmtPreSecondDimZSc(ch,:,:)=(currPmtPreSecondDim(1:currNumTrials2,:))./(currPmtPreStimRepMean(1:currNumTrials2,:)+eps);
                            currPmtPreThirdDimZSc(ch,:,:)=(currPmtPreThirdDim(1:currNumTrials3,:))./(currPmtPreStimRepMean(1:currNumTrials3,:)+eps);
                        end
                        currPmtPrSt(ch,:,:)=currPmtPreStim;
                        currPmtPsSt(ch,:,:)=currPmtPostStim;
                        currPmtStat(ch,:,:)=currPmtStationary;
                        currPmtPsCu(ch,:,:)=currPmtPostCue;
                        currPmtPFDm(ch,:,:)=currPmtPreFirstDim;
                    end
                    PmtStruct(pr).Sorted.Full.(lamLayerTag{ll}).(gratCondTag{cnd}).PostStimZSc=permute(nanmean(currPmtPostStimZSc,2),[1 3 2]);
                    PmtStruct(pr).Sorted.Full.(lamLayerTag{ll}).(gratCondTag{cnd}).StationaryZSc=permute(nanmean(currPmtStationaryZSc,2),[1 3 2]);
                    PmtStruct(pr).Sorted.Full.(lamLayerTag{ll}).(gratCondTag{cnd}).PostCueZSc=permute(nanmean(currPmtPostCueZSc,2),[1 3 2]);
                    PmtStruct(pr).Sorted.Full.(lamLayerTag{ll}).(gratCondTag{cnd}).PreFirstDimZSc=permute(nanmean(currPmtPreFirstDimZSc,2),[1 3 2]);
                    PmtStruct(pr).Sorted.Full.(lamLayerTag{ll}).(gratCondTag{cnd}).PreSecondDimZSc=permute(nanmean(currPmtPreSecondDimZSc,2),[1 3 2]);
                    PmtStruct(pr).Sorted.Full.(lamLayerTag{ll}).(gratCondTag{cnd}).PreThirdDimZSc=permute(nanmean(currPmtPreThirdDimZSc,2),[1 3 2]);
                    
                    PmtStruct(pr).Sorted.Full.(lamLayerTag{ll}).(gratCondTag{cnd}).PrStim=permute(nanmean(currPmtPrSt,2),[1 3 2]);
                    PmtStruct(pr).Sorted.Full.(lamLayerTag{ll}).(gratCondTag{cnd}).PsStim=permute(nanmean(currPmtPsSt,2),[1 3 2]);
                    PmtStruct(pr).Sorted.Full.(lamLayerTag{ll}).(gratCondTag{cnd}).Statry=permute(nanmean(currPmtStat,2),[1 3 2]);
                    PmtStruct(pr).Sorted.Full.(lamLayerTag{ll}).(gratCondTag{cnd}).PstCue=permute(nanmean(currPmtPsCu,2),[1 3 2]);
                    PmtStruct(pr).Sorted.Full.(lamLayerTag{ll}).(gratCondTag{cnd}).PrFDim=permute(nanmean(currPmtPFDm,2),[1 3 2]);
                    
                    PmtStruct(pr).Sorted.Full.(lamLayerTag{ll}).(gratCondTag{cnd}).(subjectData.penInfos(pp).grcColorOrder).PostStimZSc=permute(nanmean(currPmtPostStimZSc,2),[1 3 2]);
                    PmtStruct(pr).Sorted.Full.(lamLayerTag{ll}).(gratCondTag{cnd}).(subjectData.penInfos(pp).grcColorOrder).StationaryZSc=permute(nanmean(currPmtStationaryZSc,2),[1 3 2]);
                    PmtStruct(pr).Sorted.Full.(lamLayerTag{ll}).(gratCondTag{cnd}).(subjectData.penInfos(pp).grcColorOrder).PostCueZSc=permute(nanmean(currPmtPostCueZSc,2),[1 3 2]);
                    PmtStruct(pr).Sorted.Full.(lamLayerTag{ll}).(gratCondTag{cnd}).(subjectData.penInfos(pp).grcColorOrder).PreFirstDimZSc=permute(nanmean(currPmtPreFirstDimZSc,2),[1 3 2]);
                    PmtStruct(pr).Sorted.Full.(lamLayerTag{ll}).(gratCondTag{cnd}).(subjectData.penInfos(pp).grcColorOrder).PreSecondDimZSc=permute(nanmean(currPmtPreSecondDimZSc,2),[1 3 2]);
                    PmtStruct(pr).Sorted.Full.(lamLayerTag{ll}).(gratCondTag{cnd}).(subjectData.penInfos(pp).grcColorOrder).PreThirdDimZSc=permute(nanmean(currPmtPreThirdDimZSc,2),[1 3 2]);
                    clear -regexp ^curr % saves memory at runtime
                end
            end
            pr=pr+1;
        end
    end
    
    %% POOL CHANNELS ACROSS PENs
    for ll=1:3
        for cnd=1:3
            currPmtPostStimZSc=[];
            currPmtStationaryZSc=[];
            currPmtPostCueZSc=[];
            currPmtPreFirstDimZSc=[];
            currPmtPreSecondDimZSc=[];
            currPmtPreThirdDimZSc=[];
            currPmtPreStim=[];
            currPmtPostStim=[];
            currPmtStationary=[];
            currPmtPostCue=[];
            currPmtPreFirstDim=[];
            for pp=1:length(PmtStruct)
                currPmtPreStim=      [currPmtPreStim;         PmtStruct(pp).Sorted.Full.(lamLayerTag{ll}).(gratCondTag{cnd}).PrStim];
                currPmtPostStim=     [currPmtPostStim;        PmtStruct(pp).Sorted.Full.(lamLayerTag{ll}).(gratCondTag{cnd}).PsStim];
                currPmtStationary=   [currPmtStationary;      PmtStruct(pp).Sorted.Full.(lamLayerTag{ll}).(gratCondTag{cnd}).Statry];
                currPmtPostCue=      [currPmtPostCue;         PmtStruct(pp).Sorted.Full.(lamLayerTag{ll}).(gratCondTag{cnd}).PstCue];
                currPmtPreFirstDim=  [currPmtPreFirstDim;     PmtStruct(pp).Sorted.Full.(lamLayerTag{ll}).(gratCondTag{cnd}).PrFDim];
                currPmtPostStimZSc=     [currPmtPostStimZSc;        PmtStruct(pp).Sorted.Full.(lamLayerTag{ll}).(gratCondTag{cnd}).PostStimZSc];
                currPmtStationaryZSc=   [currPmtStationaryZSc;      PmtStruct(pp).Sorted.Full.(lamLayerTag{ll}).(gratCondTag{cnd}).StationaryZSc];
                currPmtPostCueZSc=      [currPmtPostCueZSc;         PmtStruct(pp).Sorted.Full.(lamLayerTag{ll}).(gratCondTag{cnd}).PostCueZSc];
                currPmtPreFirstDimZSc=  [currPmtPreFirstDimZSc;     PmtStruct(pp).Sorted.Full.(lamLayerTag{ll}).(gratCondTag{cnd}).PreFirstDimZSc];
                currPmtPreSecondDimZSc= [currPmtPreSecondDimZSc;    PmtStruct(pp).Sorted.Full.(lamLayerTag{ll}).(gratCondTag{cnd}).PreSecondDimZSc];
                currPmtPreThirdDimZSc=  [currPmtPreThirdDimZSc;     PmtStruct(pp).Sorted.Full.(lamLayerTag{ll}).(gratCondTag{cnd}).PreThirdDimZSc];
            end
            PoolPmtStruct.(lamLayerTag{ll}).(gratCondTag{cnd}).PostStimZSc=currPmtPostStimZSc;
            PoolPmtStruct.(lamLayerTag{ll}).(gratCondTag{cnd}).StationaryZSc=currPmtStationaryZSc;
            PoolPmtStruct.(lamLayerTag{ll}).(gratCondTag{cnd}).PostCueZSc=currPmtPostCueZSc;
            PoolPmtStruct.(lamLayerTag{ll}).(gratCondTag{cnd}).PreFirstDimZSc=currPmtPreFirstDimZSc;
            PoolPmtStruct.(lamLayerTag{ll}).(gratCondTag{cnd}).PreSecondDimZSc=currPmtPreSecondDimZSc;
            PoolPmtStruct.(lamLayerTag{ll}).(gratCondTag{cnd}).PreThirdDimZSc=currPmtPreThirdDimZSc;
            PoolPmtStruct.(lamLayerTag{ll}).(gratCondTag{cnd}).PrStim=currPmtPreStim;
            PoolPmtStruct.(lamLayerTag{ll}).(gratCondTag{cnd}).PsStim=currPmtPostStim;
            PoolPmtStruct.(lamLayerTag{ll}).(gratCondTag{cnd}).Statry=currPmtStationary;
            PoolPmtStruct.(lamLayerTag{ll}).(gratCondTag{cnd}).PstCue=currPmtPostCue;
            PoolPmtStruct.(lamLayerTag{ll}).(gratCondTag{cnd}).PrFDim=currPmtPreFirstDim;
            clear -regexp ^curr % saves memory at runtime
        end
    end
    
    %% PLOT MEAN ACROSS PENS
    if plotPmtsPenAvg
        figure('units','normalized','outerposition',[0 0 1 1]);
        for tt=1:6
            for ll=1:3
                currSpectPowRF=PoolPmtStruct.(lamLayerTag{ll}).RF.([timeWinTag{tt} 'ZSc']);
                currSpectPowOUT=nanmean(cat(4,PoolPmtStruct.(lamLayerTag{ll}).OUT1.([timeWinTag{tt} 'ZSc']),...
                    PoolPmtStruct.(lamLayerTag{ll}).OUT2.([timeWinTag{tt} 'ZSc'])),4);
                numPooledChs=size(currSpectPowRF,1);
                subplot(3,6,(ll-1)*6+tt);
                plotcmapdots([eye(3); .5 .5 .5]);
                if reportPmtMinusMean
                    plot(freqAxis,zeros(1,length(freqAxis)),'k');
                end
                for cnd=1:(1+2*double(tt<6)) % Pre3rdDim is just RF
                    plotmsem(freqAxis,PoolPmtStruct.(lamLayerTag{ll}).(gratCondTag{cnd}).([timeWinTag{tt} 'ZSc']),double([1 2 3]==cnd)); hold on;
                end
                title([timeWinTag{tt} ' - ' lamLayerTag{ll} ' Layer']);
                xlim([0 100]);
                if strcmp(sVisualArea,'V1')
                    if reportPmtMinusMean; ylim([-1 8]); else; ylim([0 8]); end
                else
                    if reportPmtMinusMean; ylim([-1 2.5]); else; ylim([0 2.5]); end
                end
                if ll==3; xlabel('Frequency [Hz]'); end
                if tt==1; ylabel('Baseline-Normalized Power');
                    if strcmp(sVisualArea,'V1')
                        text(3,7.5,['n=' num2str(numPooledChs) ';']);
                    else
                        text(3,2.3,['n=' num2str(numPooledChs) ';']);
                    end
                end
                if tt<6
                    pValuesVec=nan(1,sum(freqAxis<100));
                    for ff=1:sum(freqAxis<100)
                        pValuesVec(ff)=signrank(currSpectPowRF(:,ff),currSpectPowOUT(:,ff));
                    end
                    [~,~,~,pValuesVecFDR]=fdr_bh(pValuesVec);
                    pValuesBar05FDR=nan(1,sum(freqAxis<100));
                    if reportPmtMinusMean
                        pValuesBar05FDR(pValuesVecFDR<=.05)=-1;
                    else
                        pValuesBar05FDR(pValuesVecFDR<=.05)=0;
                    end
                    plot(freqAxis(freqAxis<100),pValuesBar05FDR,'color',[.5 .5 .5],'linewidth',2.5);
                end
                set(gca,'XTick',0:20:100);%set(gca,'XTickLabel',{'0','10','20','30','40','50','60','70','80','90','100'});
            end
        end
        legend('RF','OUT1','OUT2','p\leq 0.05');
        supertitle([sSubjectName ' ' sVisualArea ' - Avg across PENs - LFPs Spectral Power ' spectWinTag{bb} ' Band']);
    end
    
    
    %% PLOT MEAN ACROSS PENS OUT=OUT1/2+OUT2/2
    if plotPmtsPenAvgOutMean
        figure('units','normalized','outerposition',[0 0 1 1]);
        for tt=1:6
            for ll=1:3
                currSpectPowRF=PoolPmtStruct.(lamLayerTag{ll}).RF.([timeWinTag{tt} 'ZSc']);
                currSpectPowOUT=nanmean(cat(4,PoolPmtStruct.(lamLayerTag{ll}).OUT1.([timeWinTag{tt} 'ZSc']),...
                    PoolPmtStruct.(lamLayerTag{ll}).OUT2.([timeWinTag{tt} 'ZSc'])),4);
                numPooledChs=size(currSpectPowRF,1);
                subplot(3,6,(ll-1)*6+tt);
                plotcmapdots([1 0 0; 0 0 1; .5 .5 .5]);
                if reportPmtMinusMean
                    plot(freqAxis,zeros(1,length(freqAxis)),'k');
                end
                plotmsem(freqAxis,currSpectPowRF,'r');
                if tt<6
                    plotmsem(freqAxis,currSpectPowOUT,'b');
                end
                title([timeWinTag{tt} ' - ' lamLayerTag{ll} ' Layer']);
                xlim([0 100]);
                if strcmp(sVisualArea,'V1')
                    if reportPmtMinusMean; ylim([-1 8]); else; ylim([0 8]); end
                else
                    if reportPmtMinusMean; ylim([-1 2.5]); else; ylim([0 2.5]); end
                end
                if ll==3; xlabel('Frequency [Hz]'); end
                if tt==1; ylabel('Baseline-Normalized Power');
                    if strcmp(sVisualArea,'V1')
                        text(3,7.5,['n=' num2str(numPooledChs) ';']);
                    else
                        text(3,2.3,['n=' num2str(numPooledChs) ';']);
                    end
                end
                if tt<6
                    pValuesVec=nan(1,sum(freqAxis<100));
                    for ff=1:sum(freqAxis<100)
                        pValuesVec(ff)=signrank(currSpectPowRF(:,ff),currSpectPowOUT(:,ff));
                    end
                    [~,~,~,pValuesVecFDR]=fdr_bh(pValuesVec);
                    pValuesBar05FDR=nan(1,sum(freqAxis<100));
                    if reportPmtMinusMean
                        pValuesBar05FDR(pValuesVecFDR<=.05)=-1;
                    else
                        pValuesBar05FDR(pValuesVecFDR<=.05)=0;
                    end
                    plot(freqAxis(freqAxis<100),pValuesBar05FDR,'color',[0 0 0],'linewidth',2.5);
                end
                set(gca,'XTick',0:20:100);
            end
        end
        legend('RF','OUT','p\leq 0.05');
        supertitle([sSubjectName ' ' sVisualArea ' - Avg across PENs - LFPs Spectral Power']);
    end
    
    
    %% PLOT MEAN ACROSS PENS OVERLAY FIRST DIM OUT=OUT1/2+OUT2/2, PREFIRSTDIM OVERLAYED (DASHED)
    if plotPmtsPenAvgOutMeanPreFirstDimOverlayed
        figure('units','normalized','outerposition',[0 0 1 1]);
        for tt=1:3
            for ll=1:3
                currSpectPowRF=PoolPmtStruct.(lamLayerTag{ll}).RF.([timeWinTag{tt} 'ZSc']);
                currSpectPowOUT=nanmean(cat(4,PoolPmtStruct.(lamLayerTag{ll}).OUT1.([timeWinTag{tt} 'ZSc']),...
                    PoolPmtStruct.(lamLayerTag{ll}).OUT2.([timeWinTag{tt} 'ZSc'])),4);
                numPooledChs=size(currSpectPowRF,1);
                subplot(3,3,(ll-1)*3+tt);
                plotcmapdots([1 0 0; 0 0 1; .5 .5 .5]);
                if reportPmtMinusMean; plot(freqAxis,zeros(1,length(freqAxis)),'k'); end
                plotmsem(freqAxis,currSpectPowRF,'r');
                if tt<6; plotmsem(freqAxis,currSpectPowOUT,'b'); end
                currSpectPowPFDRF=PoolPmtStruct.(lamLayerTag{ll}).RF.PreFirstDimZSc;
                currSpectPowPFDOUT=nanmean(cat(4,PoolPmtStruct.(lamLayerTag{ll}).OUT1.PreFirstDimZSc,...
                    PoolPmtStruct.(lamLayerTag{ll}).OUT2.PreFirstDimZSc),4);
                plotmsem(freqAxis,currSpectPowPFDRF,'r:');
                plotmsem(freqAxis,currSpectPowPFDOUT,'b:');
                title([timeWinTag{tt} ' - ' lamLayerTag{ll} ' Layer']);
                xlim([0 100]);
                if strcmp(sVisualArea,'V1')
                    if reportPmtMinusMean
                        if strcmpi(sSubjectName,'Wyman'); ylim([-1 8]); else; ylim([-1 4]); end
                    else
                        ylim([0 8]);
                    end
                else
                    if reportPmtMinusMean; ylim([-1 2.5]); else; ylim([0 2.5]); end
                end
                if ll==3; xlabel('Frequency [Hz]'); end
                if tt==1; ylabel('Baseline-Normalized Power');
                    if strcmp(sVisualArea,'V1')
                        if strcmpi(sSubjectName,'Wyman')
                            text(3,7.5,['n=' num2str(numPooledChs) ';']);
                        else
                            text(3,3.5,['n=' num2str(numPooledChs) ';']);
                        end
                    else
                        text(3,2.3,['n=' num2str(numPooledChs) ';']);
                    end
                end
                if tt<6
                    pValuesVec=nan(1,sum(freqAxis<100));
                    for ff=1:sum(freqAxis<100)
                        pValuesVec(ff)=signrank(currSpectPowRF(:,ff),currSpectPowOUT(:,ff));
                    end
                    [~,~,~,pValuesVecFDR]=fdr_bh(pValuesVec);
                    pValuesBar05FDR=nan(1,sum(freqAxis<100));
                    if reportPmtMinusMean
                        pValuesBar05FDR(pValuesVecFDR<=.05)=-1;
                    else
                        pValuesBar05FDR(pValuesVecFDR<=.05)=0;
                    end
                    plot(freqAxis(freqAxis<100),pValuesBar05FDR,'color',[0 0 0],'linewidth',2.5);
                end
                set(gca,'XTick',0:20:100);%set(gca,'XTickLabel',{'0','10','20','30','40','50','60','70','80','90','100'});
            end
        end
        legend('RF','OUT','p\leq 0.05');
        supertitle([sSubjectName ' ' sVisualArea ' - Avg across PENs - LFPs Spectral Power']);
    end
    
    
    %% PLOT MEAN ACROSS PENS PRE-FIRST-DIM OUT=OUT1/2+OUT2/2, LAYERS POOLED, POSTSTIM OVERLAYED (DASHED)
    if plotPmtsPenAvgOutMeanPoolLayersPreFirstDimOverlayed
        figure('units','normalized','outerposition',[0 0 1 1]);
        tt=2; % Pick Stationary time window
        currSpectPowRF=[PoolPmtStruct.(lamLayerTag{1}).RF.([timeWinTag{tt} 'ZSc']);...
            PoolPmtStruct.(lamLayerTag{2}).RF.([timeWinTag{tt} 'ZSc']);...
            PoolPmtStruct.(lamLayerTag{3}).RF.([timeWinTag{tt} 'ZSc'])];
        currSpectPowOUT=[nanmean(cat(4,PoolPmtStruct.(lamLayerTag{1}).OUT1.([timeWinTag{tt} 'ZSc']),...
            PoolPmtStruct.(lamLayerTag{1}).OUT2.([timeWinTag{tt} 'ZSc'])),4);
            nanmean(cat(4,PoolPmtStruct.(lamLayerTag{2}).OUT1.([timeWinTag{tt} 'ZSc']),...
            PoolPmtStruct.(lamLayerTag{2}).OUT2.([timeWinTag{tt} 'ZSc'])),4);
            nanmean(cat(4,PoolPmtStruct.(lamLayerTag{3}).OUT1.([timeWinTag{tt} 'ZSc']),...
            PoolPmtStruct.(lamLayerTag{3}).OUT2.([timeWinTag{tt} 'ZSc'])),4)];
        numPooledChs=size(currSpectPowRF,1);
        if reportPmtMinusMean
            plot(freqAxis,zeros(1,length(freqAxis)),'k');
        end
        plotmsem(freqAxis,currSpectPowRF,'r:');
        plotmsem(freqAxis,currSpectPowOUT,'b:');
        currSpectPowPFDRF=[PoolPmtStruct.(lamLayerTag{1}).RF.PreFirstDimZSc;...
            PoolPmtStruct.(lamLayerTag{2}).RF.PreFirstDimZSc;...
            PoolPmtStruct.(lamLayerTag{3}).RF.PreFirstDimZSc];
        currSpectPowPFDOUT=[nanmean(cat(4,PoolPmtStruct.(lamLayerTag{1}).OUT1.PreFirstDimZSc,...
            PoolPmtStruct.(lamLayerTag{1}).OUT2.PreFirstDimZSc),4);...
            nanmean(cat(4,PoolPmtStruct.(lamLayerTag{2}).OUT1.PreFirstDimZSc,...
            PoolPmtStruct.(lamLayerTag{2}).OUT2.PreFirstDimZSc),4);...
            nanmean(cat(4,PoolPmtStruct.(lamLayerTag{3}).OUT1.PreFirstDimZSc,...
            PoolPmtStruct.(lamLayerTag{3}).OUT2.PreFirstDimZSc),4)];
        
        plotmsem(freqAxis,currSpectPowPFDRF,'r');
        plotmsem(freqAxis,currSpectPowPFDOUT,'b');
        title([timeWinTag{tt} ' - ' lamLayerTag{ll} ' Layer']);
        xlim([0 100]);
        if strcmp(sVisualArea,'V1')
            if reportPmtMinusMean
                if strcmpi(sSubjectName(end),'1')
                    ylim([-1 8]);
                else
                    ylim([-1 6]);
                end
            else
                ylim([0 8]);
            end
        else
            if reportPmtMinusMean
                ylim([-1 2.5]);
            else
                ylim([0 2.5]);
            end
        end
        if ll==3; xlabel('Frequency [Hz]'); end
        pValuesVec=nan(1,sum(freqAxis<=100));
        for ff=1:sum(freqAxis<100)
            pValuesVec(ff)=signrank(currSpectPowPFDRF(:,ff),currSpectPowPFDOUT(:,ff));
        end
        [~,~,~,pValuesVecFDR]=fdr_bh(pValuesVec);
        pValuesBar05FDR=nan(1,sum(freqAxis<=100));
        if reportPmtMinusMean
            pValuesBar05FDR(pValuesVecFDR<=.05)=-.5;
        else
            pValuesBar05FDR(pValuesVecFDR<=.05)=0;
        end
        plot(freqAxis(freqAxis<100),pValuesBar05FDR,'color',[0 0 0],'linewidth',2.5);
        
        set(gca,'XTick',0:20:100); %set(gca,'XTickLabel',{'0','10','20','30','40','50','60','70','80','90','100'});
        supertitle([sSubjectName ' ' sVisualArea ' - Avg across PENs - LFPs Spectral Power']);
    end
    
end