clear all

fs = 16000; % 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, 512);
l( 256 ) = 1;
r( 262 ) = 1; % ITD of 6 samples

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

% sinfrq = 1984.5; % 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 ) );

Mdist = 0.20; % m microphone distance
soundspeed = 340; % speed of sound (340 m/s)
ITDmax = Mdist / soundspeed; % maximal possible interaural time difference in seconds
I = 361; % number of coincidence detectors

fftr = fft( r ); % Fourier transform of the left microphone channel
fftl = fft( l ); % Fourier transform of the right microphone channel
M = length( l ); % number of points in the Fourier transform

% main loop

for m = 1 : ( M / 2 ) % loop through frequencies

  fm = ( m - 1 ) * fs / M; % frequency (in Hz) corresponding to frequency index 'm'

    for i = 1 : I % loop through coincidence detectors

      % compute the delay value for coincidence detector 'i'
      Taui = ( ITDmax / 2 ) * sin( ( i - 1 ) / ( I - 1 ) * pi - pi / 2 );

      % adjust phase of left signal
      Xl = fftl( m ) * exp( -j * 2 * pi * fm * Taui );

      % adjust phase of right signal
      Xr = fftr( m ) * exp( j * 2 * pi * fm * Taui );

      % magnitude of the difference of the left and right phase-delayed signals
      Delta( i, m ) = abs( Xl - Xr );

    end

end

% uncomment if you want to plot the 3D coincidence map Delta (azimuth
% vs. frequency vs. dissimilarity)
% NOTE: This could take a looooong time under Octave

% mesh( 90 - ( ( 1 : I ) - 1 ) / ( I - 1 ) * 180, ( 0 : ( M / 2 - 1 ) ) * fs / M, Delta' );

% Sum up Delta over frequency, which results in the dissimilarity
% function, which shows the similarity between the delayed signals at a
% specific coincidence detector (the higher the value, the less similar
% the signals).
D=sum(Delta');

% plot dissimilarity vs. coincidence detector index
plot(1:I,D);

% find the index 'i' at which the dissimilarity function has a minimum
% (this is the place of coincidence)
[ C, i ] = min( D );

% compute the azimuth corresponding to coincidence detector 'i'
alphai = 90 - ( i - 1 ) / ( I - 1 ) * 180
