function latencies = findthrulatencies(stereo, doplot, filt, 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
%
% This version is for testing audio thru latency; it assumes that the
% stimulus and response are both band-limited impulse trains and therefore
% just identifies "event" times by looking for local amplitude maxima.

% By Matt Wright, CCRMA, March 3, 2004

if nargin < 4
   fs = 44100;
end
if nargin < 3
   filt = [1];
end
if nargin < 2
   doplot = 0;
end


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


for channel = 1:2
%    disp(['   channel ' num2str(channel)]);
    times = [];

    mono = stereo(:,channel);
    monofilt = conv(mono, filt);
    mono_mag = abs(monofilt);

    % One extra zero sample to simplify testing below
    mono_mag = [mono_mag; 0];


    totalmax = max(mono_mag);

    % The "noise" floor is just half the global maximum; this
    % is safe because we synthesized the signals in the first place.
    noisefloor = 0.5 * totalmax;

    i = 1;
    while i <= len,
           % Look for the first sample above the noise floor
          while ((i <= len) && (mono_mag(i) < noisefloor)),
	  	 i = i + 1;
          end

	  if (i > len) 
	      break 
	  end

	  % Now find the local maximum
	  while (mono_mag(i+1) >= mono_mag(i)),
		i = i + 1;
	  end

%	  i
%	  peak = mono_mag(i)
%	  below = mono_mag(i-1)
%	  above = mono_mag(i+1)
         
          j = i;

	  % Check if there were some equal-amplitude values in a row
	  while (mono_mag(i) == mono_mag(j-1)),
	     j = j-1
          end

	  max_location = (j+i)/2;

	  if (j == i)

	       logspace = 0;

	       if logspace
	       	     ym1 = log(mono_mag(i-1));
		     y0 =  log(mono_mag(i));
		     yp1 = log(mono_mag(i+1));

		     p = exp(qint(ym1, y0, yp1));
	       else
	       	     ym1 = mono_mag(i-1);
		     y0 =  mono_mag(i);
		     yp1 = mono_mag(i+1);

		     p = qint(ym1, y0, yp1);
	       end

	     max_location = max_location + p -1;
	  end

	  % Remember this sample index as the impulse instant's time
	  times = [times max_location];

	  % Skip the rest of this event
	  i = i + 3000;
   end

   if (channel == 1)
      stimulus_times = times;
      stim_filt = monofilt;
   else
      response_times = times;
      resp_filt = monofilt;
   end
end



if length(stimulus_times) ~= length(response_times)
   warning('Not the same number of stimuli and responses');
end

if doplot
%   response_times

   figure(100);
   subplot(211);
   plot(resp_filt);
   hold on;
   stem(response_times, max(resp_filt)*ones(1, length(response_times)), 'r');
   title('filtered responses with instants marked');


   subplot(212);
   plot(stim_filt);
   hold on;
   stem(stimulus_times, max(stim_filt)*ones(1, length(stimulus_times)), 'r');
   title('Filtered stimuli with instants marked');

%   stimulus_times
end


latencies = response_times - stimulus_times;


if ((max(latencies) - min(latencies)) > 1)
   warning(['latency spread from ' num2str(min(latencies)) ' to ' num2str(max(latencies))]);
end

samps = mean(latencies);
ms = samps * (1000/fs);

disp(['   Mean latency '  num2str(samps)  ' samples ('  num2str(ms)  ' ms)']);
