clear all;

datafile = {'F:\LAMP2\LAMP2_042_190614_120612',
            'F:\LAMP3\LAMP3_010_190617_120744',
            'F:\LAMP3\LAMP3_016_190619_180636',
            'F:\LAMP3\LAMP3_019_190620_153139',
            'F:\LAMP3\LAMP3_024_190701_124743',
            'F:\LAMP37\LAMP37_009_201211_120544',  %%start 11.530 opti in intan
            'F:\LAMP37\LAMP37_011_201212_192627',  %%start 17.560 opti in intan
            'F:\LAMP37\LAMP37_012_201213_163237'} %%start 18.539 opti in intan
            
thCh = [53,16,4,16,16,1,1,1];
ALItime = [0,0,0,0,0,11.530,17.563,18.539];

datafile = {'F:\LIMP003\LIMP003_008_210530_142215'}
           
thCh = [11];

for ii = 1:length(datafile)
    cd(datafile{ii,:})
    
    basepath = pwd;
    basename = bz_BasenameFromBasepath(basepath);

    % get alighment
    info = read_Intan_RHD2000_file(basepath);
    digitalin = read_Intan_digitalin_file('digitalin.dat',info.num_board_dig_in_channels)';
    
    info = read_Intan_RHD2000_file(basepath);
    if ii < 6
        digitalin = read_Intan_digitalin_file('digitalin.dat',info.num_board_dig_in_channels)';
        lfp_fbtime(1) = find(digitalin(1,:)>(max(digitalin)+min(digitalin))/2,1);lfp_fbtime(2) = find(digitalin(1,:)>(max(digitalin)+min(digitalin))/2,1,'last');
    else 
        lfp_fbtime(1) = round(ALItime(ii)*2e4);
    end
    % get name 
    [filepath,name,ext] = fileparts(basename);
    [k y] = strtok(name,'_');
    name = [k,y(1:4)];
    % get syn
    if exist('dFoverF.mat')==0
        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));
        ACh_signal.data  = fillmissing(ACh_signal.data ,'linear');
        
        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);xlim([0 inf])
    else load('dFoverF.mat')
    end
    %% get LFPs
    THLFPs = bz_GetLFP(thCh(ii),'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,5e5]);ylim([0 50]);xlim([300 inf])
    yticks([0.5,20,50,100,250]); ylabel('Frequency (Hz)'); xlabel('Time (s)');
    set(gca,'YDir','normal');
    hold on;plot(ACh_signal.time+lfp_fbtime(1)/20000, (dFoverF)*100+35,'w-');xlim([300 inf])
    hold on;plot(tm, bands.theta/20000+10,'g-');xlim([300 inf])
%%

    datpath = 'G:\My Drive\Paper_ACh_OXT\Fig1\Data\Figure1F\ACh_theta'
    save(fullfile(datpath,[basename,'_ACh_theta.mat']),'THLFPs','lfp_fbtime','dFoverF');
end
