function spikeLocations = findSpikes(signal,handles)
% Author: Jakub Zamouøil
% Bakaláøská práce, FEL ÈVUT, 2015
% Bachelor work, FEE CTU, 2015

% SPIKE DETECTION ALGORITHM
% Version 1.0
% Inputs: filtered signal, handles with additional information
% Output: vector of indexes of samples where spikes were detected

% This program contains two almost identical functions.
% DetectPrecise uses three methods of detection while detectFast only two.

% Algorithm doesn't perform well for signals with very low AP rate

% HANDLES:
% It is recommended to set these variables:
% handles.spdet.sr = sampling frequency of your signal
% handles.spdet.fast = 0 or 1, here you choose which function to use
% Other fields are not used in this algorithm
if (~exist('handles', 'var'))
    handles = createDefaultHandles();
end
if (~isfield(handles.spdet, 'fast'))
    handles.spdet.fast = 0; % Default function is detectPrecise
end

if handles.spdet.fast
    spikeLocations = detectFast(signal,handles);
else
    spikeLocations = detectPrecise(signal,handles);
end

%% FASTER DETECTION
function spikeLocations = detectFast(signal,handles)

% Thresholds for detection
threshold = 2.47*iqr(signal); 
thresholdNEO = 3.05*iqr(signal); 

fs = handles.spdet.sr; % sampling frequency
delkaSpiku = ceil(0.0015*fs); % length of spike in samples
ctvrt = ceil((delkaSpiku)/4); % quarter of length

% padding with zeros, allocation of results vector
signal = [zeros(1,delkaSpiku), signal, zeros(1,delkaSpiku)];
spikeLocations = zeros(1,ceil(length(signal)/500));
spikeLocations(1) = 0;
index = 2;

for i = delkaSpiku+1:(length(signal)-delkaSpiku)
    
    % Is sample local maximum?
    if abs(signal(i-1)) > abs(signal(i))
        continue
    end  
    if abs(signal(i+1)) > abs(signal(i))
        continue
    end  
    
    % Is sample too close to last spike?
    if (i-delkaSpiku - spikeLocations(index-1)) < fs*0.0009
        continue
    end   
    
    % Actual spike detection
    if ((abs(signal(i))/threshold)...
            +(signal(i)^2 - signal(i-ctvrt)*signal(i+ctvrt))/thresholdNEO)...
            < 2
        continue
    end
    
    % Entry to results vector
    lokaceSpiku = i-delkaSpiku;
    spikeLocations(index) = lokaceSpiku;
    index = index+1;
end

%truncate ending zeros
spikeLocations = spikeLocations(1:find(spikeLocations,1,'last'));
end

%% MORE PRECISE DETECTION
function spikeLocations = detectPrecise(signal,handles)

% Thresholds for detection
threshold = 2.6*iqr(signal); 
thresholdNEO = 2.4*iqr(signal); 
vyskaSpiku = 6.6*iqr(signal); 

fs = handles.spdet.sr; % sampling frequency
delkaSpiku = ceil(0.0015*fs); % length of spike in samples
ctvrt = ceil((delkaSpiku)/4); % quarter of length

% padding with zeros, allocation of results vector
signal = [zeros(1,delkaSpiku), signal, zeros(1,delkaSpiku)];
spikeLocations = zeros(1,ceil(length(signal)/500));
spikeLocations(1) = 0;
index = 2;

for i = delkaSpiku+1:(length(signal)-delkaSpiku)

    % Is sample local maximum?
    if abs(signal(i-1)) > abs(signal(i))
        continue
    end  
    if abs(signal(i+1)) > abs(signal(i))
        continue
    end 
    
    % Is sample too close to last spike?
    if (i-delkaSpiku - spikeLocations(index-1)) < fs*0.0009;
        continue
    end    

    % Actual spike detection
    leftSlope = sum(diff(signal(i-ctvrt:i)));
    rightSlope = sum(diff(signal(i:i+ctvrt)));   
    if ((abs(signal(i))/threshold)...
            +((signal(i)^2 - signal(i-ctvrt)*signal(i+ctvrt))/thresholdNEO)...
            +(abs(leftSlope-rightSlope)/vyskaSpiku))...
            < 3
        continue
    end
    
    %Zápis do vektoru nalezených spikù
    lokaceSpiku = i-delkaSpiku;
    spikeLocations(index) = lokaceSpiku;
    index = index+1;
end

%truncate ending zeros
spikeLocations = spikeLocations(2:find(spikeLocations,1,'last'));
end
end