Cosas que hacemos con cierta frecuencia

Qué encontrarás en esta entrada?

• Análisis de Fourier desde Octave/Matlab para sacar las notas de una canción.
• Comparación entre algoritmos para conseguir este objetivo.

En la última entrada hablamos de cómo sacar las notas de una canción con un sencillo algoritmo en código Octave/Matlab. Allí vimos que había casos en los que funcionaba medianamente bien, y otros en los que no demasiado.

El algoritmo en cuestión se basaba en una estimación de la frecuencia en base al número de nodos. Esto puede funcionar para una perfecta onda sinusoidal, pero todos los sonidos reales (desde el mismo instante en el que son finitos) son una superposición de ondas de distintas frecuencia. Aunque nos ha interesado desde el principio quedarnos con la frecuencia principal (despreciando los armónicos que resuenen a su son), se hacía indispensable un análisis de Fourier que nos diera más información sobre el espectro de la señal.

Mientras intentaba refrescar lo que dijo nuestro amigo Fourier (y que tenía bastante olvidado), he descubierto que Octave/Matlab cuentan con una función llamada "fft" (Fast Fourier Transform). Esta función nos da la gráfica del espectro, sin embargo, presenta tres complicaciones:

• Los resultados están desordenados. Hay que ordenarlos con la función "fftshift". 

• Los resultados son complejos. Hay que hacerlos reales, para ello se puede tomar el módulo con la función "abs".

• Una vez ordenados y pasados a reales los valores de la amplitud para cada frecuencia (lo que podríamos considerar gráficamente como el "eje y"), hay que escalar bien el eje de frecuencias ("eje x"). Si no, la relación entre la amplitud de cada frecuencia y el valor de la misma (que es lo que nos interesa) no sería correcta. Para ello hay que centrar el intervalo y escalarlo todo dividiéndolo entre el tiempo entre muestras y la longitud del intervalo.

Después de todo ello,  tendríamos dos vectores: uno que hemos llamado más abajo "yt", con las amplitudes de las frecuencias, y otro "f" con las frecuencias (ya desplazadas y escaladas para que sean las correctas). La nota que buscamos es la frecuencia que resuene con mayor intensidad, o lo que es lo mismo, hay que localizar el máximo valor de "yt", y ver qué valor de "f" le corresponde.

Este método no sólo es más limpio de programar y más exacto (luego veremos que da mejores resultados), sino que además es muchísimo más rápido.

El cambio que hemos realizado es básicamente el siguiente:

Script antiguo:

    % Recuento de los nodos:
    n=0;
    for i=1:(part-1)
        if (yt(i)*yt(i+1)) < 0
            n=n+1;
        else
            if (yt(i)*yt(i+1)) == 0
                n=n+1/2;
            end
        end
    end
    n=round(n);
    % Cálculo de la frecuencia original:
    freq(j)=((n+1)*fm)/(2*(fin-ini));

Script moderno:

    % FFT:
    [Fourmax,posmax]=max(abs(fftshift(fft(yt))));
    f=((1:fin-ini+1)-ceil((fin-ini+1)/2))*fm/(fin-ini+1);
    freq(j)=abs(f(posmax));

Siendo unas 6 veces más rápido el algoritmo nuevo con respecto al viejo. 

Pero, ¿es realmente más exacto este algoritmo? ¡Vamos a ponerlo a prueba! El caso más simple es el "La" puro generado por ordenador.

Detección de notas por el método de recuento de nodos

Detección de notas por el método Fourier

El resultado es parecido, pero constatamos que los tiempos de ejecución son muy distintos según qué algoritmo utilicemos.

Ejecución con el antiguo algoritmo

Ejecución con el nuevo algoritmo

Vale, pero éste era un caso muy fácil. ¿Qué pasa si tocamos un "La" con un órgano con efecto rotatorio? Aunque escuchemos el sonido modulado, la frecuencia principal viene dada por la nota ("La"), por lo que debería ser constante.

