I² Blog

API de procesado de audio para JAVA

2010-12-22T22:01:00.003+01:00

A la hora de realizar un prototipo en JAVA sobre un proyecto que me encuentro realizando, me he encontrado con que no hay mucho material disponible escrito para dicho lenguaje de fácil uso para el procesado de señal. Por ello, yo mismo he tenido que programar métodos para realizar funciones que necesitaba de procesado. Mi proyecto versa sobre el tratamiento del audio por lo que, básicamente, el código escrito está orientado a ello. Muchas de las funciones son muy concretas, pues no quise complicarme escribiendo algo muy genérico pero, otras tantas, sí se pueden usar de forma genérica (también está la opción de modificar según la necesidad).

Podéis descargar la API como fichero con extensión .java haciendo clic aquí.

A continuación presento una descripción breve sobre las funciones incluidas, ya que la clase escrita no está documentada.

- Suavizado de Pitch

La función admite como entrada un vector de números reales en doble precisión representantes de un contorno de frecuencia de pitch en el tiempo (también podría usarse con otros propósitos). Como salida devuelve el mismo contorno suavizado. Este suavizado se basa en umbrales heurísticos supuesto que el contorno de pitch es el logaritmo en base e de la secuencia de frecuencias de pitch en hertzios. La realización de este método se encuentra motivada por el hecho de que, en muchas ocasiones, durante la detección del pitch a partir del espectro o funciones de autocorrelación, existen armónicos más energéticos que el tono fundamental, lo que se traduce en picos erróneos dentro de la secuencia. Se considera errónea una ráfaga de valores cuando dicha ráfaga es inferior a 10 cifras consecutivas y la pendiente generada en su entorno es superior a 0.5.

- Cálculo de la PSD

Calcula la PSD (módulo de la FFT pero no en dB) entre $0$ y $2\pi$ de una secuencia de muestras dada como un vector de números reales en doble precisión. La secuencia de salida tiene la misma longitud y se calcula según:

$Y(k) = \left | \sum_{n=0}^{N-1}x(n)e^{-i\frac{2\pi nk}{N}} \right |$

- Filtro de Mediana de Orden 5

Muy útil para la eliminación de ruido aditivo impulsivo. Realmente puede modificarse en su interior el orden, siempre y cuando este sea impar. La muestra k-ésima filtrada se obtiene a partir del cálculo de la mediana como:

$y(k) = mediana(x(k - (n-1)/2),...,x(k + (n+1)/2))$

- Cuantizador Musical

Aproxima una secuencia de valores de pitch frecuencial dados en forma logarítmica a semitonos de la escala musical. Esto sólo se encuentra habilitado en el rango frecuencial que va del Sol2 al Sol5.

- Cálculo de la Energía

Calcula la energía de una secuencia de muestras dadas en un vector de datos reales de doble precisión. Dicha energía se obtiene de la forma:

$E(dB) = 10log_{10}\sum_{n=0}^{N-1}x(n)^{2}$

- Interpolador de Orden 5

Este interpolador aumenta la resolución en un factor 5 de un conjunto de muestras dadas en el formato ya mentado. Este orden de interpolación es inamovible ya que el funcionamiento del método se basa en la inserción de ceros y en la aplicación de un filtro FIR simétrico especial (cuyos coeficientes sólo son válidos para este orden de interpolación).

- Enventanador de Hamming

Enventana una secuencia de muestras con Hamming. La función de la ventana de Hamming es:

$w(n) = a_{0}-a_{1}\cos(\frac{2\pi n}{N-1})$

- Cálculo de la Media

Calcula la media de un vector de datos de entrada como:

$\mu = \frac{1}{N}\sum_{n=0}^{N-1}x(n)$

- Cálculo de la Varianza

Calcula la varianza de un vector de datos de entrada como:

$\sigma^2 = \frac{1}{N-1}\sum_{n=0}^{N-1}(x(n)-\mu)^2$

- Normalizador

Normaliza un vector de muestras de entrada por el máximo del vector en valor absoluto.

- Filtrado

Método para el filtrado de un vector de datos reales en doble precisión de entrada. Devuelve otro vector del mismo tamaño al que se le ha aplicado un filtrado dados los coeficientes del correspondiente filtro FIR o IIR. Dado un filtro de orden (p,q), la salida n-ésima filtrada se obtiene como:

$y(n) = b_{0}x(n)+...+b_{q}x(n-q) - a_{1}y(n-1)-...-a_{p}y(n-p)$

- Cálculo de la Autocorrelación

Dada una secuencia de muestras de entrada, este método calcula su autocorrelación haciendo uso del estimador sesgado, el cual es de uso preferente debido a su inferior varianza, lo que se traduce en una mejor convergencia de las muestras de autocorrelación lejanas. Se devuelve un vector de salida de la misma dimensión que el de entrada, el cual contiene las muestras correspondientes a la parte positiva del eje de abscisas. El estimador sesgado de la autocorrelación para la parte positiva del eje de abscisas se define como:

$r_{x}(k) = \frac{1}{N}\sum_{n=0}^{N-1-k}x(n)x(n+k)$

- Análisis LPC

Se encarga de obtener hasta un total de 13 coeficientes del predictor lineal que modela, como un filtro todo-polos, la respuesta del tracto vocal. Estos coeficientes se obtienen a partir de la resolución de las ecuaciones de Yule-Walker a partir del cálculo de muestras de autocorrelación dado un vector de muestras reales en doble precisión de entrada.

- Detector de Pico Máximo

Esta función detecta la presencia de un pico máximo absoluto en una secuencia de valores de entrada. Es la misma función que ya presenté en otra entrada escrita en MatLab.

- Remuestreador de Orden 5

No sólo es un remuestreador de orden 5 sino que, únicamente, funciona como remuestreador para pasar de una frecuencia de 8kHz a 1.6kHz. Esto es así ya que, al principio, se aplica un filtrado FIR de fase lineal cuyos coeficientes son estáticos y almacenados para este propósito. Una vez filtrado el vector de entrada, se selecciona una de cada 5 muestras.

Detector de pitch (II)

2010-11-26T17:36:00.001+01:00

Vuelvo a presentar un nuevo detector de pitch, esta vez, basado en autocorrelación. Dicho detector está programado en MatLab y está basado en el algoritmo SIFT. Su diagrama de bloques se observa a continuación.

Primero, el audio de entrada se preprocesa remuestreando desde 44100Hz (asumimos dicha frecuencia de muestreo para el audio de entrada) a 8000Hz, habiendo ponderado previamente ambos canales estéreo para, finalmente, normalizar la señal.

Seguidamente filtramos la señal mediante un filtro tipo Butterworth paso-bajo con frecuencia de corte 800Hz (correspondiente, aproximadamente, a un Sol-5 en la escala musical) para, a continuación, aplicar un downsampling correspondiente a la nueva frecuencia máxima presente en la señal, lo que se traduce, según el criterio de Nyquist, en un remuestreo a 1600Hz. Fragmentamos en tramas el audio filtrado y, para cada trama, aplicamos el proceso descrito a partir de ahora. Se toma el frame y con él se realiza un estudio LPC (Linear Predictive Coding), con el fin de modelar la respuesta espectral referente al tracto vocal (formantes) mediante un filtro todo-polos. Una vez obtenidos los coeficientes de dicho filtro, se filtra la trama mediante otro filtro FIR cuyos coeficientes son los anteriores del modelo LPC. Con ello, intentamos eliminar la información propia del tracto vocal de la señal con el fin de mantener sólo la información de interés: el pitch. Por último se calcula la función autocorrelación de la señal filtrada y se interpola para ganar resolución temporal. El pico máximo detectado a la salida, se corresponderá con el período de pitch deseado.

La función que implementa el anterior proceso toma como datos de entrada el tamaño del frame (en segundos, con un valor típico de 40ms) y el solapamiento (también en segundos, con un valor típico de 10ms) aparte de, claro está, el propio audio cuyo pitch deseamos estimar.

