%Doppler Effect Project Part 4:
%Source from left to right in front of a binaural observer
%srclen is the length of the source trajectory in meters
%v is the velocity of the source
%obsx is the x-position along that axis
%obsy is the distance of the source trajectory away from the observer
%f0 is the starting frequency
%Mike Battaglia

function doppler4(srclen, v, obsx, obsy, f0)
fs = 8000;
bitdepth = 16;
c = 341.17;
headwidth = .17;

oxl = obsx - headwidth/2;
oxr = obsx + headwidth/2;
totaltime = srclen/v;
t = 1/fs:1/fs:totaltime;
srcx = v*t;
relxl = srcx - oxl;
relxr = srcx - oxr;
radl = sqrt(relxl.^2 +obsy^2);
radr = sqrt(relxr.^2 +obsy^2);
pressenvl = 1./radl;
pressenvr = 1./radr;
vradl = (relxl*v)./sqrt(relxl.^2 + obsy^2);
vradr = (relxr*v)./sqrt(relxr.^2 + obsy^2);
fl = f0*c./(c+vradl);
fr = f0*c./(c+vradr);

timeoffsetl = (oxl/srclen)*totaltime;
timeoffsetr = (oxr/srclen)*totaltime;

outputl = pressenvl.*sin(2.*pi.*fl.*(t-timeoffsetl));
outputr = pressenvr.*sin(2.*pi.*fr.*(t-timeoffsetr));
outputl = outputl/max(outputl);
outputr = outputr/max(outputr);
output = cat(2, outputl', outputr');
subplot(411); hold all; 
plot(t, vradl); xlabel('Time'); ylabel('Radial Velocity'); title('Radial Velocity');
plot(t, vradr);
hold off;
subplot(412); hold all;
plot(t, pressenvl); xlabel('Time'); ylabel('Pressure'); title('Sound Pressure Envelope');
plot(t, pressenvr);
hold off;
subplot(413); hold all;
plot(t, fl); xlabel('Time'); ylabel('Frequency'); title('Frequency vs. Time');
plot(t, fr);
hold off;
subplot(414); hold all;
plot(t, outputl); xlabel('Time'); ylabel('Amplitude'); title('Output Waveform');
plot(t, outputr);
hold off;

sound(output, fs, bitdepth);
wavwrite(0.99*output, fs, bitdepth, 'doppler4.wav');