load('data_fig_6.db','-mat');
close all;


%%
f=figure;
subplot(2,1,1)
sel=1:numel(data.v);
speeds=data.v(sel)';
accs=[0 speeds(3:end)-speeds(1:end-2) 0]/.04;
ths=smooth_theta2(accs,1-exp(-.1)); 

raccs=accs;
raccs(raccs<0)=0;

plot(sel/50,(speeds)/std(speeds),'r');
hold on;


plot(sel/50,accs/std(accs),'b');
plot(sel/50,ths/std(accs),'g');
xlim([0 10])


ths=8+ths/100;
da=20; abins=-5:5;
dv=5; vbins=0:12;
rv=round(speeds/dv);
ra=round(accs/da);
vfr=arrayfun(@(x) mean(ths(rv==x)),vbins);
afr=arrayfun(@(x) mean(ths(ra==x)),abins);
m8=vfr(1);
[th_f5,th_s5]=mybp(ths-m8,0.5);
[th_f1,th_s1]=mybp(ths-m8,0.1);
[th_f3,th_s3]=mybp(ths-m8,0.3);

vfr_f=arrayfun(@(x) mean(th_f5(rv==x))+m8,vbins);
vfr_s=arrayfun(@(x) mean(th_s5(rv==x))+m8,vbins);

afr_f=arrayfun(@(x) mean(th_f1(ra==x))+m8,abins);
afr_s=arrayfun(@(x) mean(th_s1(ra==x))+m8,abins);


subplot(2,2,3)
plot(vbins*dv,[vfr; vfr_s; vfr_f]');
axis([0 50 7.8 9.1])
golden(2.5)

subplot(2,2,4)
plot(abins*da,[afr; afr_s; afr_f]');
axis([-120 140 7.8 9.1])
set(gca,'YTick',6:.5:12);



[corr(speeds',ths') corr(accs',ths');
corr(speeds',th_f5') corr(accs',th_f1');
corr(speeds',th_s5') corr(accs',th_s1')]
drawnow;




function [th_f,th_s]=mybp(ths,f1)
    th_f=real(bandpass_and_hilbert(ths',f1,100,50)');
    th_s=real(bandpass_and_hilbert(ths',0,f1,50)');
end