localmax.m

function [y, x, n] = localmax(a)
% localmax - interpolated local maxima.
%
% [Y I] = localmax(X) returns indices I and values Y, each column of which
% corresponds to the set of local maxima in that column of the input X.
%
% The indices and maxima are determined by a quadratic fit to local maxima 
% in the columns of X; maxima occurring in the first and last rows of X are 
% not interpolated.
%
% Note that when the input X has more than one column, different columns of
% the input X will likely have different numbers of local maxima.  In this
% case, to indicate the absence of a local maxima an index entry of 0 and a
% value of NaN are used.  The call [Y I N] = localmin(X) also returns a row
% N containing the number of local maxima found in each column of X.

% (c) Copyright 1995 Abel Innovations.  All rights reserved.
%
% Jonathan Abel
% Created: 21-Aug-95, based on imax.m, version 1.0.
% Version: 1.0

%% verify input
%%
if (nargin < 1),
	error('Not enough input arguments.');
else,
	[r c] = size(a);
	if (r == 1),
		a = a(:);
		[r c] = size(a);
	end;
end;


%% find local maxima for each column
%%
ismax = [zeros(1,c); (a(2:r-1,:) > a(3:r,:)) & (a(2:r-1,:) > a(1:r-2,:)); zeros(1,c)];
n = sum(ismax);

x = zeros(max(n), c);
y = zeros(max(n), c);
for i = [1:c],
	index = find(ismax(:,i));

	A = (a(index+1,i) + a(index-1,i))/2 - a(index,i);
	B = (a(index+1,i) - a(index-1,i))/2;
	delta = min(ones(size(index)), max(-ones(size(index)), -0.5*B./A ));
	value = A.*(delta.^2) + B.*delta + a(index,i);

	x(1:n(i),i) = index + delta;
	y(1:n(i),i) = value;
end;