También se incluye un detector de actividad de voz basado en hangover con los parámetros por defecto fijados en la función (cambiar si fuera necesario). La salida es la secuencia de pitch, a la que se le ha aplicado el logaritmo neperiano (si se desea obtener la frecuencia en Hz's, cambiar desde el código o aplicar la exponencial al vector de salida), de tal forma que no se incluyen las estimaciones para las tramas cuyo flag proporcionado por el detector de voz se corresponda con silencio.

La gráfica adjunta muestra un ejemplo del pitch estimado para la canción Hush de Deep Purple tarareada.

IMPORTANTE: La función detPico.m se encuentra en la entrada justamente anterior.

% Detector de pitch
% Autor: Iván López Espejo

function pitch = detPitch(x,wsize,sol)

% Preprocesado.
fm1 = 44100;
fm = 8000;
x = (x(:,1) + x(:,2))/2;
x = resample(x,fm,fm1)
x = x/max(abs(x));

% Parámetros del detector de actividad de voz.
THRESHOLD = 0.01;
MIN_SPEECH_FRAME_HANGOVER = 4;
HANGOVER = 1;
meanEn = 0;
nbSpeechFrame = 0;
ind = 1;

% Diseño de filtro paso-bajo Butterworth con frecuencia de corte a 800Hz.
fc = 800;
Nf = 2; % Orden del filtro.
Wn = 2*fc/fm;
[b,a] = butter(Nf,Wn,'low');
x = filter(b,a,x); % Filtramos paso-bajo.
x = resample(x,2*fc,fm); % Downsampling.
% Análisis por frame.
wsize = round(wsize * 2 * fc);
sol = round(sol * 2 * fc);
% Definimos las constantes de la ventana.
a0 = 0.53836;
a1 = 0.46164;
Npred = 12; % Número de coeficientes del predictor lineal.
b0 = 1;
fint = 5; % Factor de interpolación.
pitch = [];
for m = 1:(wsize-sol-1):length(x)
    if wsize > length(x) - m
        N = length(x) - m + 1;
    else
        N = wsize;
    end
    % Detección de la actividad de voz.
    frameEn = dot(x(m:m+N-1),x(m:m+N-1))/N;
    if frameEn - meanEn > THRESHOLD
        VAD(ind) = 1;
        nbSpeechFrame = nbSpeechFrame + 1;
        if nbSpeechFrame > MIN_SPEECH_FRAME_HANGOVER
            hangOver = HANGOVER;
        end
    else
        if hangOver == 0
            VAD(ind) = 0;
        else
            VAD(ind) = 1;
            hangOver = hangOver - 1;
        end
    end
    ind = ind + 1; % Incremento del contador.
    % Enventanado del frame mediante Hamming.
    n = 1:N;
    win = a0 - a1*cos(2*pi*n/(N-1));
    frame = x(m:m+N-1).*win';
    % Método de la autocorrelación para modelado de los coeficientes del
    % predictor lineal asociado a la trama.
    acoef = lpc(frame,Npred);
    % Usamos los anteriores coeficientes como coeficientes de un filtro
    % FIR.
    eps = filter(acoef,b0,frame);
    reps = xcorr(eps,'biased'); % Cómputo de la autocorrelación.
    % Interpolación para la obtención de mayor resolución temporal.
    Pm = interp(reps,fint);
    Pm = Pm(floor(length(Pm)/2):length(Pm));
    % Detectamos el pico máximo en la función de autocorrelación.
    picomax = detPico(Pm);
    pause();
    if VAD(ind-1) == 1
        pitch = [pitch (log(2*fc*fint/(picomax-1)))];
    end
end
% Representamos el resultado.
plot(pitch,'-*r')
grid on
xlabel('Evolución temporal')
ylabel('Frecuencia de pitch logarítmica')
title('PITCH estimado')

Detector de actividad de voz (VAD)

2010-11-01T15:54:00.008+01:00

En esta entrada se recoge la implementación de un detector de actividad de voz a través de una red neuronal artificial entrenada mediante propagación hacia atrás.

Este tipo de redes neuronales es ideal para el reconocimiento de patrones, de modo que entrenamos el sistema con vectores de características representativos de pasajes de silencio (ruido de fondo) y pasajes de voz. Para ello caracterizamos una trama de audio mediante 12 coeficientes cepstrales MFCC (los primeros 13 coeficientes salvo el primero, ya que este está relacionado con la energía de la trama y en este caso no es representativo) y un coeficiente de entropía, computada como:

$H = -\sum_{k=1}^N p_{k}\log p_{k}$

Donde $N$ es el número de muestras del módulo del espectro de la trama en consideración y $p_{k}$ representa el valor k-ésimo de la función densidad de probabilidad calculada a partir de dicho espectro, definiéndose como:

$p_{k} = \frac{S_{k}}{\sum_{j}^N S_{j}}$

Con $S$ el módulo del espectro de la trama.

La red neuronal que he construido para solventar el problema utiliza tramas de 128 muestras (donde consideramos la señal aproximadamente estacionaria) para un audio muestreado a 8kHz y está diseñada para usar 10 capas ocultas. Dicha red está entrenada con 20 vectores de características representativos de silencio y otros 22 vectores representativos de voz. El resto de parámetros se extrae del código MatLab incluido a continuación.

Obsérvese un ejemplo de un fragmento de audio al que se le eliminan los silencios mediante esta implementación de VAD. El funcionamiento es bastante bueno. El programa va consumiendo tramas de 128 muestras a 8kHz, clasificándolas como voz o silencio gracias a la red neuronal previamente entrenada. Recomiendo sacar del código la parte del entrenamiento de la red neuronal, a fin de no estar inicializándola constantemente cada vez que se llame a la función.

Audio de entrada.

Audio de salida con silencios eliminados.

Esta herramienta que he construido está pensada para funcionar previamente a la detección del pitch, a fin de poder eliminar las partes sin voz para no obtener falsas estimaciones.

Se puede ver el código a continuación, así como descargar mi fichero con patrones de entrenamiento haciendo clic aquí.

% Detector de actividad de voz basado en red neuronal
% Autor: Iván López Espejo.

function [y,sil] = vadet(x,fm)

% Remuestreo a 8kHz.
x = resample(x,8000,fm);
x = x/max(abs(x));
% Definición de algunos parámetros.
N = 128;
% Inicialización y entrenamiento de la red neuronal.
nump = 2; % Número de clases.
% Patrones de entrenamiento.
load Patrones
P = Patrones;
% Salidas.
T = zeros(2,42);
T(1,21:42) = ones(1,22);
T(2,1:20) = ones(1,20);
S1 = 10; % Número de capas ocultas.
S2 = nump; % Número de capas de salida (= al número de clases).
[R,Q] = size(P);
epochs = 10000; % Número de iteraciones.
goal_err = 5e-5; % Error objetivo.
a = 0.7; % Define el rango de variables aleatorias.
b = -0.7;
W1 = a + (b-a) * rand(S1,R); % Pesos entre la entrada y las neuronas ocultas.
W2 = a + (b-a) * rand(S2,S1); % Pesos entre las neuronas ocultas y las de salida.
b1 = a + (b-a) * rand(S1,1); % Pesos entre la entrada y las neuronas ocultas.
b2 = a + (b-a) * rand(S2,1); % Pesos entre las neuronas ocultas y las de salida.
n1 = W1 * P;
A1 = logsig(n1);
n2 = W2 * A1;
A2 = logsig(n2);
e = A2 - T;
error = 0.5 * mean(mean(e.*e));
nntwarn off
for itr = 1:epochs
    if error <= goal_err
        break;
    else
         for i = 1:Q
            df1 = dlogsig(n1,A1(:,i));
            df2 = dlogsig(n2,A2(:,i));
            s2 = -2 * diag(df2) * e(:,i);
            s1 = diag(df1)* W2' * s2;
            W2 = W2 - 0.1 * s2 * A1(:,i)';
            b2 = b2 - 0.1 * s2;
            W1 = W1 - 0.1 * s1 * P(:,i)';
            b1 = b1 - 0.1 * s1;

            A1(:,i) = logsig(W1*P(:,i),b1);
            A2(:,i) = logsig(W2*A1(:,i),b2);
         end
            e = T - A2;
            error = 0.5 * mean(mean(e.*e));
            disp(sprintf('Iteración: %5d mse: %12.6f%',itr,error));
            mse(itr) = error;
    end
end
threshold = 0.9; % Umbral del sistema (umbral más alto = mayor precisión).
TrnOutput = real(A2>threshold);
disp('Resultado del entrenamiento: ')
disp(TrnOutput)
% Inicialización del audio de salida.
y = [];
sil = [];
for j = 1:N:length(x)-N
   coefs = calCoeficientes(x(j:j+N-1))';
   % Clasificamos el frame en voz o silencio.
   n1 = W1 * coefs;
   A1 = logsig(n1);
   n2 = W2 * A1;
   A2test = logsig(n2);
   resultado = real(A2test>threshold);
   if resultado(1) == 1 && resultado(2) == 0
       y = [y; x(j:j+N-1)];
   else
       sil = [sil; x(j:j+N-1)];
   end
end

Detector de pitch

2010-10-08T21:31:00.013+02:00

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

Corrección de ojos rojos

2010-06-14T14:03:00.000+02:00

A veces, cuando realizamos una fotografía con flash, podemos apreciar el indeseado efecto de los ojos rojos. Esto es debido a que la pupila se encuentra, normalmente, muy dilatada en el momento de realizar la fotografía con flash (ya que esta tiene lugar en un ambiente oscuro y, por tanto, la pupila ha de captar mayor luminosidad para integrar una imagen del entorno correctamente). La luz del flash entonces entra a través de la pupila y se refleja produciendo ese característico color rojo propio de los capilares sanguíneos situados en dicha zona.

Flashes modernos, previo a la realización de la fotografía, emiten un primer destello cuya finalidad es que la pupila se contraiga. Con esta misma premisa, un truco casero previo a la realización de la foto, consistiría en someterse a una fuente intensa de luz durante unos instantes a fin de conseguir la contracción de la pupila.

Si no podemos evitar el efecto, algunas cámaras digitales y programas de retoque fotográfico incluyen alguna funcionalidad para mitigarlo. De esto se trata la presente entrada, en la que desarrollamos un sistema semi-automático para intentar corregir el efecto de los ojos rojos.

A continuación detallamos los pasos de la solución particular propuesta.

Esta solución se encuentra inspirada someramente por el artículo de Jalal acerca de la detección de ojos. Observé que, ante el color rojo intenso del efecto producido por el flash sobre los ojos, las componentes cromáticas, según el modelo YCbCr, presentaban unas características bien diferenciadas y es que, los ojos en la componente Cb tendían al negro y en la componente Cr al blanco. A partir de esta idea, tratamos de aislar y localizar los ojos en la imagen. Asistamos paso a paso la solución mediante un ejemplo.

Consideremos la siguiente fotografía de partida:

Transformando al modelo de color YCbCr, obtenemos que los canales Cb y Cr son, respectivamente:

Si invertimos el canal Cr, los ojos quedarán también bastante más oscuros que el resto de la imagen, obteniendo:

Ahora, si sumamos Cb + Cr’, conseguiremos un contraste tal que los ojos serán mucho más oscuros que el resto de la imagen, de tal modo que obtenemos:

Aplicando una función de stretching sobre el anterior resultado conseguimos un mayor contraste entre los ojos y el fondo, obteniendo lo siguiente:

Justamente en este punto entra en juego el primer parámetro de caracterización de la función, y este es el nivel de umbralización de la anterior imagen con stretching. Para un valor umbral de 128, obtenemos:

Este parámetro es crítico y su nivel óptimo impredecible, pues es dependiente de la característica aleatoria de cada imagen. Si el valor es demasiado bajo, se tiende a perder potencial información sobre ojos en la imagen y, si es demasiado alto, pueden considerarse otras partes de la imagen como falsos ojos rojos, siendo ambas situaciones indeseables (podría aplicarse un algoritmo de umbralización óptimo como es el método de Otsu para la estimación del umbral tal que se maximice la varianza inter-clase). En la mayoría de los casos, un valor intermedio en un rango con bastante tolerancia debe ser más que suficiente. En el caso que nos ocupa, como se observa, el valor intermedio comentado de 128 (umbral standard) nos proporciona un resultado bastante bueno, pues ya tenemos delimitados nuestros ojos rojos.

Lo único que nos queda ya por hacer es determinar su posición dentro de la imagen. Para tal fin, usaremos proyecciones de histogramas. Para mayor comodidad de análisis con los mismos, invertiremos de nuevo la imagen anterior, a fin de observar los ojos en blanco y quedar estos caracterizados en las proyecciones de histograma como lóbulos de un cierto ancho:

A continuación calculamos las proyecciones de histograma, obteniéndose para las proyecciones horizontal y vertical, respectivamente:

Se observa en ambas imágenes dos picos, los cuales caracterizan las coordenadas del punto central de cada uno de los ojos. Para no llevar a equívocos, pues ahora precisamos determinar la posición de los picos y el ancho de los lóbulos, suavizamos las proyecciones de histograma mediante una función gaussiana 1-D, en la cual entra en juego el segundo parámetro de la rutina, y este es la varianza de la función. Conseguiremos un mayor suavizado para una mayor varianza (la longitud de la función gaussiana con la que convolucionamos las proyecciones de histograma es constante y es de 11 muestras). Un valor que funciona bien es σ^2 = 16. De hecho, este parámetro es potencialmente extraíble de la función, ya que funciona bien para el valor constante antes comentado, aunque, no obstante, he decidido dejarlo. A continuación se muestran las proyecciones de histograma suavizadas:

Como es normal, si usamos un valor bajo de varianza, obtendremos un resultado más abrupto, lo cual, nos puede llevar a equívocos a la hora de determinar el número y la posición de los ojos. Por el contrario, si el valor de varianza es excesivamente alto (o lo fuera el del número de muestras de la función gaussiana con la que convolucionamos), podríamos perder información acerca de los picos y, por tanto, del número y posición de los ojos.

Mediante un algoritmo de detección de pendiente, conseguimos contar el número de picos sobre las proyecciones de histograma (lo cual nos revela información acerca del número de ojos y las coordenadas de sus centros). Mediante otro algoritmo de detección de pendiente, conseguimos calcular el número de valles a fin de caracterizar el ancho de los lóbulos (esto nos da información acerca del radio del ojo, tanto horizontal, como vertical). Para el ejemplo que nos ocupa obtenemos los siguientes datos mediante el algoritmo que implementé en MatLab:

Posición de picos Posición de valles

Proyección Horizontal 83, 103 71, 93 y 115

Proyección Vertical 114, 242 102, 127, 230 y 255

Pero, ahora surge una duda. Normalmente, el número de picos de la proyección horizontal es igual al de la proyección vertical y, estos son iguales al número de ojos rojos en la fotografía. Sin embargo, tenemos N^2 posibles ojos rojos debido a las combinaciones posibles entre picos de ambas proyecciones (en nuestro caso particular, hasta 4, donde realmente hay 2). Para solventar esta ambigüedad, sumamos todos los niveles de gris de cada uno de los rectángulos cuyos lados quedan definidos por la posición de los valles y cuyos centros por las distintas combinaciones de los picos en la última imagen con la región de los ojos en blanco. Si la combinación de picos no se corresponde con un ojo, la suma de este rectángulo será eminentemente 0. Sin embargo, si se corresponde, la suma será un valor elevado. El punto de corte para considerar si una combinación de picos determina o no el centro geométrico de un ojo, viene definido porque, al menos, el 25% de píxels del rectángulo que comprende al hipotético ojo, sea blanco. Si es así, en la imagen RGB se pinta una cruz con centro el centro geométrico del ojo y cuya longitud de brazos viene determinada por el ancho de los lóbulos (es decir, por la posición de los valles).

Sabido esto, para la anterior imagen obtenemos:

Una vez detectada la posición de los ojos rojos, únicamente habremos de corregir su tonalidad. Aquí entra en juego el último de los parámetros. Este parámetro se basa también en establecer un umbral para la imagen Cb + Cr’, de forma que se analiza cada uno de los rectángulos que comprende aun ojo sobre dicha imagen según los valores antes obtenidos de picos y valles. Si el nivel de gris i-ésimo en el rectángulo que comprende a cada uno de los ojos es menor que el umbral, se considera ojo rojo y se castiga a los niveles RGB con factores 0.10, 0.12 y 0.12, respectivamente (es decir, se tiende a llevar al negro, que sería lo correcto para el color de la pupila). Si es superior el nivel de gris i-ésimo, se entiende que no forma parte del ojo rojo y, por tanto, se deja intacto.

Este parámetro es muy importante y participa, junto con el primero, de forma directa en la calidad final de la imagen retocada. Si el umbral es muy bajo, corremos el riesgo de no corregir el color rojo de los ojos. Sin embargo, si es muy alto, tenderemos a llevar al negro todo el rectángulo que abarca al ojo.

El resultado final para un umbral de corrección de 200, resulta ser:

La forma de utilizar la función implementada en MatLab es la siguiente y se puede descargar haciendo clic aquí:

y = correctorOjosRojos(imagen,nUmb,var,cte)

Donde en y se obtiene la matriz imagen corregida, siendo imagen la imagen de entrada con el defecto de los ojos rojos, nUmb el umbral para binarización aplicado tras el primer stretching, var la varianza de la función gaussiana de suavizado y cte el punto de corte para consideración de píxel rojo o no durante la corrección final.

Reconocimiento automático de partituras

2010-06-14T00:32:00.002+02:00

Mi amigo Jonathan y yo hemos programado un pequeño código que intenta ser un reconocedor automático de partituras basado en el procesamiento de la imagen (algo análogo al OCR para el reconocimiento de caracteres).

La idea es que una partitura escaneada pueda ser digitalizada para ser interpretada y procesada su información por un computador (por ejemplo, para conseguir una reproducción automática de la misma).

Sólo reconoce partituras sencillas, funcionando del siguiente modo (caracterizamos cada uno de los bloques funcionales):

Preprocesamiento (OMR.m)

En este primer apartado adaptamos la partitura para su tratamiento. Una vez obtenemos la imagen de la misma, supuesta su entrada con el modelo de color RGB, la transformamos a niveles de gris para seguidamente invertirla. La misión de la inversión es la de un tratamiento más sencillo mediante la contemplación de la información contenida en la partitura como niveles de gris relativamente grandes frente al fondo tendente a nivel cero. Finalmente aplicamos binarización para que los píxels de interés se mantengan en blanco frente al fondo negro.

Corrección de la inclinación (corregirInclinacion.m)

Esta etapa está incluida en el preprocesado de la imagen y como su nombre indica consiste en corregir la inclinación que pueda presentar la imagen de partida. Este paso es de vital importancia dado que los métodos de segmentación y reconocimiento desarrollados no operarían correctamente ante la presencia de este defecto.

El método seguido para llevar a cabo la corrección de la inclinación se basa principalmente en la transformada de Hough. En primer lugar se aplica la transformada de Hough con una resolución en la componente theta de 0.05º. El espacio de parámetros de salida usado es el polar. La ventaja que aporta (aparte de los problemas que nos evitamos con las líneas verticales) es que obtenemos directamente una aproximación del ángulo de inclinación de la imagen. El porqué de optar por el uso de la transformada de Hough es por las características inherentes de la imagen que vamos a emplear. Dado que estamos ante la imagen de una partitura musical que presenta numerosos pentagramas, los máximos de la transformada de Hough se corresponderán con las líneas que componen dichos pentagramas. En consecuencia uno de los máximos de la transformada de Hough nos proporcionará la ecuación de la recta en coordenadas polares asociada a una de las líneas de los pentagramas y a partir de la misma podremos obtener directamente una estimación del ángulo de inclinación. Posteriormente, por si acaso la imagen presenta aún cierta inclinación, se opta por un ajuste de tanteo. Este consiste en aplicar inclinaciones a la imagen en los márgenes de error que tenemos tras aplicar el paso de la transformada de Hough. La realimentación para determinar si las correcciones que estamos realizando están surtiendo efecto o no (y en definitiva para obtener convergencia hacia la solución) se basan en el estudio de la proyección horizontal del histograma de la imagen resultado.

Segmentación de pentagramas (OMR.m)

A continuación segmentamos la imagen separando pentagramas y título de la partitura. Esto lo logramos con ayuda de la proyección horizontal de histograma. Dado que el fondo es completamente negro, las zonas inter-pentagrama se caracterizarán en la proyección horizontal de histograma por una secuencia continua de ceros por lo que, mediante un algoritmo de búsqueda de ceros basado en comparación con el valor anterior y posterior, establecemos el punto medio entre dos pentagramas para que estos sean separados.

Segmentación de símbolos (segSimb.m)

Ahora, para cada uno de los pentagramas, así como para el título, la función programada en MatLab se encarga de segmentar cada uno de los símbolos que componen la línea con la ayuda de la proyección vertical de histograma de cada una de las líneas por separado. En el caso de un pentagrama, en su proyección vertical de histograma se podrá observar un nivel mínimo de continua correspondiente a la contribución de las cinco líneas del pentagrama. Este nivel, que conocemos como el mínimo no-cero, es detectado para cada uno de los pentagramas individualmente de modo que, con un algoritmo similar al anteriormente presentado para la segmentación de la proyección horizontal de histograma, segmentamos los símbolos de cada una de las líneas detectadas (con la diferencia del uso de un algoritmo de búsqueda del mínimo no-cero comentado). Estos límites intersimbólicos podrían aprovecharse para la segmentación directa del elemento en la imagen. No obstante, se prefirió realizar una caracterización del símbolo mediante el estudio del contorno de la proyección vertical de histograma del mismo, ya que podemos limitar su extensión mucho más fácilmente respecto de la extensión del símbolo en el interior del pentagrama (comentar en este sentido que la caracterización del símbolo mediante el estudio de seis momentos estadísticos arrojó peores resultados respecto del método finalmente implementado).

Reconocimiento de símbolos (segSimb.m y funClasif.m)

A continuación se procede a la clasificación de cada uno de los símbolos detectados según el siguiente patrón. Ya hemos dicho que realmente estamos segmentando la proyección de histograma vertical del símbolo. Pues bien, con esta información hacemos lo siguiente:

- Eliminación de la componente de continua detectada como el mínimo no-cero anteriormente comentado.

- Normalización de los valores de la proyección vertical de histograma del símbolo.

- Obtención de 24 coeficientes del módulo de la FFT de la proyección vertical de histograma del símbolo.

- Comparación con cada uno de los patrones base almacenados por clase en la base de datos. La selección de la clase a la cual pertenece el símbolo se realiza en base a la minimización del error lineal calculado como la acumulación de las diferencias en valor absoluto entre cada par de coeficientes de los módulos de las FFT’s del símbolo en estudio y del patrón de la base de datos con el cual se compara.

El estudio del símbolo en base a la FFT se basa, junto con la normalización, en la importancia de caracterizar con un número standard de coeficientes el símbolo, de modo que se contemple la propiedad de la invarianza a la escala.

La función de MatLab asume que el símbolo que se estudia se va a encontrar en la base de datos, por lo que se selecciona mediante minimización del error y no se establece ningún umbral necesario tipo tolerancia para poder clasificarse o no.

Otra limitación importante es que la función también asume que el título sólo ocupa la primera línea segmentada mediante la proyección horizontal de histograma. Esto se ha pensado así debido a la forma de segmentar un símbolo dentro de un pentagrama. Recordemos que el símbolo se lleva sobre el nivel cero (eliminación de la componente de continua introducida por las cinco líneas de pentagrama). Pues bien, el programa hace una excepción sobre el cálculo del nivel mínimo de continua para la primera línea ya que, si se selecciona un mínimo no-cero en la línea del título, podemos perder y distorsionar la información no proveyendo un resultado coherente.

Comentar que sobre los cinco sencillos ejemplos provistos se obtienen unos resultados excelentes con este método de clasificación.

Los símbolos contemplados en la base de datos son los siguientes. Notar que, para mayor precisión en la selección, por cada símbolo se incluyen varias caracterizaciones en función de si este puede aparecer reflejado, invertido (piénsese en una figura como una blanca)…

El reconocimiento de más símbolos aparte de los incluidos en la base de datos sería inmediatamente extensible, pudiendo cubrir con este método partituras de mayor complejidad si fuese necesario.

Como se observa, tampoco se ha incluido el reconocimiento de armaduras por simplicidad (sólo se admiten partituras cuya tonalidad sea de Do mayor o su relativo menor La menor) ni alteraciones accidentales. Por tanto, las partituras que se proporcionan están todas transportadas a la tonalidad de Do mayor.

Detección de tono (deteccionTono.m)

Una vez tenemos segmentada la imagen y se han reconocido las figuras musicales se procede a la detección de tono donde se requiera. El método seguido para la detección de tono está basado en el uso de la proyección horizontal del histograma. Dado el trozo de la imagen hallamos su proyección horizontal (sobre el eje y). A continuación se sigue un método en el que se determina tanto la posición de las líneas del pentagrama como las líneas virtuales (se pueden definir 13 líneas desde el Re 7 al Sol 3, 5 visibles y 8 virtuales). Lo que hacemos es recorrer la proyección horizontal hasta determinar las líneas del pentagrama. Como las líneas del pentagrama estarán representadas en blanco y ocuparán todo el trozo de la imagen se corresponderán con máximos. Así determinamos sus posiciones y las distancias de separación entre cada par de líneas. Haciendo un promedio de la separación entre pares de líneas adyacentes visibles podemos obtener una estimación de la separación (aplicando redondeo por estar en un espacio discreto) entre líneas. Con esta separación estimada podemos determinar las posiciones de las 8 líneas virtuales restantes.

Una vez determinadas todas las posiciones vamos recorriendo la proyección del histograma situándonos en las líneas estimadas y en el punto medio de ambas. En cada posición haremos la suma de todos los elementos del histograma en un entorno igual a la mitad de la separación estimada.

Finalmente en aquella posición (de las 26 posibles, líneas y puntos medios entre dos líneas) donde se obtenga una suma mayor se considera que se encuentra la figura musical y en consecuencia podremos establecer.

Los distintos ficheros y partituras de prueba se pueden descargar aquí. El fichero demo.m ejecuta un ejemplo, obteniéndose el resultado en un fichero de texto con los distintos símbolos reconocidos y las frecuencias auditivas de las notas musicales.

Limpieza de imágenes biomédicas

2010-04-20T00:36:00.004+02:00

La limpieza de imágenes biomédicas consiste en obtener una imagen de la muestra de interés lo más pura posible. Esta imagen, normalmente se obtiene aplicando operaciones matemáticas sencillas, basadas en transformaciones globales, sobre las imágenes obtenidas a través del microscopio con el portaobjetos con y sin muestra.

La función que finalmente he desarrollado es muy sencilla y en ella no se puede controlar ningún parámetro, sino que devuelve un resultado unívoco para las imágenes de entrada. Si la función del parámetro está bien definida, es interesante no pasar por alto su contribución en un problema de índole aleatoria como es este, donde no conocemos las características de las imágenes de entrada (no con total determinismo).

En primer lugar restamos píxel a píxel el fondo a la imagen de la muestra contaminada con el ruido. Esto lo hacemos así y no a la inversa debido a que la muestra contrasta con el fondo por estar a contraluz, es decir, la muestra es más oscura. Si lo hiciésemos al contrario perderíamos mucha información siendo, en algunos casos, imposible dar una solución medianamente aceptable al problema.

Una vez restadas las imágenes, obtendremos un resultado bastante limpio, como es lógico, aunque con poco contraste. Por ello, le aplicamos stretching, cuya función de transferencia adjunto en la siguiente figura, así como el código de la función implementada en MatLab.

H(l) = 255(l - mín)/(max - mín)

Donde max y mín son los niveles de gris máximo y mínimo presentes en la imagen.

% STRETCHING
% Autor: Iván López Espejo
%
% y = stretch(x)
%
% Donde en 'y' se obtiene la imagen tras aplicar la operación, siendo % 'x' la imagen de entrada.
function y = stretch(x)
% Algunos parámetros necesarios.
minimo = double(min(min(x)));
maximo = double(max(max(x)));
filas = length(x(:,1));
columnas = length(x(1,:));
% Definición de la función de transferencia y aplicación.
for j = 1:1:filas
for k = 1:1:columnas
y(j,k) = round((255/(maximo-minimo))*double(x(j,k))-minimo*(255/(maximo-minimo)));
end
end

Finalmente sólo nos queda invertir los niveles de gris de la imagen resultante a fin de obtener el objeto de interés en negro sobre el fondo blanco. El siguiente código implementa todo el sistema haciendo uso de la función anteriormente descrita. Además, las siguientes figuras recogen algunos resultados obtenidos sobre ejemplos reales.

% FUNCIÓN PARA LA LIMPIEZA DE IMÁGENES RUIDOSAS
% Autor: Iván López Espejo
function y = Biolimp2(img,fon)

% Diferencia del fondo y la imagen.
img2 = fon - img;
% Stretching del resultado.
y2 = stretch(img2);
% Inversión del resultado y muestra por pantalla.
for j = 1:1:length(y2(:,1))
for k = 1:1:length(y2(1,:))
y(j,k) = 255 - y2(j,k);
end
end
subplot(2,1,1), imshow(img), title('Imagen original')
subplot(2,1,2), imshow(y), title('Imagen limpia')

Pueden existir muchas variantes del anterior método mucho más complejas como, por ejemplo, extraer estadísticamente algún parámetro de la imagen para delimitar el tipo de crecimiento (entre exponencial y logarítmico) de una función de stretching en forma de S, que “mejorase” el contraste entre el fondo y la muestra biomédica. No obstante, dado que el análisis final será personal y cualitativo, estas pequeñas diferencias no aportarán nada o prácticamente nada, por lo que pienso que es preferible una función sencilla y eficiente a algo más complejo que será desaprovechado finalmente por el ojo humano (en principio).

Perfilador de imágenes

2010-03-29T02:02:00.000+02:00

Un perfilador de imágenes se encargaría de conseguir el efecto contrario al deseado cuando se aplica la técnica de suavizado. Para resaltar los bordes de una imagen y así perfilarla, es necesario detectarlos previamente mediante alguna técnica de derivación de la imagen (en mi caso reutilizaré el detector de bordes mediante filtrado lineal gaussiano). Una vez hecho esto, sería suficiente con realizar alguna ponderación sobre la matriz de bordes detectados, sumando el resultado a la imagen original.

La siguiente función para MatLab se encarga de ello. Funciona exactamente igual que el detector de fronteras explicado unas entradas más abajo. Sólo se ha incluido una inversión de la imagen filtrada (debido a la característica primitiva del filtrado) y la comentada ponderación y posterior adición.

A continuación se adjuntan el código y un ejemplo aplicado sobre una fotografía de la ciudad de Brisbane (Australia).

% Perfilador
% Autor: Iván López Espejo
%
% Función para perfilar imágenes (efecto contrario al suavizado).
%
% Y = perfilador(x,var,dim,f)
%
% Donde en 'Y' se almacena la matriz imagen resultado, siendo 'x' la matriz
% imagen en niveles de gris de entrada sobre la que queremos aplicar el
% perfilado. Por otro lado, 'var' es la varianza de la función gaussiana 2-D
% (a mayor varianza, mayor consideración de figuras como borde (mayor
% densidad de bordes)) y 'dim' es el tamaño, en píxels, del área del filtro
% con el que convolucionaremos (promediaremos). A mayor 'dim',
% representación de las fronteras como líneas de mayor grosor. Por último,
% 'f' es el factor de escalado de la imagen de fronteras detectada para
% ponderar su suma con la imagen original y obtener el efecto deseado.

function y = perfilador(x,var,dim,f)

% Algunos parámetros necesarios.
filas = length(x(:,1));
columnas = length(x(1,:));

% Generación de la respuesta espacial del filtro gaussiano.
H = zeros(dim,dim);
for j = 1:1:dim
    for k = 1:1:dim
        H(j,k) = exp(-(( ((j)^2) + ((k)^2) - 2*(round(dim/2))*(j) - 2*(round(dim/2))*(k) + (round(dim/2)^2) + (round(dim/2)^2) )/(2*var)));
    end
end
% Normalización del volumen del filtro en el dominio espacial.
k = (1-sum(sum(H)))/(dim^2);
H = H + k;
% Representación de la magnitud de la respuesta espacial del filtro.
surf(H)
% Filtrado mediante convolución en el dominio espacial y representación.
z = filter2(H,x);
% Inversión del resultado por característica del filtrado.
z = uint8(255*z);
for j = 1:1:filas
    for k = 1:1:columnas
        z(j,k) = 255 - z(j,k);
    end
end
% Perfilado de la imagen con la matriz de fronteras detectadas.
y = f.*z + x;
figure
subplot(2,1,1), imshow(x), title('Imagen original');
subplot(2,1,2), imshow(y), title('Imagen perfilada');

Imagen original.

Imagen perfilada con varianza unidad, dimensión de máscara 5x5 y factor de ponderación 0.2.

Detección de fronteras

2010-03-26T02:41:00.003+01:00

Casualmente he logrado hacer un detector de fronteras para imágenes en niveles de gris (extensible para imágenes en cualquier otro formato de color como RGB).

En principio estaba intentando hacer un simple filtro de suavizado (o desenfoque) gaussiano como primer paso para la realización de un filtro de Wiener de enfoque gaussiano dada una imagen de entrada desenfocada (debido a que el desenfoque producido por una cámara puede ser modelado mediante una función gaussiana).

No obstante, por error, he logrado realizar un detector de fronteras mediante un simple filtrado con el que se puede controlar tanto la resolución o densidad de bordes como la anchura de los mismos (mediante la parametrización de la varianza y tamaño del bloque espacial convolucionable).

Es una solución muy simple, rápida, eficiente y elegante basada, como digo, en un filtrado. La clave de que esto haya podido ser así, he observado que se asienta sobre la normalización del volumen que encierra el área gaussiana representada en la respuesta espacial del bloque convolucionable. Esta normalización se produce por desplazamiento en altitud de la respuesta espacial del filtro, lo que se traduce en un filtrado pasa-alta con coeficientes negativos y otro pasa-baja con coeficientes positivos. A partir de haber hecho esto (pensado como mero paso para un correcto promediado que no alterase los niveles de gris de la imagen original) pude observar sin más la detección de fronteras.

A continuación adjunto unos ejemplos y un par de respuestas espaciales de los filtros, así como el código de la función implementada en MatLab.

Imagen original 1.

Fronteras detectadas por filtrado de la imagen 1 (var = 2 dim = 10).

Imagen original 2.

Fronteras detectadas por filtrado de la imagen 2 (var = 4 dim = 12).

Respuesta espacial del filtro aplicado a la imagen 2.

Imagen original 3.

Fronteras detectadas por filtrado de la imagen 3 (var = 1 dim = 10).

Respuesta espacial del filtro aplicado a la imagen 3.

% Detector de fronteras
% Autor: Iván López Espejo
%
% Función para la detección de fronteras basada en filtrado gaussiano.
%
% Y = detBordes(x,var,dim)
%
% Donde en 'Y' se almacena la matriz imagen resultado, siendo 'x' la matriz
% imagen en niveles de gris de entrada sobre la que queremos aplicar el
% filtrado. Por otro lado, 'var' es la varianza de la función gaussiana 2-D
% (a mayor varianza, mayor consideración de figuras como borde (mayor
% densidad de bordes)) y 'dim' es el tamaño, en píxels, del área del filtro
% con el que convolucionaremos (promediaremos). A mayor 'dim',
% representación de las fronteras como líneas de mayor grosor.

function y = detBordes(x,var,dim)

% Generación de la respuesta espacial del filtro gaussiano.
H = zeros(dim,dim);
for j = 1:1:dim
    for k = 1:1:dim
        H(j,k) = exp(-(( ((j)^2) + ((k)^2) - 2*(round(dim/2))*(j) - 2*(round(dim/2))*(k) + (round(dim/2)^2) + (round(dim/2)^2) )/(2*var)));
    end
end
% Normalización del volumen del filtro en el dominio espacial.
k = (1-sum(sum(H)))/(dim^2);
H = H + k;
% Representación de la magnitud de la respuesta espacial del filtro.
surf(H)
% Filtrado mediante convolución en el dominio espacial y representación.
y = filter2(H,x);
figure
subplot(2,1,1), imshow(x), title('Imagen original');
subplot(2,1,2), imshow(y), title('Imagen filtrada');

Suavizado de bordes en imágenes

2010-03-19T04:00:00.000+01:00

En este documento explico cómo abordo el planteamiento y construcción de un sistema que he ideado para el suavizado (especialmente de bordes) de imágenes.

La idea del sistema es bien sencilla y está basada en minimización de gradiente. El método se encarga de iterar el siguiente algoritmo:

1. Para el píxel i-ésimo de la imagen (con i = 1,2,...,FxC) almacenamos en un vector el valor de las muestras diferenciales calculadas como la diferencia del nivel de gris del píxel i-ésimo con cada uno de los niveles de gris de los píxels de su vecindad.
2. Mediante comparación, para cada uno de los valores del anterior vector construido (en valor absoluto), estudiamos si alguno de ellos excede el valor diferencial máximo impuesto por el parámetro de entrada, donde este valor diferencial máximo se obtiene como:

diferencia_máxima = (1 - parámetro)(Niveles - 1)

Donde Niveles hace referencia al número de niveles de gris de la imagen (por defecto Niveles = 256) y parámetro es un valor real en el intervalo [0,1], donde 0 indicaría suavizado mínimo (ya que establecería la diferencia máxima de niveles de gris entre dos píxels contiguos en 255) y 1 haría lo propio referenciando suavizado máximo (ya que no permitiría diferencia de niveles de gris ninguna en la imagen, lo que se traduce en una imagen final de un nivel de gris continuo).
3. Si ninguno de los valores del vector, en valor absoluto, excede el diferencial máximo de niveles de gris, se vuelve al paso 1, hasta completar el recorrido por la imagen. Si alguno sí supera esta cota, calculamos el nuevo valor de gris del píxel, a partir del píxel vecino que provoca la máxima diferencia (gradiente), de la forma:

I(i) = I(i) + signo(diferencia_máxima - gradiente)

Donde signo controla el signo de la resta entre el píxel i-ésimo y uno de los vecinos.

4. Una vez se hayan recorrido todos los píxeles de la imagen con este procedimiento, se repetirá el algoritmo tantas veces como sea necesario hasta que en la n-ésima iteración no se haya producido ningún cambio en la imagen.

De lo anterior se deduce que el algoritmo ha de ser convergente en un número finito de pasos. No he procedido a su demostración matemática, aunque parece intuitivo que esto sea así. El caso que sí podría poner en peligro la convergencia del método, sería el de una imagen que no tuviese bordes, como se daría en la situación de una imagen continua proyectada sobre una superficie esférica.

Notar que el concepto de vecindad se basa en un 8-entorno, así como que el comportamiento de la implementación en los bordes pasa por considerar únicamente los píxels reales contiguos (es decir, en el caso de un píxel en una esquina, el entorno estaría compuesto por 3 píxels).

A continuación adjunto el código MatLab y un ejemplo de una fotografía a la que le aplico el sistema con diferentes coeficientes de suavizado. Notar que es mejor usar un coeficiente intermedio (aunque depende del caso) para no falsear en exceso la imagen original y dar un pequeño suavizado.

% SUAVIZADOR DE BORDES
% Autor: Iván López Espejo
%
% Esta función se encarga de suavizar una imagen (especialmente los bordes)
% mediante minimización de gradiente. Gracias a un parámetro de entrada,
% controlaremos el máximo valor permitido entre muestras diferenciales.
%
% Uso de la función:
% imagen_suavizada = suavizaBordes(x,param)
%
% Donde x es la imagen constituida por una matriz MxN (en niveles de gris)
% y param es un parámetro que puede valer entre 0 y 1, donde 0 se
% corresponde con mínimo suavizado y 1 con máximo suavizado. El resultado
% se devuelve en forma de matriz en la variable imagen_suavizada.

function y = suavizaBordes(x,param)

if param < 0 || param > 1
    disp('Error: el parámetro de suavizado es un número real entre 0 y 1 inclusives');
else
    % Algunos parámetros que usaremos.
    grayL = 256;
    y = x;
    filas = length(x(:,1));
    columnas = length(x(1,:));
    vGr = zeros(1,8);
    ind = 1;
    check = 1;
    while check ~= 0
        % Reestablecemos la variable de comprobación de no-modificación
        % para la n-ésima iteración del algoritmo. Si en la presente
        % iteración converge, no habrá modificación ni de la imagen ni de
        % la variable y habremos terminado.
        check = 0;
        for j = 1:1:filas
            for k = 1:1:columnas
                % Obtención del gradiente en el píxel j-k-ésimo.
                if j > 1 && k > 1
                    vGr(ind) = y(j,k) - y(j-1,k-1);
                    ind = ind + 1;
                end
                if j > 1
                    vGr(ind) = y(j,k) - y(j-1,k);
                    ind = ind + 1;
                end
                if j > 1 && k < columnas
                    vGr(ind) = y(j,k) - y(j-1,k+1);
                    ind = ind + 1;
                end
                if k < columnas
                    vGr(ind) = y(j,k) - y(j,k+1);
                    ind = ind + 1;
                end
                if j < filas && k < columnas
                    vGr(ind) = y(j,k) - y(j+1,k+1);
                    ind = ind + 1;
                end
                if j < filas
                    vGr(ind) = y(j,k) - y(j+1,k);
                    ind = ind + 1;
                end
                if j < filas && k > 1
                    vGr(ind) = y(j,k) - y(j+1,k-1);
                    ind = ind + 1;
                end
                if k > 1
                    vGr(ind) = y(j,k) - y(j,k-1);
                end
                vAbs = abs(vGr);
                [maxi,pos] = max(vAbs);
                valor = round((1-param) * (grayL-1));
                % En caso de que el valor diferencial de muestras exceda el
                % definido por el parámetro de entrada (respecto del peor caso),
                % indicamos que hemos hecho una modificación en el actual
                % recorrido por la imagen y establecemos el píxel de tal
                % modo que la diferencia de nivel de gris entre dicho píxel
                % y el vecino, no exceda del valor establecido por el
                % parámetro de suavizado.
                if valor < maxi
                    check = 1;
                    signo = vGr(pos)/abs(vGr(pos));
                    y(j,k) = y(j,k) + signo * (valor - maxi);
                end
                % Reinicialización de variables para un nuevo análisis de píxel.
                ind = 1;
                vGr = vGr.*0;
            end
        end
    end
    % Pintado de las imágenes original y tratada.
    figure
    subplot(2,1,1), imshow(x), title('Imagen original')
    subplot(2,1,2), imshow(y), title('Imagen suavizada')
end

Imagen original algo pixelada.

Imagen suavizada con coeficiente 0.8

Imagen suavizada con coeficiente 0.85

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

2010-03-03T18:19:00.000+01:00

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;

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

2010-03-02T01:03:00.014+01:00

Como dije, en esta entrada vamos a ver un ejemplo simplificado del sistema. En el caso de querer aplicarlo sobre instrumentos reales habría que realizar un estudio más extenso y a fondo, el cual no voy a hacer por falta de tiempo, presupuesto y apoyo, pero se me ocurriría que podría realizarse con ventanados, pues la característica del timbre del instrumento se encuentra distribuida a lo largo de todo el eje de frecuencias y no sólo se concentra en los tonos y armónicos componentes del sonido polifónico.

Supongamos que tenemos un acorde de La menor construido sobre un La 4 (La-Do-Mi). No obstante, las frecuencias de los tonos fundamentales de las notas componentes del acorde son 435Hz (para el La a 440Hz), 520Hz (para el Do a 523.25Hz) y 650Hz (para el Mi a 659.26Hz), por lo que está desafinado. (idealmente 440Hz),

En la siguiente figura se observa un conjunto de 1000 muestras del acorde:

Figura 2. 1000 muestras del acorde desafinado.

El sonido polifónico está muestreado a 44100Hz (calidad de audio) y está compuesto por, aparte de los tres tonos fundamentales, 20 armónicos cuya amplitud sufre un decaimiento exponencial por cada uno de los tonos fundamentales. A continuación, el bloque Estimador de Tonos Fundamentales calcula el módulo de la DFT del acorde entre 0 y pi. Como se aprecia en la figura 3, esta etapa implementa un algoritmo de selección de posición de máximos que se correspondan con las frecuencias de los tonos fundamentales. Por simplicidad, cada muestra se refiere a un hertzio, por lo que sobre el eje de abscisas se puede ver directamente la frecuencia fundamental de cada una de las notas componentes del acorde:

Figura 3. Espectro del acorde desafinado.

Por ello, el vector devuelto por el bloque será v = (435,520,650).

Ahora, por cada nota detectada, se procederá del siguiente modo. La etapa Filtro, conocida la frecuencia fundamental de la nota i-ésima, con i = 1,2,…,M, siendo en este caso M = 3, generará el filtro en el dominio del tiempo anteriormente descrito. La siguiente secuencia de figuras recoge 1000 muestras de cada una de las respuestas temporales de los filtros según el vector v:

Figura 4. Respuesta en el dominio temporal del filtro para la nota cuya frecuencia fundamental se sitúa en 435Hz.

Figura 5. Respuesta en el dominio temporal del filtro para la nota cuya frecuencia fundamental se sitúa en 520Hz.

Figura 6. Respuesta en el dominio temporal del filtro para la nota cuya frecuencia fundamental se sitúa en 650Hz.

Para cada uno de los filtros anteriormente generados se calcula el producto de sus respuestas en frecuencia con el espectro del acorde, lo que nos devuelve, en este caso, tres espectros (puesto que M = 3) que se corresponden con los propios de las notas componentes del sonido polifónico por separado.

“Paralelamente”, el Clasificador Lineal hace uso de las frecuencias estimadas de los tonos fundamentales para determinar hacia qué nota debe afinarse cada una de las componentes de la polifonía según lo descrito. Es decir:

Realizamos una simple clasificación lineal basada en minimización de distancia en el espacio de las notas musicales, por lo que deducimos que las notas componentes son La4-Do5-Mi5, como se puede extraer de la tabla anteriormente presentada. Las frecuencias reales de las notas son, en consecuencia:

Todos estos datos son recogidos por la etapa Afinador por Desplazamiento de Espectro. De este modo, para cada una de las notas componentes del acorde, calculamos el parámetro alpha:

Estos son los parámetros de escalado temporal motivados por la necesidad de modificar el pitch de cada una de las notas desafinadas. Pensemos en que el tono fundamental de cada una de las notas componentes de la polifonía se puede aproximar por la función cos(2pift), por lo que, para conseguir que dicho tono se sitúe ahora a la frecuencia f', precisaremos de multiplicar el anterior argumento por alpha = f'/f, de tal modo que cos(alpha2pift) = cos(2pif't).

Esta idea la hacemos corresponder con la propiedad de escalado temporal de la transformada de Fourier anteriormente expuesta. A continuación se muestra la secuencia de espectros componentes ya afinados:

Figura 7. Espectro del La 4 afinado.

Figura 8. Espectro del Do 5 afinado.

Figura 9. Espectro del Mi 5 afinado.

La última etapa del diagrama de bloques se encargaría de sumar las transformadas inversas (muestras temporales) de cada uno de los espectros de las notas componentes del acorde, lo que se traduce en el sonido polifónico ya afinado. La siguiente gráfica recoge una comparación entre las primeras 1000 muestras del acorde desafinado (en azul) y las primeras 1000 muestras del acorde ya afinado (en rojo). Observar que el filtro introduce una cierta distorsión (sólo sería cuestión de ir "afinando" el sistema :D).

Figura 10. Comparación entre acordes afinado y desafinado en el dominio del tiempo.

Podéis descargar el ejemplo comentado a continuación. Notar que por algún problema relacionado con el filtrado o con el desplazamiento de espectro, se introdujo un efecto extraño en el acorde afinado, por lo que para mayor representatividad copié y pegué una secuencia del período más limpio. No obstante, los cortes sí se notan en el ensamblaje, lo que se traduce en un pequeño piqueteo.

ACORDE DESAFINADO | ACORDE AFINADO

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

2010-02-28T01:37:00.008+01:00

Este texto recoge la descripción fundamental del procedimiento ideado para la afinación automática (asistida por computador) de sonidos polifónicos.

Motivación.

Todos conocemos de la existencia de plug-in’s y otros softwares de ordenador que nos permiten afinar en tiempo real o por postprocesado una línea melódica, ya sea voz humana u otro instrumento musical como puede ser una guitarra o un violín. No obstante aún no existe un sistema comercial de afinación polifónica. Pensemos en que en un estudio de grabación se recoge como línea de entrada una secuencia de acordes de guitarra. Por algún avatar del destino, algunos de los acordes están desafinados y no podemos volver a grabar dicha secuencia y, sin embargo, queremos respetarla porque su ausencia provocaría un vacío en la pieza musical. La solución radicaría en la afinación asistida por ordenador del acorde o acordes en cuestión. El sistema que se presenta a continuación vendría a auxiliar en caso de que este problema se diese.

Otra utilidad similar podría venir dada por la necesidad de mejora de afinación en una canción u obra que se desea remasterizar. En muchos casos de remasterización de audio, no se dispone de la disgregación de pistas necesarias como para tratar cada línea por separado. El sistema de afinación de una obra en conjunción vendría de nuevo a auxiliar este problema.

Descripción de la solución.

La clave de la solución ofertada por este sistema radica en el trabajo realizado en el dominio de la frecuencia ya que, de otro modo, a priori resultaría casi imposible la afinación de polifonía en el dominio del tiempo. Mostraremos una aproximación mediante la generación de un acorde ideal bastante limpio únicamente conformado por los tonos fundamentales y sus respectivos armónicos atenuados exponencialmente. Presentamos a continuación el diagrama de bloques, procediendo seguidamente a su desglose y descripción:

Figura 1. Diagrama de bloques del sistema afinador.

El sonido polifónico desafinado debe ser recogido en primera instancia por el bloque Estimador de Tonos Fundamentales. En esta primera etapa se determinan los tonos fundamentales que componen el sonido polifónico mediante el estudio de máximos del espectro. A tal fin, sea el acorde desafinado en el dominio del tiempo el conjunto de N muestras y(n), obtenemos el espectro del mismo mediante la operación de DFT de la forma:

Seguidamente trabajamos con el módulo de la transformada pues, recordemos que la información sonora relevante se encuentra en el módulo del espectro y no en la fase. Por tanto, calculamos las muestras del módulo (suponemos ahora mismo por simplicidad que disponemos también de N muestras de la DFT) como:

Una vez calculado habremos de obtener un conjunto representativo de posiciones de las muestras que son máximos. Estos máximos han de corresponderse con los picos de los tonos y armónicos que conforman el acorde desafinado. Como no sabemos, a priori, el número de tonos fundamentales (o no tenemos porqué saberlo), podemos curarnos en salud tomando las posiciones de un conjunto amplio de máximos aunque, con la experiencia, podremos ir adquiriendo un compromiso que maximice la eficiencia de tal modo que el programa cumpla su cometido a la perfección en el menor tiempo posible.

Mediante un algoritmo de recorrido y comparación entre todas las posiciones de los máximos obtenidos, clasificamos los tonos fundamentales y desechamos los armónicos. El desechado de los armónicos se realiza mediante comparación por cociente entre el potencial armónico y el potencial tono fundamental. De este modo, si uno de los máximos se corresponde con un armónico, el cociente, idealmente, debe ser un número natural n. Como puede haber cierta desviación, de nuevo, mediante una base de datos estocástica determinaremos el mejor épsilon que discierna entre un armónico del que no lo es. En principio, épsilon = 0.05, algo más que suficiente.

Una vez estimados todos los tonos fundamentales, la etapa devolverá un vector de tonos fundamentales (los cuales se presuponen desafinados).

A continuación, el bloque Clasificador Lineal se encargará, mediante un sistema de decisión lineal, de determinar cuál es el tono fundamental al cual se debe afinar cada uno de los sonidos componentes de la polifonía desafinada. Esto lo hace linealmente adaptando previamente la escala exponencial frecuencial a una función logarítmica. Recordemos que la percepción auditiva del ser humano no es lineal con la frecuencia, sino que la sensación de linealidad proviene del crecimiento exponencial de la frecuencia. Recurriendo a la siguiente expresión:

Donde log(2) hace referencia al logaritmo neperiano de 2. Nos interesa que la función sólo sea dependiente de n, donde n corresponde a la siguiente asociación:

Como es lógico, para terminar de caracterizar la nota necesitaremos especificar la octava en la que esta se encuentra, y de ello se encarga la variable o. Podemos dejar expresada la octava, o, en función de la frecuencia como:

En consecuencia, la función resulta:

La función anterior resulta en un decisor lineal, de modo que, por redondeo, podremos establecer la nota a la cual, respetando la percepción auditiva, se aproxima más la nota desafinada. Se deduce inmediatamente que si n = 1,2,…,12, la nota componente está perfectamente afinada y no será necesario realizar ninguna afinación por desplazamiento de espectro.

El Estimador de Tonos Fundamentales enviará el vector de tonos al Clasificador Lineal, de tal modo que, por ejemplo, imaginemos que uno de los tonos fundamentales se situase según lo estimado a 450Hz (f = 450Hz). Aplicando la anterior expresión deducida:

Como 0.3891 < 0.5, n = 10 y, por tanto, se correspondería con un La 4 (a 440Hz). Podemos modificar la tabla decisoria de forma que la afinación no se correspondiese con distancias de medio tono, sino que agrupase armonías concretas, siendo ello mera particularización de lo expuesto. Notar también que consideramos un Do# igual a un Reb, aunque ello no sea estrictamente ortodoxo en la práctica.

El bloque Filtro se encarga de filtrar todo el espectro correspondiente al acorde desafinado de entrada, únicamente respetando las componentes asociadas a cada uno de los tonos, por separado, estimados por el bloque Estimador de Tonos Fundamentales. Realizando esta separación, lo que logramos es simplificar el problema, llevándolo posteriormente a una mera afinación de un sonido aislado por desplazamiento del espectro y finalmente sumando todas las muestras afinadas en el dominio del tiempo reconstruyendo la polifonía ya arreglada.

El filtro se genera para cada caso particular según cada una de las frecuencias estimadas. Este filtro no es más que la superposición en el dominio del tiempo del tono fundamental estimado y todos los armónicos audibles con una dispersión de 2Hz, para tratar de no modificar altamente el timbre del instrumento en cuestión. Aplicando la operación de DFT anteriormente descrita, obtenemos la respuesta del filtro en el dominio de la frecuencia, de modo que la convolución en el dominio del tiempo del filtro con el acorde se transforma en el producto de los espectros del filtro y del acorde, seleccionando así una sola nota de las que componen el sonido polifónico:

Por tanto, este bloque se encargaría de, para cada una de las frecuencias fundamentales estimadas, devolver el espectro de las M notas componentes, donde, en consecuencia, i = 1,2,...,M.

Ahora, con ayuda del Clasificador Lineal, el bloque Afinador por Desplazamiento de Espectro se encargará de afinar cada una de las notas componentes a partir de los espectros proporcionados por el Filtro. Este bloque aprovecha el conocimiento de la frecuencia fundamental estimada y la real o ideal para definir el parámetro alpha = frecuencia ideal/frecuencia estimada.

Aprovechamos la propiedad de escalado temporal de la transformada de Fourier que nos dice:

En consecuencia, cada uno de los espectros devueltos por el Filtro, Xi(k), lo transformaremos en:

Donde los índices de muestas k/alpha se deben redondear al entero más cercano (el prototipo trabaja con una muestra por hertzio, aunque mejoraría mucho su funcionamiento si lo hiciese con una mayor resolución, únicamente escalando los anteriores resultados).

Una vez afinada la nota i-ésima por desplazamiento de espectro componente del sonido polifónico, esta se pasa ya convertida mediante DFT inversa al bloque Buffer Sumador en el Dominio del Tiempo:

Aquí esperan todas las notas componentes hasta que están todas listas, momento en el que son superpuestas para conformar el sonido polifónico ya afinado.

En la próxima entrada colocaré un ejemplo práctico del sistema ideado y presentaré los problemas que aún quedan por solventar.

Conversor de texto a voz

2010-01-24T21:32:00.002+01:00

Mediante la segmentación de difonemas he construido una versión primitiva de un conversor de texto a voz. Un difonema no es más que la transición entre la parte estable de un fonema y otro contiguo, incluidas las partes estacionarias. Mediante la concatenación de estas secuencias de muestras conseguimos construir frases de voz.

La mayor dificultad radica en la correcta segmentación de los difonemas. Además, es un trabajo complicado que requiere de paciencia y de asistencia en tiempo real de un usuario. Yo he conseguido realizar una aproximación más o menos válida mediante la grabación de difonemas con agregación de vocales iniciales y finales. Mediante un programa que hice, generé todos los difonemas posibles y, estos, junto a la anterior agregación de vocales, fueron grabados con una separación de 1 segundo con la ayuda de un metrónomo a una frecuencia de muestreo de 44100Hz, lo que permitió, como se explica en la presentación adjunta, una segmentación programada y asistida por ordenador.

Los códigos de MatLab se pueden ver en la presentación y, a continuación, algunos ejemplos generados con el programa creado.

"Calculadora"

"Esto es una prueba de una frase más larga"

"Explícamelo, por favor"

"Impedancia equivalente"