clear; clc;

%% read sipke2 singal.
load('LAMP2_037.mat');
sig = LAMP2_037_AChSig1; opt = LAMP2_037_OpTrTTL; sti = LAMP2_037_AChSti; opo = LAMP2_037_OpGeSti;

staced(1,1) = find(opt.values(:,1)>1,1);
staced(1,2) = find(opt.values(:,1)>1,1,'last');
statimced = opt.times(staced(1,:));   

signal = sig.values;
signal_t = sig.times - statimced(1,1);


tim_so = downsample((opt.times(staced(1,1):staced(1,2)) - statimced(1)),200);
stimulation = downsample(sti.values(staced(1,1):staced(1,2)),5);
optogenetics = downsample(opo.values(staced(1,1):staced(1,2)),5);


%% read intan data
bz_LFPfromDat(pwd);
LFPs = bz_GetLFP([0: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,:);
EMG = bz_EMGFromLFP(pwd);

%% read movement.
fps = 100; mvt = 1;
optitrack=csvread('LAMP2_041.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.
% p = polyfit(tim, signal, 1);
% signal_d = (signal - polyval(p, tim))./polyval(p, tim);
signal = interp1(signal_t,signal,tim_o);
signal = detrend(signal,'linear');

thetarange = linspace(4,12,20);
ripplerange = linspace(140,250,50);
brchan = bz_GetBestRippleChan(LFPs);
lfp = double(LFPs.data(:,LFPs.channels==brchan));
[spectrogram1,t,f] = MTSpectrogram(lfp,'window',2);
bands_all = SpectrogramBands(spectrogram1,f,'smooth',2);

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

%% REM 

% from LAMP2_041
lfp_rem(:,1) = lfp((730*1250):(850*1250));
signal_rem(:,1) = signal((730*100):(850*100));
vel_rem(:,1) = vel((730*100):(850*100));
ratios_rem(:,1) = bands_all.ratio(730:850);

lfp_rem(:,2) = lfp((730*1250):(850*1250));
signal_rem(:,2) = signal((730*100):(850*100));
vel_rem(:,2) = vel((730*100):(850*100));
ratios_rem(:,2) = bands_all.ratio(730:850);

% from LAMP2_002

lfp_rem(:,3) = lfp((744*1250):(864*1250));
signal_rem(:,3) = signal((744*100):(864*100));
vel_rem(:,3) = vel((744*100):(864*100));
ratios_rem(:,3) = bands_all.ratio(744:864);

lfp_rem(:,4) = lfp((3240*1250):(3360*1250));
signal_rem(:,4) = signal((3240*100):(3360*100));
vel_rem(:,4) = vel((3240*100):(3360*100));
ratios_rem(:,4) = bands_all.ratio(3240:3360);

lfp_rem(:,5) = lfp((3804*1250):(3924*1250));
signal_rem(:,5) = signal((3804*100):(3924*100));
vel_rem(:,5) = vel((3804*100):(3924*100));
ratios_rem(:,5) = bands_all.ratio(3804:3924);

signal_rem_avg = mean(signal_rem,2);
signal_rem_std = std(signal_rem,0,2);
vel_rem_avg = mean(vel_rem,2);
vel_rem_std = std(vel_rem,0,2);
ratios_rem_avg = mean(ratios_rem,2);
ratios_rem_std = std(ratios_rem,0,2);

REM_before.signal = REM.Signal(1:5000,:);
REM_before.vel = REM.Vel(1:5000,:);
REM_before.lfp = REM.lfp(1:62500,:);
REM_before.thetascore = REM.thetascore(1:50,:);

REM_in.signal = REM.Signal(5001:12001,:);
REM_in.vel = REM.Vel(5001:12001,:);
REM_in.lfp = REM.lfp(62501:150001,:);
REM_in.thetascore = REM.thetascore(51:121,:);




%% peri-event
opo_sti=find(optogenetics>1);
lfp1 = bz_GetLFP(brchan,'noPrompts',true);
lfp_theta = bz_Filter(lfp1, 'passband', [5 10]);

t1=signal(400000:600000);
t2=optogenetics(400000:600000);

t2_index=find(optogenetics>5.9);
t2_g=[];
 for ii=1:length(t2_index)
    if t2_index(ii+1)-t2_index(ii)>1000
        t2_g(ii,1)=t2_index(ii)
    end
 end
t2=t2_g(find(t2_g>1));
i=2
for i = 2:length(t2)
  signal_opo(:,i) = signal((t2(i)-1000):(t2(i)+1000));
end

signal_opo_before = signal_opo(1:950,1:6);
signal_opo_in = signal_opo(1050:1500,1:6);


errorbar(tim_indx-1000,value_ave,value_std,'CapSize',0)
hold on
plot(tim_indx-1000,value_ave);


%% plot
figure
subplot(6,1,[1:2])
yyaxis left;
plot(tim_so,optogenetics);
xlim([min(tim_so) max(tim_so)]); ylim([0,6]);
ylabel('Opotogenetics Stimulation');
yyaxis right;
plot(tim_o,signal);
xlim([min(tim_o) max(tim_o)]); ylim([-0.06 0.06])
ylabel('Ach signal (%)');yticks([-0.06,-0.03,0,0.03,0.06]);
box off;



subplot(6,1,[3])
plot(tim_o,vel,'b'); xlim([min(tim_o) max(tim_o)]); ylabel('Movment'); ylim([0 0.5]);
box off;
% subplot(8,1,[5])
% plot(LFPs.timestamps,lfp,'k'); xlim([min(LFPs.timestamps),max(LFPs.timestamps)]); ylabel('Raw LFP');
% box off; axis off;
subplot(6,1,[4])
plot(bands_all.ratio,'k'); ylabel('theta score');ylim([0 6]);xlim([min(tim_o) max(tim_o)]);
box off; 
% subplot(12,1,[7])
% plot(bands_all.theta,'k'); xlim([0 6102]); ylabel('theta power');
% box off; 
% subplot(12,1,[8])
% plot(bands_all.delta,'k'); xlim([0 6102]); ylabel('delta power');
% box off; 
% subplot(12,1,[9])
% plot(bands_all.ripples,'k'); xlim([0 6102]); ylabel('ripple power');
% box off; 


% subplot(12,1,[10])
% colormap('jet')
% imagesc(t,ripplerange,abs(sr)); clim([0,25000])
% yticks([140,200,250]); ylabel('Frequency (Hz)'); xticks([]);
% set(gca,'YDir','normal')?
subplot(6,1,[5:6])
colormap('jet')
imagesc(specttim,thetarange,abs(st)); clim([0,2e5])
yticks([4,8,12]); ylabel('Frequency (Hz)'); xlabel('Time (s)');
set(gca,'YDir','normal');

figure(2)
subplot(5,1,[1:2])
yyaxis left;
plot(tim,optogenetics);
xlim([5600 5800]); ylim([0,6]);
ylabel('Opotogenetics Stimulation');
yyaxis right;
plot(tim,signal);
xlim([5600 5800]); ylim([-0.06 0.06])
ylabel('Ach signal (%)');yticks([-0.06,-0.03,0,0.03,0.06]);
box off;



subplot(5,1,[3])
plot(tim_o,vel,'b'); xlim([5600 5800]); ylabel('Movment');
box off;
% subplot(8,1,[5])
% plot(LFPs.timestamps,lfp,'k'); xlim([min(LFPs.timestamps),max(LFPs.timestamps)]); ylabel('Raw LFP');
% box off; axis off;
subplot(5,1,[4])
plot(bands_all.ratio,'k'); xlim([5600 5800]); ylabel('theta score');ylim([0 7])
box off; 
% subplot(12,1,[7])
% plot(bands_all.theta,'k'); xlim([0 6102]); ylabel('theta power');
% box off; 
% subplot(12,1,[8])
% plot(bands_all.delta,'k'); xlim([0 6102]); ylabel('delta power');
% box off; 
% subplot(12,1,[9])
% plot(bands_all.ripples,'k'); xlim([0 6102]); ylabel('ripple power');
% box off; 


% subplot(12,1,[10])
% colormap('jet')
% imagesc(t,ripplerange,abs(sr)); clim([0,25000])
% yticks([140,200,250]); ylabel('Frequency (Hz)'); xticks([]);
% set(gca,'YDir','normal');


subplot(5,1,[5])
colormap('jet')
imagesc(specttim,thetarange,abs(st)); clim([0,3e5])
yticks([4,8,12]); ylabel('Frequency (Hz)'); xlabel('Time (s)');xlim([5600 5800]);
set(gca,'YDir','normal');


figure()
yyaxis left;
plot(signal_win);
yyaxis right;
plot(bands_all.theta,'k'); ylim([0,1e5])

signal_win = interp1(tim_o,signal,1:length(bands_all.theta));
bands_all.theta(bands_all.theta>1e5) = 0;
