function [pos, mag, dsig] = peakfind(tsig, tx);
%  This function uses the change in gradient to determine local maxima in
%  sampled waveform SIG.
%
%    [pos, mag, grad] = peakfind(sig, tx)
%
% TX is the abscissa index for SIG, and peak locations are interpolated
% based on a weighted average between the points where the gradient changes.
% The output values POS are the positions of the local maximum (it
% will always include the endpoints of SIG) and MAG are the corresponding
% values of the peaks.  DSIG is the gradient of SIG.

[r,c] = size(tsig);  %  Get signal/vector size paramters
intv = (tx(2)-tx(1));   %  Get sampling interval (assuming a uniform interval)
fs = 1/intv;  %  Get sampling frequency

%  Find gradient of signal
dsig = gradient(tsig,intv);
%  Find transitions from positive to negative
% Initalize frist point as a relative maximum
count = 1;
mag(count) = tsig(1);
pos(count) = tx(1);
%  Search for all subsequent postive peak points
for k=2:length(dsig)
    %  If gradient trasnistions from positive to negative (or stays at
    %  zero, indicating a flat field)
    if(dsig(k) <= 0 & dsig(k-1) >= 0)
        count = count +1;   %  Index the new peak
        [dum, idum] = max(tsig((k-1):k));  %  Find peak height and position in original figure
        mag(count) = dum(1);  %  Store the peak magnituded
        if dsig(k) == dsig(k-1)  %  if the gradient for both points is equal
            pos(count) = tx(k-1)+intv/2;  %  Store the peak position at the midpoint
        else  %  If there is a clear peak transition
            %  Interpolate with weighted average to get the estimated peak between samples
                pos(count) = (tx(k)*(tsig(k))+(tsig(k-1))*tx(k-1))/((tsig(k))+(tsig(k-1)));
        end
          
        
    end
end

%  make the end point a potential peak point.
if pos(count) ~= tx(length(tx))
    pos(count+1) = tx(length(tx));
    mag(count+1) = tsig(length(tx));
end

% Check for flat fields, choose center of field as maximum.
k=1;
nk =0;

% Select middle of flat fields
while k < length(pos)
    %  Consecutive peaks are detected, use middle point
    if (mag(k) == mag(k+1))
        count = 1;
        %  Count out consecutive points
        while count+k < length(pos) & mag(k) == mag(k+count)
            count = count + 1;
        end
        nk=nk+1;
        nmag(nk) = mag(k);
        npos(nk) = pos(k + fix(count/2));
        k = k+count;
    else
        nk=nk+1;
        nmag(nk) = mag(k);
        npos(nk) = pos(k);
        k=k+1;
    end
end
nmag(nk+1) = mag(k);
npos(nk+1) = pos(k);
mag = nmag;
pos = npos;