%%N Vs B, & Class-Specific coherence results for low and high P(Choice) and RPE
%  randomized permutation and and nanparametric stats for class specific
%  and broad vs narrow

%Kia Banaie Boroujeni....2020
%%%-----------------------------%%%

clear
close all
clc,

load('Coherence_Data.mat') % a mat file of coherence measures for each class and each condition in each area

classN={'B5','B4','B3','B2','B1','N3','N2','N1'};

foi=find(freq<80 & freq>4); %freq of interest
Fr=freq(foi);
toi=find(Time<600 & Time>-950); %time of interest
tsm=1;
fsm=7;
W1=gausswin(tsm,3); %smoothing window time
W2=gausswin(fsm,3); %smoothing window frequency
Wx=repmat(W1,1,fsm);
Wy=repmat(W2',tsm,1);
Wt=Wx.*Wy;
NColor=[1 .6 .2];
PreOnset=find(Time<0); %pre-onset time
PostOnset=find(Time>=0); %post-onset time


%Grouping RPEs and classes

N3_hRPE=cat(1,res.RES_V_dis_hRPE{6});
Brd_hRPE=cat(1,res.RES_V_dis_hRPE{1},res.RES_V_dis_hRPE{2},res.RES_V_dis_hRPE{3},res.RES_V_dis_hRPE{4},res.RES_V_dis_hRPE{5});

N3_hChP=cat(1,res.RES_V_dis_hChP{6});
Brd_hChP=cat(1,res.RES_V_dis_hChP{1},res.RES_V_dis_hChP{2},res.RES_V_dis_hChP{3},res.RES_V_dis_hChP{4},res.RES_V_dis_hChP{5});


N3_lRPE=cat(1,res.RES_V_dis_lRPE{6});
Brd_lRPE=cat(1,res.RES_V_dis_lRPE{1},res.RES_V_dis_lRPE{2},res.RES_V_dis_lRPE{3},res.RES_V_dis_lRPE{4},res.RES_V_dis_lRPE{5});

N3_lChP=cat(1,res.RES_V_dis_lChP{6});
Brd_lChP=cat(1,res.RES_V_dis_lChP{1},res.RES_V_dis_lChP{2},res.RES_V_dis_lChP{3},res.RES_V_dis_lChP{4},res.RES_V_dis_lChP{5});

Brd_area = {cat(2, Area2{1:5})};
N3_area = {cat(2, Area2{6:6})};
all_area= {cat(2, Area2{1:6})};

all_lRPE=cat(1,N3_lRPE,Brd_lRPE);
all_hRPE=cat(1,N3_hRPE,Brd_hRPE);
all_lChP=cat(1,N3_lChP,Brd_lChP);
all_hChP=cat(1,N3_hChP,Brd_hChP);

%% Plot 2D mean coh in a time window after the onset
% ACC high RPE
figure
subplot(2,2,1)
ALLID=find(strcmp(all_area{1,1},'ACC'));
N3ID=find(strcmp(N3_area{1,1},'ACC'));
A_RPE=all_hRPE(ALLID,:,:);
L_all=length(ALLID);
L_N=length(N3ID);
for i=1:1000
    perm=randperm(L_all);
    Norm_res=(A_RPE(perm(1:L_N),toi,foi)-nanmean(A_RPE(perm(1:L_N),PreOnset,foi),[2,3]))./nanstd(A_RPE(perm(1:L_N),PreOnset,foi),[],[2,3]);
    Norm_res=squeeze(nanmean(Norm_res,1))';
    mCoh=conv2(Norm_res,Wt,'same')/sum(Wt,'all');
    Rand_Coh(:,i)=nanmean(mCoh(:,PostOnset),2);
end
Norm_res=(N3_hRPE(N3ID,toi,foi)-nanmean(N3_hRPE(N3ID,PreOnset,foi),[2,3]))./nanstd(N3_hRPE(N3ID,PreOnset,foi),[],[2,3]);
Norm_res=squeeze(nanmean(Norm_res,1))';
mCoh=conv2(Norm_res,Wt,'same')/sum(Wt,'all');
N3_Coh=nanmean(mCoh(:,PostOnset),2);
plot(Fr,nanmean(Rand_Coh,2),'b')
pv=prctile(Rand_Coh,95,'all');
hold on
plot(get(gca,'xlim'),[pv pv],'k')
plot(get(gca,'xlim'),[0 0],'linestyle','--','linewidth',1,'color','k')
plot(Fr,N3_Coh,'color',NColor)
set(gca,'ylim',[-1.2 1.2],'tickdir','out','box','off')
title('ACC High RPE')
xlabel 'Freq (Hz)'
ylabel 'Mean Norm. Coh'


% LPFC low p(choice)

ALLID=find((strcmp(all_area{1,1},'8') | strcmp(all_area{1,1},'46')));
N3ID=find((strcmp(N3_area{1,1},'8') | strcmp(N3_area{1,1},'46')));
A_ChP=all_lChP(ALLID,:,:);
L_all=length(ALLID);
L_N=length(N3ID);
subplot(2,2,3)
for i=1:1000
    perm=randperm(L_all);
    Norm_res=(A_ChP(perm(1:L_N),toi,foi)-nanmean(A_ChP(perm(1:L_N),PreOnset,foi),[2,3]))./nanstd(A_ChP(perm(1:L_N),PreOnset,foi),[],[2,3]);
    Norm_res=squeeze(nanmean(Norm_res,1))';
    mCoh=conv2(Norm_res,Wt,'same')/sum(Wt,'all');
    Rand_Coh(:,i)=nanmean(mCoh(:,PostOnset),2);
end
Norm_res=(N3_lChP(N3ID,toi,foi)-nanmean(N3_lChP(N3ID,PreOnset,foi),[2,3]))./nanstd(N3_lChP(N3ID,PreOnset,foi),[],[2,3]);
Norm_res=squeeze(nanmean(Norm_res,1))';
mCoh=conv2(Norm_res,Wt,'same')/sum(Wt,'all');
N3_Coh=nanmean(mCoh(:,PostOnset),2);
plot(Fr,nanmean(Rand_Coh,2),'b')
pv=prctile(Rand_Coh,95,'all');
hold on
plot(get(gca,'xlim'),[pv pv],'k')
plot(Fr,N3_Coh,'color',NColor)
plot(get(gca,'xlim'),[0 0],'linestyle','--','linewidth',1,'color','k')
set(gca,'ylim',[-1.2 1.2],'tickdir','out','box','off')
title('LPFC Low P(choice)')
xlabel 'Freq (Hz)'
ylabel 'Mean Norm. Coh'


% ACC Low RPE
subplot(2,2,2)
ALLID=find(strcmp(all_area{1,1},'ACC'));
N3ID=find(strcmp(N3_area{1,1},'ACC'));
A_RPE=all_lRPE(ALLID,:,:);
L_all=length(ALLID);
L_N=length(N3ID);
for i=1:1000
    perm=randperm(L_all);
    Norm_res=(A_RPE(perm(1:L_N),toi,foi)-nanmean(A_RPE(perm(1:L_N),PreOnset,foi),[2,3]))./nanstd(A_RPE(perm(1:L_N),PreOnset,foi),[],[2,3]);
    Norm_res=squeeze(nanmean(Norm_res,1))';
    mCoh=conv2(Norm_res,Wt,'same')/sum(Wt,'all');
    Rand_Coh(:,i)=nanmean(mCoh(:,PostOnset),2);
end
ALLID=find(strcmp(N3_area{1,1},'ACC'));
Norm_res=(N3_lRPE(N3ID,toi,foi)-nanmean(N3_lRPE(N3ID,PreOnset,foi),[2,3]))./nanstd(N3_lRPE(N3ID,PreOnset,foi),[],[2,3]);
Norm_res=squeeze(nanmean(Norm_res,1))';
mCoh=conv2(Norm_res,Wt,'same')/sum(Wt,'all');
N3_Coh=nanmean(mCoh(:,PostOnset),2);
plot(Fr,nanmean(Rand_Coh,2),'b')
pv=prctile(Rand_Coh,95,'all');
hold on
plot(get(gca,'xlim'),[pv pv],'k')
plot(get(gca,'xlim'),[0 0],'linestyle','--','linewidth',1,'color','k')
plot(Fr,N3_Coh,'color',NColor)
set(gca,'ylim',[-1.2 1.2],'tickdir','out','box','off')
title('ACC Low RPE')
xlabel 'Freq (Hz)'
ylabel 'Mean Norm. Coh'



% LPFC high p(choice)
subplot(2,2,4)
ALLID=find((strcmp(all_area{1,1},'8') | strcmp(all_area{1,1},'46')));
N3ID=find((strcmp(N3_area{1,1},'8') | strcmp(N3_area{1,1},'46')));
A_ChP=all_hChP(ALLID,:,:);
L_all=length(ALLID);
L_N=length(N3ID);
for i=1:1000
    
    perm=randperm(L_all);
    Norm_res=(A_ChP(perm(1:L_N),toi,foi)-nanmean(A_ChP(perm(1:L_N),PreOnset,foi),[2,3]))./nanstd(A_ChP(perm(1:L_N),PreOnset,foi),[],[2,3]);
    Norm_res=squeeze(nanmean(Norm_res,1))';
    mCoh=conv2(Norm_res,Wt,'same')/sum(Wt,'all');
    Rand_Coh(:,i)=nanmean(mCoh(:,PostOnset),2);
end
Norm_res=(N3_hChP(N3ID,toi,foi)-nanmean(N3_hChP(N3ID,PreOnset,foi),[2,3]))./nanstd(N3_hChP(N3ID,PreOnset,foi),[],[2,3]);
Norm_res=squeeze(nanmean(Norm_res,1))';
mCoh=conv2(Norm_res,Wt,'same')/sum(Wt,'all');
N3_Coh=nanmean(mCoh(:,PostOnset),2);
plot(Fr,nanmean(Rand_Coh,2),'b')
pv=prctile(Rand_Coh,95,'all');
hold on
plot(get(gca,'xlim'),[pv pv],'k')
plot(Fr,N3_Coh,'color',NColor)
plot(get(gca,'xlim'),[0 0],'linestyle','--','linewidth',1,'color','k')
set(gca,'ylim',[-1.2 1.2],'tickdir','out','box','off')
title('LPFC High P(choice)')
xlabel 'Freq (Hz)'
ylabel 'Mean Norm. Coh'

set(gcf, 'position',[560   582   669   366])


%% Image plots for N3 / Condition/ Area
figure
colormap jet
subplot(3,4,1)
N3ID=find(strcmp(N3_area{1,1},'ACC'));
Norm_res=(N3_lChP(N3ID,toi,foi)-nanmean(N3_lChP(N3ID,PreOnset,foi),[2,3]))./nanstd(N3_lChP(N3ID,PreOnset,foi),[],[2,3]);
Norm_res=squeeze(nanmean(Norm_res,1))';
mCoh=conv2(Norm_res,Wt,'same')/sum(Wt,'all');

imagesc(Time,Fr,mCoh)
set(gca,'clim',[-1.2 1.2])
hold on
plot([0 0],get(gca,'ylim'),'linewidth',1,'linestyle','--','color','k')
title('ACC N3 Low ChP')
axis xy
subplot(3,4,2)
N3ID=find(strcmp(N3_area{1,1},'ACC'));
Norm_res=(N3_hChP(N3ID,toi,foi)-nanmean(N3_hChP(N3ID,PreOnset,foi),[2,3]))./nanstd(N3_hChP(N3ID,PreOnset,foi),[],[2,3]);
Norm_res=squeeze(nanmean(Norm_res,1))';
mCoh=conv2(Norm_res,Wt,'same')/sum(Wt,'all');
uimagesc(Time,Fr,mCoh)
set(gca,'clim',[-1.2 1.2])
hold on
plot([0 0],get(gca,'ylim'),'linewidth',1,'linestyle','--','color','k')
axis xy
title('ACC N3 High ChP')

subplot(3,4,5)
N3ID=find(strcmp(N3_area{1,1},'ACC'));
Norm_res=(N3_lRPE(N3ID,toi,foi)-nanmean(N3_lRPE(N3ID,PreOnset,foi),[2,3]))./nanstd(N3_lRPE(N3ID,PreOnset,foi),[],[2,3]);
Norm_res=squeeze(nanmean(Norm_res,1))';
mCoh=conv2(Norm_res,Wt,'same')/sum(Wt,'all');
uimagesc(Time,Fr,mCoh)
set(gca,'clim',[-1.2 1.2])
hold on
plot([0 0],get(gca,'ylim'),'linewidth',1,'linestyle','--','color','k')
axis xy
title('ACC N3 Low RPE')
subplot(3,4,6)
NID=find(strcmp(N3_area{1,1},'ACC'));
Norm_res=(N3_hRPE(N3ID,toi,foi)-nanmean(N3_hRPE(N3ID,PreOnset,foi),[2,3]))./nanstd(N3_hRPE(N3ID,PreOnset,foi),[],[2,3]);
Norm_res=squeeze(nanmean(Norm_res,1))';
mCoh=conv2(Norm_res,Wt,'same')/sum(Wt,'all');
uimagesc(Time,Fr,mCoh)
set(gca,'clim',[-1.2 1.2])
hold on
plot([0 0],get(gca,'ylim'),'linewidth',1,'linestyle','--','color','k')
axis xy
title('ACC N3 High RPE')
subplot(3,4,3)
N3ID=find((strcmp(N3_area{1,1},'8') | strcmp(N3_area{1,1},'46')));
Norm_res=(N3_lChP(N3ID,toi,foi)-nanmean(N3_lChP(N3ID,PreOnset,foi),[2,3]))./nanstd(N3_lChP(N3ID,PreOnset,foi),[],[2,3]);
Norm_res=squeeze(nanmean(Norm_res,1))';
mCoh=conv2(Norm_res,Wt,'same')/sum(Wt,'all');
uimagesc(Time,Fr,mCoh)
set(gca,'clim',[-1.2 1.2])
hold on
plot([0 0],get(gca,'ylim'),'linewidth',1,'linestyle','--','color','k')
axis xy
title('PFC N3 Low ChP')
subplot(3,4,4)
N3ID=find((strcmp(N3_area{1,1},'8') | strcmp(N3_area{1,1},'46')));
Norm_res=(N3_hChP(N3ID,toi,foi)-nanmean(N3_hChP(N3ID,PreOnset,foi),[2,3]))./nanstd(N3_hChP(N3ID,PreOnset,foi),[],[2,3]);
Norm_res=squeeze(nanmean(Norm_res,1))';
mCoh=conv2(Norm_res,Wt,'same')/sum(Wt,'all');
uimagesc(Time,Fr,mCoh)
set(gca,'clim',[-1.2 1.2])
hold on
plot([0 0],get(gca,'ylim'),'linewidth',1,'linestyle','--','color','k')
axis xy
title('PFC N3 High ChP')

subplot(3,4,7)
N3ID=find((strcmp(N3_area{1,1},'8') | strcmp(N3_area{1,1},'46')));
Norm_res=(N3_lRPE(N3ID,toi,foi)-nanmean(N3_lRPE(N3ID,PreOnset,foi),[2,3]))./nanstd(N3_lRPE(N3ID,PreOnset,foi),[],[2,3]);
Norm_res=squeeze(nanmean(Norm_res,1))';

mCoh=conv2(Norm_res,Wt,'same')/sum(Wt,'all');
uimagesc(Time,Fr,mCoh)
set(gca,'clim',[-1.2 1.2])
hold on
plot([0 0],get(gca,'ylim'),'linewidth',1,'linestyle','--','color','k')
axis xy
title('PFC N3 Low RPE')

subplot(3,4,8)
N3ID=find((strcmp(N3_area{1,1},'8') | strcmp(N3_area{1,1},'46')));
Norm_res=(N3_hRPE(N3ID,toi,foi)-nanmean(N3_hRPE(N3ID,PreOnset,foi),[2,3]))./nanstd(N3_hRPE(N3ID,PreOnset,foi),[],[2,3]);
Norm_res=squeeze(nanmean(Norm_res,1))';
mCoh=conv2(Norm_res,Wt,'same')/sum(Wt,'all');
uimagesc(Time,Fr,mCoh)
set(gca,'clim',[-1.2 1.2])
hold on
plot([0 0],get(gca,'ylim'),'linewidth',1,'linestyle','--','color','k')
axis xy
title('PFC N3 High RPE')


%% Coherence of classes and Statistics
% High P(choice) LPFC
figure
state=cell2mat(res.RES_V_dis_hChP');
state(isnan(sum(state,[2,3])),:,:)=[];
L1=size(state,1);
Sign=zeros(length(foi),length(toi));
for i=1:8
    subplot(8,4,4*(i-1)+1)
    Cl_Inds=find(strcmp(Area2{1, i},'46')|strcmp(Area2{1, i},'8'));
    L=sum(~isnan(sum(state(Cl_Inds,:,:),[2,3])));
    Sign=zeros(length(foi),length(toi));
    for j=1:1000
        perm = randperm(L1);
        inds = perm(1:L);
        perm_st=state(inds,:,:);
        Norm_perm =(perm_st(:,toi,foi)-nanmean(perm_st(:,PreOnset,foi),[2,3]))./nanstd(perm_st(:,PreOnset,foi),[],[2,3]);
        Norm_perm = squeeze(nanmean(Norm_perm,1))';
        mCoh=conv2(Norm_perm,Wt,'same')/sum(Wt,'all');
        Stat(j,:,:) = mCoh;
    end
    
    Norm_res=(res.RES_V_dis_hChP{1,i}(Cl_Inds,toi,foi)-nanmean(res.RES_V_dis_hChP{1,i}(Cl_Inds,PreOnset,foi),[2,3]))./nanstd(res.RES_V_dis_hChP{1,i}(Cl_Inds,PreOnset,foi),[],[2,3]);
    Norm_res=squeeze(nanmean(Norm_res,1))';
    mCoh=conv2(Norm_res,Wt,'same')/sum(Wt,'all');
    uimagesc(Time,Fr,mCoh)
    up=prctile(Stat,99.5,'all');
    Sign((mCoh>up))=1;
    hold on
    contour(Time,Fr,Sign,1,'k','linewidth',1.5)
    set(gca,'clim',[-1.2 1.2])
    
    plot([0 0],get(gca,'ylim'),'linewidth',1,'linestyle','--','color','k')
    axis xy
    colormap('jet')
    title(['high-p(ch) class-' classN{i}])
    ylabel 'Norm. Coherence'
    if i==8
        xlabel 'Time to Feature (msec)'
    end
end


% Low P(choice) LPFC
state=cell2mat(res.RES_V_dis_lChP');
state(isnan(sum(state,[2,3])),:,:)=[];
L1=size(state,1);
Sign=zeros(length(foi),length(toi));
for i=1:8
    Cl_Inds=find(strcmp(Area2{1, i},'46')|strcmp(Area2{1, i},'8'));
    L=sum(~isnan(sum(res.RES_V_dis_lChP{1,i}(Cl_Inds,:,:),[2,3])));
    Sign=zeros(length(foi),length(toi));
    for j=1:1000
        perm = randperm(L1);
        inds = perm(1:L);
        perm_st=state(inds,:,:);
        Norm_perm =(perm_st(:,toi,foi)-nanmean(perm_st(:,PreOnset,foi),[2,3]))./nanstd(perm_st(:,PreOnset,foi),[],[2,3]);
        Norm_perm = squeeze(nanmean(Norm_perm,1))';
        mCoh=conv2(Norm_perm,Wt,'same')/sum(Wt,'all');
        Stat(j,:,:) = mCoh;
    end
    subplot(8,4,4*(i-1)+2)
    Cl_Inds=find(strcmp(Area2{1, i},'46')|strcmp(Area2{1, i},'8'));
    Norm_res=(res.RES_V_dis_lChP{1,i}(Cl_Inds,toi,foi)-nanmean(res.RES_V_dis_lChP{1,i}(Cl_Inds,PreOnset,foi),[2,3]))./nanstd(res.RES_V_dis_lChP{1,i}(Cl_Inds,PreOnset,foi),[],[2,3]);
    Norm_res=squeeze(nanmean(Norm_res,1))';
    mCoh=conv2(Norm_res,Wt,'same')/sum(Wt,'all');
    uimagesc(Time,Fr,mCoh)
    up=prctile(Stat,99.5,'all');
    Sign((mCoh>up))=1;
    hold on
    contour(Time,Fr,Sign,1,'k','linewidth',1.5)
    set(gca,'clim',[-1.2 1.2])
    hold on
    plot([0 0],get(gca,'ylim'),'linewidth',1,'linestyle','--','color','k')
    axis xy
    colormap('jet')
    title(['low-p(ch) class-' classN{i}])
    ylabel 'Norm. Coherence'
    if i==8
        xlabel 'Time to Feature (msec)'
    end
    
end
% Low RPE ACC
state=cell2mat(res.RES_V_dis_lRPE');
state(isnan(sum(state,[2,3])),:,:)=[];
L1=size(state,1);
Sign=zeros(length(foi),length(toi));
for i=1:8
    Cl_Inds=find(strcmp(Area2{1, i},'ACC'));
    L=sum(~isnan(sum(res.RES_V_dis_lRPE{1,i}(Cl_Inds,:,:),[2,3])));
    Sign=zeros(length(foi),length(toi));
    for j=1:1000
        perm = randperm(L1);
        inds = perm(1:L);
        perm_st=state(inds,:,:);
        Norm_perm =(perm_st(:,toi,foi)-nanmean(perm_st(:,PreOnset,foi),[2,3]))./nanstd(perm_st(:,PreOnset,foi),[],[2,3]);
        Norm_perm = squeeze(nanmean(Norm_perm,1))';
        mCoh=conv2(Norm_perm,Wt,'same')/sum(Wt,'all');
        Stat(j,:,:) = mCoh;
    end
    subplot(8,4,4*(i-1)+3)
    
    Norm_res=(res.RES_V_dis_lRPE{1,i}(Cl_Inds,toi,foi)-nanmean(res.RES_V_dis_lRPE{1,i}(Cl_Inds,PreOnset,foi),[2,3]))./nanstd(res.RES_V_dis_lRPE{1,i}(Cl_Inds,PreOnset,foi),[],[2,3]);
    Norm_res=squeeze(nanmean(Norm_res,1))';
    mCoh=conv2(Norm_res,Wt,'same')/sum(Wt,'all');
    uimagesc(Time,Fr,mCoh)
    up=prctile(Stat,99.5,'all');
    Sign((mCoh>up))=1;
    hold on
    contour(Time,Fr,Sign,1,'k','linewidth',1.5)
    set(gca,'clim',[-1.2 1.2])
    hold on
    plot([0 0],get(gca,'ylim'),'linewidth',1,'linestyle','--','color','k')
    axis xy
    colormap('jet')
    title(['low-RPE class-' classN{i}])
    ylabel 'Norm. Coherence'
    if i==8
        xlabel 'Time to Reward (msec)'
    end
end

% High RPE ACC
state=cell2mat(res.RES_V_dis_hRPE');
state(isnan(sum(state,[2,3])),:,:)=[];
L1=size(state,1);
Sign=zeros(length(foi),length(toi));
for i=1:8
    Cl_Inds=find(strcmp(Area2{1, i},'ACC'));
    L=sum(~isnan(sum(res.RES_V_dis_hRPE{1,i}(Cl_Inds,:,:),[2,3])));
    Sign=zeros(length(foi),length(toi));
    for j=1:1000
        perm = randperm(L1);
        inds = perm(1:L);
        perm_st=state(inds,:,:);
        Norm_perm =(perm_st(:,toi,foi)-nanmean(perm_st(:,PreOnset,foi),[2,3]))./nanstd(perm_st(:,PreOnset,foi),[],[2,3]);
        Norm_perm = squeeze(nanmean(Norm_perm,1))';
        mCoh=conv2(Norm_perm,Wt,'same')/sum(Wt,'all');
        Stat(j,:,:) = mCoh;
    end
    subplot(8,4,4*(i-1)+4)
    Cl_Inds=find(strcmp(Area2{1, i},'ACC'));
    Norm_res=(res.RES_V_dis_hRPE{1,i}(Cl_Inds,toi,foi)-nanmean(res.RES_V_dis_hRPE{1,i}(Cl_Inds,PreOnset,foi),[2,3]))./nanstd(res.RES_V_dis_hRPE{1,i}(Cl_Inds,PreOnset,foi),[],[2,3]);
    Norm_res=squeeze(nanmean(Norm_res,1))';
    mCoh=conv2(Norm_res,Wt,'same')/sum(Wt,'all');
    uimagesc(Time,Fr,mCoh)
    up=prctile(Stat,99.5,'all');
    Sign((mCoh>up))=1;
    hold on
    contour(Time,Fr,Sign,1,'k','linewidth',1.5)
    set(gca,'clim',[-1.2 1.2])
    hold on
    plot([0 0],get(gca,'ylim'),'linewidth',1,'linestyle','--','color','k')
    axis xy
    colormap('jet')
    title(['high-RPE class-' classN{i}])
    ylabel 'Norm. Coherence'
    if i==8
        xlabel 'Time to Reward (msec)'
    end
end

set(gcf,'position',[561   151   596   798])
