Theta_acceleration_data/code/Fig_1b/Fig_1c.m

load('data_fig_1c.db','-mat');

ecp=bsxfun(@rdivide,bsxfun(@minus,eeghp.ecL.p,nanmean(eeghp.ecL.p)),nanmean(amp.ecL.p));
ecm=bsxfun(@rdivide,bsxfun(@minus,eeghp.ecL.m,nanmean(eeghp.ecL.m)),nanmean(amp.ecL.m));

egp=bsxfun(@rdivide,bsxfun(@minus,eeg.ecL.p,nanmean(eeg.ecL.p)),nanmean(amp.ecL.p));
egm=bsxfun(@rdivide,bsxfun(@minus,eeg.ecL.m,nanmean(eeg.ecL.m)),nanmean(amp.ecL.m));
%%
ks=20;
sel=4801:14400;
aux=ecp(sel,:);
[idx,Cp]=kmeans(aux',ks);
aux=ecm(sel,:);
[idxm,Cm]=kmeans(aux',ks);
%%
close all;
f=figure;
pp=zeros(19201,ks);
pm=pp;
for j=1:ks
   pp(:,j)=mean(ecp(:,idx==j),2);
   pm(:,j)=mean(ecm(:,idxm==j),2);
end

C=corr(pp(sel(2401:4800),:),pm(sel(2401:4800),:));
C(1:ks+1:end)=-1;
nsp=8;
fp=zeros(1,nsp);
fm=fp;
phm=fp;
for j=1:nsp
    [mC,rw]=max(C);
    [~,cl]=max(mC);
    rw=rw(cl);
  

fp(j)=rw;
fm(j)=cl;
   
   C(:,cl)=-1;
   C(rw,:)=-1;
   phm(j)=angle(mean(exp(1i*phase.ecL.m(9601,bsxfun(@eq,idxm,fm(j)))),2));
end
[~,sindex]=sort(phm,'descend');
for k=1:nsp
    j=sindex(k);
    subplot(nsp,1,k)
    rw=fp(j); cl=fm(j);
    disp([C(rw,cl),corr(pp(sel(2401:4800),rw),pm(sel(2401:4800),cl)),sum(idx==rw),sum(idxm==cl)])

    auxp=mean(ecp(:,idx==fp(j)),2); 
    plot((-9600:9600)/4800,auxp,'c','LineWidth',2);
   hold on;
   auxm=mean(ecm(:,idxm==fm(j)),2);
   plot((-9600:9600)/4800,auxm,'b','LineWidth',2);
   axis([-.5 .5 1.4*[-1 1]]);   
   set(gca,'XTick',-0.5:.5:.5,'XTickLabel',[]);
   golden(.3)
end
set(gca,'XTick',-0.5:.5:.5,'XTickLabel',-0.5:.5:.5,'FontName','Arial');