Question 2

close all;
clear all;

% define sampling rate and set up a time vector
Fs = 44100;
TotalTime = 4;
t = -TotalTime/2+1/Fs:1/Fs:TotalTime/2;

c = 345; %speed of sound

offset = 10; % distance of observer from source trajectory

% set up a major triad
freq = 300;
freq2 = freq*5/4;
freq3 = freq*3/2;

vSource = 20; % velocity of the source

distance = vSource*t; % horizontal position of source, NOT w.r.t. observer

newDistance = sqrt(offset.^2 + distance.^2); % horizontal position w.r.t. observer

envelope = 1.0./newDistance; % envelope of the sound. No need to limit this time.

% calculate the velocity of the source w.r.t the observer
theta = asin(offset./newDistance);
newVSource = vSource*cos(theta);

% initialize vector
f = zeros(1,length(t));
f2 = f;
f3 = f;

% apply doppler effect
for i = 1:length(t)
    if t(i)<=0
        f(i) = freq*c/(c-newVSource(i));
        f2(i) = freq2*c/(c-newVSource(i));
        f3(i) = freq3*c/(c-newVSource(i));
    else
        f(i) = freq*c/(c+newVSource(i));
        f2(i) = freq2*c/(c+newVSource(i));
        f3(i) = freq3*c/(c+newVSource(i));
    end
end

% turn into final sound and normalize
x = envelope.*(cos(2*pi*f.*t) + cos(2*pi*f2.*t) + cos(2*pi*f3.*t));

h = readhrtf(0,15,'H');
h = h';
M = length(h);

azimuth = theta;
azimuth = azimuth*180/pi;
azimuth = round(azimuth/5);
azimuth = 90 - azimuth*5;


x = x/max(abs(x));
x = x';
x = [x x];

N = 5*128;
iterations = ceil(length(x)/N);

y = zeros(length(x),2);

for r=0:(iterations-2)/2
    h = readhrtf(0,azimuth(N*r+1),'H')';
    y((r*N+1):(r*N+1 + N+M-2),1) = y((r*N+1):(r*N+1 + N+M-2),1) + conv(x((r*N+1):(r*N+N),1),h(:,2));
    y((r*N+1):(r*N+1 +N+M-2),2) =  y((r*N+1):(r*N+1 +N+M-2),2) + conv(x((r*N+1):(r*N+N),2),h(:,1));
end

for r = (iterations-2)/2+1:(iterations-2)
    h = readhrtf(0,azimuth(N*r+1),'H')';
    y((r*N+1):(r*N+1 + N+M-2),1) = y((r*N+1):(r*N+1 + N+M-2),1) + conv(x((r*N+1):(r*N+N),1),h(:,1));
    y((r*N+1):(r*N+1 +N+M-2),2) = y((r*N+1):(r*N+1 +N+M-2),2) + conv(x((r*N+1):(r*N+N),2),h(:,2));
end

y = 0.97*y/max(max(abs(y)));

soundsc(y,Fs);
wavwrite(y,Fs,'Prob2');