function [S, NS, LA] = loudness(x, ref, Fs)
   //   Sonoridad y nivel de sonoridad
   //   loudness(x) calcula la sonoridad y el nivel de
   //   sonoridad de un vector x por el método Mark VI
   //   de Stevens detallado en la Norma ISO 532A a
   //   partir de la señal o de su espectro en bandas 
   //   de octava.
   //
   //   Uso:   
   //         [S, NS, LA] = loudness(x, ref, Fs)
   //   donde
   //         x:   Vector fila o columna que contiene 
   //              la señal cuya sonoridad se quiere 
   //              obtener o matriz de m columnas que
   //              contiene una señal por columna, o 
   //              bien los niveles de las bandas de
   //              octava qntre 31.5 Hz y 16000 Hz
   //         ref: Nivel de presión sonora en dB 
   //              correspondiente a una onda senoidal
   //              de amplitud 1 (típicamente obtenida 
   //              a partir de una señal de calibración)
   //         Fs   Frecuencia de muestreo. 
   //         S:   Sonoridad en son
   //         NS:  Nivel de sonoridad en fon
   //         LA:  Nivel de presión sonora ponderada A

   //   Si no se indica ref se la supone 100 dB. Si no 
   //   se indica Fs se la supone 44100 Hz.
   //   loudness admite vectores fila o columna, pero 
   //   si x es una matriz de n filas y m columnas, el
   //   análisis se lleva a cabo por columnas, siendo
   //   S, NS y LA vectores fila de m componentes, una 
   //   por cada columna. 
   //   La función determina inteligentemente el tipo
   //   de entrada.
   //   La versión con un solo argumento de salida
   //   calcula un vector fila de m componentes con 
   //   las sonoridades de las columnas de x.
   //
   //   -------------------------------
   //   Autor: Federico Miyara
   //   Versión: 2.0
   //   Fecha: 2005-04-05
   //          2024-10-20  
   //          

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

   // Tamaño del argumento 
   [n, m] = size(x);
   // Si el argumento es un vector fila se convierte
   // en columna y se recalcula el tamaño del argumento
   if n == 1
      x = x';
      [n, m] = size(x);
   end

   if n > 10
      // Si x tiene más de 10 componentes, interpretar
      // x como una señal y obtener el nivel de presión
      // sonora de las bandas de octava
      prot = funcprot();
      funcprot(0);
      exec octaves.sci;
      funcprot(prot);
      Lpf = 20*log10(octaves(x, 1, Fs, 16384)) + ref;

      // Nivel de presión sonora (según Parseval)
      Lp_parseval = 10*log10(sum(10.^(Lpf/10))); 

      // Nivel de presión sonora (Cálculo directo)
      Lp = 10*log10(mean(x.^2, 1)) + ref;

      // Detecta si hay error por derrame espectral
      // comparando el nivel equivalente de presión
      // sonora calculado en forma directa a partir
      // de x y el calculado por aplicación de Parseval
      delta = abs(Lp - Lp_parseval)
      if delta > 3
         mprintf('Advertencia: error por derrame espectral > 3 dB\n')
      elseif delta > 2
         mprintf('Advertencia: error por derrame espectral > 2 dB\n')
      else delta > 1
         mprintf('Advertencia: error por derrame espectral > 1 dB\n')
      end
   else
      // Si x tiene 10 o menos componentes, interpretar
      // x como un espectro y completar con ceros las 
      // componentes faltantes de las bandas superiores
      Lpf = x;
      if length(x)< 10
         Lpf(length(x)+1:10) = 0;
      end
   end

   // Cálculo de la sonoridad en son
   // El procedimiento es el siguiente: dado un par
   // (f, Lp), se determina el nivel de presión sonora 
   // Lp1000eq tal que (1000, Lp1000eq) tenga igual
   // sonoridad que (f, Lp). Luego se aplica una fórmula
   // que aproxima a la tabla dada en ISO 532A para
   // determinar el índice de sonoridad correspondiente
   // LI.  Se repite para cada banda de octava y los
   // resultados se suman ponderados por un factor 0.3 
   // excepto el LI máximo que se pondera por un factor
   // 1, obteniéndose la sonoridad S (la diferente
   // ponderación contempla el enmascaramiento de la
   // banda más sonora sobre las otras). Hay tres
   // regiones en el plano f-Lp. Según la región en la
   // que se encuentra (f, Lp) el correspondiente
   // Lp1000eq se calcula con una fórmula diferente.

   // Cálculo de las frecuencias de bandas de octava
   // exactas Sólo se incluyen las bandas que cumplen 
   // con la condición de Nyquist. A efectos de su uso
   // posterior, se repiten las frecuencias en m bandas
   // iguales.
   bandas = floor(6.5 + log10(Fs/2000)/log10(2));  
   f = 2.^[-5:bandas-6]'*1000 * ones(1, m);

   // Vectores bandera. Cada componente h de Ri es 1 si
   // la frecuencia fh y nivel Lph están en la región i
   // y 0 si no lo está. Para cada frecuencia fh, sólo
   // uno de los Ri es 1. Se aprovecha la naturaleza de
   // operación binaria  con resultados 1 ó 0 del operador
   // relacional <. 
   R1 = log10(f)<(220-Lpf)/70;
   R2 = ((220-Lpf)/70<log10(f)) .* (log10(f)<log10(9000));
   R3 = log10(9000)<log10(f);

   // Cálculo del niveles de presión sonora equivalente 
   // para f = 1000 Hz (utilizando para ello la fórmula
   // correspondiente a cada región)
   Lp1000eq = R1 .* (6/5*Lpf + 24*log10(f) - 74) + ...
   R2 .* (Lpf + 10*log10(f) - 30) + ...
   R3 .* (Lpf - 40*log10(f) + 167.7);

   // Índices de sonoridad aplicando la aproximación de
   // Control de Ruido, Cap 1
   LI = 10.^(-3.491 + 0.2172*Lp1000eq - 5.234e-3*Lp1000eq.^2 + ...
             6.722e-5*Lp1000eq.^3 - 4.079e-7*Lp1000eq.^4 + ...
             9.524e-10*Lp1000eq.^5);

   // Máximo índice de sonoridad  
   LIm = max(LI) 

   // Sonoridad según Stevens, en son o, más propiamente,
   // son(OD)
   S = 0.3 * sum(LI) + 0.7 * LIm

   // Nivel de sonoridad según la función de transferencia
   // en fon o, más propiamente, fon (OD)
   NS = 10* log10(S)/log10(2) + 40           

   // Ponderación A
   A  =  [-39.4 -26.2 -16.1 -8.6 -3.2  0.0  1.2  1.0 -1.1 -6.6]';
   LA = 10*log10(sum(10.^((Lpf + A*ones(1,m))/10), 1));

endfunction
