freq_dep_3d_headattenuation(az_angle_rad, el_angle_rad,can_dat,fs,...

function [pspec,headdelay] = freq_dep_3d_headattenuation(az_angle_rad, el_angle_rad,can_dat,fs,...
                src,mic,aim,hrad,sos);
global microphone_number            
%
%function three_d_headattenuation

%  Computes the impulse response of a source signal with respect to 
%  both azimuth and elevation angles.  The az_angle_rad and el_angle_rad are 
%  scalars and are angles off the nominal facing direction of the talker. The procedure
%  is as follows:
%
%  1) Use the input azimuth and elevation angles to determine the four
%     frequency responses from which to interpolate; two in azimuth and
%     two in elevation
%  2) Convert the data to amplitude from dB
%  3) Interpolate in both directions to ascertain best logartithmically
%     sampled magnitude versus frequency data 
%  4) Interpolate for an appropriate frequency response at 20kHz
%  5) Determine spherical-head modelled time delay using 'head_delay'
%  6) Multiply real frequency response by e^jwt
%  7) Appropriately transform to the time domain to get impulse response
%     at the appropriate sampling rate

%  INPUTS:
%  az_angle_rad = azimuth angle of the microphone from zero radians
%                 which is at the direct front of the source. 
%  el_angle_rad = elevation angle in radians wrt zero directly
%                 in front of the mouth. 
%  can_dat      = data taken from Canadian source which is a function of
%                 azimuth, elevation and frequency.  
%  fs           = sampling rate
%  src          = (x,y,z) source location in shifted space
%  mic          = (x,y,z) mic location in shifted space
%  aim          = (x,y,z) of aiming point in shifted space
%  hrad         = radius of spherical head in meters
%  sos          = speed of sound
%  mic_no       = microphone number [0,447]
%  OUTPUT     --  the impulse response (1024 long) for the source in this
%                 direction


% Canadian measured data for the head attenuation is in 18 frequency steps
% and goes first by elevation angle from -270 to +90 (180 degrees total) in
% ten steps of screwy sizes.  Then it goes by azimuth angle
% Azimuth Datapoints in 25 15-degree steps from 0 to 360



% Conversion factors
n = 1024;               % DFT size
n2 = n / 2;
twopi = 2.0 * pi;
threepi = 3.0 * pi;
pi2 = pi/2.0;
pi32 = 3.0 * pi / 2.0;
dtr = pi / 180.0;
fres = fs / n;
fifteen_deg_in_rad = 15.0 * pi / 180.0;
nx = size(az_angle_rad);
e_angles = [-90*dtr -52*dtr -33*dtr -13*dtr 0 10*dtr 31*dtr 54*dtr 74*dtr 90*dtr];
freqs = [-50.0 160.0 200.0 250.0 315.0 400.0 500.0 630.0 800.0 1000.0 1250.0 ...
         1600.0 2000.0 2500.0 3150.0 4000.0 5000.0 6300.0 8000.0 12000.0];

% Find azimuth and elevation indices
% Note for elevation use negative angles for consistency

i_az = floor(az_angle_rad / fifteen_deg_in_rad) + 1;    % +1 for indexing

prop_angle = mod(el_angle_rad,twopi);
if prop_angle > pi2 && prop_angle < pi
    prop_angle = pi - prop_angle;
end
if prop_angle > pi && prop_angle < pi32
    prop_angle = pi - prop_angle;
end
if prop_angle > pi32
    prop_angle = prop_angle - twopi;
end

if  abs(prop_angle + pi2) < 0.00001 
    i_el = 1;
else
    k = 0;
    j = 1;
    while k == 0
        if prop_angle >= e_angles(j) && prop_angle < e_angles(j+1)
            i_el = j;
            k = 1;
        end
        j = j + 1;
    end
end

% Derive azimuthal factor

v11 = can_dat((i_az - 1)*10 + i_el,1:18);
v21 = can_dat((i_az - 1)*10 + (i_el+1), 1:18);
v12 = can_dat(i_az * 10 + i_el, 1:18);
v22 = can_dat(i_az * 10 + (i_el+1), 1:18);
az_bottom = az_angle_rad / fifteen_deg_in_rad - (i_az - 1);
az_top = i_az - az_angle_rad / fifteen_deg_in_rad;
el_range = e_angles(i_el+1) - e_angles(i_el);
el_bottom = (prop_angle - e_angles(i_el)) / el_range;
el_top = (e_angles(i_el+1) - prop_angle) / el_range;
vzbot = az_bottom * v12 + az_top * v11;
vztop = az_bottom * v22 + az_top * v21;
fout = el_bottom * vztop + el_top * vzbot;

spectrum = 10.0 .^(fout/20.0);
spectrum1 = [1.0 spectrum(1:18) 0.00001];
%plot(freqs,spectrum,'b');
%hold
%i_az = i_az
%i_el = i_el
xx2 = [fres:fres:fs/2.0];
%xx2size = size(xx2);
aspec = interp1(freqs,spectrum1,xx2,'spline');
aspec = max(aspec, 0.00001);
%plot(freqs,spectrum1,'bo',xx2,aspec,'r')
%title(['Magnitude of source spectrum for Microphone Number ' num2str(microphone_number)]);
%pause
%yy = 20.0 * log10(abs(aspec));
%plot(yy)
%pause

% Now get the delay around the head

headdelay = head_delay(hrad,src,mic,aim,sos);
headstart = round(headdelay*fs);

% Put the amplitude and the phase together

xx3 = complex(0.0,-twopi*xx2*headdelay);
phase = exp(xx3);

% put magnitude and phase together; reflect and take DFT to get impulse
% response

pspec = aspec .* phase;

% these three lines are commented out because from now on I'm returning the
% frequency response of the source transfer function upto the nyquist
% frequency

% a2 = conj(pspec);
% 
% 
% a3 = a2(n2:-1:2);
% 
% pspecn = [pspec a3];
% 
% y = real(ifft(pspecn,n));

return