clear all

fs = 44100; % sampling frequency (in Hz)

% First, generate the input signals.
% We need two signals, one for the right
% and one for the left channels.
% These are stored in the 'l' and 'r' arrays
%
% As an example, 3 signal types are provided:
% click, sine and noise (uncomment the commands
% for the signal you wish to use and comment out
% the others)
%
% In order to use signals from real microphones,
% the samples would need to somehow be recorded
% from the hardware audio device and then stored
% in the 'l' and 'r' arrays. The rest of the
% program would not need to be changed.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% click (all samples set to zero except one)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

l = r = zeros( 1, 1000);
l( 500 ) = 1;
r( 515 ) = 1; % ITD of 15 samples

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% a sinusoidal signal of frequency 'sinfrq'
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% sinfrq = 500; % sine frequency (in Hz)
% l = sin( 2 * pi * sinfrq * ( 1 : 1000 ) / fs );
% r = sin( 2 * pi * sinfrq * ( ( 1 : 1000 ) / fs + 300e-6 ) ); % ITD of 300 microseconds

%%%%%%%%%%%%%%
% random noise
%%%%%%%%%%%%%%

% l = rand( 1, 1000 );
% r( 1, 26 : 1000 ) = l( 1, 1 : ( 1000 -25 ) );

% initialize some variables needed for the correlation computation
Tint = 10e-3; % integration constant in seconds (10 ms)
Tint = Tint * fs; % convert to samples
Mdist = 0.20; % m microphone distance
soundspeed = 340; % speed of sound (340 m/s)
ITDmax = Mdist / soundspeed; % maximal possible interaural time difference in seconds
ITDmax = floor( ITDmax * fs ); % convert to samples
Psi = zeros( 1, 2 * ITDmax + 1 ); % create an array for the correlation function

% main loop

for n = ( 1 + ITDmax ) : ( 1000 - ITDmax ) % loop through the signal samples

  for Tau = -ITDmax : ITDmax % loop through all possible ITD values

    Psi( Tau + ITDmax + 1 ) = exp( -1 / Tint ) * Psi( Tau + ITDmax + 1 ) + l( n ) * r( n - Tau);

  end

end

plot( - ITDmax : ITDmax, Psi ); % plot ITD (in samples) vs. correlation
[ C, I ] = max( Psi ); % find the index 'I' at which the correlation function has a maximum
I = I - ITDmax - 1; % compute the corresponding ITD in samples
alpha = asin( 340 * I / ( Mdist * fs ) ) * 180 / pi % compute the azimuth (in degrees) corresponding to ITD 'I'
