clear; clc
warning('off');

par = {'lamp1_2019-03-21_001','_190321_161429\',6,'Acetylcholine',[0:31];
       'lamp1_2019-03-21_002-s','_190321_175559\',1.5,'Serotonin',[33:63];
       'lamp2_2019-03-21_001','_190321_171615\',4.5,'Acetylcholine',[0:35,48:63];
        };
ip = 3;
    
fps = 100; mvt = 1; otoffset = par{ip,3};
rootstr = par{ip,1};
optictrack = csvread([rootstr,'.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);
z = optictrack(:,7);
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);


sesssion = ['.\',rootstr];
channel = par{ip,5};

thetarange = linspace(0,20,50);
ripplerange = linspace(140,250,50);
LFPs = bz_GetLFP(channel,'basepath','E:\acetylcholine imaging\data\lamp\lamp2\lamp2_2019-03-21_001','noPrompts',true);
brchan = bz_GetBestRippleChan(LFPs);
lfp=double(LFPs.data(:,50));
tim = LFPs.timestamps;

num_channels = 2;
fileinfo = dir(fullfile(sesssion,'analogin.dat'));
num_samples = fileinfo.bytes/(num_channels * 2); % uint16 = 2 bytes
fid = fopen(fullfile(sesssion,'analogin.dat'), 'r'); 
v = fread(fid, [num_channels, num_samples], 'uint16'); 
fclose(fid);
v = v * 3.3/2^16; % convert to volts
v = downsample(v',16);
tmpdata.samplingRate = LFPs.samplingRate;
tmpdata.timestamps = LFPs.timestamps;
tmpdata.data = v;
sigfil = bz_Filter(tmpdata,'passband',[0,10],'filter','butter');
idx = sigfil.data(:,2) > 0.2;
% signal = detrend(sigfil.data(idx,1));
p = polyfit(sigfil.timestamps(idx), sigfil.data(idx,1), 3);
signal = sigfil.data(idx,1) - polyval(p, sigfil.timestamps(idx));
lfp = lfp(idx,:);
tim = tim(idx);

specttim = unique(floor(tim)); 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);

figure
set(gcf,'color','w')
subplot(12,1,[1:2])
plot(t(2:end)+otoffset,vel,'r'); xlim([min(tim),1100]);
box off; axis off; text(0,max(vel)*1.0,"Speed");
subplot(12,1,[3:4])
plot(tim,signal,'b'); xlim([min(tim),1100]);
box off; text(0,max(signal)*1.2,[par{ip,4}," (polynomial detrended)"]);
subplot(12,1,[5:6])
plot(tim,lfp,'k'); xlim([min(tim),1100]);
box off; axis off; text(0,max(lfp*1.2),"LFP");
subplot(12,1,[7:9])
colormap('jet')
imagesc(specttim,ripplerange,abs(sr)); clim([0,10000])
yticks([140,200,250]); ylabel('Frequency (Hz)'); xticks([]);xlim([min(tim),1100])
set(gca,'YDir','normal');
subplot(12,1,[10:12])
imagesc(specttim,thetarange,abs(st)); clim([0,2e5])
yticks([0,10,20]); ylabel('Frequency (Hz)'); xlabel('Time (s)');xlim([min(tim),1100])
set(gca,'YDir','normal');


%% plot small figure
figure
set(gcf,'color','w')
subplot(12,1,[1:2])
plot(t(2:end)+otoffset,vel,'r'); xlim([973,978]);
box off; axis off; text(0,max(vel)*1.0,"Speed");
subplot(12,1,[3:4])
plot(tim,signal,'b'); xlim([973,978]);
box off; text(0,max(signal)*1.2,[par{ip,4}," (polynomial detrended)"]);
subplot(12,1,[5:6])
plot(tim,lfp,'k'); xlim([973,978]);
box off; axis off; text(0,max(lfp*1.2),"LFP");
% subplot(12,1,[7:9])
% colormap('jet')
% imagesc(specttim,ripplerange,abs(sr)); clim([0,10000])
% yticks([140,200,250]); ylabel('Frequency (Hz)'); xticks([]);xlim([min(tim),1100])
% set(gca,'YDir','normal');
subplot(12,1,[10:12])
imagesc(specttim,thetarange,abs(st)); clim([0,2e5])
yticks([0,10,20]); ylabel('Frequency (Hz)'); xlabel('Time (s)');xlim([973,978])
set(gca,'YDir','normal');
