%% Lesson 16:  Auditory Motion
%
% In this lesson, I'll introduce a simple method that converts a sound wave
% and a vector of x and y spatial position over time into a souund wave
% that simulates what an observer would hear while standing at the origin.
%
% The method introduces the cues of interaural time delay (ITD) a doppler
% shift, and inverse square law.  It also has a bit of an interaural level
% difference (ILD) because of the effect of the inverse square law and the
% fact that one ear will be closer to the sound than the other.
%
% The trick is to calculate what the time-course of a sound wave should be
% at the location of the ear given the sound wave at the source, s(t), and
% the vectors of spatial position x(t) and y(t).
%
% We'll start with a stationary band-pass sound.  We'll then simulate what
% it should sound like when translating on a horizontal pass from left to
% right in front of your head.  This first part should look familiar.

p.Fs = 44100;     % Sampling rate (Hz)
p.dur = 4;        % Duration (seconds)
centerFreq =440;  % Center of band-pass filter (Hz)
width = 200;      % Width of band-pass filter (Hz)

% Time vector
t= linspace(0,p.dur,p.dur*p.Fs)';

% White noise
y = randn(size(t));

% Fourier Transform:
F = complex2real(fft(y),t);

% Band-pass filter in the frequency domain:
Gaussian = exp(-(F.freq-centerFreq).^2/width^2)';
plot(F.freq,Gaussian);
F.amp = F.amp.*Gaussian;

% Inverse Fourier Transform:
s = real(ifft(real2complex(F)))';

% Play the sound
sound(s,p.Fs);

%% Auditory Horizontal Motion
%
% The next step is to set up vectors x and y that determine the location of
% the sound over time.  The units will be in meters.  This will be a sound
% traveling over the four seconds from 50 meters to the left to 50 meters
% to the right, 10 meters in front of you.  (The speed of the sound source
% will therefore be a steady 25 meters/sec or about 56 mph).  The graphics
% show a simple diagram of the position relative to your (yellow) head.

x = 50*linspace(-1,1,length(s))';
y = 10*ones(size(s));

figure(1)
clf
plot(x(1:200:end),y(1:200:end),'r:');
hold on
plot(0,0,'ko','MarkerSize',20,'MarkerFaceColor','y');
axis equal
xlabel('X (meters)');
ylabel('Y (meters)');

%%
% To calculate the sound at each of the two ears, we need to define the
% speed of sound and the width of the head:

p.c = 345;  %speed of sound (meters/second)
p.a = .0875;  %with of head (meters)

%%
% Next we'll calculate the distance (squared) from the source to each ear 
% for each point in time.  For illustration, we can plot these distances as
% a function of time.

dLsq = (x-p.a).^2+y.^2;
dRsq = (x+p.a).^2+y.^2;

figure(1)
clf
plot(t,sqrt(dLsq),'r-');
hold on
plot(t,sqrt(dRsq),'g-');
set(gca,'YLim',[0,sqrt(max(dLsq))*1.1]);
xlabel('Time (s)');
ylabel('Distance to each ear (m)');
legend({'left','right'},'Location','SouthEast');

%% Time delay to each ear
%
% The sound received at each ear depends on how long it takes for the sound
% to travel from the source to each ear, as a function of time. The delay 
% (tauL and tauR) is the distance to the ears divided by the speed of
% sound:

tauL = sqrt(dLsq)/p.c;  
tauR = sqrt(dRsq)/p.c;

clf
hold on
plot(t,tauL,'r-');
plot(t,tauR,'g-');

xlabel('Time (s)');
ylabel('Time delay to each ear (m)');
legend({'left','right'},'Location','SouthEast');


%% Interaural Time Difference (ITD)
%
% The Interaural Time Difference (ITD) is caused by the difference in
% distance of the object to the two ears.  The time difference between the
% two ears will be the difference in the left and right distances divided
% by the speed of sound.  Here's a plot of the ITD for our sound as a
% function of time:

ITD = 1000*(tauR-tauL);

clf
plot(t,ITD);
xlabel('Time (s)');
ylabel('ITD (msec)');

%% 
% The ITD asymptotes at +/- 5 msec because it takes 5 msec for a sound to
% pass the width of the head.  At the beginning and end of the sound, the
% position is pretty much on one side.

%% Inverse square law
%
% The intensity of a sound drops of inversely with the square of the
% distance.  We can simulate this by dividing our sound wave by our
% distance vectors.  The two vectors are then combined to produce a
% two-channel (stereo) sound.  There will be a very slight interaural level
% difference because one ear will always be closer to the sound than the
% other. But this is negligable for the distances here (as can seen by the
% overlap in the curves in the figure).

sL = s./dLsq;  
sR = s./dRsq;

yy = [sL,sR];
yy = yy/max(yy(:));

sound(yy,p.Fs);

%% The Doppler effect
%
% The Doppler effect is the change in pitch experienced by an observer as
% an object changes its speed relative to you.  For illustration, we can
% plot the speed based on the rate of changes of the distances:

