clear all;clc   %% check the alighment 
%% STATISTIC OXYTOCIN BEHAVIRO
dataname = ['F:\LIMP013\LIMP013_009_210630_101849',
            'F:\LIMP014\LIMP014_013_211006_133329',
            'F:\LIMP014\LIMP014_014_211007_122021',
            'F:\LIMP018\LIMP018_003_220127_150338', % align with analogin
            'F:\LIMP019\LIMP019_004_220303_115500'] % align with analogin
Cxy_all = [];pxx_all_sp = [];pxx_all_df = [];rho_all = [];xcf_all =[];
Ch = [57;59;34;31];
for ii = 1:size(dataname,1)
    cd(dataname(ii,:))
    basepath = pwd;
    basename = bz_BasenameFromBasepath(basepath);
    load([basename,'.SleepState.states.mat']);
    load([basename,'.SleepScoreLFP.LFP.mat']);
    load([basename,'.accel.behavior.mat']);

        
    info = read_Intan_RHD2000_file(basepath);
    if ii<4
        digitalin = read_Intan_digitalin_file('digitalin.dat',info.num_board_dig_in_channels);
        lfp_fbtime = find(digitalin> (min(digitalin)+max(digitalin))/2,1);lfp_fbtime(1,2) = find(digitalin> (min(digitalin)+max(digitalin))/2,1,'last');
        delete digitalin
    else 
        analogin = read_Intan_analogin_file('analogin.dat',info.num_board_adc_channels);
        if ii == 3;
        lfp_fbtime = find(analogin(2,:)> (min(analogin(2,:))+max(analogin(2,:)))/2,1);lfp_fbtime(1,2) = find(analogin(2,:)> (min(analogin(2,:))+max(analogin(2,:)))/2,1,'last');
        delete analogin
        elseif ii == 4
        lfp_fbtime = find(analogin(1,:)> (min(analogin(1,:))+max(analogin(1,:)))/2,1);lfp_fbtime(1,2) = find(analogin(1,:)> (min(analogin(1,:))+max(analogin(1,:)))/2,1,'last');
        delete analogin    
        end
    end
    % sitmulation and signal
    if exist('dFoverF.mat')==0
        name = basename(1:11);
        ext = '.mat'
         if ii<4
            load([name,ext],[name,'_OpTrTTL']);
            ACh_sti = eval([name,'_OpTrTTL']);
            fbtime = ACh_sti.times(find(ACh_sti.values(:,1)>1,1));fbtime(1,2) = ACh_sti.times(find(ACh_sti.values(:,1)>1,1,'last'));
        else 
            load([name,ext],[name,'_AChSti']);
            ACh_sti = eval([name,'_AChSti']);
            fbtime = ACh_sti.times(find(ACh_sti.values(:,1)>0.5,1));fbtime(1,2) = ACh_sti.times(find(ACh_sti.values(:,1)>0.5,1,'last'));
        end
        fbIDX = round(fbtime*20000);
        load([name,ext],[name,'_AChSig1']);
        ACh_signal = eval([name,'_AChSig1']);
        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,2);
        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
    
    rippleCh = Ch(ii);
    THLFPs = bz_GetLFP(rippleCh,'basepath',basepath,'noPrompts',true);
    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);
    figure;
    colormap('jet')
    imagesc(tt,thetarange,abs(ss));clim([0,8e4]);
    yticks([0.5,20,50,100,200]); ylabel('Frequency (Hz)'); xlabel('Time (s)');
    set(gca,'YDir','normal');
    hold on; plot(tm,(smoothdata(bands.spindles,10))/1000+110,'color','r');xlim([0 inf]);
    hold on; plot([1:length(dFoverF)]/100+lfp_fbtime(1)/20000, dFoverF*100+100,'g');xlim([0 inf])
    hold on; plot(accel.timestamps/2, accel.acceleration.motion/50+150,'w');xlim([0 inf]);
    figure;plot(SleepState.idx.states);xlim([0 inf]);
    %% get NREM packet and micro arousal(accelrate > 0.5 and <40)
    statesIDX = bz_INTtoIDX(SleepState.ints);  %% need to check
    ix = linspace(1, numel(statesIDX.states), numel(accel.acceleration.motion));
    Mov = interp1(ix, accel.acceleration.motion, 1:numel(statesIDX.states));
    figure;plot(50*statesIDX.states-50);xlim([0 inf]);
    hold on; 
    %plot(Mov/50);xlim([0 inf]);
    idx_ = find(statesIDX.states ~= 2);
    %%%%%%%%%%%%%%%%%%%%%%%%%%%get microarousals 
    MovNREM1 = Mov;
    MovNREM1(idx_) = 0;
    plot(MovNREM1/10);xlim([0 inf]);
    
    MovNREM = smoothdata(MovNREM1,"gaussian",5);

    peak_mean_diff = 15;
    threshold_ = 20;
    [pks, locs, width] = findpeaks(MovNREM,'MinPeakProminence',threshold_,...
        'WidthReference','halfprom','MinPeakDistance',peak_mean_diff);
    
    figure; plot(MovNREM);
    hold on;
    plot(locs,pks,'o');
    plot(MovNREM1)

    [arousals] = FindMicroArousal(statesIDX, Mov);

    figure; plot(arousals.times,arousals.movNREM,'k');
    hold on;
    xline(arousals.timestamps(:,1),'b')
    xline(arousals.timestamps(:,2),'r')

    %%%%%%%%%%%%%%%%%get NREM packets
    nREMnoArousals = ones(size(Mov));
    
    for i = 1:length(arousals.peaks)
        idx = arousals.timestamps(i,1):arousals.timestamps(i,2);
        nREMnoArousals(idx) = 0;
    end
    idx_ = find(statesIDX.states ~= 2);
    nREMnoArousals(idx_) = 0;
    
    figure;plot(nREMnoArousals*2500);
    hold on
    plot(arousals.times,Mov,'k');
    plot(1000*statesIDX.states,'r');xlim([0 inf]);
    nNREMpackets = cell2mat(IDXtoINT(nREMnoArousals));

    %% get nREMpackets of LFP 
    nNREMpac_Data = {};
    for jj = 1:length(nNREMpackets)
        Dur = nNREMpackets(jj,:)*1250;
        Dur_ach = round((nNREMpackets(jj,:)-lfp_fbtime(1)/20000)*100);
        if Dur_ach(1)>0 & Dur_ach(2)<length(dFoverF)
            nNREMpac_Data.LFP{jj,1} = lfp(Dur(1):Dur(2));
            nNREMpac_Data.ACh{jj,1} = dFoverF(Dur_ach(1):Dur_ach(2));
        end
    end
    nNREMpac_Data.ints = nNREMpackets;
   
    %% save data


