Theta_acceleration_data/code/Fig_1/Fig_1.m

load('data_fig_1.db','-mat');
figure;

%% theta frequency
frwin=3;
for a=1:2
    rx=size(data.unwrapped_phase_right{a},1);
    freq_right{a}=([data.unwrapped_phase_right{a}(:,1+frwin:end) nan(rx,frwin)]-[nan(rx,frwin) data.unwrapped_phase_right{a}(:,1:end-frwin)])/2/pi/0.02/frwin/2;
    rx=size(data.unwrapped_phase_left{a},1);
    freq_left{a}=([data.unwrapped_phase_left{a}(:,1+frwin:end) nan(rx,frwin)]-[nan(rx,frwin) data.unwrapped_phase_left{a}(:,1:end-frwin)])/2/pi/0.02/frwin/2;    
end

%% acceleration
win=10;    %gfr=fr;


    rx=size(data.x_right,1);
    acc_right=-([nan(rx,win) data.x_right(:,1:end-win)]+[data.x_right(:,1+win:end) nan(rx,win)]-2*data.x_right)/0.02/0.02/win/win;
    rx=size(data.x_left,1);    
    acc_left=([nan(rx,win) data.x_left(:,1:end-win)]+[data.x_left(:,1+win:end) nan(rx,win)]-2*data.x_left)/0.02/0.02/win/win;    


 fct=0; 
th=8; axlim=[-120  220]; 
subplot(2,1,1)

gacc=acc_right;
gfr=freq_right;
do_plots;

set(AX(2),'XTick',[]);

subplot(2,1,2)
    gfr=freq_left; 
    gacc=acc_left;    

do_plots;
set(AX(2),'XTick',[]);