function autodetect = detectallpeaks(x,a,levpeakmark,ifir,ilas)

maxdet = 10;
scal = 0.5;
an = a/max(a);    

Val = [];
%detect local minima in the 2nd derivative by x as peaks
mx = 4;
crit = crit1(an,x,mx);
Val = crit/max(crit);

%% remove/reduce noise
Val = Val-mean(Val);
Val = filt_SG(Val,x(2)-x(1),2,2);
% Val = smooth(Val,5,'sgolay',2)';

Val = max(Val,0);
Val = peak2(Val);

%% calculate number of returned peaks

Val = Val/max(Val);
data.Val = Val;

md = min(maxdet,length(find(Val>0)));
pts = (1:md)';
[srt,idall] = sort(Val,'descend');
srt = log(srt(pts));

% straight g with the parametric form ab(1)*x+ab(2)*y=1
% ab = [pts(1) srt(1);pts(end) scal]\[1;1];
ab = [pts(1) srt(1);pts(end) srt(end)]\[1;1];
% dd includes the distances of all points to g multiplied by norm(ab,2)
dd = zeros(1,md);
for i=1:md
    dd(i) = abs(ab(1)*i+ab(2)*srt(i)-1);
end

[~,id] = max(dd);

data.autodetectall = idall(1:length(find(Val>0)));
data.bestpeak = id;

id = max(1,min(round(scal*id),md));
tmp = idall(1:id);
autodetect = [];
for i=1:length(tmp)                                                        
    if max(1,ifir-5) > tmp(i) || tmp(i) > min(length(x),ilas+5)
        autodetect(end+1) = tmp(i);
    end
end
autodetect = sort(autodetect);

function crit = crit1(Dn,x,mx)
    crit = -vand2(Dn,x(2)-x(1),mx);

% return 1 peak per non zero interval
function Val = peak2(temp)
    n = length(temp);
    Val=zeros(size(temp));
    i0 = 1;
    while i0 < n
        while i0<n && temp(i0) == 0
            i0 = i0+1;
        end
        j = i0;
        while j<n && temp(j+1)>0
            j = j+1;
        end
        [~,id] = max(temp(i0:j));
        Val(i0+id-1) = temp(i0+id-1);
        i0 = j+1;
    end

% approximate 2nd derivative
    % find polynomial ax²+bx+c to approximate D on [x-mx,x+mx]
    % a/2 approximates 2nd derivative in x
function deriv2 = vand2(D,h,mx)
    n = size(D,2);
    if ~exist('mx','var')
        mx = 9;
    end
    np = 2*mx+1;             
    M = zeros(np,3);
    
    for i=1:mx      % get M so that M*(a b c) = D(x-mx:x+mx)
        M(i,:) = [(mx+1-i)^2*h^2 -(mx+1-i)*h 1];
        M(np+1-i,:) = [(mx+1-i)^2*h^2 (mx+1-i)*h 1];
    end
    M(mx+1,:) = [0 0 1];
    
    pM = pinv(M);           %pseudoinverse  -> (a b c) = pM * D(x-mx:x+mx)
    sM = sparse(n-2*mx,n);
    for i=1:n-2*mx
        sM(i,i+[1:np]-1) = pM(1,:); %reduce pM to first column
    end
    deriv2 = D*sM';                 %calculate a

function filtD = filt_SG(D,h,mx,ord)
n = size(D,2);
np = 2*mx+1;

M = zeros(np,ord+1);
% for i=1:mx      % get M so that M*(a b c) = D(x-mx:x+mx)
%     for j = 1:ord+1
%         M(i,j) = ((mx+1-i)*h)^(ord+1-j);
%         M(np+1-i,j) = ((mx+1-i)*h)^(ord+1-j);
%     end
% end
% M(mx+1,ord+1) = 1;

for i=-mx:mx
       for j=0:ord
              M(i+mx+1, j+1) = (i)^j;
       end
end

% F = M' * M;
pM = M * (M' * M)^(-1) * M';
filtD = D;
for n = mx+2:n-mx-1
    filtD(n) = D(n-mx:n+mx)*pM(:,mx+1);
end
for i = 1:mx+1
    filtD(i) = D(1:np)*pM(:,i);
end
for i = 0:mx
    filtD(n-mx+i) = D(n-np+1:n)*pM(:,mx+1+i);
end

function m = ismaxbool(t)
m = zeros(size(t));
    for i = 2:length(t)-1
        if t(i)>t(i-1) && t(i)>t(i+1)
            m(i) = 1;
        end
    end