cd ('G:\My Drive\Paper_ACh_OXT\Fig2\nonNREMpackets_data');

clear all;
%%
files = dir('*.mat');

%% get ACh(1:5).OXT(6:9) and LFP 
NREM_spindles_all = [];
nREM_ACh_all = [];
for ii = 6:9
    
    load(files(ii).name)
    
    % exclude data shorter than 6 seconds
    idx = find(diff(nNREM_Data.nNREMepsode.ints,1,2)>100);
    idx = idx(idx<length(nNREM_Data.nNREMepsode.LFP));
    nREM_ACh = nNREM_Data.nNREMepsode.ACh(idx,1);
    nREM_LFP = nNREM_Data.nNREMepsode.LFP(idx,1);
    
    % get spectrum NREM packects
    THLFPs.samplingRate = 1250;
    thetarange = linspace(1,100,250);
    NREM_spindles = [];
    for jj = 1:length(nREM_LFP)
        if ~isempty(nREM_LFP{jj,1})>0
            lfp = nREM_LFP{jj,1};
            [spectrogrammt,tm,fm] = MTSpectrogram(lfp,'window',2);
            bands = SpectrogramBands(spectrogrammt,fm);
            NREM_spindles{jj,1} = bands.spindles;
        end
    end
    ii
    NREM_spindles_all = [NREM_spindles_all;NREM_spindles ];
    nREM_ACh_all = [nREM_ACh_all;nREM_ACh];
    
end
%%
pxx_all_sp = [];pxx_all_df =[]; Cxy_all = [];rho_all =[];xcf_all = [];
for kk = 1: length(nREM_ACh_all)
    fs = 10;
    ups = 10;
    dns = 1;
    y = resample(NREM_spindles_all{kk,1},ups,dns);
    x = resample(nREM_ACh_all{kk,1},dns,ups);
    y = bandpass(y,[0.01 9],fs);
    x = bandpass(x,[0.01 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(10),[],[],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] = xcorr(y(1:length(x)),x,250,'normalized');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] = xcorr(y,x(1:length(y)),250,'normalized');plot(lags,xcf)
    end
    close all
    Cxy_all = [Cxy_all,Cxy];
    rho_all = [rho_all,rho];
     xcf_all = [xcf_all,xcf];
end

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