%{

Matt Montag
University of Miami
EEN502 Project 2

Experimenting with HRTF and Loudspeaker Arrays
Using the UC Davis CIPIC HRTF database, subject 17.

%}
clear all; clc;
c = 345;                        % speed of sound - m/s
fs = 44100;                      % sample rate - 1/s

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% PART A %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% prob 2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

f0 = 440;                       % source frequency - hz
A = 1;                          % source amplitude
W = 2*pi*f0;                    % omega - rad/s
T = 8;                          % total time - seconds
l = 100;                        % length of path - meters
Vs = l/T;                       % speed of source - m/s
y = 5;							% distance of source path - meters
t = -T/2:1/fs:T/2;		        % shifted array of time
x = t .* Vs;                     % array of source position
d = sqrt(x.^2 + y.^2);           % array of source distance - m
Va = Vs .* cos((atan2(y,x)));    % array of apparent velocity
d(d<1) = 1;                      % prevent divide by zero
env = 1./d;                      % array of amplitude envelope
Fa = f0.*(c./(c+Va));            % array of apparent frequency - hz
Wa = 2*pi*Fa;                    % array of apparent frequncy - rad/s
wav = env.*square(Wa.* t);         % array of waveform samples
wav = wav + env.*square(Wa*0.7.*t);
wav = wav + env.*0.3.*randn(1,length(wav));

signal = .95 * wav ./ max(abs(wav));

load hrir_final.mat;
hrir_index = -1;
segments = 100; % segmentation of HRTF convolution
% zero-pad signal so that it becomes divisible by # of segments
signal = [ signal zeros(1,segments-mod(length(signal),segments))]';
segment_dur = length(signal)/segments;

% allocate stereo signal for convolution sum length = A + B - 1
output = zeros(fs*T + length(hrir_l) - 1,2);

for segment=0:segments-2
	% We need to get the azimuth as a function of time.
	% So we take inverse tan of (y/x) and the way we obtain x is
	% by sampling from the x vector a number of times equal to Segments.
	azim = 90-rad2deg(atan2(y, x(1+floor((segment/segments)*length(x)))));
	%clc;
	% convert azimuth-first to elevation-first coordinate
	if azim > 270
		azim = azim - 360;
		elev = 0;
	elseif azim > 90
		azim = 180 - azim;
		elev = 180;
	else 
		azim = azim;
		elev = 0;
	end;
	
	% check if this position is closer to a different HRIR sample point
	if hrir_index ~= getNearestUCDindex(azim,elev);
		hrir_index = getNearestUCDindex(azim,elev);
		hl = getNearestUCDpulse(azim,elev,hrir_l);
		hr = getNearestUCDpulse(azim,elev,hrir_r);
	end;

	sampstart = segment * length(signal)/segments + 1;
	sampend = (segment+1) * length(signal)/segments;
	
	output(sampstart:sampend+length(hrir_l)-1,1) = ...
		output(sampstart:sampend+length(hrir_l)-1,1) + ...
		conv(hl,signal(sampstart:sampend));
	output(sampstart:sampend+length(hrir_r)-1,2) = ...
		output(sampstart:sampend+length(hrir_r)-1,2) + ...
		conv(hr,signal(sampstart:sampend));
	
end;

plot(output);
output = .9.*output/max(max(abs(output)));
wavwrite(output,fs,16,'partA_prob2.wav');
soundsc(output,fs);