clear; clc;

sessionname = 'ham37_04_002_2019-03-18';

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 
load([sessionname,'.mat']);
ACh_signal = ham37_04_002_2019_03_18_untitled;
time = ham37_04_002_2019_03_18_optitrac;
a = find(ACh_signal.times > time.times(1) & ACh_signal.times < time.times(2));
b = find(ACh_signal.times > time.times(end-1) & ACh_signal.times < time.times(end));

start(1,1)= a(1);start(1,2)=b(1);

ACh_signal_f.data =downsample(ACh_signal.values(start(1,1):start(1,2)),200);
ACh_signal_f.times = downsample(ACh_signal.times(start(1,1):start(1,2)),200)-ACh_signal.times(start(1,1));
dff = (ACh_signal_f.data-median(ACh_signal_f.data))./median(ACh_signal_f.data);
ACh_d = msbackadj(ACh_signal_f.times,dff,'WindowSize',400);
% vel = interp1(t(1:(end)-1),vel,ACh_signal_f.times);
figure
plot(ACh_signal_f.times,ACh_d);hold on;ylim([-0.5 1]);
order = 3;
framelen = 1001;
ACh_g = sgolayfilt(ACh_d,order,framelen);
plot(ACh_signal_f.times,ACh_g,'linewidth',1.5);ylim([-0.5 1]);
hold on;
yyaxis right;ylim([0 2]);
plot(t(1:(end-1)),vel);

figure();

subplot(2,1,1);
plot(ACh_signal_f.times(10000:60000),ACh_d(10000:60000));ylim([-0.1 0.3]);xlim([100 600])
hold on;
order = 3;
framelen = 1001;
ACh_g = sgolayfilt(ACh_d,order,framelen);
plot(ACh_signal_f.times(10000:60000),ACh_g(10000:60000),'linewidth',1.5);ylim([-0.1 0.3]);xlim([100 600])
hold on;
yyaxis right;ylim([-0.1 0.3]);xlim([100 600])
plot(t((10000:60000)),vel(10000:60000),'linewidth',1.5);
xlabel('Time(s)');
ylabel('Acetylcholine')

subplot(2,1,2);
plot(ACh_signal_f.times(470000:520000),ACh_d(470000:520000));ylim([-0.1 0.3]);xlim([4700 5200])
hold on;
order = 3;
framelen = 1001;
ACh_g = sgolayfilt(ACh_d,order,framelen);
plot(ACh_signal_f.times(470000:520000),ACh_g(470000:520000),'linewidth',1.5);ylim([-0.1 0.3]);xlim([4700 5200])
hold on;
yyaxis right;ylim([-0.1 0.3]);xlim([4700 5200])
plot(t((470000:520000)),vel(470000:520000),'linewidth',1.5);
xlabel('Time(s)');
ylabel('Acetylcholine')

%% statistic
statistic = [];
QW = find(vel<0.1);QW_before = QW(QW<360000);QW_after = QW(QW>360000);
AE = find(vel> 0.4);AE_before = AE(AE<360000);AE_after = AE(AE>360000);
ACh.QB = ACh_d(QW_before(1):QW_before(end)); 
ACh.QA = ACh_d(QW_after(1):QW_after(end)); 
ACh.AB = ACh_d(AE_before(1):AE_before(end)); 
ACh.AA = ACh_d(AE_after(1):AE_after(end)); 



