%% Ripple detection
clearvars -except subnum sub Summary
addpath \\mycloudpr4100\ivan\Ivan_analysis\Ivan\human-20190516T170723Z-001\human\FMAToolbox\Helpers
addpath  \\mycloudpr4100\ivan\Ivan_analysis\Ivan\human-20190516T170723Z-001\human\FMAToolbox\Analyses
addpath  \\mycloudpr4100\ivan\single_unit\human\wcst\awake_swr
addpath C:\HaoxinAnalyses\EMOP\Data\codes

subnum = 0;
for sub = {'IR29','IR35','IR51','IR54','IR57' ,'IR75','IR84'}    
    clearvars -except subnum sub Summary
    subnum=subnum+1;
    sub=char(sub)
    cd(['C:\HaoxinAnalyses\EMOP\Data\' sub])

    load([sub '_HaoxinROI.mat'])
    try
        MTL_roi.EC_Stan;
    catch
        MTL_roi.EC_Stan=[];
    end
    try
        MTL_roi.PRC_Stan;
    catch
        MTL_roi.PRC_Stan=[];
    end
    try
        MTL_roi.PHC_Stan;
    catch
        MTL_roi.PHC_Stan=[];
    end
    try
        MTL_roi.Sub_Stan;
    catch
        MTL_roi.Sub_Stan=[];
    end
    try
        MTL_roi.AMY_Stan;
    catch
        MTL_roi.AMY_Stan=[];
    end
    
    chan_ROI = unique([MTL_roi.CA1_Stan MTL_roi.CA3_Stan]);
    
    
    
    load([sub '.mat'])
    load([sub '_evnt_clean.mat'])
    load([sub '_evnt.mat'])
    load([sub '_log_EMOPSData.mat'])
    
    
    Fs=2000;
    raw_eeg_encoding = extractTrials(data_ecog,event_study_onsets(ok_epochs_study),chan_ROI,[-1 5]*Fs);
    raw_eeg_retrieval = extractTrials(data_ecog,event_test_onsets(ok_epochs_test),chan_ROI,[-1 5]*Fs);
    
    trials_time_axis_encoding = ( -Fs+1 : 5*Fs )'; trials_time_axis_encoding = repmat(trials_time_axis_encoding,1,size(raw_eeg_encoding,3));
    trials_time_axis_retrieval = ( -Fs+1 : 5*Fs )'; trials_time_axis_retrieval = repmat(trials_time_axis_retrieval,1,size(raw_eeg_retrieval,3));
    
    trials_ind_mat_encoding = repmat(ok_epochs_study,size(raw_eeg_encoding,2),1);
    trials_ind_mat_retrieval = repmat(ok_epochs_test,size(raw_eeg_retrieval,2),1);
    
    N_chanel = length(chan_ROI);
    for nch = 1:N_chanel
        tic
        %% Encoding
        lfp_task_range = reshape(raw_eeg_encoding(nch,:,:),1,[]);
        trials_time_axis = reshape(trials_time_axis_encoding,1,[]);
        trials_ind_mat = reshape(trials_ind_mat_encoding ,1,[]);
        [ripples_encoding] = detect_awake_ripple_haoxin(lfp_task_range, trials_time_axis, Fs);
        ripple_encoding_ch(nch).ripples=ripples_encoding; ripple_peaks_ind = round(ripples_encoding(:,2).*2000);
%         
%         for j = 1:numel(ripple_peaks_ind)
%             if ripple_peaks_ind(j) > 2000 & (ripple_peaks_ind(j) + 2000) < length(lfp_task_range)
%                 ripple_encoding_ch(nch).ripple_traces(j,:) = lfp_task_range((ripple_peaks_ind(j)-2000):(ripple_peaks_ind(j)+2000));
%                 ripple_encoding_ch(nch).trial_time(j) = trials_time_axis(ripple_peaks_ind(j));
%                 ripple_encoding_ch(nch).trial(j) = trials_ind_mat(ripple_peaks_ind(j));
%             else
%                 ripple_encoding_ch(nch).ripple_traces(j,1:4001) = NaN;
%                 ripple_encoding_ch(nch).trial_time(j) = NaN;
%                 ripple_encoding_ch(nch).trial(j) = NaN;
%             end
%         end

         %% Simulation analyses 
        Nfft = length(lfp_task_range) ; 
        psd = fft(lfp_task_range, Nfft); 
        psd = abs(psd/Nfft).^2 ; 
        psd = psd(1:Nfft/2+1) ; 
        psd = smoothdata(psd,'movmedian',5); 
        fVec = Fs*(0:(Nfft/2))/Nfft; 
        simulated_lfp = mean(lfp_task_range) + randn(size(lfp_task_range))*std(lfp_task_range) ;
        psd_simulated = fft(simulated_lfp,Nfft); 
        psd_simulated = psd_simulated(1:Nfft/2+1) ; 
        psd_simulated_new = psd_simulated.*(psd/psd(1)) ; 
        xf_filtered = [psd_simulated_new fliplr(psd_simulated_new)];
        simulated_LFP = real(ifft(xf_filtered)) ; 

        [ripples_encoding_sim] = detect_awake_ripple_haoxin(simulated_LFP, trials_time_axis, Fs);
        ripple_encoding_sim_ch(nch).ripples=ripples_encoding_sim; ripple_peaks_sim_ind = round(ripples_encoding_sim(:,2).*2000);
        
        
        Summary{subnum}.ripple_encoding_sim_ch=ripple_encoding_sim_ch;
        Summary{subnum}.ripple_encoding_ch=ripple_encoding_ch;

        
        toc
    end
    
%     
%     save(['R00_Output_RippleDetection_' sub],'ripple_encoding_ch','ripple_retrieval_ch','chan_ROI','MTL_roi',...
%         'trials_time_axis_encoding','trials_time_axis_retrieval')
%     
end


rmpath  \\mycloudpr4100\ivan\Ivan_analysis\Ivan\human-20190516T170723Z-001\human\FMAToolbox\Helpers
rmpath \\mycloudpr4100\ivan\Ivan_analysis\Ivan\human-20190516T170723Z-001\human\FMAToolbox\Analyses



for nsub = 1:7
    SWRnum=[];SWRSIMnum=[];
    for nch = 1: length( Summary{nsub}.ripple_encoding_sim_ch   )
    SWRnum(nch) = length(Summary{nsub}.ripple_encoding_ch(nch).ripples)+length(Summary{nsub}.ripple_retrieval_ch(nch).ripples)  ; 
    SWRSIMnum(nch) = length(Summary{nsub}.ripple_encoding_sim_ch(nch).ripples) +length(Summary{nsub}.ripple_retrieval_sim_ch(nch).ripples)  ; 
    end
    [SWRnum_sub(nsub),nchselect ] = max( SWRnum  );
    SWRnum_SIM_sub(nsub) = SWRSIMnum(nchselect); 
end

figure, plot( [SWRnum_sub; SWRnum_SIM_sub],'LineWidth',2)
set(gca, 'LineWidth',2)
xlim([0.5 2.5])
legend({'S1' , 'S2' , 'S3' ,'S4', 'S5', 'S6','S7'})



