localmax2.m

function [xout,yout,imaxout] = localmax2(X,Y);
% LOCALMAX2 - find local maxima by parabola fitting
% 		[xm,ym]=LOCALMAX2(x,y) returns maxima for every triplet of points (along 1st dim). 
% 		xm(i,j)=NaN means that there is no local maximum of y(:,j)=F{x(:,j)} between x(i,j) and x(i+2,j)
%		x may be the same size as y or as y(:,1) 
%  added to find the three highest points of array

% August 2001, ashcherbina@ucsd.edu
% April 2004, Visbeck visbeck@ldeo.columbia.edu

% fix array sizes
if size(X,1)==size(Y,1) & size(X,2)==1
   X=repmat(X,1,size(Y,2));
end

% find maximum
[dum,imax]=max(Y);
xout=imax+nan;
yout=imax+nan;
imaxout=imax+nan;

% select valid indices for fit
iok=find(imax>1 & imax<(size(Y,1)));

% restrict to valid values
imaxout(iok)=imax(iok);
imaxn=imax(iok);

if length(iok)>0
  disp([' localmax2: found ',int2str(length(iok)),' valid values'])
else
 disp([' localmax2 found no valid values'])
 return
end

for i=1:length(imaxn)
 j=iok(i);
 x1(i)=X(imaxn(i)-1,j);
 x2(i)=X(imaxn(i),j);
 x3(i)=X(imaxn(i)+1,j);
 y1(i)=Y(imaxn(i)-1,j);
 y2(i)=Y(imaxn(i),j);
 y3(i)=Y(imaxn(i)+1,j);
end




% compute fit to parabola 

a=(x3.*y2+x1.*y3-x1.*y2-y3.*x2-y1.*x3+y1.*x2)./...
   (x3.*x2.^2-x1.*x2.^2+x1.*x3.^2-x1.^2.*x3+x1.^2.*x2-x3.^2.*x2);
b=-(-x2.^2.*y3+x2.^2.*y1-y2.*x1.^2+y3.*x1.^2-x3.^2.*y1+y2.*x3.^2)./...
   ((-x3+x2).*(x2.*x3-x2.*x1+x1.^2-x3.*x1) );
c=(x2.^2.*y1.*x3-x2.^2.*x1.*y3-x3.^2.*y1.*x2+...
   y3.*x1.^2.*x2+x3.^2.*x1.*y2-y2.*x1.^2.*x3)./...
   ((-x3+x2).*(x2.*x3-x2.*x1+x1.^2-x3.*x1));

% check for maximum in range
xn=x1+nan;
ii=find(a<0);
xn(ii)=-b(ii)./a(ii)/2;

% compute max taget strength if desired
xout(iok)=xn;
if nargout>1
   yn=xn.^2.*a+xn.*b+c;
   yout(iok)=yn;
end