function [y, ir, HRTF, f, t, HRTF0, f0] = hrtf(theta, x, Fs, N)
    // Head-related transfer function
    //
    // Usage:
    //       [y, ir, HRTF, f, t HRTF0, f0] = hrtf(theta, x, Fs, N)
    // where
    //       theta: Azimuth angle (clock-wise)
    //       x:     Monophonic anechoic input signal. It may be
    //              a vector, in which case its sample rate is
    //              specified by Fs, or it may be a path to a
    //              wav file from which the sample rate is also
    //              obtained from the file. In any case the 
    //              sample rate is converted to 44100 Hz. In the
    //              latter case Fs is ignored
    //       Fs:    Sample rate in Hz. Default: 44100 
    //       N:     Size of the analysis window. Default: 4096
    //       y:     Binaural signal (row 1 contains left channel
    //              response to x; row2, right channel)
    //       t:     Time of binaural signal in s 
    //       ir:    Binaural impulse response of the HRTF (row 
    //              1 contains left channel impulse response;
    //              row 2, right channel)    
    //       HRTF:  Binaural frequency response of the HRTF, in
    //              dB (row  1 contains left channel frequency
    //              response, row 2, right channel). It has N/2
    //              amplitude values corresponding to equally 
    //              spaced frequencies between 0 and Fs
    //       f:     Frequency vector in Hz 
    //       HRTF0: Original Shaw HRTF data, presented as a
    //              matrix whose columns correspond to 
    //              frequencies and whose rows correspond to
    //              azimuth values
    //       f0:    Frequencies of original Shaw data 
    //
    // NOTE: The HRTF is interpolated between the curves given
    //       at 15° intervals 
    //       
    // References
    // Blauert, J. (1997). Spatial Hearing: The Psychophysics
    //    of Human Sound Localization (MIT Press, Cambridge, MA).
    // Shaw, E . A. G. (1974) Transformation of sound pressure 
    //    level from the free field to the eardrum in the
    //    horizontal plane. Journal of the Acoustical Society of
    //    America. 56, 1848-1861.
    // Shaw, E. A. G., and Vaillancourt, M. M. (1985). 
    //    Transformation of sound-pressure level from the free
    //    field to the eardrum presented in numerical form. 
    //    Journal of the Acoustical Society of America 
    //    78:1120-1122.
    //
    // ---------------------------
    // Author: Federico Miyara
    // Date:   2024-06-07
    
    // Missing data
    if argn(2)<4
        N = 4096;
    end
    if argn(2)<3
        Fs = 44100;
    end
    if argn(2)<2
        x = zeros(1,2*N);
        x(1) = 1;
    end
             
    // For testing purposes
    if 1==2
        N = 4096
        Fs = 44100
        theta = 32
        caso = 2
        select caso 
        case 1
            x = zeros(1,2*N);
            x(1) = 1;
        case 2
            [x, Fsx] = wavread("s:\curiosamente.wav");
            x = intdec(x, Fs/Fsx); 
        end
    end

    // If x is a path to a wav file, read the file
    // and convert sample rate to 44100
    if type(x)==10
        [x, Fsx] = wavread(x);
        x = intdec(x, 44100/Fsx);
    end
    
    // Frequencies in Hz
    f0 =  [ 0.2   0.25  0.3   0.32  0.4  0.5   0.6    0.63  0.7   0.8 ..
            0.9   1.0   1.2   1.25  1.4  1.6   1.8    2.0   2.3   2.5 ..
            2.7   2.9   3.0   3.2   3.5  4.0   4.5    5.0   5.5   6.0 ..
            6.3   6.5   7.0   7.5   8.0  8.5   9.0    9.5  10.0  10.5 ..
           11.0  11.5  12.0]*1000;

    // Azimuths in °
    theta0 =  [  0  15  30  45  60  75  90 105 120 135 150 165 ...
    180 195 210 225 240 255 270 285 300 315 330 345]';

    // HRTFs in dB
    H1 = [  0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0 
            0.4   0.4   0.5   0.5   0.8   1.1   1.4   1.4   1.4   1.3
            0.7   0.8   1.0   1.1   1.7   2.2   2.6   2.6   2.6   2.5
            1.0   1.2   1.4   1.5   2.4   3.2   3.6   3.6   3.5   3.5
            1.3   1.5   1.8   1.9   2.8   3.8   4.3   4.3   4.3   4.4
            1.5   1.7   2.1   2.2   3.2   4.3   4.8   4.8   4.8   4.9
            1.6   1.8   2.2   2.3   3.4   4.5   5.0   5.0   5.0   5.1
            1.6   1.8   2.1   2.2   3.3   4.3   4.8   4.9   4.9   5. 
            1.4   1.6   1.9   2.0   3.0   3.9   4.4   4.5   4.6   4.7
            1.1   1.3   1.5   1.6   2.4   3.1   3.5   3.6   3.8   4.1
            0.7   0.8   0.9   0.9   1.4   2.0   2.4   2.5   2.7   3.2
            0.2   0.1   0.0   0.0   0.3   0.7   1.1   1.2   1.4   1.9
           -0.2  -0.5  -0.7  -0.7  -0.7  -0.6  -0.4  -0.3  -0.1   0.4
           -0.5  -0.9  -1.2  -1.3  -1.6  -1.7  -1.6  -1.6  -1.4  -1.2
           -0.8  -1.3  -1.7  -1.8  -2.2  -2.5  -2.6  -2.6  -2.7  -2.8
           -1.1  -1.5  -1.8  -1.9  -2.2  -2.5  -3.   -3.1  -3.5  -3.9
           -1.4  -1.6  -1.7  -1.7  -1.9  -2.0  -2.5  -2.6  -2.9  -3.2
           -1.6  -1.6  -1.5  -1.5  -1.4  -1.4  -1.7  -1.8  -1.9  -2.0
           -1.6  -1.5  -1.4  -1.4  -1.2  -1.0  -1.2  -1.3  -1.5  -1.5
           -1.5  -1.4  -1.4  -1.4  -1.4  -1.4  -1.6  -1.7  -2.0  -2.3
           -1.3  -1.3  -1.4  -1.4  -1.8  -2.2  -2.5  -2.6  -2.9  -3.4
           -1.0  -1.1  -1.3  -1.4  -1.9  -2.4  -2.9  -3.0  -3.4  -4.1
           -0.7  -0.8  -1.0  -1.1  -1.5  -2.0  -2.4  -2.5  -2.8  -3.2
           -0.4  -0.4  -0.5  -0.5  -0.8  -1.1  -1.3  -1.3  -1.4  -1.6];

    H2 = [  0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0 
            1.4   1.4   1.5   1.5   1.4   1.3   1.2   1.2   1.3   1.5 
            2.6   2.7   2.9   2.8   2.7   2.5   2.3   2.3   2.4   2.8 
            3.7   3.8   4.0   3.9   3.7   3.4   3.1   3.0   3.2   3.6 
            4.5   4.7   4.9   4.9   4.7   4.2   3.7   3.3   3.5   4.0 
            5.1   5.3   5.7   5.7   5.5   4.6   3.6   2.9   3.0   3.5 
            5.3   5.6   6.1   6.1   5.7   4.5   3.1   2.0   1.7   2.1 
            5.2   5.5   6.1   6.    5.5   4.1   2.2   0.8  -0.2  -0.2 
            5.0   5.4   6.0   5.9   5.4   3.9   2.0   0.2  -1.3  -1.5 
            4.5   5.0   5.9   5.9   5.4   4.0   2.2   0.4  -1.6  -1.9 
            3.7   4.3   5.2   5.2   5.    3.8   2.2   0.6  -0.9  -1.5 
            2.6   3.3   4.1   4.1   3.9   2.8   1.2  -0.2  -1.2  -1.5 
            1.1   1.8   2.3   2.3   2.1   1.3  -0.2  -1.6  -2.5  -2.6 
           -0.8  -0.3   0.2   0.2   0.2  -0.4  -2.   -3.4  -4.3  -4.4 
           -2.9  -2.8  -2.3  -2.2  -2.0  -2.4  -3.8  -5.4  -6.5  -6.6 
           -4.3  -4.5  -4.5  -4.5  -4.5  -4.9  -6.2  -7.8  -9.0  -9.1 
           -3.4  -3.6  -3.9  -4.1  -5.2  -7.1  -8.9 -10.0 -10.5 -10.6
           -2.0  -2.0  -2.1  -2.4  -3.8  -5.8  -7.7  -9.0  -9.7 -10.0 
           -1.4  -1.2  -1.0  -1.1  -1.7  -3.1  -5.1  -6.9  -8.0  -8.3 
           -2.4  -2.5  -2.9  -3.1  -4.0  -5.9  -6.8  -7.0  -6.8  -7.0  
           -3.9  -4.4  -5.9  -6.3  -7.1  -7.3  -6.6  -6.0  -6.2  -6.6 
           -4.8  -5.4  -6.3  -6.3  -6.1  -5.3  -4.5  -4.3  -5.1  -5.8 
           -3.6  -3.9  -4.0  -3.9  -3.6  -3.1  -2.8  -2.7  -3.4  -3.9 
           -1.8  -1.9  -1.8  -1.8  -1.6  -1.4  -1.3  -1.3  -1.6  -1.9];

    H3 = [  0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0  
            1.8   2.0   2.1   2.2   2.2   2.1   2.1   2.4   2.6   2.6 
            3.1   3.5   3.7   3.9   3.7   3.5   3.8   4.5   5.    5.1 
            4.1   4.6   4.8   4.8   4.3   4.0   4.6   5.8   6.8   7.4 
            4.6   4.9   4.8   4.6   3.8   3.3   4.4   6.2   7.7   8.8 
            4.0   4.1   4.0   3.5   2.5   1.8   3.0   5.6   7.9   9.4 
            2.4   2.4   2.3   1.9   0.7  -0.3   0.8   4.0   7.0   9.3 
            0.0   0.1   0.1  -0.1  -1.2  -2.4  -1.7   1.7   4.8   7.7 
           -1.4  -1.2  -1.2  -1.4  -2.3  -4.3  -4.2  -1.0   2.0   4.8 
           -1.8  -1.8  -1.8  -1.9  -2.7  -5.2  -6.   -3.8  -0.8   1.7 
           -1.4  -1.3  -1.3  -1.4  -2.5  -4.9  -5.9  -4.9  -3.0  -1.0 
           -1.6  -1.5  -1.4  -1.5  -2.4  -4.3  -5.3  -4.9  -4.2  -3.3 
           -2.5  -2.5  -2.5  -2.6  -3.4  -5.0  -5.4  -4.5  -5.0  -5.1 
           -4.3  -4.2  -4.2  -4.4  -5.3  -6.7  -6.6  -5.8  -6.4  -6.7 
           -6.6  -6.7  -6.7  -6.9  -7.8  -9.1  -9.2  -8.7  -9.0  -9.4 
           -9.1  -9.1  -9.2  -9.5 -10.6 -12.3 -12.2 -11.4 -11.6 -11.9
          -10.6 -10.8 -10.9 -11.3 -12.5 -13.8 -13.6 -12.8 -12.3 -12.5
          -10.2 -10.6 -10.8 -11.3 -12.5 -13.8 -13.6 -12.6 -12.  -12.2
           -8.4  -8.8  -9.1  -9.9 -11.2 -13.2 -13.  -11.9 -11.5 -11.4
           -7.2  -7.7  -8.1  -8.9 -10.7 -12.9 -12.7 -11.7 -11.3 -11.1
           -6.9  -7.0  -7.2  -8.0 -10.0 -11.9 -11.9 -11.1 -11.0 -10.9
           -6.0  -6.0  -6.1  -6.7  -8.1  -9.3  -9.0  -8.6  -8.8  -8.9 
           -4.2  -4.4  -4.5  -4.8  -5.2  -5.8  -5.7  -5.6  -5.6  -5.5 
           -2.1  -2.3  -2.4  -2.4  -2.5  -2.6  -2.6  -2.6  -2.6  -2.6];

    H4 = [  0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0  
            2.6   2.6   2.4   2.3   2.1   2.1   2.1   2.2   2.3   2.4 
            5.1   5.0   4.8   4.4   4.0   3.8   3.7   3.7   3.7   3.6 
            7.4   7.3   7.0   6.3   5.6   5.1   4.7   4.4   4.4   4.3 
            8.9   8.9   8.5   7.9   7.0   6.3   5.6   5.0   4.7   4.4 
            9.9  10.0   9.7   9.1   8.1   7.2   6.4   5.9   5.4   5.2 
           10.1  10.3  10.0   9.5   8.7   7.7   7.0   6.4   6.2   6.4 
            8.8   9.2   9.0   8.6   7.9   7.2   6.7   6.4   6.4   6.6 
            6.3   7.0   7.2   6.9   6.3   5.9   5.5   5.5   5.7   6.0
            3.3   4.0   4.6   4.3   3.8   3.4   3.3   3.6   4.0   4.5 
            0.1   0.7   1.8   1.6   1.0   0.7   0.8   1.1   1.6   2.1 
           -2.7  -2.2  -0.8  -1.0  -1.4  -1.8  -1.8  -1.5  -1.0  -0.7 
           -4.8  -4.5  -3.2  -3.2  -3.8  -4.2  -4.4  -4.4  -4.0  -3.7 
           -6.6  -6.4  -5.7  -5.6  -5.9  -6.4  -6.8  -7.0  -6.9  -6.7 
           -9.3  -9.0  -8.1  -8.0  -8.2  -8.8  -9.3  -9.6  -9.7  -9.6 
          -11.7 -11.4 -10.6 -10.2 -10.5 -11.0 -11.9 -12.5 -12.9 -13.0
          -12.4 -12.1 -11.5 -11.2 -11.5 -12.3 -13.3 -14.1 -14.6 -14.8
          -12.0 -11.7 -10.7 -10.4 -10.9 -11.9 -13.0 -13.9 -14.5 -14.8
          -11.1 -10.6  -9.3  -8.9  -9.0  -9.6 -10.7 -12.0 -13.1 -13.7
          -10.8 -10.4  -8.8  -8.0  -8.0  -8.9 -10.0 -11.0 -11.8 -12.1
          -10.6 -10.2  -8.8  -8.2  -8.3  -9.0 -10.0 -11.0 -11.3 -11.3
           -8.7  -8.4  -7.8  -7.6  -7.7  -8.7  -9.8 -10.7 -10.8 -10.4
           -5.5  -5.5  -5.4  -5.2  -5.0  -5.4  -6.3  -7.0  -7.3  -7.0 
           -2.6  -2.6  -2.6  -2.5  -2.4  -2.4  -2.7  -3.0  -3.1  -3.0];

    H5 = [  0.0   0.0   0.0  
            2.3   2.0   1.2 
            3.3   2.8   1.6 
            3.7   3.0   1.6 
            4.0   3.3   2.3 
            5.1   5.0   5.0  
            6.8   7.4   8.2 
            7.3   8.0   9.0  
            6.7   7.3   8.0  
            5.1   5.6   6.2 
            2.6   3.1   3.7 
           -0.2  -0.2   0.5 
           -3.3  -2.9  -2.5 
           -6.3  -5.9  -5.2 
           -9.4  -9.0  -8.5 
          -12.9 -12.4 -11.6
          -14.8 -14.3 -13.9
          -15.0 -14.8 -14.5
          -14.0 -14.2 -14.2
          -12.3 -12.2 -12.0 
          -11.0 -10.0  -9.0  
           -9.6  -8.1  -6.5 
           -6.4  -5.1  -4.0  
           -2.8  -2.3  -1.6];

    // HRTF matrix, each row corresponds to an
    // azimuth and each column to a frequency
    HRTF0 = [H1 H2 H3 H4 H5];
    
    // Replace theta with congruent angle between 
    // 0° and 360° 
    theta = theta - floor(theta/360)*360;
    
    // Set flag if left side is the exposed one
    if theta<180
        left = 1; 
    else
        left = 0;
    end
    
    // Left ear HRTF
    // -------------
    // Azimuth respect to left side
    thetaL = theta
    // Indices of azimuths before and after the
    // azimuth for the left ear
    ntheta1L = floor(thetaL/15) + 1;
    ntheta2L = floor(thetaL/15) + 2; 
    // Interpolation of the HRTF for the input
    // azimuth
    hrtfL = interp1(theta0, HRTF0, thetaL);
    
    // Right ear HRTF
    // --------------
    // Azimuth respect to right side
    thetaR = 360 - theta;
    // Interpolation of the HRTF for the input
    // azimuth
    hrtfR = interp1(theta0, HRTF0, thetaR);

    // Frequency vector between 0 and Fs/2
    f = (0:N/2-1)*Fs/N;
    
    // Mean distance between ears in m
    L = 0.17;
    // Speed of sound
    c =345;
    // Estimated interaural time delay (right - left)
    // (Blauert, 1997)
    ITD = L/2*(theta*%pi/180 + sin(theta*%pi/180))/c;
    // Interaural phase delay 
    phiITD = ITD * 2*%pi*f;
    // Delay
    Td = 0.01;
    // Phase due to delay
    phiTd = Td * 2*%pi*f; 
    
    // Interpolated frequency responses in dB
    // As the Shaw HRTF response is given from 
    // 200 Hz to 12000 Hz only, a slope starting 
    // from 0 dB and reaching the value at 200 Hz
    // is assumed between 0 Hz and 200 Hz. This 
    // means very little attenuation at low
    // frequency. 
    // For the exposed ear, constant HRTF is 
    // assumed from 12000 Hz to Fs/2, while for 
    // the unexposed ear a slope from the value
    // at 12000 Hz to -100 dB is assumed between 
    // 12000 Hz and 12500 Hz, i.e., very large 
    // attenuation at high frequency
    // Left ear
    if left==1
       HL = interp1([0,f0,f($)], [0,hrtfL,hrtfL($)], f);
    else
       HL = interp1([0,f0,12500,f($)], [0,hrtfL,-100,-100], f);
    end
        
    // Right ear
    if left==0  
       HR = interp1([0,f0,f($)], [0,hrtfR,hrtfR($)], f);
    else
       HR = interp1([0,f0,12500,f($)], [0,hrtfR,-100,-100], f);
    end
     
    // Stereo unilateral frequency response in dB
    HRTF = [HL; HR];

    // Interpolated linear-scale frequency response
    if left==1
        // Left ear
        hL = 10.^(HL/20).*exp(-%i*phiTd);
        // Right ear 
        hR = 10.^(HR/20).*exp(-%i*(phiTd + phiITD));
    else
        // Left ear
        hL = 10.^(HL/20).*exp(-%i*(phiTd + phiITD));
        // Right ear 
        hR = 10.^(HR/20).*exp(-%i*phiTd);
    end
    
    // Symmetrically-extended frequency response
    // Left ear 
    HHL = [hL, 0, conj(hL($:-1:2))];
    // Right ear 
    HHR = [hR, 0, conj(hR($:-1:2))];
        
    // Impulse response
    // Right ear
    irL = ifft(HHL);
    // Left ear
    irR = ifft(HHR);

    // Binaural impulse response
    ir = [irL; irR];
    
    // Binaural signal
    y = [convol(irL, x); convol(irR, x)];

    // Binaural signal time vector
    t = (0:size(y,2)-1)/Fs;
   
    
    // For testing purposes
    if 1==2
    
    // Play mono and binaural signal
    playsnd(x,Fs,1)
    playsnd(y,Fs,1)
    
    // Mono signal time vector
    tx = (0:size(x,2)-1)/Fs;
        
    // Impulse response time vector
    tir = (0:N-1)/Fs; 

    // Extended frequency vector between 0 and Fs 
    FF = (0:N-1)*Fs/N;

    // Indices of azimuths before and after the
    // azimuth for the right ear
    ntheta1R = floor(thetaR/15) + 1;
    ntheta2R = floor(thetaR/15) + 2; 

    
    // Plot original Shaw data and interpolated HRTF
    // for left and right ears
    scf(1);
    clf(1);
    subplot(2,1,1)
    plot(f0, HRTF0(ntheta1L,:), f0, hrtfL, f0, HRTF0(ntheta2L,:))
    xlabel("f  [Hz]")
    ylabel("HRTF_L  [dB]")
    legend(string(theta1L)+"° ",string(thetaL)+"° ",string(theta2L)+"° ");
    subplot(2,1,2)
    plot(f0, HRTF0(ntheta1R,:), f0, hrtfR, f0, HRTF0(ntheta2R,:))
    xlabel("f  [Hz]")
    ylabel("HRTF_R  [dB]")
    legend(string(theta1R)+"° ",string(thetaR)+"° ",string(theta2R)+"° ");
    
    // Plot left and right interpolated, bilateral
    // HRTF  
    scf(2);
    clf(2);
    subplot(2,1,1)
    plot(FF, 20*log10(abs(HHL)), f0, hrtfL, "or")
    xlabel("f  [Hz]")
    ylabel("hrtf_L  [dB]")
    gca().data_bounds = [0 min([hrtfR,hrtfL]); 44100 max([hrtfR,hrtfL])];
    subplot(2,1,2)
    plot(FF, 20*log10(abs(HHR)), f0, hrtfR, "or")
    xlabel("f  [Hz]")
    ylabel("hrtf_R  [dB]")
    gca().data_bounds = [0 min([ hrtfR, hrtfL]); 44100 max([hrtfR,hrtfL])];

    // Plot left and right impulse responses 
    scf(3);
    clf(3);
    subplot(2,1,1)
    plot(tir, irL)
    ylabel("irL")
    xlabel("t   [s]")
    gca().data_bounds = [0 -0.3; 0.1 1.1]; 
    subplot(2,1,2)
    plot(tir, irR)
    xlabel("t   [s]")
    ylabel("irR")
    gca().data_bounds = [0 -0.3; 0.1 1.1]; 
    
    
   // Axis limits
    maxt = max(t)
    maxy = max(y)
    miny = min(y)    
    
    // Plot mono signal and left and right 
    // binaural signal
    scf(4);
    clf(4);
    subplot(3,1,1)
    plot(tx, x)
    ylabel("x")
    xlabel("tx   [s]")
    gca().data_bounds = [0 miny; maxt maxy]; 
    subplot(3,1,2)
    plot(t, y(1,:))
    xlabel("t   [s]")
    ylabel("yL")
    gca().data_bounds = [0 miny; maxt maxy]; 
    subplot(3,1,3)
    plot(t, y(2,:))
    xlabel("t   [s]")
    ylabel("yR")
    gca().data_bounds = [0 miny; maxt maxy]; 
    
    end

endfunction
