function [ startles_beg, startles_end ] = cleandetect( t, s, Fs, thresh, gapthresh )
% CLEANDETECT 
% @param t time
% @param s signal
% @param Fs optional sampling frequency (default to 2*pi/mean dt)
% @param thresh optional threshhold to define significant startle (default to 0.003)
% @param gapthresh optional threshhold to define significant gaps -- second-order diff (default to 31)
% @return vectors startles_beg and startles_end containing the indices of the beginning and end of the startles

% Startle Detection Script
% by Jennifer Healey - Dec 13 1999
% adapted by Nada Amin - Mar 6 2005

if ~exist('Fs')
  Fs = mean(diff(t));
end

if ~exist('thresh')
  thresh = 0.003;
end

if ~exist('gapthresh')
  gapthresh = 31;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%Filter out the high frequency noise%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

ls = length(s);
ls_half = ls/2;

[b,a] = ellip(4, 0.1, 40, 4*2/Fs);
sf = filter(b,a,s);

sf_diff = diff(sf);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%Find Significant Startles%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

overthresh = find(sf_diff>thresh);
gaps = diff(overthresh);
% big_gaps contain the indices of beginning of startles

big_gaps = find(gaps>gapthresh);

ngroups = length(big_gaps);

iend = overthresh(big_gaps);
ibegin = iend+1;

%%%%%%%%%%%%%%%
%%Fine Tuning%%
%%%%%%%%%%%%%%%

% find zero crossing close to threshholds
overzero = find(sf_diff>0);
zerogaps = diff(overzero);
z_gaps = find(zerogaps>1);

idown = overzero(z_gaps);
iup = idown+1;

% find crossing closest to ibegin
new_begin = zeros(1, ngroups);
for i=1:ngroups,
  temp=find(iup<ibegin(i));
  choice=temp(length(temp));
  new_begin(i) = iup(choice);
end

% find the end of the startle as the maximum between 
% startle beginning
new_end = zeros(1, ngroups-1);
for i=1:ngroups-1,
  startit = new_begin(i);
  endit = new_begin(i+1);
  [val,loc] = max(s(startit:endit));
  new_end(i) = startit+loc;
end

startles_beg = new_begin(2:ngroups);
startles_end = new_end;
