function latencies = findlatencies(stereo, doplot, uselocalmaximum, fs)
% Given a two-column array with the stimulus in the first channel
% and the response in the second, return the list of latencies
% between each stimulus and its response

% By Matt Wright, CCRMA, March 3, 2004

if nargin < 4
   fs = 44100;
end

if nargin < 3
   uselocalmaximum = 1;
end

if nargin < 2
   doplot = 0;
end

stimulus = stereo(:,1);
response = stereo(:,2);
len = length(response);



% Define noise floor in terms of the max. amplitude in first 10000 samples
% Since the response is synthesized this should be pretty safe
response_noisefloor = 10. * max(abs(response(1:10000)))


response_times = [];
i = 1;
while i <= len,
   % Look for the first sample above the noise floor
   while ((i <= len) && (abs(response(i)) < response_noisefloor)),
	 i = i + 1;
   end

   if (i > len) 
      break 
   end

   % Remember this sample index as the next response time
   response_times = [response_times i];

   % Skip the rest of this event by looking for 1000 consecutive 
   % samples below the noise floor again

   Nsilent = 0;
   while i <= len,
       if (abs(response(i)) < response_noisefloor)
	  Nsilent = Nsilent + 1;
	  if Nsilent >= 1000
	     break;
	  end
       else
	  Nsilent = 0;
       end
       i=i+1;
   end
end

if doplot
   response_times

   figure(1);
   subplot(211);
   plot(response);
   hold on;
   stem(response_times, ones(1, length(response_times)), 'r');
   title('responses with instants marked');

end


% Assume latency will always be within 100 ms
max_latency = fs * 0.1;

stimulus_times = [];
for response_index = response_times


   interval = [response_index - max_latency: response_index];
   possible_stimulus = stimulus(interval);


   if uselocalmaximum
       [amp, offset_index] = max(possible_stimulus);
   else
       %% Same technique as response: look for first magnitude above noise floor

       if possible_stimulus(1) > response_noisefloor
	  warning('Stimulus 100ms before response above noise floor');
       end

       for i = 1:max_latency
          if (abs(possible_stimulus(i)) > response_noisefloor)
	     offset_index = i;
	     break;
          end
       end

       if (i == max_latency)
           warning('entire 100ms before response below noise floor: zero');
	   offset_index = max_latency;
       end
   end

   true_index = offset_index + response_index - max_latency;


   stimulus_times = [stimulus_times true_index];
end


if doplot
   subplot(212);
   plot(stimulus);
   hold on;
   stem(stimulus_times, ones(1, length(stimulus_times)), 'r');
   title('stimuli with instants marked');

   stimulus_times
end


latencies = response_times - stimulus_times;
