Detector de pitch

Como pequeña parte de mi Proyecto Fin de Carrera, en el cual ya estoy sumergido y que desvelaré más adelante, he requerido programar en MatLab un detector del pitch para una secuencia de audio de entrada. Con esta herramienta pretendo modelar la secuencia de notas musicales desconocidas contenidas en un audio de entrada.

El detector de pitch programado en MatLab se basa en la idea de máxima semejanza espectral. 

Se construye en primer lugar una base de datos de espectros ideales de las distintas notas plausibles de detectar mediante convolución de una ventana espectral con un tren de impulsos que simboliza las posiciones del tono fundamental y sus armónicos.

El usuario entonces podrá controlar como parámetros el tamaño de la ventana y el solapamiento con el entorno, determinando así el tamaño de los fragmentos temporales que aproximaremos como estacionarios en el audio de entrada. A cada uno de estos fragmentos enventanados (mediante hamming por defecto) se le calcula el módulo de la FFT y se realiza un producto escalar de vectores entre dicho módulo y cada uno de los espectros ideales de la base de datos. El espectro ideal, tal que multiplicado escalarmente por el módulo de la FFT del fragmento enventanado del audio de entrada devuelva el máximo valor, será finalmente el correspondiente con la nota detectada, habiéndose minimizado así el error cuadrático:

$E = ||H(\omega)-G(\omega)||^{2} = ||H(\omega)||^{2} + ||G(\omega)||^{2}-2H(\omega)G(\omega)$

Donde $H(\omega)$ es el espectro ideal y $G(\omega)$ es el espectro real.

Se han realizado distintas pruebas y los resultados óptimos se obtienen utilizando un único armónico aparte del tono fundamental en la construcción del espectro ideal. También se han fijado otros parámetros como un umbral energético para establecer si existe nota musical o el fragmento se corresponde con silencio.

El código para MatLab se encuentra a continuación. El éxito en las pruebas realizadas es del 100% y su funcionamiento óptimo. La ayuda se encuentra detallada en la propia función.

Nota: He limitado la detección de notas en el rango Sol2 - Sol5 (3 octavas), correspondiente a mi registro vocal (ya se explicará el porqué más adelante).

% DETECTOR DE PITCH
%
% Autor: Iván López Espejo
%
% y = detectorPitch(x,wsize,sol,fm)
%
% Donde 'x' es el audio mono a analizar de entrada, 'wsize' el tamaño de la
% ventana o fragmentos en los que se segmentará el audio de entrada (valor
% típico en el rango 50ms-250ms) dado en segundos, 'sol' el solapamiento
% entre fragmentos adyacentes también dado en segundos (valor típico en el
% rango 10ms-50ms) y 'fm' la frecuencia de muestreo del audio de entrada
% dada en Hz's.

function y = detectorPitch(x,wsize,sol,fm)

% Lo damos en término de muestras.
wsize = round( wsize * fm );
sol = round( sol * fm );

% Obtención de los espectros ideales de referencia.
% k = floor((fm/(2*freq)) - 1);
k = 1;
l = 1;
h = [0.2 0.55 0.8 0.95 1 0.95 0.8 0.55 0.2];
o = 2;
for n = 8:12
        freq = 440*(2^((o-4) + ((n-10)/12)));
        vec = zeros(1,fm/2);
        for arm = 1:(k+1)
            vec(round(freq*arm)) = 1;
        end
        vec = conv(vec,h);
        vec = vec(5:length(vec)-4);
        sideal(l,:) = vec;
        l = l + 1;
end
for o = 3:4
    for n = 1:12
        freq = 440*(2^((o-4) + ((n-10)/12)));
        vec = zeros(1,fm/2);
        for arm = 1:(k+1)
            vec(round(freq*arm)) = 1;
        end
        vec = conv(vec,h);
        vec = vec(5:length(vec)-4);
        sideal(l,:) = vec;
        l = l + 1;
    end
end
o = 5;
for n = 1:8
        freq = 440*(2^((o-4) + ((n-10)/12)));
        vec = zeros(1,fm/2);
        for arm = 1:(k+1)
            vec(round(freq*arm)) = 1;
        end
        vec = conv(vec,h);
        vec = vec(5:length(vec)-4);
        sideal(l,:) = vec;
        l = l + 1;
end

% Calculamos el espectro de cada trozo enventanando según hamming.
% Definimos las constantes de la ventana.
a0 = 0.53836;
a1 = 0.46164;
% Umbral de detección de nota.
umbral = 2;
l = 1;
for j = 1:(wsize-sol-1):length(x)
    if wsize > length(x) - j
        N = length(x) - j + 1;
    else
        N = wsize;
    end
    n = 1:N;
    vent = a0 - a1*cos(2*pi*n/(N-1));
    trozo = x(j:j+N-1).*vent';
    % Calculamos la FFT del trozo.
    spc = abs(fft(trozo,fm));
    spc = spc(1:fm/2);
    spc = spc/max(spc);
    for n = 1:1:length(sideal(:,1))
        vecAux(n) = sum(spc'.*sideal(n,:));
    end
    [maxim,pos] = max(vecAux);
    y(l) = 0;
    if maxim > umbral
        y(l) = pos;
    end
    l = l + 1;
end

% Bloque para devolver la secuencia de notas.
ind = 1;
for j = 1:length(y)-1
    if y(j) ~= y(j+1)
       y2(ind) = y(j);
       ind = ind + 1;
    end
end
y2(ind) = y(length(y));
vNotas = ['Sol2 ';'Sol#2';'La2  ';'La#2 ';'Si2  ';'Do3  ';'Do#3 ';'Re3  ';'Re#3 ';'Mi3  ';'Fa3  ';'Fa#3 ';'Sol3 ';'Sol#3';'La3  ';'La#3 ';'Si3  ';'Do4  ';'Do#4 ';'Re4  ';'Re#4 ';'Mi4  ';'Fa4  ';'Fa#4 ';'Sol4 ';'Sol#4';'La4  ';'La#4 ';'Si4  ';'Do5  ';'Do#5 ';'Re5  ';'Re#5 ';'Mi5  ';'Fa5  ';'Fa#5 ';'Sol5 '];
ind = 1;
for j = 1:length(y2)
    if y2(j) ~= 0
        y3(ind,:) = vNotas(y2(j),:);
        ind = ind + 1;
    end
end
% Imprimimos la secuencia de notas detectada por pantalla.
y3