Detección de notas por el método de recuento de nodos

Detección de notas por el método Fourier

Como vemos, el método de los nodos se despista con las modulaciones mientras que el de Fourier resiste como debe.

Vamos un paso más allá: hagamos un arpegio en "Do mayor", es decir, las notas "Do", "Mi", "Sol" y "Do".

Detección de notas por el método de recuento de nodos

Detección de notas por el método Fourier

Viendo cómo reacciona el método de los nodos comprobamos que ni se acerca a identificar uno de los acordes más reconocibles que hay. El método de Fourier, sin embargo, lo hace a la perfección (e, insisto, muchísimo más rápido).

La última prueba es una escala que recorra todos los semitonos: "Do", "Do#", "Re", "Re#", "Mi", "Fa", "Fa#", "Sol", "Sol#", "La" y "Si".

Detección de notas por el método de recuento de nodos

Detección de notas por el método Fourier

Este es el único ejemplo donde falla el algoritmo de Fourier en una de las notas ("Re"), aunque todas las demás las saca exactas. El algoritmo de los nodos, sin embargo, parece mostrar notas casi al azar.

Por último, destacando la mucho mayor eficiencia del nuevo algoritmo, podemos procesar archivos mucho mayores y con un mayor número de particiones en un tiempo razonable. Valga de ejemplo el siguiente.

Se trata de una canción entera de unos dos minutos y medio grabada con una frecuencia de muestreo de 44.000Hz (6.613.200 valores), subdividida en 67 intervalos de 100.000 valores cada uno (2,27 segundos). Se tarda en procesar un tiempo similar a la duración de la canción, es decir, algo más de dos minutos. Además de ser más fiable este resultado, ésto no era procesable en un tiempo razonable con el algoritmo antiguo en un equipo como el mío.

Os dejo el código actualizado. ¡Espero que os sea útil!

function [freq]=notasFFT(p,part)
%
% Este script pretende dar una estimación de la nota que suena en un
%fragmento de audio WAVE. La sintaxtis sería:
%
% >>notasFFT('archivo.wav')
% >>notasFFT('archivo.wav',1000)
%
% Donde 'archivo.wav' es el archivo a analizar y 1000 es el ancho de
%la partición. Una partición con 1000 registros significa que de los
%miles de registros que se han muestreado a una tasa de muestreo dada
%se cogen bloques de 1000 para ser analizados. La relación entre el
%número de registros y la duración del intervalo viene dada por la
%tasa de muestreo, la cual se lee directamente del archivo WAV.
%
% Si no se especifica el tamaño de la partición (part), por defecto
%el script hace una división en 10 fragmentos.
%
% Creado por Astaroth (O.R.G.) el 17/10/2015.
% Actualizado el 21/10/2015.
%
% http://astarothsworld.blogspot.com
tic
clc


% Tabla de relaciones "Notas vs Frecuencias".
r=2^(1/12);

LA=440;
%LAs=LA*r;
%SI=LA*r^2;
%DO=LA*r^(3);
%DOs=LA*r^(4);
%RE=LA*r^(5);
%REs=LA*r^(6);
%MI=LA*r^(7);
%FA=LA*r^(8);
%FAs=LA*r^(9);
%SOL=LA*r^(10);
%SOLs=LA*r^(11);

Colores=[0,0,0;    %La
0,0,1;        %La#
0,1,0;        %Si
1,0,0;        %Do
1,0,1;        %Do#
1,1,0;        %Re
0,0,0;        %Re#
0,0,0.5;    %Mi
0,0.5,0;    %Fa
0.5,0,0;    %Fa#
0.5,0,0.5;    %Sol
0.5,0.5,0];    %Sol#
       

NOTAS=['La';'La#';'Si';'Do';'Do#';'Re';'Re#';'Mi';'Fa';'Fa#';'Sol';'Sol#';'La';'La#'];

close all

% Adquisición de datos:
[y,fm]=wavread(p);
y=y(:,1);
long=max(size(y));
Y=max(abs(y));

