%The latest startle detection program
%Programmer: Jennifer Healey, Dec 13, 1999
%Last Modified: by me, now
%
% Usage: this is a script so you dont' have to pass a vector.  
% You can easily make it into a function if you prefer.
%
% The signals is "s"  the sampling frequency is "Fs"
%
% the results are s_freq, s_mags and s_dur
%

% Edit all here or comment out for use in a giant script
Fs=31;

%get the test data
cd /u/fenn/SignalTEST
load testgsr
cd /u/fenn/matlab/StartleDetection
s=testgsr;


%assign s

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

lgsr=length(s);
lgsr2=lgsr/2;
t=(1:lgsr)/Fs;

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

sf_prime=diff(sf);

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

%There is ringing in the signal so the first 35 points are excluded

l35=length(sf_prime);

sf_prime35=sf_prime(35:l35);

%Set a threshhold to define significant startle

thresh=0.01;

vector=sf_prime35;

figure(2)
plot(vector)

overthresh=find(vector>thresh);

%overthresh is the values at which the segment is over the threshold

overthresh35=overthresh+35;

%the true values of the segment

gaps=diff(overthresh35);

big_gaps=find(gaps>31)

%big_gaps returns the indices of gaps that exceed 31
%
% eg - big_gaps=[60 92 132 168....]
%
%  gaps(60)=245; gaps(92)=205 ...
%
% overthresh35(58:62)=  [346   347   348   593   594]
%
% so overthresh(61) is where the startle starts (ish)
%
% is overthresh (60) where the startle ends?



%check the results

iend=[];
ibegin=[];

for i=1:length(big_gaps)
  iend=[iend overthresh35(big_gaps(i))];
  ibegin=[ibegin overthresh35(big_gaps(i)+1)];
end;

%iend
%ibegin

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

% The idea being this is to find the zero crossing closet to where it goes 
% over threshold 

% find all zero crossings


overzero=find(sf_prime>0);
zerogaps=diff(overzero);
z_gaps=find(zerogaps>1);

iup=[];
idown=[];

for i=1:length(z_gaps)
  idown=[idown overzero(z_gaps(i))];
  iup=[iup overzero(z_gaps(i)+1)];
end;

% find up crossing closest to ibegin

new_begin=[];
for i=1:length(ibegin)
  temp=find(iup<ibegin(i));
  choice=temp(length(temp));
  new_begin(i)=iup(choice);
end;

% find down crossing closest to iend

new_end=[];
for i=1:length(ibegin)
  temp=find(idown>iend(i));
  choice=temp(1);
  new_end(i)=idown(choice);
end;

% plot the results

hold off
figure(2)
t=1/31:1/31:length(s)/31;
plot(t,s)
hold on
for i=1:(new_end/31)
  plot((new_end(i)/31),s(new_end(i)),'go')
  plot((new_begin(i)/31),s(new_begin(i)),'rx')
end;



s_mag=[];   %initialize a vector of startle magnitudes
s_dur=[];    %initialize a vector of startle durations

for i=1:(length(iend)-1)
  s_dur(i)=new_end(i+1)-new_begin(i);
  s_mag(i)=s(new_end(i+1))-s(new_begin(i));
end;

s_freq=length(ibegin);



cd /u/fenn/matlab/StartleDetection/Archive

hold off

