PeakPosFFT (spectrum, axis, max_pos)

%PeakPosFFT
%
%This function determines the position of a peak. The maximum of the
%discrete signal is determined, whereafter the signal is converted into the
%Fourier domain and the exact peak position is determined.
%
%Syntax:
%   [pos_peak, window_max, window_as, pos_intensity, IsAccurate] = PeakPosFFT (spectrum, axis, max_pos)
%
%Input parameters:
%   spectrum: double vector: the spectrum
%   axis: double vector: the corresponding x-axis
%   max_pos: double: position of the local maximum in the spectrum (has to
%      be determined before PeakPosFFT is called)
%
%Output parameters:
%   pos_peak: double: exact position of the peak
%   window_max: double vector: intensity values of the window around the
%       peak maximum
%   window_as: double vector: axis of the window around the peak maximum
%   pos_intensity: double: intensity of the peak at the maximum
%   IsAccurate: boolean: indicates whether the calculated peak position is
%       not more than one channel from the given spectral channel
%
%Example:
%   test = GSSpcRead ('c:\test.spc');
%   m = max (test.data);
%   p = find (test.data==m);
%   position_of_max = PeakPosFFT (test.data, test.xaxis, p(1))

%This software package is dual licensed. You can use it according to the terms
%of either the GPLv3 or the BSD license.
%
%C2004 Didier Hutsebaut, Raman Spectroscopy Research Group, department
%of Analytical Chemistry, Ghent University
%Small update by Kris De Gussem, C2007, Raman Spectroscopy Research Group,
%department of Analytical Chemistry, Ghent University
%
%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
%
%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_peak, window_max, window_as, pos_intensity,IsAccurate] = PeakPosFFT (spectrum, axis, max_pos)
tmpspec = spectrum(max_pos-1:max_pos+1);
m=max(tmpspec);
i=find(tmpspec == m);
max_pos = max_pos+i-2; %just alter the maxpos if one of the pixels next to this pixel is higher in intensity
clear i tmpspec m

if (max_pos - 3) < 1
    max_pos = 4;
    warning ('Biodata:msg', 'Position of peakmax is altered. Given peak position will be inaccurate');
    disp (['max_pos: ' str(max_pos)]);
    disp(max_pos-3:1:max_pos+4);
elseif (max_pos >length (axis)-3)
    max_pos = length (axis) - 3;
    warning ('Biodata:msg', 'Position of peakmax is altered. Given peak position will be inaccurate');
    disp (['max_pos: ' str(max_pos)]);
    disp(max_pos-3:1:max_pos+4);
end

window_max = spectrum(1,max_pos-3:1:max_pos+4);
window_as = axis(1,max_pos-3:1:max_pos+4);
fft_window = interpft(window_max,64*8);
fft_as = interpft(window_as,64*8);

% pos_max = find(fft_window == max(fft_window));
%    --> this might give troubles if we just use this command, if there is
%    no peak: the command will return a result and by taking 3 pixels to
%    the left and 4 to the right, the FFT_spectrum is heavily distorted and
%    maxima of other peaks may be present
%    ==> therefore, we use the interval [2 pixels to the left, 2 pixels to
%    the right]: in this way we can be quite sure that the maximum we
%    obtain is the maximum of the peak
lowlim=axis(max_pos-2);
highlim=axis(max_pos+2);
mypos = (fft_as >= lowlim) & (fft_as <= highlim);
% mypos2 = fft_as(mypos);
mywin2 = fft_window(mypos);
[m,pos_max] = max(mywin2); %#ok<ASGLU>
% we have the index of the trimmed axis, transform it into the index of the
% complete axis
tmp=1:length(fft_as);
tmp = tmp(mypos);
pos_max = tmp(pos_max); 


if ((pos_max+8 > 512) || (pos_max-7 < 1))
    pos_peak = max_pos;
    pos_intensity = spectrum(max_pos);
    warning ('Biodata:msg', ['Position of maximum could not be returned. Returning initial peak position (' num2str(max_pos)  ')']);
else
    window_max2 = fft_window(pos_max-7:pos_max+8); % 10
    window_as2 = fft_as(pos_max-7:pos_max+8);% -7,:1:8

    % POLYNOOM FITTEN DOOR WINDOW OM PIEKPOSITIE TE BEPALEN
    [P1,S1, MU] = polyfit(window_as2,window_max2,2); %#ok<ASGLU> % vroeger 2
    [dP1,dS1] = polyder(P1);
    options = optimset('Display','final','TolX',0.01,'MaxFunEvals',1000); % Niet te laag RS zelf maar precies tot op 1 cijfer na komma !!
    f = inline('polyval(dP1,window_as2,dS1,MU)','window_as2','dP1','dS1', 'MU');
    pos_peak = fzero(f,mean(window_as2),options,dP1,dS1,MU);
    pos_intensity = polyval(P1,pos_peak,S1,MU);
end

IsAccurate = (pos_peak >= axis(1,max_pos-1)) && (pos_peak <= axis(1,max_pos+1));