clear all;clc;

basepath = pwd;
files = dir('*.mat');
%% ach syn and caculation
ii = 2
[filepath,name,ext] = fileparts(files(ii).name);
% syn
load([name,ext],[name,'_OpTrTTL']);
Opti = eval([name,'_OpTrTTL']);
Opti.data = downsample(Opti.values,10);
Opti.time = downsample(Opti.times,10);
time_st = round(Opti.time(find(Opti.data>1,1))*2e4);time_nd = round(Opti.time(find(Opti.data>1,1,'last'))*2e4);
% signal
load([name,ext],[name,'_AChSig1']);
ACh_signal = eval([name,'_AChSig1']);
ACh_signal.data = downsample(ACh_signal.values(time_st:time_nd),200);
ACh_signal.time = downsample(ACh_signal.times(time_st:time_nd),200)- Opti.time(find(Opti.data>1,1));

p = [];
p = polyfit(ACh_signal.time,ACh_signal.data,5);
Ybaseline = polyval(p,ACh_signal.time);

figure
plot(ACh_signal.time, ACh_signal.data); hold on;
plot(ACh_signal.time,Ybaseline);
dFoverF = smoothdata((ACh_signal.data - Ybaseline)./Ybaseline,'movmean',10);
hold on;
plot(ACh_signal.time, dFoverF);

%% optogenetics
load([name,ext],[name,'_OpGeSti']);
Opto = eval([name,'_OpGeSti']);
Opto.data = downsample(Opto.values,5);
Opto.time = downsample(Opto.times,5);
Index = Opto.time(Opto.data>1);
pulse = Index(diff(Index)>10)-10;pulse(:,2) = pulse(:,1)+15;pulse = pulse-Opti.time(find(Opti.data>1,1));
ACh_Opto = [];
OPindex = round(pulse*100);
for jj = 1:length(pulse )
    ACh_Opto(jj,:) = dFoverF(OPindex(jj,1):OPindex(jj,2));
    ACh_Opto(jj,496:549)= NaN; ACh_Opto(jj,991:1016) = NaN;
    ACh_Opto(jj,:) = fillmissing(ACh_Opto(jj,:),'linear');
end
h = figure;
options.color_area = [128 193 219]./255;    % Blue theme
options.color_line = [ 52 148 186]./255;
%options.color_area = [243 169 114]./255;    % Orange theme
%options.color_line = [236 112  22]./255;
options.alpha      = 0.5;
options.line_width = 2;
options.error      = 'sem';
plot_areaerrorbar(ACh_Opto, options);set(gcf,'Position',  [100, 100, 400, 300])
xlim([0 1500]);title(name)
savefig(h,[name,'.fig']);
print(h,name,'-dsvg');
save([name,'_Peri','.mat'],'ACh_Opto');

