Comparative_MEA_dataset/Codes/MATLAB spike detection/detect_spikes.m

function [spikes,thr,thrmax, noise_std_detect, noise_std_sorted, 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 (only amp thresholding) version was taken by Meeri from
%some ?? Wave Clus PROGRAM Get_spikes ?? or so
%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 error window)        
               
x = double(x);

% Depends:
% toolbox             :  \signal\signal\ellip.m
% other               :  filtfilt
% other               :  fix_filter
% other               :  int_spikes
%
% 090814 korjattu ohi-indeksointi lis??m?ll? tarkistukset spike storing
% osioon

sr=handles.par.sr;
w_pre=handles.par.w_pre;
w_post=handles.par.w_post;
ref=handles.par.ref;
detect = handles.par.detection;
stdmin = handles.par.stdmin;
stdmax = handles.par.stdmax;
fmin_detect = handles.par.detect_fmin;
fmax_detect = handles.par.detect_fmax;
fmin_sort = handles.par.sort_fmin;
fmax_sort = handles.par.sort_fmax;
index1 = [];
spikes = [];
thr = [];
thrmax = [];
noise_std_detect = [];
noise_std_sorted = [];

% HIGH-PASS FILTER OF THE DATA

xf=zeros(length(x),1);
if exist('ellip')                         %Checks for the signal processing toolbox
    [b,a]=ellip(2,0.1,40,[fmin_detect fmax_detect]*2/sr);
    xf_detect=filtfilt(b,a,x);
    [b,a]=ellip(2,0.1,40,[fmin_sort fmax_sort]*2/sr);
    xf=filtfilt(b,a,x);
else
    xf=fix_filter(x);                   %Does a bandpass filtering between [300 3000] without the toolbox.
    xf_detect = xf;
end
lx=length(xf);

%clear x; MOVED TO THE END 

noise_std_detect = median(abs(xf_detect))/0.6745;
noise_std_sorted = median(abs(xf))/0.6745;
thr = stdmin * noise_std_detect;        %thr for detection is based on detected settings.
thrmax = stdmax * noise_std_sorted;     %thrmax for artifact removal is based on sorted settings.

% LOCATE SPIKE TIMES
switch detect
    case 'pos'
        nspk = 0;
        xaux = find(xf_detect(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_detect(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;
        % palauttaa ei-nollat ja lis?? niihin w_pre+1
        % absoluuttisista arvoista
        % jotka xf_detectiss? v?lill? w_pre+2:end-w_post-2
        xaux = find(abs(xf_detect(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; % heitt?? virheilmon: index sis?lt?? spiketimet, maxi ja iaux tyhji?
                xaux0 = index1(nspk);
            end
        end
end
%some values for SWTTEO
in.M = xf_detect;
in.SaRa = sr;
%Parameters
params.method = 'numspikes';
params.numspikes = numel(index1);
params.filter = 0;
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)
% tehd??n vain jos spikej? indexiss?!
if ~isempty(index)
    ls=w_pre+w_post;
    spikes=zeros(nspk,ls+4);
    xf=[xf zeros(1,w_post)];              % xf ja zeros ei yhteensopivia
    parfor i=1:nspk                          %Eliminates artifacts
        if max(abs( xf(index(i)-w_pre:index(i)+w_post) )) < thrmax
            % jos piikin indeksi datan vika tai tokavika piste niin aiheuttaa
            % ohi indeksoinnin, eik? piikkimuotokaan ole kokonainen
            % eli eliminoidaan artifactina
            if ~(index(i)+w_post+2 > length(xf))
                % lis?t??n nollaborderia, koska piikki aivan datan alussa
                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));
                    % aiheuttaa ongelman jos index(i)-w_pre-1 < 1
                    % tai
                    % jos index(i)+w_post+2 > length(xf)
                    % eli jos piikki on datassa alle w_pren p??ss? alusta
                    % tai aivan datan lopussa
                    % --> t?ll?in ei saada kokonaista piikkimuotoa, joten
                end
            end
        end
    end
    aux = find(spikes(:,w_pre)==0);       %erases indexes that were artifacts
    spikes(aux,:)=[];
    index(aux)=[];
    
    clear x; %moved from line 57

end
        
switch handles.par.interpolation
    case 'n'
        %eliminates borders that were introduced for interpolation 
        spikes(:,end-1:end)=[];       % laittaa kaikille vikan ja tokavikan []
        spikes(:,1:2)=[];             % laittaa kaikille ekan ja tokan []
    case 'y'
        %Does interpolation
        spikes = int_spikes(spikes,handles);   
end