clear all;clc;


basepath = pwd;
files = dir('*.mat');
%% ach syn and caculation
ii = 1
[filepath,name,ext] = fileparts(files(ii).name);
% syn
load([name,ext],[name,'_OpTrTTL']);
Opti = eval([name,'_OpTrTTL']);
Opti.data = downsample(Opti.values,10);
Opti.time = downsample(Opti.times,10);
time_st = round(Opti.time(find(Opti.data>1,1))*2e4);time_nd = round(Opti.time(find(Opti.data>1,1,'last'))*2e4);
% signal
load([name,ext],[name,'_AChSig1']);
ACh_signal = eval([name,'_AChSig1']);
ACh_signal.data = downsample(ACh_signal.values(time_st:time_nd),200);
ACh_signal.time = downsample(ACh_signal.times(time_st:time_nd),200)- Opti.time(find(Opti.data>1,1));

p = [];
p = polyfit(ACh_signal.time,ACh_signal.data,10);
Ybaseline = polyval(p,ACh_signal.time);

figure
plot(ACh_signal.time, ACh_signal.data); hold on;
plot(ACh_signal.time,Ybaseline);
dFoverF = smoothdata((ACh_signal.data - Ybaseline)./Ybaseline,'movmean',10);
% dfff(:,1) = [1043:3000:731366];
% dfff(:,2) = dfff(:,1)+60;
% dFoverF(dfff(1:2,:))= NaN;
% dFoverF = fillmissing(dFoverF,'linear');
hold on;
plot(ACh_signal.time, dFoverF);xlim([0 inf])

%% MOVEMENT

sessionname = [name];
mvt = 1;fps = 120;
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);
tim = ACh_signal.time;

hold on; plot(t(1:end-1),vel-0.6);

%% STIMUALTION
load([name,ext],[name,'_40hz']);
door = eval([name,'_40hz']);

door.data = downsample(door.values(time_st:time_nd),200);
door.time = downsample(door.times(time_st:time_nd),200)- Opti.time(find(Opti.data>1,1));

Index = door.time(door.data>1);
pulse = Index(diff(Index)>3)-15;pulse(:,2) = pulse(:,1)+20;

ACh_Opto = [];
OPindex = round(pulse*100);
for jj = 1:length(pulse)
    ACh_Opto(jj,:) = dFoverF(OPindex(jj,1):OPindex(jj,2));
%     ACh_Opto(jj,496:549)= NaN; ACh_Opto(jj,991:1016) = NaN;
%     ACh_Opto(jj,:) = fillmissing(ACh_Opto(jj,:),'linear');
end
h = figure;
subplot(2,1,1)
options.color_area = [128 193 219]./255;    % Blue theme
options.color_line = [ 52 148 186]./255;
%options.color_area = [243 169 114]./255;    % Orange theme
%options.color_line = [236 112  22]./255;
options.alpha      = 0.5;
options.line_width = 2;
options.error      = 'sem';
plot_areaerrorbar(ACh_Opto, options);set(gcf,'Position',  [100, 100, 400, 300])
xlim([200 1800]);title('ach_door')
savefig(h,[name,'.fig']);
print(h,name,'-dsvg');
save([name,'_Peri','.mat'],'ACh_Opto');


vel_Opto = [];
OPindex = round(pulse*120);
for jj = 1:length(pulse)
    vel_Opto(jj,:) = vel(OPindex(jj,1):OPindex(jj,2));
%     ACh_Opto(jj,496:549)= NaN; ACh_Opto(jj,991:1016) = NaN;
%     ACh_Opto(jj,:) = fillmissing(ACh_Opto(jj,:),'linear');
end
subplot(2,1,2)
% options.color_area = [128 193 219]./255;    % Blue theme
% options.color_line = [ 52 148 186]./255;
options.color_area = [243 169 114]./255;    % Orange theme
options.color_line = [236 112  22]./255;
options.alpha      = 0.5;
options.line_width = 2;
options.error      = 'sem';
plot_areaerrorbar(vel_Opto, options);set(gcf,'Position',  [100, 100, 400, 300])
xlim([200*1.2 1800*1.2]);title('ach_door')
savefig(h,[name,'.fig']);
print(h,name,'-dsvg');
save([name,'_Peri','.mat'],'ACh_Opto');








