clear; clc;

sessionname = 'LAMP34_001';

mvt = 1;
fps = 100;

%% position
optictrack = csvread([sessionname,'.csv'],7,1);
idx = find(optictrack(:,2)==0);
optictrack(idx,2:8) = nan(size(idx,1),7);
optictrack = fillmissing(optictrack,'linear');
t = optictrack(:,1);
x = optictrack(:,6);
y = optictrack(:,8);
vel = sqrt((x(2:end) - x(1:end-1)).^2 + (y(2:end) - y(1:end-1)).^2)*fps;
vel = smoothdata(vel,'movmean',fps*mvt);
plot(vel);

%%processing 
srdown = 200;
load([sessionname,'.mat']);
ACh_signal = LAMP34_001_AChSig1
time = LAMP29_002_OpTrTTL;
a = time.times(find(time.values>3));

start(1,1)= a(1)*20000;start(1,2)=a(end)*20000;

ACh_signal_f.data = downsample(ACh_signal.values(start(1,1):start(1,2)),20000/srdown);
ACh_signal_f.times = downsample(ACh_signal.times(start(1,1):start(1,2)),20000/srdown);
ACh_signal_f.times = ACh_signal_f.times-ACh_signal_f.times(1);
p = polyfit(ACh_signal_f.times,ACh_signal_f.data,1);
Ybaseline = polyval(p,ACh_signal_f.times);
dFoverF = smoothdata((ACh_signal_f.data - Ybaseline)./Ybaseline,'movmean',5*srdown);
dff = (ACh_signal_f.data - Ybaseline)./Ybaseline;

Ndown = 2;
dff_downsampled = nanmean(reshape([dff(:); nan(mod(-numel(dff),Ndown),1)],Ndown,[]));
tic
dff_downsampled = dff_downsampled';
ACh_d = msbackadj(ACh_signal_f.times(1:Ndown:end),dff_downsampled(:),'WindowSize',100,'SMOOTHMETHOD','rlowess','REGRESSIONMETHOD','spline');
toc



% vel = interp1(t(1:(end)-1),vel,ACh_signal_f.times);
% vel = smoothdata(vel,'movmean',fps*mvt*2);

figure
% plot(ACh_signal_f.times,ACh_signal_f.data);hold on
plot(ACh_signal_f.times,dFoverF);hold on
plot(t(1:(end)-1),vel);
ylim([-0.2,0.4])

%%
dFoverF = interp1(ACh_signal_f.times,dFoverF,vel);

[R,P] = corr(dFoverF((200*100):(600*100)),vel((200*100):(600*100)))
%%
figure
scatter(dFoverF,vel)