function [y, yenv, fo] = bands(x, typ, Fs, tau, nd)
   // Given a frequency band type (for instance one third
   // octave) this function provides filtered versions of
   // an input signal for all the bands covering the audio
   // spectrum (20 Hz to 20000 Hz), as well as the envelope
   // of each filtered signal and the corresponding center
   // frequencies. 
   //
   // Usage:
   //       [y, yenv, fo] = bands(x, typ, Fs, tau, nd)
   //
   // where 
   //       y:     filtered otput; a matrix with as many 
   //              rows as there are bands, each one of 
   //              them containing the signal filtered by 
   //              one band filter
   //       yenv:  energy envelopes; a matrix with an 
   //              envelope per row
   //       fo:    Column vector containing the geometric 
   //              centre of each band; in the case of 
   //              custom bands, two-column vector whose
   //              first column is the lower limit and 
   //              whose second column is the upper limit
   //       x:     vector containing input signal
   //       typ:   equal speech intelligibility (-1), 
   //              critical (0), octave (1), one-third 
   //              octave (3), etc.; also any custom set 
   //              of frequencies which define adjacent 
   //              frequency bands, arranged as a vector. 
   //       Fs:    sampling rate in Hz
   //       tau:   time constant for energy envelope in s
   //       nd:    no-downsampling flag; if 1, signal is
   //              not downsampled at low-frequency bands,
   //              increasing accuracy at the expense of
   //              computation time  
   //
   // Fs, typ, tau and nd are optional, being their default 
   // values 44100, 1, 0.035 and 0.
   // The envelope of a signal is the square root of its
   // square filtered with a low-pass filter of time 
   // constant tau.
   // ---------------------------------------------------
   // Author: Federico Miyara
   // Date:   2019-03-16
   //         2020-01-14 
   //         2021-08-29
   //         2025-01-12
   
   // Dealing with input arguments

   [lhs, rhs] = argn()

   if rhs<5 
      // If no-downsampling flag is omitted, set it
      // to 0 (downsamplig enabled, gaining speed) 
      nd = 0;
   end

   if rhs<4 
      // If tau is not given, set it to 0.035, i.e.,
      // "impulsive" time weighting  
      tau = 0.035;
   end

   if rhs<3 
      // If sampling rate is not given, set it to
      // 44100 
      Fs = 44100;
   end

   if rhs<2 
      // If band type is not given, set it to 1
      // (octave bands) 
      typ = 1;
   end

   // Band frequencies
   if isscalar(typ)
      select typ
      case -1
         finf = [ 250  375  505  645  795  955 1130 1315 1515 1720 ..
                 1930 2140 2355 2600 2900 3255 3680 4200 4860 5720];
         fsup = [ 375  505  645  795  955 1130 1315 1515 1720 1930 ..
                 2140 2355 2600 2900 3255 3680 4200 4860 5720 7000];
      case 0
         // If it is a critical band 
         exec zbark.sci;
         exec zkrab.sci;
         f = zkrab(1:24)';
         finf = zkrab(0.5:23.5)';
         fsup = zkrab(1.5:24.5)';
      else
         // If it is a fractional octave band
         nmin = floor(typ*log2(20/1000));
         nmax = ceil(typ*log2(20000/1000));
         fo = 1000 * 2.^((nmin:nmax)/typ)';
         finf = fo / 2^(1/typ/2);
         fsup = fo * 2^(1/typ/2);
         // Adjust first and last band so to slightly
         // overlap with the audio range [20, 20000]
         imin = max (find(round(finf)<20));
         imax = min (find(round(fsup)>=20000));
         fo = fo(imin:imax);
         finf = finf(imin:imax);
         fsup = fsup(imin:imax);
      end
   else
      // typ is a frequency vector
      // Reorder in ascending order 
      typ = gsort(typ);
      typ = typ($:-1:1);
      // Convert to column  
      if isrow(typ)
         typ = typ';
      end
      // Define finf, fsup and fo
      finf = typ(1:$-1);
      fsup = typ(2:$);
      fo = sqrt(finf.*fsup);
   end
   //[finf,fo, fsup] 

   // Pre-allocate y and yenv
   n = length(fo);
   Lx = length(x); 
   y = zeros(n, Lx);
   yenv = zeros(n, Lx);

   // Numerator and denominator coefficients of low-pass
   // filter used to get the energy envelope (in descending
   // power order). They were derived applying the bilinear
   // transform to a first order low-pass Laplace transfer 
   // function with time constant tau.
   NUM = [1, 1];
   DEN = [1 + 2*tau*Fs, 1 - 2*tau*Fs];

   // Load function
   exec bandfilter.sci;

   for k=1:n
      // Use custom syntax of bandfilter where fo is a 
      // two-component vector containing finf and fsup 
      fn = [finf(k), fsup(k)];
      disp([finf(k),fsup(k)])
      y(k, :) = bandfilter(x, fn, 1, Fs, nd);
      yenv(k, :) = sqrt(filter(NUM, DEN, y(k, :).^2));

   end

endfunction