end

    %% find cohenrence
    fs = 10;
    ups = 10;
    dns = 1;
    y = resample(bands.delta,ups/2,dns);
    %y = resample(bands.spindles,ups/2,dns);
    x = resample(dFoverF_SWS.data,dns,ups);
    y = bandpass(y,[0.001 9],fs);
    x = bandpass(x,[0.001 9],fs);

    window = 40;
    noverlap = window/2;
    [pxx,f] = pwelch(y,window,noverlap,500,fs);
    figure; plot(f,normalize(10*log10(pxx)));
    pxx_all_sp = [pxx_all_sp,10*log10(pxx)];

    [pxx,f] = pwelch(x,window,noverlap,500,fs);
    figure; plot(f,normalize(10*log10(pxx)));
    pxx_all_df = [pxx_all_df,10*log10(pxx)];

    f_all=f;

    
    if length(x)<length(y)
        [Cxy,f] = mscohere(x,y(1:length(x)),hamming(500),[],[],fs,'mimo');
        figure;plot(f,Cxy);
        rho = corr(x,y(1:length(x)));
        mdl = fitlm(x,y(1:length(x)));
        figure;plot(mdl)
        figure;[xcf,lags] = crosscorr(y(1:length(x)),x,'NumLags',600,'NumSTD',3);plot(lags,xcf)
    elseif length(x)>length(y)
        [Cxy,f] = mscohere(x(1:length(y)),y,hamming(500),[],[],fs,'mimo');
        figure;plot(f,Cxy)
        rho = corr(x(1:length(y)),y);
        mdl = fitlm(x(1:length(y)),y);
        figure;plot(mdl)
        figure;[xcf,lags] = crosscorr(y,x(1:length(y)),'NumLags',600,'NumSTD',3);plot(lags,xcf)
    end
    close all
    Cxy_all = [Cxy_all,Cxy];
    rho_all = [rho_all,rho];
    xcf_all = [xcf_all,xcf];

figure;
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      = 'std';
options.x_axis = f;
plot_areaerrorbar((Cxy_all(:,[1,3,4]))', options);set(gcf,'Position',  [100, 100, 400, 300])
xlim([0 max(f)]);title('MScoherence')

figure;
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      = 'std';
options.x_axis =  f_all;
plot_areaerrorbar((normalize(pxx_all_sp)'), options);set(gcf,'Position',  [100, 100, 400, 300])
xlim([0 max(f)]);title('spindle_psd')

figure;
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      = 'std';
options.x_axis =  f_all;
plot_areaerrorbar((normalize(pxx_all_df)'), options);set(gcf,'Position',  [100, 100, 400, 300])
xlim([0 max(f)]);title('dF_psd')

figure;
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';
options.x_axis =  lags;
plot_areaerrorbar((xcf_all)', options);set(gcf,'Position',  [100, 100, 400, 300])
title('cross-corr')