Comparative_MEA_dataset/Codes/MATLAB spike detection/amp_detect.m

function [spikes,thr,thrmax, noise_std_detect, index] = amp_detect(x,handles)
% Detect spikes with amplitude thresholding. Uses median estimation.
% Detection is done with filters set by fmin_detect and fmax_detect. Spikes
% are stored for sorting using fmin_sort and fmax_sort. This trick can
% eliminate noise in the detection but keeps the spikes shapes for sorting.

%The initial version involves only amp thresholding 
%Stationary Wavelet Transform Teager Energy Operator block is added by Andrey.
%Workflow: %1. Amplitude thresholding detects spikes in the prefiltered signal with TCF=4.5
           %2. The number of the detected spikes is fed into SWTTEO algorithm
           %3. The spike detected by both algorithms are considered as True Positives (within an error window)        
               
x = double(x);

% Depends:
% toolbox             :  \signal\signal\ellip.m
% other               :  filtfilt
% other               :  fix_filter
% other               :  int_spikes

sr=handles.sr;
w_pre=handles.w_pre;
w_post=handles.w_post;
ref=handles.ref;
detect = handles.detection;
stdmin = handles.stdmin;
stdmax = handles.stdmax;
fmin_detect = handles.detect_fmin;
fmax_detect = handles.detect_fmax;
index1 = [];
spikes = [];
thr = [];
thrmax = [];
noise_std_detect = [];


%% HIGH-PASS FILTER OF THE DATA
%Checks for the signal processing toolbox
[b,a]=ellip(2,0.1,40,[fmin_detect fmax_detect]*2/sr);
xf=filtfilt(b,a,x);

noise_std_detect = median(abs(xf))/0.6745;

thr = stdmin * noise_std_detect;        %thr for detection
thrmax = stdmax * noise_std_detect;     %thrmax for artifact removal 
%% LOCATE SPIKE TIMES
switch detect
    case 'pos'
        nspk = 0;
        xaux = find(xf(w_pre+2:end-w_post-2) > thr) +w_pre+1;
        xaux0 = 0;
        for i=1:length(xaux)
            if xaux(i) >= xaux0 + ref
                [maxi, iaux]=max((xf(xaux(i):xaux(i)+floor(ref/2)-1)));    %introduces alignment
                nspk = nspk + 1;
                index1(nspk) = iaux + xaux(i) -1;
                xaux0 = index1(nspk);
            end
        end
    case 'neg'
        nspk = 0;
        xaux = find(xf(w_pre+2:end-w_post-2) < -thr) +w_pre+1;
        xaux0 = 0;
        for i=1:length(xaux)
            if xaux(i) >= xaux0 + ref
                [maxi, iaux]=min((xf(xaux(i):xaux(i)+floor(ref/2)-1)));    %introduces alignment
                nspk = nspk + 1;
                index1(nspk) = iaux + xaux(i) -1;
                xaux0 = index1(nspk);
            end
        end
    case 'both'
        nspk = 0;
        xaux = find(abs(xf(w_pre+2:end-w_post-2)) > thr) +w_pre+1;
        xaux0 = 0;
        for i=1:length(xaux)
            if xaux(i) >= xaux0 + ref
                [maxi, iaux]=max(abs(xf(xaux(i):xaux(i)+floor(ref/2)-1)));    %introduces alignment
                nspk = nspk + 1;
                if isempty(iaux)
                    iaux = 0;
                end
                index1(nspk) = iaux + xaux(i) -1; 
                xaux0 = index1(nspk);
            end
        end
end
%% some values for SWTTEO
in.M = xf;
in.SaRa = sr;
%Parameters
params.method = 'numspikes'; %detection with predefined number of spikes
params.numspikes = numel(index1); %same number of spikes as detected with the amplitude thresholding
params.filter = 0; %no addition filtering for SWTTEO algorithm
%Wavelet decomposition level selection based on the sampling frequency
if sr==10000
    params.wavLevel=2;
elseif sr==12500
    params.wavLevel=3;
elseif sr==25000
    params.wavLevel=4;
elseif sr==50000
    params.wavLevel=5;
else
    error('Is it smth else than MCS or Axion? I know only sampling frequencies 10000, 12500, 25000 and 50000!')
end

%% Getting result for the SET NUMSPIKES from SWTTEO algorithm
index2= SWTTEO(in,params);

%% Selecting the spikes detected by both algorithms within some []error interval
errortime=ceil(2e-4*sr);%(in samples)

check=zeros(1, params.numspikes);
for i=1:params.numspikes
    check(i)=any(abs(index2(i)-index1)<=errortime);
end

index=index2.*check;
index=index(index>0);
index=sort(index);
 
nspk=numel(index);
%% SPIKE STORING (with or without interpolation) 
% [SPIKE STORING OUTPUT IS NOT IN USE CURRENTLY, ONLY ARTIFACT REMOVAL]  

if ~isempty(index)
    ls=w_pre+w_post;
    spikes=zeros(nspk,ls+4);
    xf=[xf zeros(1,w_post)];              
    parfor i=1:nspk                          %Eliminates artifacts
        if max(abs( xf(index(i)-w_pre:index(i)+w_post) )) < thrmax
            if ~(index(i)+w_post+2 > length(xf))
                if index(i)-w_pre-1 < 1 
                    spikes(i,:)=[0 xf(abs(index(i)-w_pre:index(i)+w_post+2))];
                else
                    spikes(i,:)=xf(abs(index(i)-w_pre-1:index(i)+w_post+2));
                end
            end
        end
    end
    aux = find(spikes(:,w_pre)==0);       %erases indexes that were artifacts
    spikes(aux,:)=[];
    index(aux)=[];
    
    clear x; 

end
        
switch handles.interpolation
    %spike storing is currently not in use
    case 'n'
        %eliminates borders that were introduced for interpolation 
        spikes(:,end-1:end)=[];       
        spikes(:,1:2)=[];             
    case 'y'
        %Does interpolation
        spikes = int_spikes(spikes,handles);   
end