%% bahvior state and syn
info = read_Intan_RHD2000_file(basepath);
digitalin = read_Intan_digitalin_file('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;

load(files(4).name);
SleepState.ints.REMstate = SleepState.ints.REMstate - lfptime(1);
SleepState.ints.NREMstate = SleepState.ints.NREMstate - lfptime(1);
SleepState.ints.WAKEstate = SleepState.ints.WAKEstate - lfptime(1);
REM.opto = [];NREM.opto = [];WAKE.opto = [];
for jj = 1:length(SleepState.ints.REMstate)
    REM.optoin = find(pulse(:,1)>SleepState.ints.REMstate(jj,1)& pulse(:,2)< SleepState.ints.REMstate(jj,2))
    REM.opto = [REM.opto;REM.optoin];
end
for jj = 1:length(SleepState.ints.NREMstate)    
    NREM.optoin = find(pulse(:,1)>SleepState.ints.NREMstate(jj,1)& pulse(:,2)< SleepState.ints.NREMstate(jj,2))
    NREM.opto = [NREM.opto;NREM.optoin]
end
for jj = 1:length(SleepState.ints.WAKEstate)
    WAKE.optoin = find(pulse(:,1)>SleepState.ints.WAKEstate(jj,1)& pulse(:,2)< SleepState.ints.WAKEstate(jj,2))
    WAKE.opto = [WAKE.opto;WAKE.optoin]
end
ACh_REMopto = ACh_Opto(REM.opto,:);
ACh_NREMopto = ACh_Opto(NREM.opto,:);
ACh_WAKEopto = ACh_Opto(WAKE.opto,:);

g = figure(2);
options.color_area = [128 193 219]./255;    % Blue theme
options.color_line = [ 52 148 186]./255;
options.alpha      = 0.5;
options.line_width = 2;
options.error      = 'sem';
plot_areaerrorbar(ACh_REMopto, options);set(gcf,'Position',  [100, 100, 400, 300])
xlim([0 1500]);title('REM')
savefig(g,'REM.fig');
print(g,'REM','-dsvg');
save([name,'REM_Peri','.mat'],'ACh_REMopto');

o = figure(3);
options.color_area = [128 193 219]./255;    % Blue theme
options.color_line = [ 52 148 186]./255;
options.alpha      = 0.5;
options.line_width = 2;
options.error      = 'sem';
plot_areaerrorbar(ACh_NREMopto, options);set(gcf,'Position',  [100, 100, 400, 300])
xlim([0 1500]);title('NREM')
savefig(o,'NREM.fig');
print(o,'NREM','-dsvg');
save([name,'NREM_Peri','.mat'],'ACh_NREMopto');

k = figure(4);
options.color_area = [128 193 219]./255;    % Blue theme
options.color_line = [ 52 148 186]./255;
options.alpha      = 0.5;
options.line_width = 2;
options.error      = 'sem';
plot_areaerrorbar(ACh_WAKEopto, options);set(gcf,'Position',  [100, 100, 400, 300])
xlim([0 1500]);title('WAKE')
savefig(k,'WAKE.fig');
print(k,'WAKE','-dsvg');
save([name,'WAKE_Peri','.mat'],'ACh_WAKEopto');

%% LFP analysis
analogin = read_Intan_analogin_file('analogin.dat',info.num_board_adc_channels);
analogin = downsample(analogin,16);
load(files(3).name);
thetaCh = SleepScoreLFP.THchanID;
rippleCh = SleepScoreLFP.SWchanID;

THLFPs = bz_GetLFP(12,'basepath',basepath,'noPrompts',true);
thetarange = linspace(0,250,1000);
lfp = double(THLFPs.data);
lfp(abs(lfp)>5000) = nan;
lfp = fillmissing(lfp,'linear');
[ss,ft,tt] = spectrogram(lfp,2*THLFPs.samplingRate,0.6*THLFPs.samplingRate,thetarange,THLFPs.samplingRate);
[spectrogrammt,tm,fm] = MTSpectrogram(lfp,'window',1);
bands = SpectrogramBands(spectrogrammt,fm);

figure;
colormap('jet')
imagesc(tt,thetarange,abs(ss));clim([0,1e5])
yticks([4,8,12]); ylabel('Frequency (Hz)'); xlabel('Time (s)');
set(gca,'YDir','normal');
hold on; plot(THLFPs.timestamps,analogin,'color','r')

%%%%%%%%%%%%%%%%%REM TIME
Index = THLFPs.timestamps(analogin>1);
pulse = Index(diff(Index)>10)-20;pulse(:,2) = pulse(:,1)+30;

REM_time = pulse(REM.opto,:);
lfp_remC = [];  lfp_remS = [];lfp_remP = [];
for ii = 1:length(REM_time)
    lfp_remCin = lfp(THLFPs.timestamps > REM_time(ii,1)-5 & THLFPs.timestamps < REM_time(ii,1));
    lfp_remC = [lfp_remC;lfp_remCin];
    lfp_remSin = lfp(THLFPs.timestamps > REM_time(ii,1)+3 & THLFPs.timestamps < REM_time(ii,2)+2 ); %% final 2 second
    lfp_remS = [lfp_remS;lfp_remSin];
    lfp_remPin = lfp(THLFPs.timestamps > REM_time(ii,2)+2 & THLFPs.timestamps < REM_time(ii,2)+5);
    lfp_remP = [lfp_remP;lfp_remPin];
end
%%%%%%%%%%%%%%%%%NREM TIME
NREM_time = pulse(NREM.opto,:);
lfp_NremC = [];  lfp_NremS = [];lfp_NremP = [];
for ii = 1:length(NREM_time)
    lfp_NremCin = lfp(THLFPs.timestamps > NREM_time(ii,1)-5 & THLFPs.timestamps < NREM_time(ii,1));
    lfp_NremC = [lfp_NremC;lfp_NremCin];
    lfp_NremSin = lfp(THLFPs.timestamps > NREM_time(ii,1)+3 & THLFPs.timestamps < NREM_time(ii,2));
    lfp_NremS = [lfp_NremS;lfp_NremSin];
    lfp_NremPin = lfp(THLFPs.timestamps > NREM_time(ii,2) & THLFPs.timestamps < NREM_time(ii,2)+5);
    lfp_NremP = [lfp_NremP;lfp_NremPin];
end
%%%%%%%%%%%%%%%%%WAKE TIME
WAKE_time = pulse(WAKE.opto,:);
lfp_WAKEC = [];lfp_WAKES = [];lfp_WAKEP = [];
for ii = 1:length(WAKE_time)
    lfp_WAKECin = lfp(THLFPs.timestamps > WAKE_time(ii,1)-5 & THLFPs.timestamps < WAKE_time(ii,1));
    lfp_WAKEC = [lfp_WAKEC;lfp_WAKECin];
    lfp_WAKESin = lfp(THLFPs.timestamps > WAKE_time(ii,1)+2.5 & THLFPs.timestamps < WAKE_time(ii,2)+0.5);
    lfp_WAKES = [lfp_WAKES;lfp_WAKESin];
    lfp_WAKEPin = lfp(THLFPs.timestamps > WAKE_time(ii,2)+2 & THLFPs.timestamps < WAKE_time(ii,2)+5);
    lfp_WAKEP = [lfp_WAKEP;lfp_WAKEPin];
end
%%%%%%%%%%%PSD 
window = 2048;
noverlap = window/2;
fs = 1250
[pxx,f] = pwelch(lfp_WAKEC,window,noverlap,500,fs)
figure; plot(f,10*log10(pxx));xlim([0 30])
[pxx,f] = pwelch(lfp_WAKES,window,noverlap,500,fs)
hold on; plot(f,10*log10(pxx));xlim([0 30])
[pxx,f] = pwelch(lfp_WAKEP,window,noverlap,500,fs)
hold on; plot(f,10*log10(pxx));xlim([0 30])

[pxx,f] = pwelch(lfp_NremC,window,noverlap,500,fs)
figure; plot(f,10*log10(pxx));xlim([0 30])
[pxx,f] = pwelch(lfp_NremS,window,noverlap,500,fs)
hold on; plot(f,10*log10(pxx));xlim([0 30])
[pxx,f] = pwelch(lfp_NremP,window,noverlap,500,fs)
hold on; plot(f,10*log10(pxx));xlim([0 30])

[pxx,f] = pwelch(lfp_remC,window,noverlap,500,fs)
figure; plot(f,10*log10(pxx));xlim([0 30])
[pxx,f] = pwelch(lfp_remS,window,noverlap,500,fs)
hold on; plot(f,10*log10(pxx));xlim([0 30])
[pxx,f] = pwelch(lfp_remP,window,noverlap,500,fs)
hold on; plot(f,10*log10(pxx));xlim([0 30])

%% 






