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;
optspk2 = time.times(find(time.values>0.01,1));

%% 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;

ACh_signal.times = ACh_signal.times - optspk2 + staintan(1,1);
digital = downsample(digital',16);

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

ACh_signal_f.data =downsample(ACh_signal.values,2000);
ACh_signal_f.times = downsample(ACh_signal.times,2000);
dff = (ACh_signal_f.data-median(ACh_signal_f.data))./median(ACh_signal_f.data);
ACh_d = msbackadj(ACh_signal_f.times,dff,'WindowSize',40);

% vel = interp1(t(1:(end)-1),vel,ACh_signal_f.times);

clear st ft tt spectrogram;
thetarange = linspace(0,15,50);
lfp = double(LFPs.data(:,LFPs.channels==28));
lfp(abs(lfp)>5000) = nan;
lfp = fillmissing(lfp,'linear');
[st,ft,tt] = spectrogram(lfp,2*LFPs.samplingRate,0.6*LFPs.samplingRate,thetarange,LFPs.samplingRate);
[spectrogrammt,t,f] = MTSpectrogram(lfp,'window',3);
bands = SpectrogramBands(spectrogrammt,f);

%% plot
figure
subplot(3,1,1)

order = 3;
framelen = 101;
ACh_g = sgolayfilt(ACh_d,order,framelen);
idxs = ACh_signal_f.times>700 & ACh_signal_f.times<1000;
plot(ACh_signal_f.times(idxs),ACh_d(idxs));hold on;
%box off; axis off;
plot(ACh_signal_f.times(idxs),ACh_g(idxs),'linewidth',1.5); 
ylim([-0.05 0.07]);xlim([700 1000]); %box off; axis off;

subplot(3,1,2);
idxs = t>700 & t<1000;
plot(t(idxs),imgaussfilt(bands.ratios.hippocampus(idxs),5));
%box off; axis off;
ylim([0,6])

subplot(3,1,3);
imagesc(tt,thetarange,abs(st)); clim([0,3e5]);
ylim([0,15]);yticks([0,5,10,15]); 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');


