function [spikes,thr,index] = Ex23_amp_detect(x)
% 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.

% code written by Rodrigo Quian Quiroga, 2004 
% WaveClus spike detection algorithm
% http://www.vis.caltech.edu/~rodri/Wave_clus/Wave_clus_home.htm (publicly available)
%
% adapted for Neuroscience course at FEL CVUT Prague by Eduard Bakstein,
% eduard.bakstein@fel.cvut.cz, 2013
%
%

par.w_pre=20;                       %number of pre-event data points stored
par.w_post=44;                      %number of post-event data points stored
par.detection = 'pos';              %type of threshold
par.stdmin = 3;                     %minimum threshold (number of standard devioations)
par.stdmax = 30;                    %maximum threshold (number of standard devioations)
par.detect_fmin = 300;              %high pass filter for detection (default 300)
par.detect_fmax = 3000;             %low pass filter for detection (default 3000)
par.interpolation = 'y';            %interpolation for alignment
par.sr=24000;                          % sampling freq
par.sort_fmin = 300;                %high pass filter for sorting (default 300)
par.sort_fmax = 3000;               %low pass filter for sorting (default 3000)
%min_ref_per = 1;                             %minimum refractory period (in ms)
%s.ref = floor(min_ref_per*s.sr/1000);                  %number of counts corresponding the minimum refractory period
par.ref = 36;
par.int_factor = 4;                 %interpolation factor




sr=par.sr;
w_pre=par.w_pre;
w_post=par.w_post;
ref=par.ref;
detect = par.detection;
stdmin = par.stdmin;
stdmax = par.stdmax;
fmin_detect = par.detect_fmin;
fmax_detect = par.detect_fmax;
fmin_sort = par.sort_fmin;
fmax_sort = par.sort_fmax;

% BAND-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);
% figure; plot(lx);
clear x;

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;
                index(nspk) = iaux + xaux(i) -1;
                xaux0 = index(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;
                index(nspk) = iaux + xaux(i) -1;
                xaux0 = index(nspk);
            end
        end
    case 'both'
        nspk = 0;
        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;
                index(nspk) = iaux + xaux(i) -1;
                xaux0 = index(nspk);
            end
        end
end

% SPIKE STORING (with or without interpolation)
ls=w_pre+w_post;
spikes=zeros(nspk,ls+4);
xf=[xf zeros(1,w_post)];
for i=1:nspk                          %Eliminates artifacts
    if max(abs( xf(index(i)-w_pre:index(i)+w_post) )) < thrmax               
        spikes(i,:)=xf(index(i)-w_pre-1:index(i)+w_post+2);
    end
end
% aux = find(spikes(:,w_pre)==0);       %erases indexes that were artifacts
% spikes(aux,:)=[];
% index(aux)=[];
        
switch par.interpolation
    case 'n'
        spikes(:,end-1:end)=[];       %eliminates borders that were introduced for interpolation 
        spikes(:,1:2)=[];
    case 'y'
        %Does interpolation
        spikes = int_spikes(spikes,par);   
end

    function [spikes1] = int_spikes(spikes,par) 
        %Interpolates with cubic splines to improve alignment.

        w_pre=par.w_pre;
        w_post=par.w_post;
        ls = w_pre + w_post;
        detect = par.detection;
        int_factor = par.int_factor;
        nspk=size(spikes,1);

        s=1:size(spikes,2);
        ints=1/int_factor:1/int_factor:size(spikes,2);

        intspikes=zeros(1,length(ints));
        spikes1=zeros(nspk,ls);
        switch detect
            case 'pos'
                for i=1:nspk
                    intspikes(:) = spline(s,spikes(i,:),ints);
                    [maxi iaux]=max((intspikes(w_pre*int_factor:w_pre*int_factor+8))); 
                    iaux = iaux + w_pre*int_factor -1;
                    spikes1(i,w_pre:-1:1) = intspikes(iaux:-int_factor:iaux-w_pre*int_factor+int_factor);
                    spikes1(i,w_pre+1:ls) = intspikes(iaux+int_factor:int_factor:iaux+w_post*int_factor);

                end
            case 'neg'
                for i=1:nspk
                    intspikes(:) = spline(s,spikes(i,:),ints);
                    [maxi iaux]=min((intspikes(w_pre*int_factor:w_pre*int_factor+8))); 
                    iaux = iaux + w_pre*int_factor -1;
                    spikes1(i,w_pre:-1:1) = intspikes(iaux:-int_factor:iaux-w_pre*int_factor+int_factor);
                    spikes1(i,w_pre+1:ls) = intspikes(iaux+int_factor:int_factor:iaux+w_post*int_factor);
                end
            case 'both'
                for i=1:nspk
                    intspikes(:) = spline(s,spikes(i,:),ints);
                    [maxi iaux]=max(abs(intspikes(w_pre*int_factor:w_pre*int_factor+8))); 
                    iaux = iaux + w_pre*int_factor -1;
                    spikes1(i,w_pre:-1:1) = intspikes(iaux:-int_factor:iaux-w_pre*int_factor+int_factor);
                    spikes1(i,w_pre+1:ls) = intspikes(iaux+int_factor:int_factor:iaux+w_post*int_factor);
                end
        end
    end
    
 end
