GetPeakPos2 (yaxis, xaxis)

%GetPeakPos2
%
%This function looks automatically for peak maxima by using
%the first and second derivative of the signal.
%
%Syntax:
%   pos = GetPeakPos2 (yaxis, xaxis)
%
%With:
%   yaxis: the spectrum
%   xaxis: the corresponding xaxis
%
%Examples:
%   peaks = GetPeakPos2 (spectrum, xaxis)

%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 = GetPeakPos2 (yaxis, xaxis)
yaxis = savgol(yaxis,9, 2, 0);

%hier selectie van pieken op basis van afgeleiden
%signaal = yaxis
SD1 = savgol( yaxis, 9, 3, 1);
SD2 = -1 .* savgol( yaxis, 9, 3, 2); %we zoeken max, dus negatieve 2de afgeleides
SD3 = savgol( yaxis, 9, 3, 3);

%zoek posities met tweede afgeleide > threshold (> 0): zoek berg in functie
pos2 = find (SD2 > 0);

%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
pos = intersect (pos1, pos2);

%doorsnede van pos1 en pos2: geschikte maxima
pos3 = find ((SD3(i).*SD3(i+1)) <= 0); %als product van 2 opeenvolgende getallen -1: dan wisseling teken

tmp = sortrows (abs(SD2(:)))';
tmp2 = ceil(length (SD2) * 0.8);
threshhold = tmp(tmp2);
pos2b = find (SD2 >= threshhold);
pos = intersect (pos, pos2b);


pos123 = [];
for i=1:length (pos)
    tf = pos(i);
    if ((isempty(find(pos3 == tf-2, 1))==false)||(isempty(find(pos3 == tf-1, 1))==false)||(isempty(find(pos3 == tf, 1))==false)||(isempty(find(pos3 == tf+1, 1))==false)||(isempty(find(pos3 == tf+2, 1))==false))
        pos123(length(pos123)+1) = tf; %#ok<AGROW>
    end
end


pos = pos123;

%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
%pos = xaxis (pos); %xas begint niet noodzakelijk op datapunt 1

%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