clear; clc;

%% read sipke2 singal.
load('LAMP2_003_mat.mat');
sig = LAMP2_003_signal; opt = LAMP2_003_optitrck; LAMP2_003_stim; opo = LAMP2_003_Opotogen;

staced(1,1) = find(opt.values(:,1)>0.005,1);
staced(1,2) = find(opt.values(:,1)>0.005,1,'last');
statimced = opt.times(staced(1,:));
   
tim = downsample((opt.times(staced(1,1):staced(1,2)) - statimced(1)),200);
signal = downsample(sig.values(staced(1,1):staced(1,2)),200);
stimulation = downsample(sti.values(staced(1,1):staced(1,2)),200);
optogenetics = downsample(opo.values(staced(1,1):staced(1,2)),200);

%% 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.
fps = 100; mvt = 1;
optitrack=csvread('LAMP2_002.csv',7,1);
idx = find(optitrack(:,2)==0);
optitrack(idx,2:8) = nan(size(idx,1),7);
optitrack = fillmissing(optitrack,'linear');
tim_o = optitrack(:,1);
x = optitrack(:,6);
y = optitrack(:,8);
z = optitrack(:,7);
vel = [0;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.
% signal = detrend(signal,'linear');
p = polyfit(tim, signal, 1);
signal_d = (signal - polyval(p, tim))./polyval(p, tim);

thetarange = linspace(0,20,50);
ripplerange = linspace(140,250,50);

brchan = bz_GetBestRippleChan(LFPs);
lfp = double(LFPs.data(:,LFPs.channels==brchan));

specttim = unique(floor(LFPs.timestamps)); specttim = specttim(1:end-1);
[st,ft,tt] = spectrogram(lfp,1*LFPs.samplingRate,0,thetarange,LFPs.samplingRate);
[sr,fr,tr] = spectrogram(lfp,1*LFPs.samplingRate,0,ripplerange,LFPs.samplingRate);

%% plot
subplot(12,1,[1:3])
yyaxis left;
plot(tim,optogenetics);
xlim([min(tim),max(tim)]); ylim([0,6]);
ylabel('Opotogenetics Stimulation');
yyaxis right;
plot(tim,signal_d);
xlim([min(tim),max(tim)]); ylim([-0.1 0.1])
ylabel('Ach signal');
box off;axis off

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,5e5])
yticks([0,10,20]); ylabel('Frequency (Hz)'); xlabel('Time (s)')
set(gca,'YDir','normal');


