clear; clc;
warning('off');

sesslist = {'LAMP2_006_190415_154337';
            'LAMP2_010_190425_164518';
            'LAMP2_042_190614_120612';
            'LAMP3_012_190617_161856';
            'LAMP3_016_190619_180636';
            'LAMP3_020_190629_152350';
            'LAMP3_024_190701_124743'};

%% parameters
ip = 7; nchan = 64;
session = sesslist{ip,1};
outsr = 100;
datapath = 'C:\Users\Liang\Dropbox\AchProject\Fig2 data\Theta_Ach\';

%% load achfile and lfpfile
% ach 
achfile = dir([datapath,session,filesep,'LAMP*.mat']);
if ismember(ip,[3,4,5,6,7])
    load([achfile.folder,filesep,achfile.name],[achfile.name(1:end-4),'_AChSig1']);
    achdata = eval([achfile.name(1:end-4),'_AChSig1']);
    load([achfile.folder,filesep,achfile.name],[achfile.name(1:end-4),'_OpTrTTL']);
    achopt = eval([achfile.name(1:end-4),'_OpTrTTL']);
elseif ismember(ip,[1,2])
    load([achfile.folder,filesep,achfile.name],[achfile.name(1:end-4),'_AchSig']);
    achdata = eval([achfile.name(1:end-4),'_AchSig']);
    load([achfile.folder,filesep,achfile.name],[achfile.name(1:end-4),'_optitrac']);
    achopt = eval([achfile.name(1:end-4),'_optitrac']);
end


achtime(1) = achopt.times(find(achopt.values >= 1,1));
achtime(2) = achopt.times(find(achopt.values >= 1,1,'last'));

% lfp
load([datapath,session,filesep,'LFPs.mat'],'LFPs');
info = read_Intan_RHD2000_file([datapath,session],'*.rhd');
digitalin = read_Intan_digitalin_file([datapath,session,filesep,'digitalin.dat'],info.num_board_dig_in_channels)';

lfptime(1) = find(digitalin >= 1,1) / info.sample_rate;
lfptime(2) = find(digitalin >= 1,1,'last') / info.sample_rate;

% sync timestamp
syncachts = achdata.times - achtime(1);
synclfpts = LFPs.timestamps - lfptime(1);

%% processing ach signal
edgepts = 90*outsr; endpts = 180*outsr;
achsig = downsample(achdata.values,round(1/(achdata.interval*outsr)));
achts = downsample(syncachts,round(1/(achdata.interval*outsr)));

achsig = achsig(edgepts:end-endpts);
achts = achts(edgepts:end-endpts);

p = polyfit(achts,achsig,5);
Ybaseline = polyval(p,achts);
dFoverF = smoothdata((achsig - Ybaseline)./Ybaseline,'movmean',5*outsr);

fig1 = figure;
plot(achts,achsig); hold on;
plot(achts,Ybaseline,'LineWidth',2);
plot(achts, dFoverF);
title('Ach dFoverF');
savefig(fig1,[datapath,session,filesep,'dFoverF']);

%% LPFs power
[spectrogram,chrt,freq] = MTSpectrogram(double(LFPs.data),'window',1);
bands = SpectrogramBands(spectrogram,freq);
thetapower = interp1(chrt-lfptime(1),pow2db(bands.theta),achts);
thetapower_sm = smoothdata(thetapower,'movmean',5*outsr);
% thetapower_sm(thetapower_sm>5e4) = nan;

fig2 = figure;
plot(achts,thetapower_sm);
title('Theta power');
savefig(fig2,[datapath,session,filesep,'Theta_Power']);


syncts = achts;
%% scatter plot
win = 5;
lfpids = length(syncts);
nwins = floor(length(syncts)/(win*outsr));
thetapower_sm_bined = nanmean(reshape(thetapower_sm(1:nwins*win*outsr),win*outsr,[]),1)';
dFoverF_bined = nanmean(reshape(dFoverF(1:nwins*win*outsr),win*outsr,[]),1)';

fig3 = figure;
scatter(dFoverF_bined,thetapower_sm_bined,'filled');
title('dFoverF vs. Theta power');
savefig(fig3,[datapath,session,filesep,'dFoverF_vs_Theta_Power']);

scatdata = [dFoverF_bined,thetapower_sm_bined];
save([datapath,session,filesep,'scatdata.mat'],'scatdata');

clear fig1 fig2 fig3;

% %% perievent
% [pks,locs,w,~] = findpeaks(thetapower_sm,syncts,'MinPeakWidth',5);
% 
% plot(syncts,thetapower_sm); hold on;
% for iloc = 1:length(locs)
%     plot([1,1]*locs(iloc),[0,1e5]);
% end
% 
% eventwidth = 30;
% achevents = zeros(0);
% for iloc = 5:length(locs)-5
%     tspts = [-eventwidth*outsr:eventwidth*outsr]./outsr + locs(iloc);
%     achevents(iloc,:) = interp1(syncts,achsig,tspts);
% end
% 
% figure;
% subplot(1,2,1)
% imagesc([-eventwidth*outsr:eventwidth*outsr]./outsr,5:length(locs)-5,zscore(achevents,[],2));
% subplot(1,2,2)
% plot([-eventwidth*outsr:eventwidth*outsr]./outsr,nanmean(achevents,1));
% 
% %% 
% plot(syncts,(thetapower_sm-nanmean(thetapower_sm))./nanstd(thetapower_sm)+3); hold on;
% plot(syncts,zscore(dFoverF)-3);