function bst = burst_detect(spktimes, sr, isi_th, spk_num, bstd_th)

% Input
% -------------------------------------------------------------------------
% spktimes: spike times (array)
% sr: sampling rate
% isi_th: maximum inter-spike-interval to determine if two spikes belong
%       to the same burst
% spk_num: minimum spike number that a burst should contain for data
%       analysis consideration.
% bstd_th: minimum burst duration
%
% Output
% -------------------------------------------------------------------------
% bst: burst onset and ofset times (2-column matrix)
  
% --------------------------------------------------------------
% main function
%---------------------------------------------------------------

% initialization
bonset = spktimes(1);       % burst onset time

n = 1;                      % spike index
count = 0;                  % number of spikes 
bst = [];               % burst onset and offset times and spike number 

while n < length(spktimes)
    
    n = n + 1;
    count = count + 1;
    
    if spktimes(n,1) - spktimes(n-1,1) > isi_th*sr
        
        boffset = spktimes(n-1,1);
        bst = [bst;[bonset,boffset,count]];
        
        % initialize variables
        bonset = spktimes(n);
        count = 0;
        
    else
        
        if n == length(spktimes)
            
            boffset = spktimes(n,1);
            bst = [bst;[bonset,boffset,count]];
        
        end
        
    end
    
end

% delete bursts that contain very few spikes or with short duration 
tempi = bst(:,3)<spk_num | bst(:,2)-bst(:,1)<bstd_th;
bst(tempi,:) = [];
