%MMI 361 Project 1
%Used for Parts 1, 2, and 3

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

for i = 1:length(R)
    dist(i) = sqrt( R(i)*R(i) + d*d );
    for j = 1:num_tones
        if (R(i) == 0) & (d == 0)
            f_obs(j,i) = f0(j);                                   %observed instantaneous frequency
            loudness(i) = 1;                                      %loudness of received signal
        else
            vso(i) = vs*(R(i)/dist(i));                           %velocity of obj, respect to medium
            f_obs(j,i) = f0(j) * ( v / ( v - vso(i) ) );          %observed instantaneous frequency
            loudness(i) = 1 / dist(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([1,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;
    A = loudness(LF_index);                                       %sets A to loudness
    freq = f_obs(:,LF_index);                                     %sets freq to observed frequency(s)
    for i = 1:num_tones
        x(index) = x(index) + A * sin(2*pi*freq(i)*t(index));
    end;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Tried useing this for single tone to build the sound  one period at a time
%The tone did not slide from one pitch to the next, so I had to use the for statment above
%even though it has artifacts
%    for k = index:period*fs+index   %adds observed frequency one period at a time
%        x(k) = A * sin(2*pi*freq*t2(k-index+1));
%        index2 = index2 + 1;
%        if index2+index > length(t)
%            break;
%        end;
%    end;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    index = index+1;
end;

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

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

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

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

soundsc( x, fs );

