function [cbw] = critband(f, m)
    // Computes the bandwidth of the critical band centered at 
    // frequency f. Applies to scalar, vector or matrix 
    //
    // Usage:
    //       cbw = critband(f, m)
    // where
    //       f:   frequency in Hz
    //       m:   method
    //            0 or none: Zwicker formula
    //            1:         Derivative of Zwicker formula for 
    //                       critical band rate
    //            2:         Internally consistent with Zwicker's
    //                       formula for critical band rate 
    //                       +/- 0.5 Bark
    //            3:         Buus formula
    //            4:         Modified Buus formula 1
    //            5:         Modified Buus formula 2 
    //            6:         Traunmüller formula
    //            7:         Schroeder formula 
    //       cbw: critical bandwidth in Hz computed by method 
    //            specifidby optional argument m
    //
    // See also fcritband, zbark, zkrab
    //
    // References: 
    // Zwicker & Fastl, "Psychoacoustics, Facts and Models"
    // Buus, "Auditory Masking". In Crocker, Malcolm (ed), 
    //      "Encyclopedia of Acoustics". John Wiley, 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     
    // ----------------------------------------------------
    // Author: Federico Miyara
    // Date: 2006-07-26
    //       2008-06-17
    //       2019-03-18

    // If only one argument is provided, set method to default
    if nargin == 1
        m = 0;
    end

    select m
    case 0
        // Zwicker formula for critical band width
        cbw = 25 + 75 * (1 + 1.4*(f/1000).^2).^0.69;
    case 1
        // Derivative of Zwicker formula for critical band rate 
        cbw = (zkrab(zbark(f) + 0.0001) - f) / 0.0001;
    case 2
        // Internally consistent with Zwicker formula for critical 
        // band rate 
        cbw = zkrab(zbark(f) + 0.5) - zkrab(zbark(f) - 0.5);
    case 3
        // Buus formula
        cbw = (f + 1750).^1.81 / 9800;
    case 4
        // Modified Buus formula 1
        cbw = (f + 2070).^1.8322 / 13704 
    case 5
        // Modified Buus formula 1
        a = 2090;
        b = 1.8225;
        c = 12900;
        cbw = ((f + a).^(1 - b) + (1-b)/c*0.5).^(1/(1 - b)) - ..
              ((f + a).^(1 - b) - (1-b)/c*0.5).^(1/(1 - b));     
    case 6
        // Traunmüller formula
        z = 26.81 * f ./ (1960 + f) - 0.53;
        cbw = 52547.6 ./ (z.^2 - 52.56*z + 690.3884);
    case 7
        // Schroeder formula
        z = 7*asinh(f/650);
        cbw = 650 * (sinh((z + 0.5)/7) - sinh((z - 0.5)/7));
    end
endfunction 
