function [f] = zkrab(z, m)
    // Computes the frequency f corresponding to critical band 
    // rate z according to various optional models. In the case
    // of Zwicker and Terhardt's equation, the method of 
    // Newton-Raphson is applied to the original. In all other 
    // cases an explicit inverse formula is used. Applies to 
    // scalar, vector or matrix.
    //
    // Usage:
    //       f = zkrab(z, m)
    // where
    //       z: scalar, vector or matrix critical band rate in Bark
    //       f: corresponding frequency in Hz
    //       m: method
    //          0 or none: Zwicker formula
    //          1:         Buus formula
    //          2:         Modified Buus formula 1
    //          3:         Modified Buus formula 2 
    //          4:         Traunmüller formula
    //          5:         Schroeder formula
    //          6:         Wang-Sekey-Gersho formula 
    //
    // NOTE: The play on words ("krab" is the reverse of "bark" and 
    //       a paronym of "crab") is intentional!
    //
    // See also: zbark, critband, fcritband
    //
    // References: 
    // Zwicker & Fastl, "Psychoacoustics, Facts and Models".
    //       Springer 1999
    // Buus, Søren, "Auditory Masking", in Malcolm Crocker (ed),
    //      “Encyclopedia of Acoustics”. John Wiley & Sons. 1997
    // Traunmüller, H. "Analytical Expressions for the Tonotopic 
    //      Sensory Scale". JASA V. 88, No. 1 July 1990 pp 97-100    
    // Schroeder, Manfred R. (1977) Recognition of complex 
    //      acoustic signals. in Life Sciences Research Report 5 
    //      (Dahlem Konferenzen). edited by T. H. Bullock (Abakon 
    //      Verlag, Berlin) pp 323-328 
    // Wang, S.; Sekey, A.; Gersho, A. (1992) "An Objective Measure
    //      for predicting Subjetive Quality of Speech Coders" 
    //      IEEE Journal Sel. Areas Commun. Vol. 10, pp 819-829
    //
    // ---------------------------------------------------------------
    // Author: Federico Miyara
    // Date: 2006-07-22
    //       2008-06-17
    //       2012-06-21
    //       2018-03-18

    // If only one argument is provided, set method to default (Zwicker)
    if argn(2) == 1
        m = 0;
    end

    select m
    case 0
        // Numerical inversion of Zwicker formula
        // Initial point for Newton-Raphson
        f1 = zeros(z);
        z1 = 13 * atan(f1/1316) + 3.5 * atan((f1/7500).^2);

        // Iterate a number of times increasing with z
        // as the derivative decreases to keep accuracy
        for i=1:1 + max(max(ceil(z/3)))
            // Derivative of z at f1
            dzdf = (13/1316) ./ (1 + (f1/1316).^2) + ...
                   (3.5/7500^2)*2*f1./(1 + (f1/7500).^4);
            // New value of f1
            f1 = (z - z1)./dzdf + f1;
            // z1 corresponding to f1
            z1 = 13 * atan(f1/1316) + 3.5 * atan((f1/7500).^2);
        end
        // Last f1 is the solution
        f = f1;
    case 1
        // Buus formula
        a = 1750;
        b = 1.81;
        c = 9800;
        k = 28.32;
        f = (c / (b - 1) ./ (k - z)).^(1/(b - 1)) - a;
    case 2
        // Modified Buus formula 1
        a = 2070;
        b = 1.8322;
        c = 13704;
        k = 28.43;
        f = (c / (b - 1) ./ (k - z)).^(1/(b - 1)) - a;
    case 3
        // Modified Buus formula 2
        a = 2090;
        b = 1.8225;
        c = 12900;
        k = 29.15;
        f = (c / (b - 1) ./ (k - z)).^(1/(b - 1)) - a;
    case 4
        // Traunmüller formula
        f = 1960*(z + 0.53)./(26.28 - z);
    case 5
        // Shroeder formula
        f = 650 * sinh(z/7);
    case 6
        // Wang-Sekey-Gersho formula 
        f = 600 * sinh(z/6);
    end
    
endfunction


