%MMI 361 Project 1
%used for part 4

clear all;
close all;
clc;

fs = 8000;                                                        %sampeling frequency
num_tones = input( 'Enter number of tones: ' );
for i = 1:num_tones
    f0(i) = input( 'Enter frequency f0: ' );                      %Hz, source frequency
end;
v = 344;                                                          %m/s, speed of sound in air
vs = input( 'Enter velocity of the object: ' );                   %m/s, speed of source

r1 = input( 'Enter distance of path in front of observer: ' );    %meters
r2 = input( 'Enter distance of path behind observer: ' );         %meters

r = horzcat(r1:-.5:0, -1:-.5:-r2);                                %path of obj in .5m intervals
d = input( 'Enter distance of observer from object path: ' );     %meters

L = r - 0.085;
R = r + 0.085;

for i = 1:length(R)
    distL(i) = sqrt( L(i)*L(i) + d*d );
    distR(i) = sqrt( R(i)*R(i) + d*d );
    for j = 1:num_tones
        if ( (L(i) == 0) | (R(i) == 0) ) & (d == 0)
            f_obsL(j,i) = f0(j);                                  %observed instantaneous frequency
            f_obsR(j,i) = f0(j);
            loudnessL(i) = 1;                                     %loudness of received signal
            loudnessR(i) = 1;
        else
            vsoL(i) = vs*(L(i)/distL(i));                         %velocity of obj, respect to medium
            vsoR(i) = vs*(R(i)/distR(i));                         %velocity of obj, respect to medium
            f_obsL(j,i) = f0(j) * ( v / ( v - vsoL(i) ) );        %observed instantaneous frequency
            f_obsR(j,i) = f0(j) * ( v / ( v - vsoR(i) ) );        %observed instantaneous frequency
            loudnessL(i) = 1 / distL(i);                          %loudness of received signal
            loudnessR(i) = 1 / distR(i);                          %loudness of received signal
        end;
    end;
end;

time = (length(R)/vs)/2;                                          %sec, lenght of sound
t = 0:1/fs:time - 1/fs;
x = zeros([2,length(t)]);

index = 1;
while index < length(t)
    LF_index = round(index/length(t)*length(R));                  %index for loudness & f_obs
    if LF_index == 0
        LF_index = 1;
    end;
    AL = loudnessL(LF_index);                                     %sets A to loudness
    AR = loudnessR(LF_index);                                     %sets A to loudness
    freqL = f_obsL(:,LF_index);                                   %sets freq to observed frequency
    freqR = f_obsR(:,LF_index);                                   %sets freq to observed frequency
    for i = 1:num_tones
        x(1,index) = x(1,index) + AL * sin(2*pi*freqL(i)*t(index));
        x(2,index) = x(2,index) + AR * sin(2*pi*freqR(i)*t(index));
    end;
    index = index+1;
end;

disp( 'Perceived frequency at the begining:' );
for i = 1:num_tones
    disp( ['   left: ', num2str(f_obsL(i,1)), ' Hz right: ', num2str(f_obsR(i,1))] );
end;
disp( 'Perceived frequency at the end:' );
for i = 1:num_tones
    disp( ['   left: ', num2str(f_obsL(i,end)), ' Hz right: ', num2str(f_obsR(i,end))] );
end;
disp( 'Perceived Frequency at the closest point:' );
for i = 1:num_tones
    disp( ['   left: ', num2str(f0(i)), ' Hz right: ', num2str(f0(i))] );
end;

x = .99*x/max(max(abs(x)));

figure(1);
subplot( 3,2,1 );
plot( vsoL(1:2:end) );
title( 'Velocity of source relative to Observer - Left' );
xlabel( 'Distance traveled - in meters' );
ylabel( 'Velocity' );
subplot( 3,2,2 );
plot( vsoR(1:2:end) );
title( 'Velocity of source relative to Observer - Right' );
xlabel( 'Distance traveled - in meters' );
ylabel( 'Velocity' );
subplot( 3,2,3 );
plot ( t,x(1,:));
title( 'Signal Envelope - Left' );
xlabel( 'Time - in Seconds' );
ylabel( 'Amplitude' );
subplot( 3,2,4 );
plot ( t,x(2,:));
title( 'Signal Envelope - Right' );
xlabel( 'Time - in Seconds' );
ylabel( 'Amplitude' );
subplot( 3,2,5);
hold;
for i = 1:num_tones
    plot( f_obsL(i,1:2:end) );
end;
title( ['Perceived Frequency - F_0 = ', int2str(f0), ' - Left'] );
xlabel( 'Distance traveled - in meters' );
ylabel( 'Frequency' );
subplot( 3,2,6);
hold;
for i = 1:num_tones
    plot( f_obsR(i,1:2:end) );
end;
title( ['Perceived Frequency - F_0 = ', int2str(f0), ' - Right'] );
xlabel( 'Distance traveled - in meters' );
ylabel( 'Frequency' );

x = x';
wavwrite(x, fs, 'doppler1.wav' );

soundsc( x, fs );

