function [Y, X] = critenv(x, tau, Fs, Fsenv) 
   // Obtención de las envolventes energéticas de las bandas
   // críticas de una señal
   //
   // Uso:
   //       [Y, X] = critenv(x, tau, Fs, Fsenv)
   // donde
   //       x:     Vector columna con la señal a analizar,
   //       tau:   Constante de tiempo en s del detector 
   //              de envolvente (0,01 por defecto)
   //       Fs:    Tasa de muestreo en Hz (44100 por defecto) 
   //       Fsenv: Tasa de muestreo en Hz de las envolventes 
   //              (100 por defecto)
   //       Y:     Matriz de hasta 24 columnas, cada una de 
   //              las cuales contiene las envolventes de las
   //              bandas críticas ordenadas en bandas
   //              crecientes
   //       X:     Matriz de hasta 24 columnas con la señal 
   //              filtrada con filtros de ancho de banda 
   //              crítica, ordenados en bandas crecientes.
   // 
   // Se utilizan filtros multitasa (multirate) para filtrar
   // las sucesivas bandas críticas. Por debajo de Fs/5 se
   // submuestrea para evitar inestabilidades debidas a que
   // los polos de la transferencia Z están muy cerca de la
   // circunferencia unitaria y cualquier error de precisión
   // aritmética en los cálculos puede llevar a una respuesta
   // divergente. Luego se obtienen las envolventes, para lo
   // cual se requiere la función envelope().
   //
   // NOTA: Esta función puede reemplazarse por bands(), que 
   //       permite obtener las envolventes espectrales de 
   //       bandas de fracción de octava, bandas críticas
   //       y bandas personalizadas
   // 
   // -----------------------------------------------
   // Autor: Federico Miyara
   // Fecha: 2006-08-14
   //        2021-08-29

   // Argumentos por defecto si se omiten
   rhs = argn(2);
   if rhs < 4
      Fsenv = 100;
   end
   if rhs < 3
      Fs = 44100;
   end
   if rhs < 2
      tau = 0.01;
   end

   // Valores tabulados de los anchos de banda críticos
   // sucesivos según Zwicker
   BC = [ 100  100  100  100  110  120  140  150  160 ...
          190  210  240  280  320  380  450  550  700 ...
          900 1100 1300 1800 2500 3500];

   // Valores tabulados de límites de las bandas críticas
   // según Zwicker
   fC = [   0  100  200  300  400  510  630  770  920 ...
         1080 1270 1480 1720 2000 2320 2700 3150 3700 ...
         4400 5300 6400 7700 9500 12000 15500];

   // Búsqueda de la banda crítica más alta compatible con
   // la tasa de Nyquist
   kmax = max(find(fC < Fs/2)) - 1;

   // Frecuencia de muestreo para la iteración inicial
   fs = Fs;

   Lx = length(x);

   // Filtrado por bandas críticas y determinación de 
   // las envolventes. Se parte de la banda más alta y
   // se va bajando
   for k = kmax:-1:1
      // Frecuencias inferior y superior
      finf = fC(k);
      fsup = fC(k+1);
      // Frecuencias inferior y superior normalizadas 
      // (w = 1 corresponde a f = fs/2) 
      winf = finf / fs;
      wsup = fsup / fs;
      // Si la frecuencia inferior es muy pequeña se 
      // submuestrea para a) evitar problemas numéricos
      // y b) ahorrar tiempo de cálculo 
      if k > 1
         if winf < 0.2
            fs = fs/2;
            winf = finf / (fs/2);
            wsup = fsup / (fs/2);
            x = intdec(x, 0.5)
            //x = resample(x,1,2);
         end
      end

      // Diseño del filtro de banda crítica de Butterworth
      // de orden 4
      if k >1
         // Si no es la banda crítica más baja se usa un
         // pasabanda
         hz = iir(4, "bp", "butt", [wsup, winf], [])
      else
         // Si es la banda crítica más baja se usa un 
         // pasabajos 
         hz = iir(4, "lp", "butt", wsup, [])
      end
      // Numerador
      B = hz(2);
      // Denominador
      A = hz(3);

      // Filtrado
      y = filter(B, A, x);

      // Si hay 2 argumentos de salida se guarda la banda
      // crítica filtrada
      lhs = argn(1);
      if lhs == 2
         // Se remuestrea (si es necesario) la banda
         // filtrada para tener la tasa de muestreo original
         if fs < Fs
            // Se obtiene una representación racional con
            // tolerancia 1e-12
            [p, q] = rat(Fs/fs, 1e-12);
            // Remuestreo a la tasa original
            XX = intdec(y, p/q);
         else
            XX = y;
         end
         // Se asigna el vector filtrado a la columna k-ésima
         X(:,k) = XX(1:Lx); 
      end

      // Obtención de la envolvente energética
      y = envelope(y, tau, fs);
      // Submuestreo de la envolvente hasta el valor indicado
      // en Fsenv o 100 Hz si se omitió
      // Obtención de la añproximación racional
      [p, q] = rat(Fsenv/fs, 1e-12);
      // Remuestreo
      YY = intdec(y, p/q);
      Y(:,k) = YY(1:floor(Lx/Fs*Fsenv));
   end
   
endfunction
