function [y, f] = fftmeandB(x, Fs, calwav, caldB, figura, N)
   // Análisis de espectro FFT promedio de señales de audio
   //
   // Este script obtiene el análisis de espectro por FFT
   // de una señal x muestreada a tasa Fs con ventana de
   // Blackman.
   // 
   // Uso:
   //       [y, f] = fftmeandB(x, calwav, caldb, figura, N)
   //
   // donde
   //       y:      Vector que contiene los valores del
   //               espectro en dB
   //       f:      Vector que contiene los valores de
   //               frecuencia 
   //       x:      vector que contiene la señal
   //       Fs:     Tasa de muestreo
   //       calwav: nivel del tono de calibración en dB con
   //               respecto al valor máximo (0 dB = +/- 1)
   //       caldB:  nivel acústico en dB correspondiente al
   //               tono de calibración 
   //       figura: indica si se dibuja un gráfico y en tal
   //               caso de qué tipo
   //       N:      tamaño de la ventana de análisis en
   //               muestras 
   //
   // El vector y tiene N/2 componentes, cada una de ellas 
   // corresponde a una frecuencia k*Fs/N, k = 1, ..., N/2. 
   // Las  componentes  son  los  valores medios  de las 
   // amplitudes (no los valores RNS) de las componentes 
   // espectrales homólogas para sucesivios cuadros de N
   // muestras de la señal.
   //
   // Para el procedimiento de calibración se supone que 
   // se ha grabado un tono de calibración equivalente a
   // un nivel de presión sonora caldB (en dB) y que el
   // valor eficaz de su versión numérica (limitada al 
   // intervalo entre -1 y 1) pasado a dB referidos al
   // máximo es calwav. Los valores por defecto son 94 dB
   // y 0 dB respectivamente, lo que significa que una
   // señal de máximo valor eficaz (0 dB digital) es
   // equivalente a un nivel de presión sonora de 94 dB.
   // El argumento figura deberá ser 'log', 'lin', o 
   // 'nof' según se desee que el eje de frecuencias sea
   // logarítmico, lineal o que directamente no haya
   // figura. El valor por defecto es 'log'. El valor N
   // por defecto es 32768.
   //
   // NOTAS: 
   // 1) Para evitar problemas de memoria la señal
   //    x no debería superar una extensión aproximada
   //    de 1 minuto.
   // 2) Si el tono de calibración es de máxima amplitud,
   //    su nivel respecto a la referencia será de 
   //    -3.0103 dB
   //
   // ----------------------
   // Autor: Federico Miyara
   // Fecha: 2008-04-24
   //        2008-08-02
   //        2010-04-14
   //        2011-03-20
   //        2011-10-15
   //        2012-12-24
   //        2016-01-02
   //        2024-12-10

   // ARGUMENTOS --------------------------------------

   rhs = argn(2);
   if rhs < 6 
      // Longitud de la FFT
      N = 32768;
   end
   if rhs < 5
      // Tipo de figura
      figura = 'log'; 
   end
   if rhs < 4
      // Calibración acústica
      caldB = 94;
   end
   if rhs < 3
      // Calibración digital   
      calwav = 0;
   end

   // Constante de calibración
   cal = caldB - calwav; 
   
   // Separación entre ventanas           
   L = N/2;
   // Si L < N entonces las ventanas se solapan. Si L > N 
   // las ventanas están separadas. Por simplicidad L no
   // es un argumento. Si se lo requiere distinto de N/2 
   // se debe modificar la función.

   // Asegurarse de que x es un vector columna
   x = x(:);

   if 1==2
      Fs = 44100
      caldB =0
      calwav = 0
      figura = 'nof'
      N = 4096
      x = y;
      cal = caldB - calwav;                                 
      L = N/2;  
   end

   // -------------------------------------------------

   // Determina tamaño de la señal

   // Cantidad de muestras
   muest = length(x);
   // Periodo de muestreo
   T = 1 / Fs;
   // Cantidad de cuadros de análisis
   q = floor((muest - N) / L) + 1;   
   // Cantidad de muestras a analizar
   n = (q - 1) * L + N;             

   // El proceso que sigue toma los cuadros de longitud
   // N de "x" espaciados L puntos y los coloca como
   // columnas de una matriz "cuadros" de q columnas
   // (una por cuadro) y N filas (una por cada muestr
   //a del correspondiente cuadro). Las últimas N - L
   // muestras de cada cuadro se repiten al comienzo
   // del siguiente. Ello proporciona un solapamiento
   // que permite una mejor localización temporal a
   // costa de mayor cantidad de cómputos. Tomando 
   // L = N los cuadros son adyacentes. 

   cuadros = zeros(N,q);
   for k = 1:q
      cuadros(:,k) = x(L*(k-1)+1: L*(k-1)+N);
   end;

   // A continuación se multiplica cada columna de 
   // "cuadros" por la ventana de Blackman de N puntos.
   // Se obtiene una matriz b cada una de cuyas
   // columnas es un cuadro enventanado por Blackman.
   // Luego se le aplica una FFT de N puntos, que 
   // devuelve en cada columna la FFT de la columna
   // original correspondiente 

   wb = blackman(N);
   //wb = hamming(N);
   //wb = hanning(N);
   //wb = boxcar(N);

   b = zeros(N,q);
   for k = 1:q
      b(1:N,k) = wb' .* cuadros(1:N,k);
   end;

   // Obtención de la FFT
   z = fft(b, -1, 1);

   // En el siguiente paso se obtiene el espectro
   // promedio del segmento. Para ello se obtiene
   // primero una matriz con los valores absolutos de
   // todos los coeficientes. Luego se aplica la
   // función mean, especificando 2 como segundo
   // argumento para que promedie por filas, inclusive
   // si la señal fuera tan corta que incluye sólo un
   // cuadro (recordemos que los cuadros están por
   // columnas). Promediar una fila significa promediar
   // los valores del espectro de los diferentes
   // cuadros a la frecuencia que corresponde a esa
   // fila. Resulta un vector columna de espectro
   // promedio del segmento. Para ahorrar tiempo de
   // cálculo el promedio se efectúa sólo sobre la
   // primera mitad de los coeficientes, ya que los
   // otros son idénticos (simétricos). 
   fftprom = mean(abs(z(1:N/2,:)),2);

   // Alternativamente, se puede hacer el promedio
   // energéticamente (al final se volverá a expresar
   // en amplitud aplicando sqrt) 
   fftprom2 = mean(abs(z(1:N/2,:)).^2, 2);

   // Finalmente se multiplican por 1/N el primer
   // coeficiente y por 2/N los coeficientes 
   // siguientes, para obtener un espectro en términos
   // de cosenos (espectro real).
   c = fftprom * 2 / N;
   c(1) = c(1) / 2;

   c2 = sqrt(fftprom2) * 2 / N;
   c2(1) = c2(1) / 2;

   // Ahora es necesario corregir por efecto de la
   // ventana y por calibración. Dada una señal de 
   // régimen estacionario de frecuencia mayor que 
   // 2*Fs/N, la ventana de Blackman atenúa su energía
   // en un valor teórico dado por su propio valor
   // eficaz:
   bef = sqrt(0.42^2 + 0.5^2/2 + 0.08^2/2);
   Lb = 20*log10(bef);

   // Considerando el valor máximo alcanzable por una
   // onda digitalizada (+/- 1) como 0 dB, el nivel RMS
   // para una senoide de máxima amplitud es -3,01 dB. 
   // El valor obtenido por FFT en la frecuencia más
   // próxima  a  la  de  la  senoide con ventana de
   // Blackman está entre -8,64 dB y -7,54 dB, según 
   // si la frecuencia de la senoide cae justo en el 
   // medio entre dos líneas espectrales o si coincide
   // con una línea espectral. En promedio es -8,1 dB.
   // Para obtener el valor  RMS  referido al valor
   // máximo se debe, por lo tanto,  sumar 5.1 dB. Con
   // respecto a esta misma referencia el tono de
   // calibración grabado, que corresponde a un valor
   // acústico igual a caldB, tiene un nivel igual a 
   // calwav. Para obtener los valores acústicos se
   // debe sumar, entonces:
   // 
   // caldB - calwav + 5.1 dB
   //
   // es decir, multiplicar c y c2 por 
   // 
   // 10^((caldB - calwav + 5,1)/20)
   //
   c = c * 10^((caldB - calwav + 5.1)/20);
   c2 = c2 * 10^((caldB - calwav + 5.1)/20);

   // Espectro en dB
   y = 20 * log10(c2);

   // Obtención de un vector f que contiene los
   // valores de frecuencia real
   f = [1 : N/2] * Fs/N;

   select figura
   case 'log'
      plot("ln", f, 20*log10(c2))
      gca().grid = [1 1];
      //axis([10 20000 -30 110])
   case 'lin'
      plot(f, 20*log10(c2))
      gca().grid = [1 1];
      //axis([10 20000 -30 110])
   case 'nof'
   end

endfunction

function w = blackman(N)
    // Periodic Blackman window
    // Usage:
    //       w = blackman(N)
    // Where 
    //       N: Length of Blackman window
    //       w: Vector containing the periodic Blackman window
    //
    // NOTE: The symmetric Blackman window can be obtained as
    //       ws = [blackman(N-1), 0]
    //     
    // ---------------------------------- 
    w = 0.42 - 0.5*cos(2*%pi*(0:N-1)/N) + 0.08*cos(4*%pi*(0:N-1)/N);
endfunction
