clear; clc;

%% read sipke2 singal.
load('LAMP2_002_mat.mat');
ACh_signal = LAMP2_002_signal; time = LAMP2_002_optitrck; sti = LAMP2_002_stim; opo = LAMP2_002_Opotogen;


%% read stimulation
opto_sti=find(opo.values(:) > 5);

%% read intan data
LFPs = bz_GetLFP([0:35,48:63],'noPrompts',true);

info = read_Intan_RHD2000_file(pwd);
num_channels = info.num_board_dig_in_channels;
digital_ch = read_Intan_digitalin_file('digitalin.dat',num_channels);
digital = digital_ch(:,1);

staintan(1,1) = find(digital,1)/info.sample_rate;
staintan(1,2) = find(digital,1,'last')/info.sample_rate;

digital = downsample(digital',16);
idx = LFPs.timestamps >= staintan(1,1) & LFPs.timestamps <= staintan(1,2);
LFPs.timestamps = LFPs.timestamps(idx)-staintan(1,1);
LFPs.data = LFPs.data(idx,:);

%% read movement.
mvt = 1;
fps = 100;
optictrack = csvread('LAMP2_002_mat.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);


%% processing.
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);

opo.data =downsample(opo.values(start(1,1):start(1,2)),200);
opo.times = downsample(opo.times(start(1,1):start(1,2)),200)-opo.times(start(1,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);

thetarange = linspace(0,20,50);
lfp = double(LFPs.data(:,LFPs.channels==29));
specttim = unique(floor(LFPs.timestamps)); specttim = specttim(1:end-1);
[st,ft,tt] = spectrogram(lfp,1*LFPs.samplingRate,0,thetarange,LFPs.samplingRate);
[spectrogram,t,f] = MTSpectrogram(lfp);
bands = SpectrogramBands(spectrogram,f);

%% optogenetics peri-event
opo_idx = find(opo.data> 5.9);
opo_time = [];
i = 1;
for i=1:length(opo_idx)
    if opo_idx(i+1,1)-opo_idx(i,1)>1000
         opo_time(i,1) = opo_idx(i);
    end
end

opo_time=opo_time(find(opo_time>1));
opto_sti_ns(:,1)= opo_time-1000;opto_sti_ns(:,2)= opo_time-500;
opto_sti_ss(:,1)= opo_time-500;opto_sti_ss(:,2)= opo_time;

ACh_opo = [];
i=1;
for i = 1:length(opo_time)
    ACh_opo(:,i)=ACh_d((opo_time(i)-1500):(opo_time(i)+500))
end
ACh_opo = ACh_opo(:,2:8);

idx = 1:length(ACh_opo);
order = 3;
framelen =501;
ACh_opos = sgolayfilt(ACh_opo,order,framelen);
ACh_opos_ms (:,1)=mean(ACh_opos,2);
ACh_opos_ms (:,2)=std(ACh_opos,0,2);
X = idx; Y_p =  ACh_opos_ms (:,1)+ACh_opos_ms (:,2);Y_n =  ACh_opos_ms (:,1)-ACh_opos_ms (:,2);
XX = [X flip(X)];YY = Y_p;YY(2002:4002,1) = flip(Y_n);
fill(XX,YY,'r');fig = gcf;alpha(fig,.5)
hold on;plot(ACh_opos_ms (:,1),'LineWidth',2);
axis off;axis on


%% statistic
ACh_oposs(:,1) = (mean(ACh_opos(501:1000,1:7)))';
ACh_oposs(:,2) = (mean(ACh_opos(1001:1500,1:7)))'; [h,p,ci,stats] = ttest2(ACh_oposs(:,1),ACh_oposs(:,2))
g = boxplot(ACh_oposs);
hold on;
a(1:size(ACh_oposs,1)) = 1;
scatter(a,ACh_oposs(:,1),'ro','filled')
b(1:size(ACh_oposs,1)) = 2;
scatter(b,ACh_oposs(:,2),'ro','filled')
for i = 1:length(a)line([a(i),b(i)], [ACh_oposs(i,1) ACh_oposs(i,2)] )
end
set(g,'color','b');ylim([-0.01 0.07]);yticks([-0.01 0.03 0.07]);yticklabels({'-1','3','7'});
xticks([1 2]);xticklabels({'Control','Sti'});ylabel('Ach(?f/f)')
set(gca,'FontSize',16,'LineWidth',1);box off





%% plot
figure
subplot(3,1,1)
plot(ACh_signal_f.times(70800:102000),ACh_d(70800:102000));hold on;ylim([-0.1 0.15]);xlim([708 1020]);
ylabel('ACh(?f/f)','FontSize',18);yticks([-0.1 ,0, 0.15]);box off;

order = 3;
framelen = 1001;
ACh_g = sgolayfilt(ACh_d,order,framelen);
plot(ACh_signal_f.times(70800:102000),ACh_g(70800:102000),'linewidth',1.5);ylim([-0.1 0.15]);xlim([708 1020]);box off;

subplot(3,1,2);
idx = 1: length(bands.theta);
plot(idx(280:404),bands.ratios.hippocampus(280:404));xlim([280 390]);%% box off; axis off;
ylabel('Theta/delta','FontSize',18);box off;


subplot(3,1,3);
imagesc(specttim,thetarange,abs(st)); clim([0,3e5])
yticks([0,10,20]); ylabel('Frequency (Hz)'); xlabel('Time (s)');xlim([700 1000])
set(gca,'YDir','normal')




%% plot


subplot(12,1,[4])
plot(tim_o,vel,'b'); xlim([min(tim_o),max(tim_o)]); ylabel('Movment');
box off; axis off;
subplot(12,1,[5:6])
plot(LFPs.timestamps,lfp,'k'); xlim([min(LFPs.timestamps),max(LFPs.timestamps)]); ylabel('Raw LFP');
box off; axis off;
subplot(12,1,[7:9])
colormap('jet')
imagesc(specttim,ripplerange,abs(sr)); clim([0,25000])
yticks([140,200,250]); ylabel('Frequency (Hz)'); xticks([]);
set(gca,'YDir','normal');
subplot(12,1,[10:12])
imagesc(specttim,thetarange,abs(st)); clim([0,3e5])
yticks([0,10,20]); ylabel('Frequency (Hz)'); xlabel('Time (s)')
set(gca,'YDir','normal');


