function [E, f, LTH, EE] = mask(Lm, fm, N, m, delta, scale)
    // Combined masking pattern of a group of tones
    // The function computes the hearing threshold in the presence 
    // of given tones (fm, Lm). 
    // 
    // Usage:
    //        [E, f, LTH, EE] = mask(Lm, fm, N, m, delta, scale)
    // 
    // where
    //        E:     combined masking pattern in dB
    //        f:     frequencies where E is specified in Hz
    //        LTH:   Threshold level in silence in dB
    //        EE:    masking patterns of individual tones in dB,
    //               where rows contains the masking curve 
    //               corresponding to each masking tone 
    //        Lm:    Sound pressure level of tones in dB, may be
    //               scalar or vector 
    //        fm:    Frequencies of tones in Hz, may be scalar
    //               or vector 
    //        N:     Number of frequencies
    //        m:     Method or model to be used
    //               1: Terhardt
    //               2: Schroeder
    //               3: Schroeder modified to include effect of Lm  
    //               4: ISO/IEC MPEG Psychoacoustic Model 2
    //               5: ISO/IEC MPEG Psychoacoustic Model 1
    //        delta: correction term accounting for type of masker
    //               and maskee; see below
    //        scale: 'log' for logarithmic distribution of points
    //               'lin' for linear distribution of points
    //
    // Arguments scale, delta, m, N and Fs are optional. Their 
    // default values are scale = 'log', delta = 18, m = 1, 
    // N = 4096, Fs = 44100 respectively.
    // Lm and fm may be scalars or vectors of the same size.
    //
    // The underlying model for each masker is:
    //
    // E = Lm - delta + F
    //
    // where
    //       Lm:    Sound pressure level of masking signal 
    //       delta: correction term related to the differential
    //              threshold of the basilar membrane excitation;
    //              it depends on level, frequency and type of
    //               masking and maskee signal (see table below) 
    //       F:     spreading curve; it varies across empirical
    //              models, usually has a constant slope for
    //              z < zm and a variable slope for z > zm (where
    //              z is the critical band rate of masked signal
    //              and zm the critical band rate of masker)
    //
    // Some delta examples:
    // 
    // Situation                         delta [dB]
    // --------------------------------------------
    // Noise band masking tone                4
    // Tone masking tone                     18
    // Tone masking noise band               24
    // Noise band masking noise band         26
    //
    // References:
    //
    // Terhardt, Stoll, Seewann. "Algorithm for Extraction of
    //      Pitch and Pitch Salience From Complex Tonal Signals".
    //      JASA, Vol 71, Nº 3 (March  1982), pp. 679-688
    // Schroeder, Atal, Hall. "Optimizing Digital Speech Coder
    //      by Exploiting Masking Properties of the Human Ear".
    //      JASA vol 66, Nº 6 (Dec 1979), pp. 1647-1652
    // Painter, Spanias. "Perceptual Coding of Digital Audio". 
    //      Proceedings of the IEEE, VOL. 88, NO. 4, April 2000
    // ISO/IEC. "Information technology – Coding of moving pictures 
    //      and associated audio for digital storage media at up 
    //      to about 1,5 Mbit/s - Part 3: Audio. ISO/IEC 11172-3
    //      International Standard,  1993.
    // Johnston, "Transform Coding of Audio Signals Using Perceptual 
    //      Noise Criteria". IEEE Journal on Selected Areas in
    //      Communications, vol.  6, pp. 314-323 (1988 Feb.)] 
    //
    // NOTE: Requires zbark.sci
    //
    // -------------------------------------------------------------
    // Federico Miyara
    // 2011-08-01
    // 2012-06-21
    // 2014-11-30
    // 2015-10-28
    // 2020-10-02

    // Dealing with missing arguments
    rhs = argn(2);
    if rhs <6
        scale = 'log';
    end
    if rhs<5
        delta = 18;
    end
    if rhs<4
        m = 1;
    end
    if rhs<3
        N = 4096;
    end

    // Frequency vector
    Fmax = 20000;
    Fmin = 10;
    if scale=='lin'
        f = Fmin:(Fmax-Fmin)/(N-1):Fmax;
    else
        f = 10.^(log10(Fmin):log10(Fmax/Fmin)/(N-1):log10(Fmax));
    end

    // Standard hearing threshold in silence according to
    // Terhardt-Stoll-Seewann, JASA, Vol 71, No 3, p 683
    LTH = 3.64*(f/1000).^(-0.8) - 6.5*exp(-0.6*(f/1000 - 3.3).^2) + ...
          0.001*(f/1000).^4;

    // Critical band rates corresponding to frequency vector
    exec zbark.sci;
    z = zbark(f);

    // Critical band rate of masker
    zm = zbark(fm);

    H = length(Lm); 

    // Allocating space for intermediate variable 
    EE = zeros(H,N);

    for h = 1:H
        // Critical band rate increment from masking tone to signal
        dz = z - zm(h);

        // Masking spreading models
        select m
        case 1
            // Terhardt spreading model
            // We subtract 18 dB from the masker's level and 
            // starting from that point two different slopes at  
            // either side
            // Computation of masking curves slope in dB/bark
            // NOTE: 1) There are two different slopes when test 
            //          tone frecuency is smaller than masker's 
            //       2) double converts the logical type result 
            //          of batch comparisons into nunerical type
            smu = double(f<=fm(h)) * 27 + ...
            double(f>fm(h)) .* (-24 - 230./fm(h) + 0.2*Lm(h));
            F = smu.*dz;
        case 2
            // Schroeder spreading model 
            // NOTES: 1) High precision constant 15.8111389 is to
            //           get F = 0 at dz = 0 (maximum normalized
            //           spreading = 1)
            //        2) This curve is independent of Lnu 
            F = 15.8111389 + 7.5 * (dz + 0.474) - ...
                17.5 * sqrt(1 + (dz + 0.474).^2);
        case 3
            // Schroeder model modified to account for Lm
            // Derivative of f with respect to z (from classical 
            // analytical model due to Zwicker)
            df_dz = (13/1316) ./ (1 + (fm(h)/1316).^2) + ...
            (3.5/7500^2)*2*fm(h) ./ (1 + (fm(h)/7500).^4);
            // Adjustment due to level Lm and f derivative
            ILm_fm = min([2, 5*10.^((Lm(h) - 96)/10 * df_dz)])
            F = 15.8111389 - ILm_fm + 7.5 * (dz + 0.474) - ...
            (17.5 - ILm_fm) .* sqrt(1 + (dz + 0.474).^2);
        case 4
            // ISO/IEC MPEG Psychoacoustic Model 2 (modified 
            // from Schroeder model)
            // Does not account for masker level
            F = 15.8111389 + 7.5 * (1.05*dz + 0.474) - ...
            17.5 * sqrt(1 + (1.05*dz + 0.474).^2) + ...
            8 * min([0, (1.05*dz - 0.5).^2 - 2*(1.05*dz - 0.5)]);
        case 5
            // ISO/IEC MPEG Psychoacoustic Model 1 
            // Takes into account the masker level
            F = (-17*dz + 0.15*Lm(h)*(dz - 1).*double(dz>1)) .* ...
            double(dz>=0) + ... 
            (-(6 + 0.4*Lm(h))*abs(dz) - ...
            (11 + 0.4*LM(h))*(abs(dz)-1).*double(abs(dz)>1)) .* ...
            double(dz<0);
        end

        // Masking by tones   
        EE(h,:) = Lm(h) - delta + F;
    end

    // Compressive masking addition exponent
    alfa = 0.33;

    // Combined masking
    // NOTE: It is necessary to take the real part in order that 
    //       the max function operates correctly, since at f=0
    //       EE is NaN and NaN^(1/alfa) is a complex NaN (may be
    //       a bug of this version of Matlab)
    if H == 1
        EEE = real(10*log10((10.^(EE/10*alfa)).^(1/alfa))); 
    else
        EEE = real(10*log10(sum(10.^(EE/10*alfa), 1).^(1/alfa))); 
    end
    // Masked threshold including hearing threshold in silence
    // NOTE: Here we used the energy superposition of threshold 
    //       in silence and masked threshold. Alternatively, we
    //       could get the maximum between threshold in silence 
    //       and masked threshold, i.e., E = max([EEE; LTH]);
    E = 10*log10(10.^(EEE/10) + 10.^(LTH/10));

endfunction

