function [z] = zbark(f, m)
    // Computes critical band rate z corresponding to frequency f
    // according to various optional models, such as the Zwicker 
    // and Terhardt equation and others. Applies to scalar, vector 
    // or matrix. 
    //
    // Usage: 
    //       z = zbark(f, m)
    // where
    //       f: scalar, vector or matrix frequency in Hz
    //       z: corresponding critical band rate in Bark
    //       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 
    //
    // See also zkrab, 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-26
    //       2008-06-17
    //       2020-03-07
    //       2024-10-27   

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

    select m
    case 0
        // Zwicker formula
        z = 13 * atan(f/1316) + 3.5 * atan((f/7500).^2);
    case 1
        // Buus formula
        a = 1750;
        b = 1.81;
        c = 9800;
        k = 28.32;
        z = k - c/(b - 1) ./ (f + a).^(b - 1);
    case 2
        // Modified Buus formula 1
        a = 2070;
        b = 1.8322;
        c = 13704;
        k = 28.43;
        z = k - c/(b - 1) ./ (f + a).^(b - 1);
    case 3
        // Modified Buus formula 2
        a = 2090;
        b = 1.8225;
        c = 12900;
        k = 29.15;
        z = k - c/(b - 1) ./ (f + a).^(b - 1);
    case 4
        // Traunmüller formula
        z = 26.81*f./(1960 + f) - 0.53;
    case 5
        // Schroeder formula
        z = 7 * asinh(f/650);
     case 6
        // Wang, Sekey & Gersho, 1992
        z = 6 * asinh(f/600);
    end

endfunction
