%Doppler Effect Project Part 2:
%Source from left to right in front of 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 doppler2(srclen, v, obsx, obsy, f0)
fs = 8000;
bitdepth = 16;
c = 341.17;
totaltime = srclen/v;
t = 1/fs:1/fs:totaltime;
srcx = v*t;
relx = srcx - obsx;
rad = sqrt(relx.^2 +obsy^2);
pressenv = 1./rad;
vrad = (relx*v)./sqrt(relx.^2 + obsy^2);
f = f0*c./(c+vrad);
timeoffset = (obsx/srclen)*totaltime;
output = pressenv.*sin(2.*pi*(f).*(t-timeoffset));
output = output/max(output);
subplot(411); plot(t, vrad, 'LineWidth', 1); xlabel('Time'); ylabel('Radial Velocity'); title('Radial Velocity'); axis([min(t) max(t) min(vrad) max(vrad)]);
subplot(412); plot(t, pressenv, 'LineWidth', 1); xlabel('Time'); ylabel('Pressure'); title('Sound Pressure Envelope'); axis([min(t) max(t) min(pressenv) max(pressenv)]);
subplot(413); plot(t, f, 'LineWidth', 1); xlabel('Time'); ylabel('Frequency'); title('Frequency vs. Time'); axis([min(t) max(t) min(f) max(f)]);
subplot(414); plot(t, output, 'LineWidth', 1); xlabel('Time'); ylabel('Amplitude'); title('Output Waveform'); axis([min(t) max(t) min(output) max(output)]);

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