function [ startles_beg, startles_end ] = msdetect( s, thresh, gapthresh, t, Fs )
% MSDETECT
% @param s signal
% @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)
% @param t optional time
% @param Fs optional sampling frequency (2*pi/mean dt by default)
% @return vectors startles_beg and startles_end containing the indices of the beginning and end of the startles

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

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

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

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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%Filter out the high frequency noise%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if exist('Fs')
    [b,a] = ellip(4, 0.1, 40, 4*2/Fs);
    sf = filter(b,a,s);
    sf_diff = diff(sf);
else
    sf_diff = diff(s);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%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);

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

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

% find zero crossings
overzero = find(sf_diff>0);
zerogaps = diff(overzero);
z_gaps = find(zerogaps>1);

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

% find the crossing closest to ibegin
new_begin = [];
for i=1:length(ibegin),
    temp=find(iup<ibegin(i));
    ltemp = length(temp);
    if ltemp>0
        choice=temp(ltemp);
        new_begin = [new_begin iup(choice)];
    end
end

n = length(new_begin);
% find the end of the startle as the maximum between 
% startle beginnings
new_end = zeros(1, n);
for i=1:n,
    startit = new_begin(i);
    if i<n
        endit = new_begin(i+1);
    else
        endit = length(s);
    end
    [val,loc] = max(s(startit:endit));
    new_end(i) = startit+loc;
end

startles_beg = new_begin;
startles_end = new_end;
