%burstanalysis.m
%Jessica Parker, October 10, 2024
%
%Function performs burst analysis, measures spike times, and burst properties such as BP, BD, IBI, and spike frequency.

function [bc tpks lowAmpOsc] = burstanalysis(V,tint,spth,ibith)

lnt = length(V);
tt = 0:tint:tint*(lnt-1);

sps = find(V>spth); %find indexes of voltages greater than spike threshold
lowAmpOsc = 0;

if length(sps)>1
    [Vpks pks] = findpeaks([-10.0; V(sps)]); %find indexes of sps corresponding to peaks of spikes
    pks = pks-1;
    tpks = [];
    
    remv = [];
    ampthrsh = 7.0; %Spikes whose amplitude is below this threshold are removed
    bb = 1;
    for aa=1:length(pks)
        if aa>1 && aa<length(pks)
            if Vpks(aa)-min(V(sps(pks(aa-1)):sps(pks(aa+1)))) < ampthrsh
                remv(bb) = aa;
                bb = bb+1;
            end
        elseif aa==1
            if Vpks(aa)-min(V(sps(pks(aa)):sps(pks(aa+1)))) < ampthrsh
                remv(bb) = aa;
                bb = bb+1;
            end
        else
            if Vpks(aa)-min(V(sps(pks(aa-1)):sps(pks(aa)))) < ampthrsh
                remv(bb) = aa;
                bb = bb+1;
            end
        end
    end
    lowAmpOsc = length(remv);
    
    pks(remv) = [];
    tpks = tt(sps(pks)); %find times of the peaks of spikes
    ipks = sps(pks); %find indexes of times of the peaks of spikes
    isi = tpks(2:end)-tpks(1:end-1); %interspike intervals (ISI)
    
    if std(isi)/mean(isi) > 0.05   %if CoV of ISI is bigger than 0.05, look
        bb = 1;                    %for bursting
        spkcnt = 0;
        tbb = [];
        
        if (tpks(1)>ibith) && (bb==1)
            tbb(1) = tpks(1);
            Vbpk(1) = V(ipks(1));
            spkcnt = 1;
            spf(1) = 0.0;
            bb = bb+1;
        end
        for aa=1:length(ipks)-1 %running through all ISIs
          if (tpks(aa+1)-tpks(aa)>ibith && min(V(ipks(aa):ipks(aa+1)))<-50.0) %if ISI is bigger than IBI threshold, and minV during ISI < -50
            tbb(bb) = tpks(aa+1); %then time of burst beggining is the time
            Vbpk(bb) = V(ipks(aa+1));
				%of the spike at the end of ISI
            if bb>1                   %end of previous burst duration is the
              teb(bb-1) = tpks(aa); %the time of the spike at the
              spb(bb-1) = aa+1-spkcnt; %Spikes per burst
              spf(bb-1) = spf(bb-1)/(spb(bb-1)-1); %Spike frequency, mean 1/ISI
            end                       %beginning of the ISI
            spf(bb) = 0;
            spkcnt = aa+1;
            bb = bb+1; %counts number of bursts
          else
            if bb>1
              spf(bb-1) = spf(bb-1) + 1.0/(tpks(aa+1)-tpks(aa));
            end
          end
        end
        
        if length(tbb)>1 && mean(spb) > 1.5
            bp = tbb(2:end)-tbb(1:end-1); %burst period
            tbb(end) = []; %take out last time of burst beginning due to no IBI
            Vbpk(end) = [];
            spf(end) = [];
            bd = teb-tbb; %burst duration
            ibi = bp-bd; %interburst interval
            
            bc = [tbb; bp; bd; ibi; spb; spf]; %will become bc1 or bc2
        else
            bc = [];
        end
    else %if CoV of ISI <= 0.05 but length of sps is at least 2, then spiking occurs
        bc = [];
    end
else %silence
    bc = [];
    tpks = [];
end

end