function peaks = find_peaks( ts, thresh, interval )
% function peaks = find_peaks( ts, thresh, interval )
%
% Return a list of peak indices.
%
% Nigel Stepp
%
% 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 );