% Control de argumentos de entrada:
if nargin <2
    % Ancho del intervalo:
    part=round(long/10)-1;    % Ancho del intervalo.

    % Ancho del intervalo (ligado al tempo):
%    tempo=150;
%    part=fm/(150/60);
end

% Representación incial:
x=linspace(1,long,long);
subplot(2,1,1)
plot(x,y)
hold on
xlabel('Registros muestreados')
ylabel('Amplitud de la onda')
plot([0,long],[0,0],'g')
axis([0,long,(-Y-Y/4),(Y+Y/4)])

subplot(2,1,2)
hold on
xlabel('Tiempo (s)')
ylabel('Notas')
axis([0,long/fm,0,14])

% Análisis de los datos:
tram=fix(long/part);
for j=1:(tram)
    ini=(j-1)*part+1;
    fin=min(j*part,long);
    yt(1:(fin-ini+1))=y(ini:fin);
    % FFT:
    [Fourmax,posmax]=max(abs(fftshift(fft(yt))));
    f=((1:fin-ini+1)-ceil((fin-ini+1)/2))*fm/(fin-ini+1);
    freq(j)=abs(f(posmax));

    % Cálculo de una frecuencia equivalente, reducida a un intervalo concreto:
    freqred=freq(j);
    while ( freqred < LA*r )
        freqred=freqred*2;
    end
    % Incluyo la octava superior por posible superación del intervalo por redondeo:
    while ( freqred > (LA*r^(13)) )
        freqred=freqred/2;
    end
    % Cálculo del exponente (distancia en semitonos desde LA):
    rcal=log(freqred/LA)/log(r);
    rcal=round(rcal+1);
    % Solución al rebosamiento por encima:
    if rcal >= 13
        rcal=rcal-12;
    end
    % Identificación de la nota:
    Nota(j,:)=NOTAS(rcal,:);

    % Representaciones gráficas de los resultados:
    subplot(2,1,1)
    plot(ini*ones(2),linspace(-Y,Y,2),'r-')
    text(ini+(fin-ini)/2,Y+Y/8,[num2str(freq(j)),' Hz'],'HorizontalAlignment','center')
    title(['Análisis del archivo: "',p,'" - Frecuencias'])

    subplot(2,1,2)
    plot([ini/fm,ini/fm],[0,12],'r--')
    plot([ini/fm,fin/fm],[rcal,rcal],'Color',Colores(rcal,:))
    text(ini/fm+(fin-ini)/(2*fm),13,Nota(j,:),'HorizontalAlignment','center')
    title(['Análisis del archivo: "',p,'" - Notas'])
    notes=gca;
    set(notes,'YTick',1:12)
    set(notes,'YTickLabel',{'La','La#','Si','Do','Do#','Re','Re#','Mi','Fa','Fa#','Sol','Sol#'})

end
% Presentación de los resultados:
disp ('Este script intentará determinar las notas de un fragmento de audio.')
disp('')
disp('')
disp('Datos del audio original:')
disp('')
disp([' - Número de muestras: ',num2str(long),'.'])
disp([' - Frecuencia de muestreo: ',num2str(fm),' Hz.'])
disp([' - Duración del audio: ',num2str(long/fm),' s.'])
disp('')
disp('')
disp('Datos del análisis:')
disp('')
disp([' - Número de tramos: ',num2str(tram),'.'])
disp([' - Ancho del tramo: ',num2str(part),' muestras (',num2str(part/fm),' s).'])
disp(' - Análisis de los tramos:')
for j=1:(tram)
    disp(['    - Tramo ',num2str(j),': ',num2str(freq(j)),' Hz (  ',Nota(j,:),').'])
end
disp('')
disp('***Fin del programa***')
toc

Referencias:

• Práctica "Señales y Análisis de Fourier", de Belén Pérez-Verdú.
• Documentación Octave: Procesamiento de señales.
• Artículo "FFT in Octave/Matlab" de Reference Voltage.