function Result = var_cv(raster, params)

% Inputs:
% raster: nx1 cell: the timing of spikes as array in n trials 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
%         keep_indx: mean matching index

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

raster = reshape(raster(:,1),size(raster,1),1);

n = length(raster);

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

if(~isfield(params,'min_spike'))
    params.min_spike = 0;
end

% find average spike count of each neuron
valid_index = ones(n,1);

for i = 1:n
    raster_local = raster{i};
    spike_mat = cell2mat(cellfun(@(x) ...
        histcounts(x,t_start:bin_width:t_end),...
        raster_local,'UniformOutput',false));
    if mean(spike_mat(:)) < params.min_spike
        valid_index(i) = 0;
    end
end

Result.valid_index = valid_index;
Result.keep_perc = sum(valid_index)/length(valid_index)*100;
raster = raster(valid_index == 1);
n = length(raster);

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

if(~isfield(params,'keep_indx'))
    params.keep_indx = repmat(1:n,length(Result.times),1)';
end

% method_params
params.cv2_params.bin_width = bin_width;
i_max = (bin_width/slide) - 1;

Result.K = zeros(size(Result.times));
Result.cv2 = zeros(size(Result.times));

Result.K_SE = zeros(size(Result.times));
Result.cv2_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);
    
    params.cv2_params.keepindx =...
        params.keep_indx(:,ismember(Result.times, times));

    cv2 = cv2_method(raster,times,params.cv2_params);
    Result.K(ismember(Result.times, times)) = cv2.shape_K;
    Result.cv2(ismember(Result.times, times)) = cv2.cv2_est;
  
    Result.cv2_all_neu(:,ismember(Result.times, times)) = cv2.cv2_est_tot; %%% new
    
    Result.K_SE(ismember(Result.times, times)) = cv2.shape_K_SE;
    Result.cv2_SE(ismember(Result.times, times)) = cv2.cv2_est_SE;
    
    Result.cvn_dist(1,ismember(Result.times, times),:) = cv2.cv_n_dist;
    
end