
%% 
clear all;

dataname = {'F:\LIMP014\LIMP014_012_211005_112625',  %% optitrack  old data
            'F:\LIMP014\LIMP014_014_211007_122021',  %% optitrack   old data
            'F:\LIMP016\LIMP016_004_220101_144210',  %% analogCh  = 2    old data
            'F:\LIMP018\LIMP018_003_220127_150338',
            'F:\LAMP2\LAMP2_002_190405_154550',            
            'F:\LAMP37\LAMP37_006_201201_192319',
            'F:\LIMP011\LIMP011_010_210705_120340',
            'F:\LIMP013\LIMP013_009_210630_101849'};  %% exclude data from 5700 to 11250    old data

RipCh = [57,57,53, 26 ,29, 1, 35, 2];
GammaCh = [17,17,46,122, 29, 0, 29, 27];


for ii = 1:size(dataname,1)
    cd(dataname{ii,:})
    basepath = pwd;
    basename = bz_BasenameFromBasepath(basepath);
    load([basename,'.SleepState.states.mat']);
    
    if ii>0 & ii<3
        info = read_Intan_RHD2000_file(basepath);
        digitalin = read_Intan_digitalin_file('digitalin.dat',info.num_board_dig_in_channels)';
        lfp_fbtime = find(digitalin(1,:)>(max(digitalin)+min(digitalin))/2,1);
        if isempty(lfp_fbtime) == 1
            lfp_fbtime = [0]
        else lfp_fbtime(1,2) = find(digitalin(1,:)>0.6,1,'last');
        end
        lfp_fbtime = lfp_fbtime*3;
    elseif ii > 2 & ii < 5
        info = read_Intan_RHD2000_file(basepath);
        analogin= read_Intan_analogin_file('analogin.dat',info.num_board_adc_channels);
        if ii == 3
            DD = analogin(2,:);
        elseif ii > 3
            DD = analogin(1,:);
        end
        lfp_fbtime = find(DD>(max(DD)+min(DD))/2,1);
       
        if isempty(lfp_fbtime) == 1
            lfp_fbtime = [0]
        else lfp_fbtime(1,2) = find(DD>0.6,1,'last');
        end 
    else
        info = read_Intan_RHD2000_file(basepath);
        digitalin = read_Intan_digitalin_file('digitalin.dat',info.num_board_dig_in_channels)';
        lfp_fbtime = find(digitalin(1,:)>(max(digitalin)+min(digitalin))/2,1);
        lfp_fbtime(1,2) = find(DD>0.6,1,'last');
    end
    % sitmulation and signal
    if exist('dFoverF.mat')==0
        if ii == 5
            name = basename(1:9);
        elseif ii == 6
            name = basename(1:10);
        else
            name = basename(1:11);
        end
  
        ext = '.mat'; 
        if ii ~= 3 & ii ~= 4 & ii ~= 5
            load([name,ext],[name,'_OpTrTTL']);
            ACh_sti = eval([name,'_OpTrTTL']);

            load([name,ext],[name,'_AChSig1']);
            ACh_signal = eval([name,'_AChSig1']);

        elseif ii == 5
            load([name,ext],[name,'_optitrck']);
            ACh_sti = eval([name,'_optitrck']);
            load([name,ext],[name,'_signal']);
            ACh_signal = eval([name,'_signal']);

        else             
            load([name,ext],[name,'_AChSti']);
            ACh_sti = eval([name,'_AChSti']);

            load([name,ext],[name,'_AChSig1']);
            ACh_signal = eval([name,'_AChSig1']);
        end

        if ii == 5
            otr = 0.03;
        else
            otr = 0.5;
        end
        fbtime = ACh_sti.times(find(ACh_sti.values(:,1)>otr,1));fbtime(1,2) = ACh_sti.times(find(ACh_sti.values(:,1)>otr,1,'last'));
        fbIDX = round(fbtime*20000);
        
        ACh_signal.data = downsample(ACh_signal.values((fbIDX(1,1)):(fbIDX(1,2))),200);
        ACh_signal.time = downsample(ACh_signal.times((fbIDX(1,1)):(fbIDX(1,2))),200)- fbtime(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);xlim([0 inf]);
    else load('dFoverF.mat'); figure;plot(dFoverF);
    end
    
%% get LFP
    THLFPs = bz_GetLFP(GammaCh(ii),'basepath',basepath,'noPrompts',true);
    THLFPs.data = double(THLFPs.data);
    THLFPs.data(abs(THLFPs.data)>5000) = nan;
    THLFPs.data = fillmissing(THLFPs.data,'linear');

    RipLFPs = bz_GetLFP(RipCh(ii),'basepath',basepath,'noPrompts',true);
    RipLFPs.data = double(RipLFPs.data);
    RipLFPs.data(abs(RipLFPs.data)>5000) = nan;
    RipLFPs.data = fillmissing(RipLFPs.data,'linear');
%% chcek data;
    thetarange = linspace(0.5,250,1000);
    lfp = double(THLFPs.data);
    lfp(abs(lfp)>8000) = 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',2);
    bands = SpectrogramBands(spectrogrammt,fm);
    
    % for figure B
    figure;
    colormap('jet')
    imagesc(tt,thetarange,abs(ss));clim([0,3e5]);
    yticks([0.5,20,50,100,200]); ylabel('Frequency (Hz)'); xlabel('Time (s)');
    set(gca,'YDir','normal');
    hold on;plot([1:length(dFoverF)]/100+lfp_fbtime(1)/20000, dFoverF*500+100,'w');xlim([0 inf]);
    hold on; plot(tm, bands.highGamma/100+50)
%% save data

    datpath = 'G:\My Drive\Papers_Buzsaki\Paper_ACh_OXT\RevisionFigures\REM_Gamma-analysis';
    save(fullfile(datpath,[basename,'ACh_LFP_REM.mat']),'RipLFPs','THLFPs','lfp_fbtime', 'dFoverF')

end

