function peaks = find_peaks( ts, thresh, interval )
% function peaks = find_peaks( ts, thresh, interval )
%
% Return a list of peak indices.
%
% Nigel Stepp <stepp@atistar.net>
%

% The first thing to do is to recall that the derivative
% of the data at each peak will be zero (or cross zero)
ts_diff = diff( ts );

% We need to find the zero crossings for peaks. These are precisely where
% ts_diff goes from being positive to negative.
zero_crossings = (ts_diff(1:end-1) > 0) & (ts_diff(2:end) <= 0);

% We have to add one because ts_diff is missing a data point.
peaks = find(zero_crossings == 1) + 1;

% This is the end of the naive peak picker, which will do
% very well for simple sinusoids.  Only 3 lines long!

% The valleys occur when ts_diff goes from negative to positive
%zero_crossings = (ts_diff(1:end-1) < 0) & (ts_diff(2:end) >= 0);
%valleys = find(zero_crossings == 1) + 1;

% For a peak to count as a peak, it has to meet some additional
% criteria. Otherwise, any local maximum would be a peak.
%
% I will put two restrictions on peaks:

% The absolute peak height must be some percentage of the max
mn = min(ts);
ts = ts - mn;

mx = max(ts);
ts = ts/mx;

for i=1:length(peaks)
	if ts(peaks(i)) < thresh
		peaks(i) = 0;
	end;
end;
peaks = peaks( peaks ~= 0 );

% Within a given window, the peak must be the maximum point
for i=1:length(peaks)
	ind = peaks(i);
	ind_lo = ind-interval;
	ind_hi = ind+interval;
	
	if ind_lo < 1
		ind_lo = 1;
	end;
	if ind_hi > length(ts)
		ind_hi = length(ts);
	end;
	maximum = max( ts(ind_lo:ind_hi) );
	
	if ts(ind) < maximum
		peaks(i) = 0;
	end;
end;
peaks = peaks( peaks ~= 0 );