GetPeakPos (yaxis, xaxis, treshhold)

%GetPeakPos
%
%This function calculates the peak maxima of a signal. The function looks
%for peaks based on the first and second derivative.
%
%Syntax:
%   pos = GetPeakPos (yaxis, xaxis, treshhold)
%
%With:
%   yaxis: the spectrum
%   xaxis: the corresponding xaxis
%   treshhold: a tyhreshold for the second derivative
%
%Examples:
%   peaks = GetPeakPos (spectrum, xaxis, 1e-4) 

%This software package is dual licensed. You can use it according to the terms
%of either the GPLv3 or the BSD license.
%
%The ShowPeaks toolbox for MATLAB: routines for peaks detection
%C 2004-2008, Kris De Gussem, Raman Spectroscopy Research Group, Department
%of analytical chemistry, Ghent University
%C 2009 Kris De Gussem
%
%This file is part of ShowPeaks. ShowPeaks is used a.o. by the Biodata toolb ox.
%
%ShowPeaks is free software: you can redistribute it and/or modify
%it under the terms of the GNU General Public License as published by
%the Free Software Foundation, either version 3 of the License, or
%(at your option) any later version.
%
%ShowPeaks is distributed in the hope that it will be useful,
%but WITHOUT ANY WARRANTY; without even the implied warranty of
%MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%GNU General Public License for more details.
%
%You should have received a copy of the GNU General Public License
%along with ShowPeaks.  If not, see <http://www.gnu.org/licenses/>.

%Copyright (c) 2004-2009, Kris De Gussem
%All rights reserved.
%
%Redistribution and use in source and binary forms, with or without 
%modification, are permitted provided that the following conditions are 
%met:
%
%    * Redistributions of source code must retain the above copyright 
%      notice, this list of conditions and the following disclaimer.
%    * Redistributions in binary form must reproduce the above copyright 
%      notice, this list of conditions and the following disclaimer in 
%      the documentation and/or other materials provided with the distribution
%    * Neither the name of Raman Spectroscopy Research Group, Department of
%	  analytical chemistry, Ghent University nor the names 
%      of its contributors may be used to endorse or promote products derived 
%      from this software without specific prior written permission.
%      
%THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" 
%AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE 
%IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE 
%ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE 
%LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR 
%CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF 
%SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS 
%INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN 
%CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) 
%ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE 
%POSSIBILITY OF SUCH DAMAGE.

function pos = GetPeakPos (yaxis, xaxis, treshhold)
%hier selectie van pieken op basis van afgeleiden
%signaal = yaxis
tmp = isnan(xaxis);
if isempty (tmp) == false;
    warning ('Calibration:msg', 'Possible error in calibration, xaxis contains NaN''s');
    yaxis (tmp) = [];
    xaxis (tmp) = [];
end
SD1 = savgol( yaxis, 13, 3, 1);
SD2 = -1 .* savgol( yaxis, 13, 3, 2); %we zoeken max, dus negatieve 2de afgeleides

%zoek posities met tweede afgeleide > threshold (> 0): zoek berg in functie
pos2 = find (SD2 >= treshhold);
if isempty (pos2)
    error ('Calibration:msg', 'Treshhold value for automatic peak detection is to high. Please check.');
end

%zoek posities met wisseling van teken van eerste afgeleide: eerste
%afgeleide = O, wil zeggen extremum
i =  1: length (SD1)-1;
pos1 = find ((SD1(i).*SD1(i+1)) <= 0); %als product van 2 opeenvolgende getallen -1: dan wisseling teken
%doorsnede van pos1 en pos2: geschikte maxima
pos = intersect (pos1, pos2);

%nog opletten: voor positie van maximum gelden twee opeenvolgende punten op
%x-as: maximum = wisseling van teken van eerste afgeleide
%   ==> kijk voor elk punt uit pos of het erop volgende punt groter is, dan
%       is dit het maximum dat we zoeken
for i = 1:length (pos);
    p = pos(i); % krijg de positie van het getal
    selspecpos = max(p-2, 1):min (p+2, length (yaxis));
    selspec = yaxis (selspecpos);
    v = max (selspec);
    tmp = find (selspec == v); % positie in geselecteerd stuk
    pos(i) = selspecpos(tmp(1));
end

%artefactjes opheffen van polynoombenadering: soms kleine piekjes als
%dubbel geteld: eerste afgeleide zeer klein en wisseling van teken 
i = 1;
while i < length (pos)
    if pos (i)+1 == pos(i+1);
        pos(i+1) = [];
        i = i-1;
    end
    i = i+1;
end
while pos (length (pos)) >= length (xaxis)-1;
    pos (length (pos)) = []; % anders problemen met interpolatie,...
end