set(0, 'DefaultLineLineWidth', 2);
set(0,'defaultTextFontSize', 10);
set(0,'defaultAxesFontSize',10);
set(0,'DefaultFigureColormap',jet)
set(0,'defaultfigurecolor',[1 1 1])
set(0,'DefaultAxesTitleFontWeight', 'normal')
set(0, 'defaulttextinterpreter', 'none') % or 'tex'

%% Load the data
clear all;clc
filename = dir('*_OXT_Arti_ripples.mat');
for ii = 7:length(filename)
    %% load data
    basename = filename(ii).name;
    load(basename);
    %%
    fs = Sp2.FluoData.sampling_rate;
    LaserShutter = double(Sp2.FluoData.laser_times); %%camsignal is the laser signal, digital_2 is cam signal
    
    %%
    OXT_dff = double(Sp2.FluoData.dFF_OXT);
    figure;plot(Sp2.FluoData.time_seconds,OXT_dff+0.1);hold on
    peakLocation = round(0.5*(LaserShutter(:,1)+LaserShutter(:,2))*fs);
    samplesToRemove = 16;
    for i = 1:length(peakLocation)
        % Calculate the range of samples to remove
        rangeStart = max(1, peakLocation(i) - samplesToRemove);
        rangeEnd = min(length(OXT_dff), peakLocation(i) + round(1.5*samplesToRemove));
        
        % Remove the samples within the defined range
        OXT_dff(rangeStart:rangeEnd) = NaN;
    end
    OXT_dff = fillmissing(OXT_dff, 'nearest');
    
    plot(Sp2.FluoData.time_seconds,OXT_dff);
    
%% random chosse time points
    data = Sp2.FluoData.time_seconds; % replace with your 1xN data
    numPoints = length(peakLocation); % number of random time points to select
    numIterations = 1000; % number of times to select random time points
    
    % Initialize a cell array to hold the results
    dFoverF_NREM_Session = [];
    % Loop over the number of iterations
    for i = 1:numIterations
        %% alighed dff with cam signal
        cam_signal = double(Sp2.FluoData.cam_signal); 
        % Define your threshold
        threshold = 0.5;                                                        
        % Calculate the difference between consecutive elements
        diff_signal = diff(cam_signal > threshold);
        % Find the start times (where diff_signal is 1)
        start_times = find(diff_signal == 1,1) + 1;
        % Find the end times (where diff_signal is -1)
        end_times = find(diff_signal == -1,1,'last');
        dFoverF_ds_aligned = OXT_dff(start_times:end_times);

        % Generate random indices
        indices = randi([start_times end_times], 1, numPoints);
    
        % Select random time points and store in the cell array
        randomTimePoints = data(indices);

        randomTimePoints_aligned = randomTimePoints - Sp2.FluoData.time_seconds(start_times);
        randomTimePoints_aligned = randomTimePoints_aligned' ;

        Npp = 1;
        NREM_flag = 3;
        WAKE_flag = 1;
        REM_flag = 5;
        pre = 5; % in seconds
        post = 5;  % in seconds
        smoothing_sigma = 15; % smoothing sigma in samples
        lfp_time = Sp2.lfp_time;
        randomTimePoints_synched= round(randomTimePoints_aligned + lfp_time(1));
        randomTimePoints_synched = randomTimePoints_synched(randomTimePoints_synched>0);
        SleepState = Sp2.SleepState;
        randomTimePoints_sleepstate = SleepState.idx.states(randomTimePoints_synched);
        stim_locs = randomTimePoints_synched*fs;
        %% get data of NREM
        pre = 5;
        post = 5;
        eventAveragedDF = []; % Initialize array for event-averaged dFoverF
        Ns = length(OXT_dff);
        windowSize = fs; % Define the window size around each event time
        dFoverFWindow = [];
        stim_NREM  = stim_locs(randomTimePoints_sleepstate==NREM_flag & circshift(randomTimePoints_sleepstate,5)==NREM_flag & circshift(randomTimePoints_sleepstate,-5)==NREM_flag);
        dFoverF_NREM = zeros(length(stim_NREM), pre*windowSize+post*windowSize+1);
        %stim_locs = round(0.5*(pulse_times_aligned(:,1)+pulse_times_aligned(:,2))*130);
        dFoverF_NREM=[];
        for k = 1:Npp:length(stim_NREM)-1
            eventTime = stim_NREM(k); % Get the current event time
            % Extract dFoverF values within the window around the event time
            windowStart = max(1, eventTime - pre*windowSize);
            windowEnd = min(length(dFoverF_ds_aligned), eventTime + post*windowSize);
            if windowEnd<length(dFoverF_ds_aligned) 
                dFoverF_NREM(k,:) = dFoverF_ds_aligned(windowStart:windowEnd) - mean(dFoverF_ds_aligned(windowStart+fs*(pre-3):windowStart+fs*pre));
            end
        end

        dFoverF_NREM_Session(i,:) = mean(dFoverF_NREM,1);
    end
    dFoverF_NREM_all{ii} = zscore(mean(dFoverF_NREM_Session,1));
end


