%Doppler Effect Project Part 1:
%Source "passing through" observer
%srclen is the length of the source trajectory in meters
%obsx is the x-position along that axis
%v is the velocity of the source
%f0 is the starting frequency
%Mike Battaglia

function doppler1(srclen, v, obsx, f0)
fs = 8000;
bitdepth = 16;
c = 341.17;
totaltime = srclen/v;
t = 1/fs:1/fs:totaltime;
srcx = v*t;
relx = srcx - obsx;
rad = abs(relx);
pressenv = 1./rad;
vrad = v*sign(relx);       % same as (reldisp*v)./sqrt(reldisp.^2 + 0^2), which is the actual formula, but this is easier to compute
f = f0*c./(c+vrad);

output = pressenv.*sin(2.*pi.*f.*t);

%clip the output signal
output = min(output, 1);
output = max(output, -1);

subplot(411); plot(t, vrad, 'LineWidth', 2); 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) -1 15]);
subplot(413); plot(t, f, 'LineWidth', 2); 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) -1 1]);

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