function Result = emp_method(raster, params)

% Inputs:
% spi_trials: nxm cell: the timing of spikes as array in n trials
%             and m tries with 0 as the stim onset
% params: struct with following att.s:
%         t_start:   time of start of recording in ms
%         t_end:     time of end of recording in ms
%         bin_width: bin width in ms
%         slide:     sliding window length in ms

% Citation: Saleh F, Fakharian M & Ghazizadeh A "Stimulus presentation can enhance spiking irregularity across subcortical and cortical regions." 2021

cv_count = 100;
trial_count = size(raster,1);

% params
t_start = params.t_start;
t_end = params.t_end;
bin_width = params.bin_width;
slide = params.slide;

% finding times
Result.times = t_start+(bin_width/2):slide:t_end-(bin_width/2);

i_max = (bin_width/slide) - 1;

Result.VEC = zeros(size(Result.times));
Result.EVC = zeros(size(Result.times));
Result.nSI = zeros(size(Result.times));
Result.nRV = zeros(size(Result.times));

Result.VEC_SE = zeros(size(Result.times));
Result.EVC_SE = zeros(size(Result.times));
Result.nSI_SE = zeros(size(Result.times));
Result.nRV_SE = zeros(size(Result.times));

for i = 0:i_max
    
    perccount(i,i_max)
    times = t_start+(bin_width/2)+i*slide:bin_width:t_end-(bin_width/2);
    
    % compute Var and E
    E = zeros(trial_count,length(times));
    Var = zeros(trial_count,length(times));
    for trial_num=1:trial_count
        temp_raster = raster(trial_num,:)';
        spi_cell_org = cellfun(@(x) ...
            histcounts(x,(times(1)-bin_width/2):bin_width:(times(end)...
            + bin_width/2)),temp_raster,'UniformOutput',false);
        spi_mat_org = cell2mat(spi_cell_org);
        E(trial_num,:) = mean(spi_mat_org,1);
        Var(trial_num,:) = var(spi_mat_org,[],1);
    end
    
    VEC = zeros(size(times));
    EVC = zeros(size(times));
    nRV = zeros(size(times));
    nSI = zeros(size(times));
    
    VEC_SE = zeros(size(times));
    EVC_SE = zeros(size(times));
    nRV_SE = zeros(size(times));
    nSI_SE = zeros(size(times));
    
    parfor t = 1:length(times)
        expc_count = E(:,t);
        m = bootstrp(cv_count,@var,expc_count);
        VEC(t) = mean(m);
        VEC_SE(t) = std(m)/sqrt(cv_count);
        
        expc_var = Var(:,t);
        m = bootstrp(cv_count,@mean,expc_var);
        EVC(t) = mean(m);
        EVC_SE(t) = std(m)/sqrt(cv_count);
        
        m = bootstrp(cv_count,@(x,y) mean(x)/mean(y),expc_var,expc_count);
        nSI(t) = mean(m);
        nSI_SE(t) = std(m)/sqrt(cv_count);
        
        m = bootstrp(cv_count,@(x,y) var(x)/mean(y),...
            expc_count,expc_count)/(bin_width*1e-3);
        nRV(t) = mean(m);
        nRV_SE(t) = std(m)/sqrt(cv_count);
    end
    
    Result.VEC(ismember(Result.times, times)) = VEC;
    Result.nRV(ismember(Result.times, times)) = nRV;
    Result.EVC(ismember(Result.times, times)) = EVC;
    Result.nSI(ismember(Result.times, times)) = nSI;
    
    Result.VEC_SE(ismember(Result.times, times)) = VEC_SE;
    Result.nRV_SE(ismember(Result.times, times)) = nRV_SE;
    Result.EVC_SE(ismember(Result.times, times)) = EVC_SE;
    Result.nSI_SE(ismember(Result.times, times)) = nSI_SE;
end

end