dt = 1/(t(2)-t(1));
speedL = diff(sqrt(dLsq))*dt;
speedR = diff(sqrt(dRsq))*dt;
clf
hold on
plot(t(1:end-1),speedL,'r-');
plot(t(1:end-1),speedR,'g-');

xlabel('Time (s)');
ylabel('Speed relative to each ear (m/s)');
legend({'left','right'},'Location','SouthEast');
%%
% For the first half, the object is coming toward you at 25 m/sec(so the
% speed is negative).  This will increase the pitch of the sound because
% the sound waves will effectively bunch up in front of the source. In the
% second half, as the object retreats, the pitch will be lower than the
% actual sound because the sound waves will be relatively spaced out behind
% the object.  
%
%% Interpolating to simulate the sound at your head
%
% Putting this all together might seem like a calculus nightmare.
% Fortunately it's really easy.  All we need to do is re-sample the sound
% using linear interpolation based on our calculations of the time delay,
% after attenuating based on the inverse square law.  The interpolation
% step uses 'interp1' which basically re-evaluates our sound wave with a
% different vector of time.  
%
% zero out the auditory stimulus matrix
yy = zeros(length(t),2);

%% 
% Here's the interesting part: Interpolate back in time to determine sound
% for each ear and attenuate the sound based on inverse-square law
yy(:,1) = interp1(t,sL,t-tauR);
yy(:,2) = interp1(t,sR,t-tauL);

%clear out the NaN's (at the beginning of time where sound hasn't come to
%the head yet)
yy(isnan(yy))= 0;

yy = yy/max(yy(:));
sound(yy,p.Fs);

%% 'makeAuditoryMotion.m'
%
% I've provided a function 'makeAuditoryMotion' that performs the above
% calculations.  

yy = makeAuditoryMotion(p,s,x,y);
yy = yy/max(yy(:));
sound(yy,p.Fs);

%% A cheesy animation
%
% For fun, we can make a movie in the figure that moves our sound source
% over time while the sound is playing. 

n= 1;  % number of times to loop.
figure(1)
clf

hold on

% Draw the 'Head'
plot(0,0,'ko','MarkerSize',10,'MarkerFaceColor','y');
axis equal

% Draw the path 
plot(x(1:200:end),y(1:200:end),'r:');
h = plot(x(1),y(1),'ko','MarkerFaceColor','r','MarkerSize',5);
xlabel('X (meters)');
ylabel('Y (meters)');

%loop n times

for i=1:n
    %play the sound
    sound(yy,p.Fs);
    %animate the object
    tic
    while toc<p.dur
        %Find index that matches closest to the current time
        id = ceil(toc*p.Fs);
        set(h,'XData',x(id));
        set(h,'Ydata',y(id));
        drawnow
    end
end

%% animateAuditoryMotion
%
% I've also provided a function 'animateAuditoryMotion' that runs
% 'makeAuditoryMotion' and plays it along with our cheesy animation.  It's
% essentially the same code as above.

animateAuditoryMotion(p,s,x,y,2);


%% Recording a sound
%
% Matlab's function 'audiorecorder' allows you to record your own sounds
% with a built-in microphone or an external one.  It's easy to use.  Here
% we'll sample a little over 4 seconds of sound and clip it so it's exactly
% four seconds so that it's compatible with our existing x and y vectors.

r = audiorecorder(p.Fs, 16, 1);

record(r);     % speak into microphone...
pause(4.2)
stop(r);

s =  getaudiodata(r);
s = s(1:length(t));

sound(s,p.Fs);   % listen to complete recording

%%
% Now we'll make it move!

animateAuditoryMotion(p,s,x,y,2);

%% Other sounds, other paths
%
% Since this is all done with simulation and not with calculus, we can
% simulate any sound with any path.  Here are some examples:

load handel  
s = y; 
p.Fs = Fs;
p.dur = length(s)/p.Fs;  %duration of sound (second)

t= linspace(0,p.dur,length(s))';

%position of sound over time (in meters, with observer at [0,0])
x = 20*sin(2*pi*t/p.dur);
y = 5*cos(2*pi*t/p.dur)+10;

%translate sound wave and position into stereo sound wave
[yy,p] = animateAuditoryMotion(p,s,x,y,1);

%% 

p.dur = 4;
p.Fs = 44100;

t= linspace(0,p.dur,p.dur*p.Fs)';
s = floor(rand(size(t))+.001);

x = linspace(-.5,.5,length(t))';
y = .05*ones(size(t));

[yy,p] = animateAuditoryMotion(p,s,x,y,1);

%% Exercises
%
% 1) Combining sounds is easy.  Just add two sound waves together - but be
% sure that the sum stays between -1 and 1. Generate the sound of low
% frequency moving from left to right and high frequency moving from right
% to left.  Can you attentively separate the two sounds?  What if they are
% the same frequency?  
%
% 2) A sonic boom occurs when a noisy object moves at the speed of sound,
% where the sound waves all pile up in front of the object.  Can you
% simulate a sonic boom?  

% 3) What happens to a sound that is receeding from you faster than the
% speed of sound?
%
% 4)