Sistema para la Afinación Automática de Sonidos Polifónicos (III)

Un ejemplo de una burda implementación del sistema en MatLab. El primer compendio de código responde al texto central que hay que ejecutar y que usa otras funciones programadas que adjunto justo a continuación.

% A³ v0.1 : Afinador Automático de Acordes.
% Autor: Iván López Espejo

clear;
% Generación de un acorde desafinado con armónicos.
chord = [];
% Definición de frecuencias fundamentales y de la frecuencia de muestreo.
frec1 = 435;
frec2 = 520;
frec3 = 650;
frecM = 44100;
for j = 1:1:frecM
    chord(j) = cos(2*pi*frec1*j/frecM);
    chord(j) = chord(j) + cos(2*pi*frec2*j/frecM);
    chord(j) = chord(j) + cos(2*pi*frec3*j/frecM);
    % Decremento exponencial de la amplitud de los armónicos.
    for k = 1:1:20
        chord(j) = chord(j) + exp(-k)*cos(2*pi*(k+1)*frec1*j/frecM);
        chord(j) = chord(j) + exp(-k)*cos(2*pi*(k+1)*frec2*j/frecM);
        chord(j) = chord(j) + exp(-k)*cos(2*pi*(k+1)*frec3*j/frecM);  
    end
end
% Calculamos los tonos fundamentales que componen el sonido polifónico.
spectrum = fft(chord);
vF = estTon(chord);
% Filtramos el acorde de forma que obtenemos tantos espectros como tonos
% fundamentales haya y afinamos por separado cada uno de ellos.
sonidoFinal = zeros(1,length(chord));
for k = 1:1:length(vF)
    % Espectro filtrado con el k-ésimo tono y sus armónicos.
    H = chordFilter(spectrum,vF(k));
    % Estimamos hacia qué frecuencia debemos afinar con la función lineal
    % clasificadora.
    num = clasif(vF(k));
    % Afinación del espectro con establecimiento de la variable de control
    % a 0 por proporción del espectro.
    notaAf = afin(H,num,vF(k),0);
    % Agregamos los sonidos resultantes afinados en el dominio del tiempo.
    sonidoFinal = sonidoFinal + notaAf;
end
% Finalmente reproducimos los sonidos polifónicos afinado y desafinado en el dominio del
% tiempo.
soundsc(chord,frecM)
soundsc(sonidoFinal,frecM)

% Función que se encarga de afinar un determinado conjunto de todo y
% armónicos mediante desplazamiento del espectro.
% Autor: Iván López Espejo.

function y = afin(nota,num,fA,ToF)

% Cálculo de la octava a la que corresponde la frecuencia dada.
octava = floor(log2(fA/32.7) + 1);
% Cálculo de la frecuencia dado el número de nota.
fN = 440*exp(((octava - 4) + ((num-10)/12))*log(2));
% Parámetro alpha de escalado del tiempo.
alpha = fN/fA;
% Espectro de la nota. Si ToF (variable de control) está a 1, al comienzo
% se proporciona el sonido en el dominio del tiempo, si está a 0, en el
% dominio de la frecuencia.
if ToF == 1
    spec = fft(nota);
else
    spec = nota;
end
% Obtención del espectro afinado haciendo uso de la propiedad de escalado
% temporal de la transformada de Fourier.
specAf = [];
for k = 1:1:length(spec)
    if round(k/alpha) > 0 && round(k/alpha) < length(spec) + 1
        specAf(k) = spec(round(k/alpha))*abs(1/alpha);
    end
end
specAf = [specAf(1:length(spec)/2) specAf(length(spec)/2:-1:1)];
% Devolución del resultado.
trans = ifft(specAf);
sonFin = abs(trans);
for k = 1:1:length(trans);
    signo = real(trans(k))/abs(real(trans(k)));
    sonFin(k) = sonFin(k)*signo;
end
y = sonFin;

% Función que filtra del espectro todo sonido a excepción del tono
% fundamental y armónicos asociados de interés.
% Autor: Iván López Espejo.

function y = chordFilter(spec,freq)

% Generación del filtro en el dominio del tiempo como la señal a la
% frecuencia fundamental con una cierta dispersión y sus armónicos
% audibles.
tFilter = [];
frecM = 44100;
for j = 1:1:frecM
    tFilter(j) = cos(2*pi*freq*j/frecM);
    tFilter(j) = tFilter(j) + cos(2*pi*(freq-1)*j/frecM);
    tFilter(j) = tFilter(j) + cos(2*pi*(freq+1)*j/frecM);
    for k = 2:1:ceil(22050/freq)-1
        tFilter(j) = tFilter(j) + cos(2*pi*k*freq*j/frecM);
        tFilter(j) = tFilter(j) + cos(2*pi*k*(freq-1)*j/frecM);
        tFilter(j) = tFilter(j) + cos(2*pi*k*(freq+1)*j/frecM);
    end
end
% Filtro en el dominio de la frecuencia (normalizado).
fFilter = fft(tFilter);
fFilter = fFilter/max(fFilter);
% Aplicación del filtro sobre el espectro dado.
y = fFilter.*spec;

% Decisor lineal a partir de la transformación logarítmica de la escala
% musical exponencial.
% Autor: Iván López Espejo.

function y = clasif(freq)

% Cálculo de la octava a partir de la frecuencia.
octava = floor(log2(freq/32.7) + 1);
% Cálculo de la nota dentro de la octava a partir de la frecuencia (Do = 1,
% Do# = 2, ... , Si = 12).
n = 10 + 12*(log(freq/440) - log(2)*(octava - 4))/log(2);
if n == round(n)
    disp('La nota está perfectamente afinada.')
    y = 0;
else
    if round(n) == 13
        y = 1;
    else
        if round(n) == 0
            y = 12;
        else
            y = round(n);
        end
    end
end

% Función auxiliar para la estimación de los tonos fundamentales
% componentes del sonido polifónico.
% Autor: Iván López Espejo.

function y = estTon(x)

% Cálculo del espectro del sonido polifónico dado en el dominio del tiempo.
spectrum = abs(fft(x));
% Nos quedamos con la representación del espectro entre 0 y pi.
spectrum = spectrum(1:length(spectrum)/2);
% Obtenemos los 10 máximos del espectro.
vM = [];
for k = 1:1:10
    [maxim,pos] = max(spectrum);
    vM(k) = pos;
    spectrum(pos) = min(spectrum);
end
% Clasificamos en tonos fundamentales o desechamos armónicos.
vF(1) = vM(1);
% Tolerancia: valor de desviación máximo entre el natural ideal que
% multiplica la frecuencia fundamental y el real.
epsilon = 0.05;
mark = 0;
for k = 2:1:length(vM)
    for l = 1:1:k-1
        dv = vM(k)/vM(l);
        if abs(round(dv) - dv) > epsilon
            mark = mark + 1;
        end
    end
    if mark == k-1
        vF = [vF vM(k)];
    end
    mark = 0;
end
% Devolución del vector de tonos fundamentales detectados.
y = vF;