Formato de publicación de la escuela técnica
Abstract
En este trabajo se presenta una solución para el problema de Separación Ciega de Fuentes, especialmente para el caso de fuentes de audio unísonas y moduladas, tanto en tiempo como en frecuencia mediante un vibrato. A partir de una base teórica sobre el problema y sobre el método para solucionarlo, la Factorización No Negativa de Matrices, se expone un modelo de forma teórica y práctica y su posible ejecución en Matlab ® mediante un algoritmo. El caso en el que se centra este trabajo, el de dos fuentes moduladas y unísonas, representa un gran desafío ya que con NMF no se pueden asumir ciertas suposiciones de importante relevancia para una correcta separación. En la búsqueda de una solución eficiente, vamos a trabajar con tensores, haciendo uso de la Factorización No Negativa de Tensores, que es una ampliación de NMF a tensores. Dividiremos la STFT de la señal de audio original en parches solapados y calcularemos la 2D-DFT a cada parche, obteniendo así un tensor de 4 dimensiones. A esto es a lo que se ha llamado Transformada de Destino Recurrente, que se incluye en el Modelo de Destino Recurrente. Sobre este tensor aplicaremos el algoritmo multiplicativo de separación, lo que marca la diferencia frente a NMF, que lo aplica sobre la STFT. Obtenida la separación, resulta necesario evaluar la calidad de la misma, para lo que se ha usado la herramienta BSS Eval. Después se ha hecho un estudio sobre la dependencia del algoritmo a los parámetros alfa y beta. Se han buscado valores óptimos de separación en función de estos parámetros. Por último, se presentan diversas simulaciones con las que se ha buscado comprobar las ventajas de este modelo y sus posibles carencias, así como sus posibles líneas futuras y las conclusiones que hemos obtenido tras el trabajo realizado.
Full text
Proyecto Fin de Carrera Ingeniería de Telecomunicación Formato de Publicación de la Escuela Técnica Superior de Ingeniería Autor: F. Javier Payán Somet Tutor: Juan José Murillo Fuentes Dep. Teoría de la Señal y Comunicaciones Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2013 Trabajo Fin de Grado Grado en Ingeniería de las Tecnologías de Telecomunicación Separación de audio con modulación Autor: Antonio Márquez Tristán Tutor: Iván Durán Díaz Dep. Teoría de la Señal y Comunicaciones Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2018
Trabajo Fin de Grado Grado en Ingeniería de las Tecnologías de Telecomunicación Separación de audio con modulación Autor: Antonio Márquez Tristán Tutor: Iván Durán Díaz Profesor Titular Dep. Teoría de la Señal y Comunicaciones Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2018
Trabajo Fin de Grado: Separación de audio con modulación Autor: Antonio Márquez Tristán Tutor: Iván Durán Díaz El tribunal nombrado para juzgar el trabajo arriba indicado, compuesto por los siguientes profesores: Presidente: Vocal/es: Secretario: acuerdan otorgarle la calificación de: El Secretario del Tribunal Fecha:
Agradecimientos Este trabajo representa el fin de una de las etapas más emocionantes de mi vida. Tras un inicio universitario muy duro, en el cual estuve a punto de tirar la toalla, me causa tremenda alegría verme redactando estas líneas. Han sido muchas las personas que me han apoyado desde el principio y han hecho que haya podido vencer todos los momentos poco gratificantes que esta carrera me ha brindado, incluso cuando los resultados no eran muy halagüeños, ellos han confiado en mis posibilidades, en ocasiones, incluso más que yo mismo. Llegados a este punto solo me queda agradecerles su confianza, ya que para mí es todo un orgullo poner fin a esta etapa y verme totalmente capacitado para enfrentarme a la siguiente. De todas esas personas tengo que destacar a mis padres, Manolo y Antonia, que son para mí un modelo a seguir y el mayor apoyo que he tenido durante estos años. A mis hermanos, Manolo y Andrés, que también han confiado ciegamente en mí y que han financiado buena parte de mi vida de ocio durante la carrera. Y a mis amigos de toda la vida: Joaquín, Pepe, Jesuli, Sergio y Jero. Sumamos etapas siendo siempre un apoyo los unos para los otros aún cuando por motivos académicos hemos podido vernos menos que en otras épocas. A todas las buenas personas que he conocido en la escuela, sin ellos toda esta lucha hubiera sido mucho más dura. En especial, a mi grupo de amigos "Los Peppers", sin vosotros mis futuros recuerdos de la universidad hubieran estado faltos de buenos ratos, de cante y palmas después de esas largas jornadas de estudio. A Mari Ángeles, mi compañera y la que ha tenido que aguantar todos mis malos momentos durante la ejecución de este proyecto, sin tu incondicional apoyo y tu cariño durante estos meses, no hubiera sido lo mismo. Y por último y no menos importante, agradecer a mi tutor Iván Durán Díaz que me haya dado la posibilidad de hacer este proyecto que tan interesante me ha resultado, y que me haya ayudado a poder realizarlo con éxito. Al señor Fabian-Robert Stöter, por resolverme varias dudas sobre su trabajo y por compartir amablemente material imprescindible para la resolución de este proyecto. También es justo agradecer a todos los profesores que han compartido sus conocimientos conmigo y mis compañeros durante estos años. Antonio Márquez Tristán Sevilla, 2018 I
Resumen En este trabajo se presenta una solución para el problema de Separación Ciega de Fuentes, especialmente para el caso de fuentes de audio unísonas y moduladas, tanto en tiempo como en frecuencia mediante un vibrato. A partir de una base teórica sobre el problema y sobre el método para solucionarlo, la Factorización No Negativa de Matrices, se expone un modelo de forma teórica y práctica y su posible ejecución en Matlab ® mediante un algoritmo. El caso en el que se centra este trabajo, el de dos fuentes moduladas y unísonas, representa un gran desafío ya que con NMF no se pueden asumir ciertas suposiciones de importante relevancia para una correcta separación. En la búsqueda de una solución eficiente, vamos a trabajar con tensores, haciendo uso de la Factorización No Negativa de Tensores, que es una ampliación de NMF a tensores. Dividiremos la STFT de la señal de audio original en parches solapados y calcularemos la 2D-DFT a cada parche, obteniendo así un tensor de 4 dimensiones. A esto es a lo que se ha llamado Transformada de Destino Recurrente, que se incluye en el Modelo de Destino Recurrente. Sobre este tensor aplicaremos el algoritmo multiplicativo de separación, lo que marca la diferencia frente a NMF, que lo aplica sobre la STFT. Obtenida la separación, resulta necesario evaluar la calidad de la misma, para lo que se ha usado la herramienta BSS Eval. Después se ha hecho un estudio sobre la dependencia del algoritmo a los parámetros alfa y beta. Se han buscado valores óptimos de separación en función de estos parámetros. Por último, se presentan diversas simulaciones con las que se ha buscado comprobar las ventajas de este modelo y sus posibles carencias, así como sus posibles líneas futuras y las conclusiones que hemos obtenido tras el trabajo realizado. III
XÍndice 3.3.1 Definición de Alfa-Beta divergencia 30 3.3.2 Propiedades 31 3.3.3 Justifcación del estudio 32 4 Simulaciones 35 4.1 Datos de entrada 35 4.2 Algoritmo paso a paso 36 4.3 Evaluación de los resultados 37 4.4 Simulación 1 38 4.5 Simulación 2: Estudio de las Alfa-Beta divergencias 40 4.5.1 Análisis de los resultados 40 4.6 Simulación 3 42 4.7 Simulación 4 44 4.8 Simulación 5 46 4.9 Simulación 6 47 4.10 Simulación 7 49 5 Conclusiones y Líneas Futuras 51 5.1 Trabajo realizado y conclusiones 51 5.2 Líneas futuras 52 Índice de Figuras 53 Índice de Tablas 55 Índice de Algoritmos 57 Bibliografía 59
Notación RCuerpo de los números reales CCuerpo de los números complejos kvkNorma del vector v hv,wiProducto escalar de los vectores vyw |A|Determinante de la matriz cuadrada A det(A)Determinante de la matriz (cuadrada) A A>Transpuesto de A A−1Inversa de la matriz A A†Matriz pseudoinversa de la matriz A AHTranspuesto y conjugado de A A∗Conjugado c.t.p. En casi todos los puntos c.q.d. Como queríamos demostrar Como queríamos demostrar Fin de la solución e.o.c. En cualquier otro caso enúmero e ejxExponencial compleja ej2πxExponencial compleja con 2π e−jxExponencial compleja negativa e−j2πxExponencial compleja negativa con 2π IRe Parte real IIm Parte imaginaria sen Función seno tg Función tangente arctg Función arco tangente sinyxFunción seno de xelevado a y cosyxFunción coseno de xelevado a y Sa Función sampling sgn Función signo rect Función rectángulo Sinc Función sinc ∂y ∂xDerivada parcial de yrespecto a x x◦Notación de grado, xgrados. Pr(A)Probabilidad del suceso A E[X]Valor esperado de la variable aleatoria X σ2 XVarianza de la variable aleatoria X ∼fX(x) Distribuido siguiendo la función densidad de probabilidad fX(x) NmX,σ2 X Distribución gaussiana para la variable aleatoria X, de media mXy varianza σ2 X XI
XII Notación InMatriz identidad de dimensión n diag(x)Matriz diagonal a partir del vector x diag(A)Vector diagonal de la matriz A SNR Signal-to-noise ratio MSE Minimum square error :Tal que def =Igual por definición kxkNorma-2 del vector x |A|Cardinal, número de elementos del conjunto A xi,i=1,2,...,nElementos i,de1an, del vector x dxDiferencial de x ⩽Menor o igual ⩾Mayor o igual \Backslash ⇔Si y sólo si x=a+3= ↑ a=1 4Igual con explicación a bFracción con estilo pequeño, a/b ∆Incremento b·10aFormato científico −→ xTiende, con x OOrden TM Trade Mark E[x]Esperanza matemática de x CxMatriz de covarianza de x RxMatriz de correlación de x σ2 xVarianza de x
1 Introducción 1.1 Motivación del Proyecto La separación de señales de audio ha sido un campo muy prolífico de estudio en los últimos 30 años. Aunque para muchas situaciones existen algoritmos conocidos que proporcionan muy buenos resultados, todavía hay casos límite, con condiciones extremas, en los que existe mucho margen de mejora y se debe seguir investigando. El caso en concreto que nos ocupa, es la separación de dos fuentes mono-canal, unísonas y moduladas tanto en amplitud como en frecuencia, esto no es más que dos instrumentos tocando la misma nota mientras ejecutan un vibrato. En este trabajo, nos centraremos exclusivamente en el área del procesamiento digital de señales enfocada a la separación de señales de audio, que se engloba dentro del problema de Separación Ciega de Fuentes (Blind Source Separation, BSS) y empleamos la Factorización de Matrices No Negativas (Nonnegative Matrix Factorization, NMF), especialmente su extensión al uso de tensores (NTF), para resolver el problema. Las técnicas de BSS son aplicables a multitud de problemas como: el tratamiento de imágenes médicas, de vídeo, comunicaciones, etcétera. La separación de fuentes de audio continúa siendo un campo de investigación muy activo desde que se empezara a investigar a principios de la década de 1980. Se han desarrollado diversos métodos de separación, que explotan diferentes características de las señales, entre ellos se puede usar NMF, que factoriza la matriz de un espectrograma en el producto de dos matrices, una llamada de frecuencia y otra de activación, haciendo posible diseñar fácilmente algoritmos eficientes que buscan minimizar la diferencia entre la matriz o tensor original y el producto matricial o tensorial de sus componentes. Para ello, se busca la divergencia y el tipo de actualización óptimos para cada algoritmo. Al mismo tiempo, aporta una reducción de rango, necesaria para descomponer mezclas en sus componentes asociadas a las fuentes. Aplicando los conceptos de NMF a los tensores se pudieron desarrollar modelos más complejos, útiles en muchas aplicaciones, como la separación multi-canal [ 43 ]. Algunos de los casos particulares de NMF, como el convolutivo o NMF invariante en el tiempo, también se han aplicado a los algoritmos de NTF. Estos enfoques, aplicados a la descomposición de mezclas de instrumentos musicales, funcionan cuando determinadas suposiciones son ciertas. Una es que los armónicos espectrales solo se solapan parcialmente. Sin embargo, cuando dos fuentes comparten la misma frecuencia fundamental, la mayoría de los armónicos se solapan, reduciendo así el porcentaje de éxito de los algoritmos basados en NMF en el aprendizaje de matrices únicas. Otra suposición es que todas las matrices temporales y espectrales semánticamente corresponden a notas musicales, formando un diccionario de átomos con sentido musical. Esto no se cumple para instrumentos con fluctuaciones variables en el tiempo. Estos efectos se pueden encontrar en instrumentos como los de cuerda o los de viento-metal cuando tocan con vibrato. En el caso en el que dos instrumentos tocan con vibrato la misma nota, las dos suposiciones anteriores pueden no cumplirse, lo que convierte este escenario en un desafío [ 54 ]. En vez de aumentar el número de plantillas por fuente, Hennequin propone [ 26 ] usar matrices de activación dependientes en frecuencia mediante el uso de un modelo basado en fuente/filtro. Como el vibrato no solo causa modulaciones en frecuencia (FM), si no que también causa modulaciones de amplitud (AM), esto recibe el nombre de espectros de modulación, que pueden ser usados para identificar el patrón de modulación. Estos espectros, a veces, son calculados aplicando la transformada de Fourier a un espectro de magnitud. El spectrograma de 1
2Capítulo 1. Introducción modulación ya ha captado mucha atención en el campo del reconocimiento [ 23 ][ 30 ] y clasificación de voz [ 31 ] [ 39 ]. Barker y Virtanen [ 7 ] fueron los primeros en proponer una modulación representada con tensores para una separación de fuentes monocanal. Esto permite aplicar la factorización al tensor usando la conocida descomposición CANDECOMP/PARAFAC (CP) [25]. 1.2 Objetivo del trabajo El objetivo de este trabajo es el estudio de la separación ciega de fuentes de audio para el caso concreto de una mezcla mono-canal de dos instrumentos que realizan un vibrato mientras tocan ambos la misma nota, para cuya resolución, haremos uso del método conocido como Factorización No Negativa de Tensores (NTF). Partiendo de los resultados presentados por Stöter et al. en [ 53 ], se va a hacer un estudio de los valores de los parámetros alfa y beta para los que se obtienen mejores resultados en la separación. El método de descomposición tensorial empleado en [ 53 ] explota las similitudes en frecuencia. También nos permite hacer uso de las dependencias entre las modulaciones de los intervalos vecinos. Esto tiene ciertas coincidencias con el modelo HR-NMF, del cual se habla en la Sección 2.2.4 y que tiene en cuenta las dependencias en el plano tiempo-frecuencia. El método propuesto en [ 53 ], relaja algunas suposiciones tomadas en HR-NMF con la intención de simplificar el proceso de estimación. Por último, se ha hecho un estudio sobre cómo afectan distintos valores de los parámetros alfa y beta, en función del tipo de fuente, a la calidad de los resultados obtenidos por el método propuesto en [53]. En resumen, el objetivo de este proyecto es estudiar una solución al problema de BSS para un caso muy concreto donde algunas suposiciones esenciales de NMF no se cumplen, estimando los parámetros NMF basándonos en la divergencia AB y utilizando un algoritmo MU, que sea capaz de combinar otros algoritmos multiplicativos existentes, de manera que pueda ser aplicado a distintos casos, ajustando el valor de solo dos parámetros. 1.3 Estructura del proyecto Este documento está dividido en cinco capítulos, cada uno de los cuales consta a su vez de diferentes secciones y subsecciones. Capítulo 1: corresponde a la introducción del trabajo. Incluye las motivaciones que nos han llevado a hacerlo, el objetivo para el que se ha realizado y una explicación de la estructura de la memoria. Capítulo 2: en este capítulo se expone el problema a resolver, se definen de manera amplia conceptos generales como la Separación Ciega de Fuentes y su aplicación a las señales de audio, por ser el motivo de nuestro estudio. Finalmente, se ahonda en el método de la Factorización No Negativa de Matrices, en el que está basado el algoritmo que se ha empleado para resolver el problema de separación. Capítulo 3: se presenta el método utilizado para resolver el problema, propuesto en [ 53 ]. En primer lugar, se presenta el modelo de descomposición del espectrograma expuesto en [ 53 ], Modelo de Destino Recurrente, después el proceso matemático que nos servirá para calcular la separación, la Transformada de Destino Recurrente, también se detalla el algoritmo MU y cómo estimar sus parámetros. Por último, hay una sección dedicada a explicar el funcionamiento de las divergencias AB y su sentido en este trabajo. Capítulo 4: dedicado a las simulaciones realizadas en Matlab ® y a los resultados obtenidos por éstas. En este capítulo se comenta cómo se ha implementado el algoritmo. Capítulo 5: finalmente, se presentan las conclusiones y se proponen líneas futuras de trabajo.
2 Técnicas de Separación. Separación de Señales de Audio. En el campo de la separación de señales, existen numerosas técnicas para separar las señales que se encuentran mezcladas en un conjunto de observaciones. Estas técnicas se agrupan bajo el nombre de Separación Ciega de Fuentes (Blind Source Separation, BSS) y una frecuentemente usada, cuando se cumplen ciertas hipótesis, es la Factorización de Matrices No Negativas (Non-Negative Matrix Factorization, NMF), que es una técnica general de descomposición de observaciones de valores no negativos. 2.1 Separación Ciega de Fuentes (Blind Source Separation, BSS) En el campo de la física y la ingeniería se llama Separación Ciega de Fuentes (BSS) a la recuperación de señales no observadas o fuentes que se encuentran mezcladas en un conjunto de señales conocidas. El hecho de que las señales de origen no se conozcan y no haya información disponible sobre la mezcla, es el motivo por el que se usa el adjetivo ciega en este tipo de separación [ 9 ]. Las técnicas de BSS han sido desarrolladas durante las últimas dos décadas; muchos algoritmos han sido desarrollados y aplicados en una amplia gama de aplicaciones que incluyen ingeniería biomédica, imágenes médicas, reconocimiento de voz, imágenes astronómicas y sistemas de comunicación [41]. El problema de la separación de fuentes fue formulado alrededor de 1982 por Bernard Ans, Jeanny Hérault y Christian Jutten, en el marco del modelado neuronal, para la decodificación del movimiento de vertebrados. El problema también se ha planteado de forma independiente en el marco de las comunicaciones [ 6 ]. Las primeras contribuciones a conferencias de procesamiento de señales y de redes neuronales, aparecieron alrededor de 1985. Inmediatamente, estos documentos llamaron la atención de los investigadores enfocados en el procesamiento de señales, principalmente en Francia y más tarde en Europa. En la comunidad de redes neuronales, el interés surgió mucho más tarde, en 1995, pero de forma muy masiva. Inicialmente, se investigó la separación de fuentes para mezclas lineales instantáneas (sin memoria). A principios de la década de 1990, buena parte de los estudios en el campo ya estaban centrados en las mezclas convolutivas. Finalmente, las mezclas no lineales, excepto unos pocos estudios aislados, se abordaron a finales de la década de 1990 [ 14 ]. 2.1.1 BSS de mezclas lineales e instantáneas El modelo de mezclas lineales e instantáneas es el más simple y el problema más fácil de resolver cuando se tiene un número suficiente de sensores. El problema tratado en este trabajo es un problema de mezclas lineales e instantáneas, pero no contamos con un número suficiente de canales para realizar la separación, lo que complica notablemente el proceso. Dadas J señales desconocidas s1(t),s2(t),...,sJ(t) , a las que llamaremos fuentes en un modelo de mezclas lineales e instantáneas, cada una de las observaciones x1(t),x2(t),...,xM(t) , pueden escribirse como 3
4Capítulo 2. Técnicas de Separación. Separación de Señales de Audio. Figura 2.1 Modelo BSS lineal instantáneo [60]. combinación lineal de las fuentes. El modelo se muestra en la Figura 2.1. xi(t) = J ∑ j=1 ai jsj(t)i=1,...,M(2.1) Los escalares ai j son los coeficientes de mezcla. Agrupando las fuentes en el vector de fuentes s(t) = [s1(t),s2(t),...,sJ(t)]T y las observaciones en el vector de observaciones x(t) = [x1(t),x2(t),...,xM(t)T , podemos escribir la ecuación (2.1) de forma matricial como x(t) = As(t)(2.2) donde A es la matriz de mezcla de dimensión M×J , cuyos elementos son ai j . Tanto s(t) como x(t) pueden contener valores complejos. Para que la mezcla sea invertible, y recuperar el vector de fuentes s(t) , es necesario que el número de filas de A sea mayor o igual que el de columnas, es decir, que el número de observaciones sea mayor igual que el de fuentes ( M≥J ). Este modelo es generativo, es decir, describe cómo los datos observados son generados mediante un proceso de mezcla de las componentes sj. La idea básica de la separación ciega de fuentes es estimar las señales originales a través de una matriz de separación W , de dimensión J×M , siendo la matriz de mezcla A y el vector de fuentes s(t) desconocidos [60]. y(t) = Wx(t)(2.3) Siendo y(t)=[y1,...,yJ] el vector de señales de salida que tratan de estimar las fuentes s(t) . En la Figura 2.1, se representa el proceso de mezcla y separación. Si los sistemas A y W pueden representarse por matrices constantes en el tiempo, estamos ante un problema de BSS con mezcla lineal e instantánea. Puesto que las fuentes y la mezcla son desconocidas, para resolver el problema es necesario utilizar cierta información a priori, en forma de hipótesis. En función de las hipótesis empleadas, se obtienen diferentes criterios de BSS, algunos de los cuales se detallan a continuación. Por otra parte, los algoritmos empleados para optimizar estos criterios pueden ser de procesamiento por bloques o bien algoritmos adaptativos. • Análisis de Componentes Principales (Principal Component Analysis, PCA): transformación del vector de datos x(t) en un vector de señales incorreladas, es decir, asumir como hipótesis que las fuentes son incorreladas. Estas señales incorreladas son llamadas componentes principales y se obtienen mediante descomposición en autovectores y autovalores o descomposición en valores singulares. Un uso común del PCA es la reducción de la dimensión de la matriz de datos y así, solo un conjunto de componentes principales se mantienen para preservar la máxima varianza de los datos. En PCA, la unicidad se consigue imponiendo ortogonalidad en la matriz de transformación [62].
2.1 Separación Ciega de Fuentes (Blind Source Separation, BSS) 5 • Análisis de Componentes Independientes (Independent Component Analysis, ICA): es una generalización del Análisis de Componentes Principales. Siguiendo las definiciones de los pioneros, Jutten y Hérault en 1991 y Common en 1994, podemos suponer que, tanto las variables de mezcla como las componentes independientes tienen una media cero: si esto no es cierto, las variables observables xi siempre se pueden centrar restando la media de la muestra, lo que hace que el modelo tenga una media cero. El punto de partida del modelo ICA es la suposición de que las componentes si son estadísticamente independientes. También se asume que las componentes independientes tienen una distribución no Gaussiana, ya que si hay más de una Gaussiana no hay forma de separar usando independencia, debido a que la mezcla de dos variables Gaussianas independientes puede dar lugar a variables Gaussianas independientes. En el modelo básico no se suponen conocidas estas distribuciones (si se conocen, se simplifica el problema considerablemente). Para simplificar, se puede suponer cuadrada la matriz de mezcla, que es desconocida. En audio, es conocido que las fuentes son independientes y no Gaussianas, así que su aplicación en este campo está muy extendida [28]. Dentro del ámbito más general del Análisis de Variables Latentes, se dice que las componentes independientes (las fuentes) son variables latentes, en el sentido en que no pueden ser directamente observadas. • Análisis de Componentes Escasas (Sparse Component Analysis, SCA): se asume que las fuentes son escasas, es decir, que las fuentes sean cero con frecuencia. El principio básico del SCA consiste en cuatro pasos: 1. Aplicar una transformada lineal de dispersión a la mezcla. Una de las transformadas que se suele usar para la dispersión es la STFT. La transformada se usa para dispersar la representación de las fuentes, así la representación de cada fuente tiene sólo algunos coeficientes significativos. 2. Estimar la matriz de mezcla del gráfico de dispersión. Además del uso del gradiente natural del modelo ICA, un enfoque común en la actualidad es confiar en las técnicas de agrupamiento (clustering), con variantes de K-medias ponderadas. Para que estas técnicas funcionen de forma eficiente, la hipótesis clave es asumir que, como máximo, una fuente contribuye significativamente a cada punto del gráfico de dispersión. En el caso de las fuentes de audio, normalmente se asume que, en el dominio tiempo-frecuencia, la actividad de cada fuente muestra cierta persistencia local dentro de las pequeñas regiones de la distribución tiempo-frecuencia donde son "visibles". 3. Consiste en estimar la representación de las fuentes basándose en la suposición de dispersión. En un escenario libre se ruido, Bofill y Zibulevsky [ 47 ] propusieron una estimación que se puede interpretar como de máxima probabilidad, asumiendo que los coeficientes de las fuentes tengan una distribución laplaciana. 4. Reconstrucción de las fuentes invirtiendo la transformada de dispersión. SCA es muy útil a la hora de aislar ruidos y distorsiones, ya que normalmente estos suelen tener un nivel bajo en proporción a la señal de interés. 2.1.2 BSS para mezclas convolutivas Cuando las fuentes contribuyen a la mezcla con numerosas versiones retardadas se consideran mezclas convolutivas [ 14 ]. Esto puede ocurrir en diversas aplicaciones como el audio. Las mezclas de audio en entornos reales, debido a la reverberación, se consideran siempre mezclas convolutivas (además de variantes en el tiempo). La diferencia entre el modelo de mezcla lineal convolutivo y el instantáneo es que, versiones retrasadas de las fuentes contribuyen a la salida del modelo en momentos dados. Modelo En el modelo convolutivo, la matriz de mezcla se sustituye por un sistema MIMO (múltiples entradas y múltiples salidas) lineal e invariante en el tiempo (LTI) con respuesta impulsiva (A(n))n∈Z . Las señales
6Capítulo 2. Técnicas de Separación. Separación de Señales de Audio. de observación son, por tanto, determinadas por las fuentes conforme al siguiente modelo de convolución multicanal: ∀n∈Zx(n) = ∑ k∈Z A(k)s(n−k).(2.4) Figura 2.2 Modelo de mezcla convolutiva [14]. Una estructura como la mostrada en la Figura 2.2 puede ser invertida con un sistema MIMO-LTI. Recuperar las fuentes es equivalente a encontrar un sistema MIMO-LTI inverso, llamado separador. Si su respuesta impulsiva se denota por (B(n))n∈Z, las salidas separadas son dadas por: ∀n∈Zy(n) = ∑ k∈Z B(k)x(n−k).(2.5) Debido al contexto convolutivo, debemos usar la transformada Z de los sistemas LTI. Para los sistemas de mezcla y separación, con respuestas impulsivas (A(n))n∈Zy(B(n))n∈Zrespectivamente se define: A[z]4 =∑ k∈Z A(k)z−kyB[z]4 =∑ k∈Z B(k)z−k(2.6) Es conveniente introducir el sistema que combina mezcla y separación. Se obtiene de las ecuaciones 2.5 y 2.6, se aprecia que la salida global en el separador recibe: ∀n∈Zy(n) = ∑ k∈Z G(k)s(n−k)(2.7) donde la respuesta impulsiva y la transformada Z del sistema global (G(n))n∈Z son dadas por las ecuaciones: ∀n∈ZG(n) = ∑ k∈Z G(n−k)A(k)yG[z] = B[z]A[z].(2.8) Evolución del modelo y de las técnicas de separación A continuación, vamos a hablar de algunos modelos específicos dentro de la BSS para mezclas convolutivas en el campo del audio [ 55 ]. Antes de introducir dichos modelos, es conveniente aclarar que, el modelo general tiene limitaciones intrínsecas, especialmente para el audio. Primero, el modelado del sistema como respuestas impulsivas entre la localización de cada fuente y la localización de cada micrófono implícitamente asume que, cada fuente emite sonido desde un único punto en el espacio, previniendo así el modelado de fuentes espacialmente difusas. Segundo, a no ser que se conozca información adicional, las fuentes se pueden recuperar, a lo sumo, hasta un filtrado indeterminado. Tercero, el sistema lineal A(t) puede ser invertido solo en determinados escenarios, en los que el número de fuentes es menor que el de micrófonos (J≤I). Debido a estas limitaciones, muchos investigadores propusieron enfocar este problema en el dominio del tiempo-frecuencia por medio de la Transformada Localizada de Fourier (STFT) compleja. En 1998, Cardoso [8] propuso reformular el proceso de mezcla como n(t) = J ∑ j=1 cj(n)(2.9) de forma que el problema de separación de fuente se convirtiera en un problema basado en extraer la contribución cj(t) = [cj1(t),...,cjI(t)]T de cada fuente a la mezcla. Con el tiempo, cj(t) fue llamado imagen
2.1 Separación Ciega de Fuentes (Blind Source Separation, BSS) 7 espacial de la fuente j-ésima [ 58 ]. Con esta reformulación se evitó la indeterminación provocada por el filtrado, uniendo aj(t)ysj(t)en una sola cantidad cj(t) = aj∗sj(n)(2.10) y el modelo general (2.9) se volvió aplicable a fuentes espacialmente difusas, que no puede expresarse como (2.10). Al mismo tiempo, numerosos investigadores, propusieron pasar el problema al dominio del tiempofrecuencia, mediante medias de la STFT compleja. Se reformuló el proceso de mezcla en cada cuadro temporal ny en cada intervalo de frecuencia f, de forma que se expresó como: x(n,f) = J ∑ j=1 cj(n,f),(2.11) En el dominio tiempo-frecuencia, el vector de fuentes se define como s(n,f) = [s1((n,f),...,sJ(n,f)] y el vector de observaciones como x(n,f) = [x1(n,f),...,xm(n,f)] . El modelo de mezcla convolutivo se aproxima bajo la suposición de banda estrecha, por la multiplicación de valores complejos en cada intervalo de frecuencias cj(n,f) = aj(f)sj(n,f),(2.12) donde la transformada de Fourier aj(f) de aj(t) es el llamado vector de mezcla de la fuente j-ésima o en la forma matricial x(n,f) = A(f)s(n,f), donde A(f) = [a1(f),...,aJ(f)] es la llamada matriz de mezcla. La separación de fuentes se reformuló de varias formas, entre ellas, como un problema similar al de agrupación (clustering), por lo que el sonido en un intervalo de tiempo-frecuencia dado debe asignarse a la única o pocas fuentes activas en ese intervalo, y así la separación se hizo viable en escenarios indeterminados, con más fuentes que micrófonos ( J≤I ) [ 61 ]. Otra de estas reformulaciones fue resolver, para cada frecuencia, el problema de la separación, y posteriormente, resolver el de las permutaciones. Mientras que las primeras técnicas de separación de fuentes se basaban en la diversidad espacial, es decir, en la suposición de que las fuentes tienen diferentes direcciones de llegada, el cambio al dominio del tiempo-frecuencia habilitó la explotación de la diversidad espectral, es decir, la suposición de que sus STFTs seguían distintas distribuciones. Esto posibilitó trabajar con mezclas mono-canal y mezclas de fuentes con la misma dirección de llegada. En los últimos años se han propuesto importantes mejoras en las técnicas de separación de fuentes de audio cada vez más adecuadas a las propiedades de las fuentes sonoras y a las especificaciones de las mezclas acústicas: numerosos modelos y sofisticados algoritmos se han desarrollado para incorporar información adicional sobre las fuentes o el entorno de la mezcla para guiar el proceso de separación. Estos modelos rompen un poco con las restricciones propias de BSS, por lo que se engloban bajo el término modelos de separación guiada de fuentes. Dentro de estos algoritmos, aquellos que emplean información sobre el comportamiento general de las fuentes de audio y/o del proceso acústico de mezcla, por ejemplo, "las fuentes están escasamente distribuidas" o "la mezcla fue realizada en exterior", se consideran algoritmos suavemente guiados. Mientras que los algoritmos que aprovechan información específica sobre la mezcla para la separación, como las posiciones de las fuentes o el género musical, se consideran algoritmos fuertemente guiados [55]. Antes de introducir algunos tipos de guía en los algoritmos, es necesario aclarar algunos conceptos comunes de los algoritmos ciegos y guiados. La separación se basa en dos paradigmas de modelado alternativos: la no gaussianidad o no estacionariedad, donde la no estacionariedad se puede manifestar en el tiempo, en frecuencia o en ambos [ 10 ]. Estos paradigmas son perfectamente intercambiables: eligiendo uno de ellos no se restringe el tipo de información que se puede incluir como guía o los escenarios prácticos que pueden ser considerados.
14 Capítulo 2. Técnicas de Separación. Separación de Señales de Audio. Otro avance importante, ha sido el interés de muchos investigadores por explotar la información codificada mediante redundancia y patrones repetitivos en escalas de tiempo muy largas, para optimizar así el uso de la información disponible sobre la duración total de la señal. Huang et al. [ 27 ], usaron el Análisis Robusto de Componentes Principales (RCPA), el cual descompone un espectrograma de entrada como la suma de una matriz de rango bajo y una matriz dispersa, para separar fuentes de batería y melodía, de fuentes de acompañamiento tonal repetitivo. La búsqueda de patrones repetitivos en la música también ha sido explotado por Rafii et al. [ 49 ] mediante la identificación de segmentos repetidos (de un máximo de 40s), modelando y extrayendo a través de un enmascarado en tiempo-frecuencia. 2.2 Factorización No Negativa de Matrices (NMF) y Factorización No Negativa de Tensores (NTF) El método desarrollado por [ 53 ] utiliza NTF, es por eso que en esta sección se va a desarrollar la técnica usada, tomando como referente para todo la sección el libro [13]. 2.2.1 Introducción La Factorización No Negativa de Matrices (Non-Negative Matrix Factorization, NMF), consiste en la descomposición de una matriz como producto de dos o más matrices. La única restricción que exige este método es que todos los coeficientes de las matrices han de ser positivos. Las primeras referencias que se tienen sobre NMF son de Paatero y Tapper en unos trabajos publicados en 1991 [ 45 ], donde se expone el método como una variante de la Factorización Positiva de Matrices (PMF), aunque fue con los trabajos de Lee y Seung publicados en Nature and NIPS [ 37 ] [ 36 ] cuando ganó popularidad, ya que éstos aportaron los primeros algoritmos de aplicación. En la actualidad, NMF es uno de los métodos más usados en BSS. En este problema se ha usado la Factorización No Negativa de Tensores (NTF), método análogo a NMF aplicado a tensores, entendiendo los tensores como matrices de N dimensiones o conjuntos de datos (datasets) indexados por N índices, donde N puede tomar valores mayores que 2 [ 19 ]. Para N=1, un tensor equivale a un escalar y para N=2 a una matriz, en nuestro trabajo usaremos tensores de N=4. Figura 2.6 Tensor de N=3 [13].
2.2 Factorización No Negativa de Matrices (NMF) y Factorización No Negativa de Tensores (NTF) 15 La principal diferencia de NMF respecto a otros métodos de factorización, es la no negatividad de sus coeficientes, la cual es muy importante en la percepción. Muchos datos del mundo real son no negativos y las componentes ocultas solo tienen significado físico cuando son positivas. Esto ocurre en varios campos como el tratamiento de imagen y vídeo, economía y por supuesto en el que nos ocupa, el tratamiento de señales de audio. En este campo, la no negatividad cobra una gran importancia, ya que suele realizarse la separación de audio en el dominio tiempo-frecuencia, usando generalmente la magnitud de las componentes transformadas. NMF es un modelo aditivo, en el que un valor cero representa la ausencia de componentes de la magnitud con la que se esté tratando y un número positivo representa la presencia de alguna componente, lo que permite que cada una de las partes que conforman la suma pueda ser considerada como parte de los datos originales. Gracias a esto, podemos mantener un buen equilibrio entre la interpretabilidad de los datos y la fidelidad estadística de los mismos, hecho que hace al método óptimo para nuestro trabajo. De este tipo de factorización existen varias versiones, podemos hablar de NMF simétrica, convolutiva o multicapa entre otras. Estas diferentes versiones permiten simplificar los modelos en diferentes casos. En nuestro trabajo nos centraremos en el modelo básico, que es el más común y en NMF de Alta Resolución. 2.2.2 Modelo NMF básico El problema básico de NMF se puede expresar de la siguiente manera: dada una matriz de coeficientes no negativos Y∈RJ×T + ( yu≥0 o equivalentemente Y≥0 ) y un rango reducido J(J≤m´ ın(I,T)) , el objetivo es encontrar dos matrices no negativas A= [a1,a2,...,aJ]∈RI×J + y X=BT= [b1,b2,...,bJ]T∈RJ×T + tales que factoricen Ylo mejor posible, eso es: Y=AX +E=ABT+E(2.27) donde la matriz E∈RI×T representa el error aproximado en la descomposición. Las matrices A y X pueden tener diferentes sentidos físicos, dependiendo de la aplicación. En los problemas de BSS, A representa la matriz de mezcla y Xlas señales fuente. En NMF estándar, solo asumimos la no negatividad de las matrices A y X . Al contrario que en los métodos para BSS basados en el Análisis de Componentes Independientes (ICA), aquí no se asume la independencia de las fuentes, en cambio, se introducen otras suposiciones y restricciones para A y/o X posteriormente. Esta simetría en las suposiciones, conduce a una simetría en la factorización: podríamos simplemente escribir YT≈XTAT , esto hace que a menudo el significado de "fuente" y "mezcla" en NMF sea algo arbitrario. El modelo NMF también puede ser representado como una forma especial del modelo bilineal, donde los vectores son no negativos (ver Figura 2.7): Y= J ∑ j=1 aj◦bj+E= J ∑ j=1 ajbT j+E(2.28) donde el símbolo ◦ representa el producto externo de dos vectores. Por lo tanto, podemos construir una representación aproximada de la matriz de datos no negativos Y , como una suma de matrices no negativas de rango unidad ajbT j . El caso en el que esta descomposición sea exacta ( E=0 ), se llama Factorización No Negativa de Rango (Nonnegative Rank Factorization, NRF), este caso en la realidad es muy complejo de conseguir, por lo que en este trabajo se considera la descomposición como una aproximación a la naturaleza, pero no exacta. Aunque NMF se puede aplicar a los problemas de BSS para fuentes y matrices de mezcla no negativas, su aplicación no está limitada a la BSS, de hecho, puede ser usada en diversas aplicaciones. En varias de estas otras aplicaciones se requieren restricciones adicionales para los elementos de las matrices A y/o X , como suavidad, dispersión, simetría y ortogonalidad.
16 Capítulo 2. Técnicas de Separación. Separación de Señales de Audio. Figura 2.7 Modelo NMF bilineal. La aproximación de la matriz de datos no negativos Y∈RI×T + se representa con una suma o una combinación lineal de matrices no negativas de rango unidad Y(j)=aj◦bj= ajbT j∈RI×T +[13]. 2.2.3 Casos particulares de NMF Como se ha expuesto en el inicio de este capítulo, para este tipo de factorización existen varios casos particulares derivados del modelo básico, aunque no se han usado en este trabajo se van a exponer brevemente para tener una idea más amplia del alcance de esta factorización. NMF simétrica Para el caso particular en el que A=B∈RI×J + , la descomposición se denomina NMF simétrica, y puede expresarse como: Y=AAT+E(2.29) Si existe la simetría exacta (cuando E=0 ), se dice que la matriz no negativa Y∈RI×I + es completamente positiva (CP). NMF semi-ortogonal Se define igual que el modelo básico: Y=AX +E=ABT+E,(2.30) la diferencia radica en que, además de la restricción de no negatividad de las matrices A y X , se añade la de ortogonalidad: ATA=IjoXXT=Ij. Semi-NMF En algunas aplicaciones, los datos de entrada observados no tienen signo: Y=Y±∈RI×T . Esto nos permite relajar las restricciones con respecto a la no negatividad de las matrices. Así, Semi-NMF se puede expresar como: Y±=A±X++E,or Y±=A+X±+E,(2.31) Tri-NMF También conocida como NMF de tres factores. Es un caso particular de NMF multicapa, en el que entra en juego una nueva matriz, quedando el modelo de la siguiente forma: Y=ASX +E,(2.32) donde las restricciones de no negatividad pueden ser impuestas a todas o solo a las matrices de factorización elegidas: A∈RI×J,S∈RJ×R, y/o X∈RR×T. Si no se añaden restricciones adicionales en la factorización, este modelo se puede reducir al estándar con la transformación A←AS o X←SX . Sin embargo, Tri-NMF no es equivalente al modelo básico si aplicamos restricciones o condiciones especiales, así aparecen varios modelos como: Tri-NMF Ortogonal, Tri-NMF No Suave, Filtrado NMF o la Descomposición CGR/CUR. NMF con offset El objetivo es eliminar el valor de referencia o el nivel de continua de la matriz Y , usando un modelo NMF ligeramente modificado: Y=AX +a0lT+E,(2.33)
2.2 Factorización No Negativa de Matrices (NMF) y Factorización No Negativa de Tensores (NTF) 17 donde l∈RTes un vector todo unos y a0∈RI +es un vector escogido para que la matriz Xtenga la tierra a cero. El término Y0=a0lT denota el offset, que junto a la restricción de no negatividad, a menudo asegura la poca dispersión de las matrices factorizadas. El papel principal de este término es absorber los valores constantes de la matriz de datos. Figura 2.8 Esquema NMF con offset [13]. NMF multicapa En este caso, la matriz A se remplaza por un conjunto de matrices en cascada (capas). El modelo se describe como (ver Figura 2.9): Y=A(1)A(2)···A(L)X+E,(2.34) Como el modelo es lineal, todas las matrices pueden ser fusionadas en una sola matriz A , si no se han impuesto restricciones especiales a las matrices que conforman las capas. Este modelo se puede utilizar para mejorar considerablemente el rendimiento del modelo NMF estándar gracias a la estructura distribuida en capas y al alivio del problema de los mínimos locales. Figura 2.9 Esquema NMF multicapa [13]. NMF simultánea En NMF Simultánea se tienen dos o más matrices de entrada de datos vinculadas (llamadas Y1 e Y2 ) y el objetivo es descomponerlas en matrices de factorización no negativas de forma que una de las matrices de factorización sea común a ambas, por ejemplo: Y1=A1X+E1, Y2=A2X+E2,(2.35) NMF proyectiva Un modelo NMF Proyectivo puede formularse como la estimación de una matriz dispersa y no negativa W∈RI×J +que satisfaga la ecuación matricial: Y=WWTY+E,(2.36) En una forma general no simétrica, NMF proyectiva implica la estimación de 2 matrices no negativas: A∈RI×J +yB∈RI×J +en el modelo (ver Figura 2.10): Y=ABTY+E(2.37)
18 Capítulo 2. Técnicas de Separación. Separación de Señales de Audio. Figura 2.10 Esquema NMF Proyectiva [13]. NMF convexa En NMF convexa se asume que los vectores base A= [a1,a2,...,aJ] tienen como restricción ser combinaciones convexas de la matriz de datos de entrada Y= [y1,y2,...,yT]. Es decir: aj= T ∑ t=1 wt jyt=YwjoA=YW,(2.38) donde W∈RT×J + y X=BT∈RJ×T + . El modelo NMF Convexo puede ser escrito de forma matricial como: Y=YWX +E(2.39) aplicando el operador de transposición obtenemos: YT=XTWTYT+ET(2.40) En la Figura 2.11 se puede apreciar que, NMF convexa se puede representar de una forma similar a NMF Proyectiva. Figura 2.11 Esquema NMF Convexa [13]. Kernel NMF Considere un mapeo yt→φ(yt) o Y→φ(Y) = [φ(y1),φ(y2),...,φ(yT)] , así Kernel NMF puede definirse como: φ(Y)∼ =φ(Y)WBT.(2.41) NMF convolutiva Este caso es una generalización de NMF básica, donde se trabaja con versiones de la matriz X desplazadas horizontalmente. Matemáticamente, se puede expresar el modelo como: Y= P−1 ∑ p=0 Ap p→ X+E,(2.42) donde X=0→ X , representa la matriz de fuentes primaria, y los términos p→ X representan los vectores de la matriz X desplazados p columnas. Las columnas desplazadas hacia fuera son fijadas a cero, tal como puede
2.2 Factorización No Negativa de Matrices (NMF) y Factorización No Negativa de Tensores (NTF) 19 verse en el siguiente ejemplo: X=135 2461→ X=013 0242→ X=001 0021← X=350 460(2.43) En la Figura 2.12 queda reflejado este modelo donde el operador Sp=T1 denota el desplazamiento horizontal. Figura 2.12 Esquema NMF Convolutiva [13]. NMF superpuesta Se trata de una extensión del caso convolutivo, mientras que en éste se realiza un desplazamiento horizontal de las columnas, en el caso de NMF superpuesta, se realizan diferentes transformaciones, como desplazamientos verticales, muy útiles por ejemplo, a la hora de trabajar con espectrogramas. La expresión matemática de este modelo varía en función de las transformaciones realizadas sobre las filas y columnas de la matriz X, por ejemplo, podría expresarse como: Y≈ P ∑ p=0 (→p X)TAT p= P ∑ p=0 (XTp)TAT p= P ∑ p=0 TT pXTAT p,(2.44) Figura 2.13 Esquema NMF Superpuesta [13].
20 Capítulo 2. Técnicas de Separación. Separación de Señales de Audio. 2.2.4 NMF de alta resolución (High Resolution NMF, HR-NMF) Debido a que éste es el modelo que más se asemeja al que se desarrolla en [ 53 ], se va a introducir en esta sección. Este modelo fue presentado en 2011 por Roland Badeau en [ 3 ], se generalizó para mezclas multicanal en [ 5 ] y se demostró que proporciona un rendimiento considerablemente mejor para la separación de fuentes que los modelos anteriores en [ 48 ]. Aunque algunas aproximaciones variacionales fueron introducidas en [ 4 ] para reducir su complejidad, estos algoritmos son, a menudo, muy exigentes para aplicaciones prácticas. Según el trabajo de Roland Badeau y A.Dremeau [ 3 ], HR-NMF es un modelo que permite superar las limitaciones de resolución espectral que tiene el modelo NMF, teniendo en cuenta tanto la fase, como las correlaciones locales en cada banda de frecuencia. Este modelo, se estima implementando de forma recursiva un algoritmo EM [16], que se aplica de forma satisfactoria a los problemas de separación de fuentes. A continuación, vamos a introducir el modelo presentado en [ 3 ]. El modelo de mezcla x(f,t) , se define para todas las frecuencias 1≤f≤F y tiempos 1≤t≤T como la suma de K componentes ocultas ck(f,t) más un ruido blanco n(f,t)∼N(0,σ2): x(f,t) = n(f,t)+ K ∑ k=1 ck(f,t)(2.45) donde •ck(f,t) = ∑P(k,f) p=1a(p,k,f)ck(f,t−p)+bk(f,t) se obtiene filtrando de forma autorregresiva una señal no estacionaria bk(f,t)( y P(k,f)∈Ntal que a(P(k,f),k,f)6=0), •bk(f,t)∼N(0,υk(f,t)) donde υk(f,t)se define como υk(f,t) = w(k,f)h(h,t),(2.46) con w(k,f)≥0yh(k,t)≥0, •Los procesos nyb1...bKson mutuamente independientes. Dado que N denota tanto la distribución normal real como la circular compleja, el modelo (2.45) puede tomar valores reales o complejos. Además, para instantes anteriores al inicial se asume ck(f,t)∼N(0,1) y no se dispone de las observaciones x(f,t). Los parámetros a estimar son σ2,a(p,k,f),w(k,f)yh(k,t). Este modelo en el dominio tiempo-frecuencia, ha servido para generalizar algunos modelos muy usados en varios sectores del procesado de señal: • Si σ2=0 y ∀k,f,P(k,f) = 0 , el modelo (2.45) se convierte en x(f,t) = ∑K k=1bk(f,t) , así x(f,t)∼ N(0,ˆ Vf t ) , donde ˆ V se define por NMF como ˆ V=WH con Wf k =w(k,f) y Hkt =h(k,t) . La estimación de máxima verosimilitud de W y H es entonces equivalente a la minimización de la divergencia de Itakura-Saito entre la matriz modelo ˆ V y el espectrograma V (donde Vft =|x(f,t)|2 ), por ello este modelo toma el nombre de IS-NMF [20]. • Para valores conocidos de k y f , si ∀t,h(k,t) = 1 , entonces ck(f,t) es un proceso autorregresivo de orden P(k,f). • Para valores conocidos de k y f , si P(k,f)≥1 y ∀t≥P(k,f) + 1 , h(k,t) = 0 , entonces ck(f,t) se puede escribir como ck(f,t) = ∑P(k,f) p=1αpzt p donde z1...zP(k,f) son las raíces del polinomio zP(k,f)− ∑P(k,f) p=1a(p,k,f)zP(k,f)−p . Esto corresponde al Modelo Exponencial Sinusoidal (ESM), frecuentemente usado en análisis espectral HR de series temporales [2]. Por todas estas razones, nos referimos al modelo (2.45) como HR-NMF.
2.2 Factorización No Negativa de Matrices (NMF) y Factorización No Negativa de Tensores (NTF) 21 2.2.5 Estimación de los parámetros NMF Estimación basada en la medida Para estimar las matrices de factorización A y X en el estándar de NMF, es necesario considerar alguna medida de similitud para cuantificar la diferencia entre la matriz de datos Y y su aproximación NMF ˆ Y=AX . La elección de la medida de similitud (también llamada distancia o divergencia), depende mayormente de la distribución de probabilidad de las señales estimadas o de las componentes y la estructura de los datos o del ruido. Una vez elegida la distancia, la función de coste será la función dada por dicha distancia, y el objetivo, obtener un algoritmo que permita minimizar dicho coste. A continuación, se expondrán algunas de las distancias más generalizadas, incluyendo un sencillo algoritmo de aplicación. Norma de Frobenius Esta norma, que toma nombre del matemático alemán Ferdinand Georg Frobenius es en la que se basa la medida más simple y comúnmente usada: DF(YkAX) = 1 2kY−AX k2 F(2.47) Se debe resaltar que la función de costes anterior es convexa con respecto, tanto a los elementos de la matriz A, como a los de la matriz X, no para ambos, si intentamos optimizar los dos se pierde la convexidad. La minimización de la función de costes dada por la norma de Frobenius, da lugar al algoritmo de minimización de los mínimos cuadrados (en inglés, "Alternating Least Squares", ALS), que es uno de los métodos de optimización más conocidos. A continuación, se describe este algoritmo con unos sencillos pasos: 1. Inicializar Ade forma aleatoria o usando una estrategia determinista específica. 2. Estimar Xde la ecuación matricial ATAX =ATYresolviendo m´ ın XDF(YkAX) = 1 2kY−AX k2 Fsiendo Afija. (2.48) 3. Fijar a cero, o a algún valor positivo pequeño, todos los elementos negativos de X. 4. Estimar Ade la ecuación matricial XXTAT=XYTresolviendo m´ ın ADF(YkAX) = 1 2kYT−XTATk2 Fsiendo Xfija. (2.49) 5. Fijar a cero. o a algún valor positivo pequeño, ε, todos los elementos negativos de A. Observando estos pasos, se hace evidente la sencillez de este algoritmo, que podemos resumir mediante las siguientes ecuaciones: X→ {ε,(ATA)−1ATY}= [A†Y]+, A→ {ε,YXT(XXT)−1}= [YX†]+,(2.50) donde A† es la inversa de Moore-Penrose de A y ε es una constante pequeña (típicamente 10−16 ), que se usa para forzar a las entradas a ser positivas. Varias restricciones adicionales pueden imponerse sobre AyX . Por último, decir que se ha considerado oportuno exponer aquí el algoritmo ALS debido a que es tomado como un enfoque básico en algunos métodos actuales, siendo utilizado frecuentemente en inicializaciones previas a la aplicación de otros algoritmos más complejos. Tiene la ventaja de que su implementación es bastante sencilla, aunque no garantiza la convergencia hacia mínimos globales y sus soluciones no son muy precisas.
22 Capítulo 2. Técnicas de Separación. Separación de Señales de Audio. Divergencia de Kullback-Leibler (KL) Otra función de costes popular en NMF es la divergencia Kullback-Leibler, también llamada divergencia de la información o divergencia-I, es un caso especial de la llamada divergencia de Bregman y se define como: DKL(YkAX) = ∑ it yit ln yit [AX]it −yit + [AX]it .(2.51) Esta medida fue introducida en 1951 por Solomon Kullback yRichard Leibler, como una divergencia dirigida entre dos distribuciones [ 35 ]. Actualmente es muy usada en estadística y está muy relacionada con el método de ajuste de distribuciones por máxima verosimilitud. Divergencia de Itakura-Saito (IS) La divergencia IS, al igual que la divergencia KL, es una extensión de la divergencia de Bregman, muy usada en la estimación de los parámetros de NMF y que puede definirse de la siguiente forma: DIS(YkAX) = ∑ it ln [AX]it yit +yit [AX]it −1.(2.52) Fue propuesta por Fumitada Itakura yShuzo Saito, cuando trabajaban en la compañía NTT (Nipppon Telegraph and Telephone) [ 29 ]. En la actualidad es una de las medidas más usadas en las técnicas de separación ciega de fuentes de audio. 2.2.6 Otros aspectos a considerar en NMF Inicialización de parámetros La solución y la convergencia dada por los algoritmos NMF, normalmente dependen mucho de las condiciones iniciales, es decir, sus valores iniciales supuestos, especialmente en un contexto multivariable. Por ello, es importante tener formas eficientes y consistentes de inicializar las matrices A y/o X . En otras palabras, la eficiencia de la mayoría de las estrategias NMF se ve claramente afectada por la selección de las matrices iniciales [ 13 ]. Inicializaciones pobres nos llevan a convergencias lentas y, en algunos casos, incluso a soluciones incorrectas o irrelevantes. Por otro lado, una determinada inicialización no se comporta igual para distintas matrices de entrada de datos, mientras que para unas puede aportar buenos resultados para otras puede resultar pobre, por ejemplo, este problema puede volverse bastante complejo cuando se trata con matrices factorizadas a las que se han impuesto ciertas limitaciones. Resulta útil, para evaluar la eficacia de la estrategia de inicialización y del propio algoritmo, realizar análisis de incertidumbre como simulaciones de Monte Carlo. A continuación, a modo de ejemplo, se presentan una serie de pasos, que proporcionarían una estrategia de inicialización robusta: 1. Generar de forma iterativa un número R de matrices iniciales A y X (normalmente con 10-20 iteraciones es suficiente), mediante una inicialización aleatoria o cualquier algoritmo sencillo como el ALS. 2. Ejecutar algún algoritmo NMF específico para cada par de matrices iniciales generalizadas y con un número fijado de iteraciones (típicamente 10-20). Como resultado, se obtienen R estimaciones de las matrices A(r)yXr. 3. Seleccionar de entre las estimaciones anteriores, aquella pareja que proporcione el menor valor para la función de coste, denotadas como A(rmin)yX(rmin). Criterios de parada Hay numerosos criterios de parada para los algoritmos iterativos usados en NMF, a continuación vamos a exponer algunos: • Que la función de costes alcance el valor cero (o cercano a cero) o un valor por debajo de un umbral dado ε, por ejemplo, D(k) F(Ykˆ Y(k)) =kˆ Y(k)−ˆ Y(k+1)k2 F≤ε,(2.53) o |D(k) F−D(k−1) F D(k) F ≤ε.(2.54)
2.2 Factorización No Negativa de Matrices (NMF) y Factorización No Negativa de Tensores (NTF) 23 • Que haya cambios insignificantes o que no haya cambios en las iteraciones sucesivas de las matrices A yX. •El número de iteraciones alcanza o supera el número máximo de iteraciones predefinido. En la práctica, las iteraciones continúan hasta que se satisface alguna combinación de criterios de parada. Ambigüedades Como se vio en el apartado 2.2.5, generalmente la estimación en NMF se realiza mediante la minimización de una o varias funciones objetivos. Sin embargo, en general, estas minimizaciones no garantizan una solución única. Incluso la función cuadrática con respecto a ambos conjuntos de argumentos {A} y {X} puede tener varios mínimos locales, lo que hace que los algoritmos NMF puedan sufrir indeterminaciones rotacionales (ambigüedades). Debido a estas indeterminaciones o ambigüedades, es posible que la solución a NMF no sea única, por ejemplo, considerando la siguiente ecuación cuadrática: DF(YkAX) =kY−AX k2 F=kY−AR−1−RX k2 F=kY−˜ A˜ Xk2 F,(2.55) donde la matriz rotacional R , ha de escogerse de manera que las matrices transformadas ˜ A6=A y ˜ X6=X sean no negativas. Es importante resaltar, que la inversa de una matriz no negativa es no negativa, si y solo si, esta es una matriz de permutación generalizada. Una matriz de permutación generalizada es aquella que tiene un solo elemento positivo distinto de cero en cada fila y en cada columna. Si asumimos que R≥0 y R−1≥0 , las cuales son condiciones suficientes para la no negatividad de las matrices AR−1 y RX , entonces R tiene que ser una matriz de permutación generalizada, es decir, R se puede expresar como el producto de una matriz diagonal no singular definida positiva y una matriz de permutación. Si las matrices originales X y A están suficientemente dispersas solo una matriz de permutación P=R puede satisfacer la restricción de no negatividad de cualquier matriz de transformación y obtener una NMF única. Para ilustrar la indeterminación rotacional descrita, se expone el siguiente ejemplo, dadas las siguientes matrices de mezcla y fuente: A= 3 2 7 2 ,X= x1(t) x2(t) su producto da como resultado: Y= y1(t) y2(t) =AX 3x1(t)+2x2(t) 7x1(t)+2x2(t) No obstante, es evidente que exiten otras descomposiciones no negativas que dan el mismo resultado, por ejemplo: Y= y1(t) y2(t) =˜ A˜ X= 0 1 4 1 x1(t) 3x1(t)+x2(t) , donde: ˜ A= 0 1 4 1 ,˜ X= x1(t) 3x1(t)+x2(t) son nuevos componentes no negativos que no provienen de las indeterminaciones de permutación o escalado. Sin embargo, incorporando alguna medida de dispersión o suavidad a la función objetivo, es suficiente para resolver el problema NMF de forma única. Cuando no hay información a priori disponible, debemos normalizar las columnas de A y/o las filas de X para ayudar a mitigar los efectos de las indeterminaciones de rotación. Dicha normalización se suele hacer escalando las columnas ajde A= [a1,...,aJ]de la siguiente manera: A←ADA,donde DA=diagka1k−1 p,ka2k−1 p,...,kaJk−1 p,p∈[0,∞) Varios experimentos demuestran que, los mejores resultados se obtienen para p=1 . Por otro lado, para evitar las indeterminaciones rotacionales, las filas de X deben ser dispersas o tener nivel de referencia cero, como
30 Capítulo 3. Modelo de Destino Recurrente (Common Fate Model, CFM) convencional de Factorización No Negativa de Matrices (NMF, ver Sección 2.2) excepto que se aplica a la CFT en vez de a la STFT, así que la factorización que se va a utilizar es (3.3). Figura 3.3 Modelo de Destino Recurrente, CFM [53]. Básicamente, lo que se busca es estimar los parámetros {Aj,Hj} para que el módulo de la CFT, elevado a la potencia α , esté lo más cerca posible de ∑jPα j , con alguna función de coste particular como criterio de ajuste de datos, llamada β−divergencia que incluye casos especiales como la Euclídea, Kullback-Leibler e Itakura-Saito [ 21 ]. Como es común en los modelos no negativos, cada parámetro es actualizado por turno, mientras que otros se dejan fijos. Se proporcionan las actualizaciones multiplicativas en el Algoritmo 1. Después de varias iteraciones, los parámetros pueden usarse en (3.2) para separar las fuentes. Algoritmo 1: Ajuste de los parámetros NMF de la CFM no negativa(3.3) [53] Con υα=|x|αy usando siempre los parámetros actualizados para calcular ˆ Pα(a,b,f,t) = ∑J j=1Aj(a,b,f)Hj(t), iterar: Aj(a,b,f)←Aj(a,b,f)∑tυα(a,b,f,t)ˆ Pα(a,b,f,t)·(β−2)Hj(t) ∑tˆ Pα(a,b,f,t)·(β−1)Hj(t) Hj(t)←Hj(t)∑a,b,fυα(a,b,f,t)ˆ Pα(a,b,f,t)·(β−2)Aj(a,b,f) ∑a,b,fˆ Pα(a,b,f,t)·(β−1)Aj(a,b,f) 3.3 Estudio de las Alfa-Beta divergencias Una vez estudiado e implementado el modelo, en la parte práctica de este trabajo, se hizo un estudio de cómo afectaban los valores de los parámetros α y β a los resultados de la separación, es decir, cómo afectan a la Relación Señal a Ruido (SDR), a la Relación Señal a Artefacto (SAR) y a la Relación Señal a Interferencia (SIR). No vamos a entrar en los detalles prácticos en esta sección ya que se exponen en el siguiente capítulo, pero vamos a introducir el concepto de Alfa-Beta divergencia en la separación de fuentes de audio. Los parámetros αyβestán relacionados con el concepto de Alfa-Beta divergencia que fue propuesto en [12]. 3.3.1 Definición de Alfa-Beta divergencia Sean las matrices positivas P y Q de dimensiones I×T y, pit y qit las entradas de dichas matrices, se define la divergencia alfa-beta o divergencia AB, como la medida de similitud entre ambas matrices, de la siguiente manera: D(α,β) AB (PkQ) = −1 αβ ∑ it pα it qβ it −α α+βpα+β it −β α+βqα+β it para α,β,α+β6=0(3.4)
3.3 Estudio de las Alfa-Beta divergencias 31 O de forma equivalente: Dα,λ−α AB (PkQ) = 1 (α−λ)α∑ it pλ it qλ−α it −α λpλ it −λ−α λqλ it para α6=0,α6=λ,λ=α+β6=0. (3.5) Para evitar indeterminaciones o singularidades para algunos valores de los parámetros, la divergencia AB puede extenderse por continuidad (aplicando la fórmula de l’Hôpital) para cubrir todos los valores de α,β∈R, así la divergencia AB puede expresarse de una forma más explícita: D(α,β) AB (PkQ) = ∑ it d(α,β) AB (pit ,qit ),(3.6) donde d(α,β) AB (pit qit ) = −1 αβ pα it qβ it −α (α+β)pα+β it −β (α+β)qα+β it para α,β,α+β6=0 1 α2pα it ln pα it qα it −pα it +qα it para α6=0,β=0 1 α2ln qα it pα it +qα it pα it −1−1para α=−β6=0 1 β2qβ it ln qβ it pβ it −qβ it +pβ it para α=0,β6=0 1 2(ln pit −lnqit)2para α,β=0 (3.7) Sustituyendo los parámetros α y β con los valores adecuados, pueden obtenerse otras distancias conocidas, como la divergencia de Kullback-Leibler (para α=1 y β=0 ), o la divergencia de Itakura-Saito (para α=1 yβ=−1), entre otras D(1,0) AB (PkQ) = DKL (PkQ) = pit ln pit qit −pit +qit D(1,−1) AB (PkQ) = DIS (PkQ) = ln pit qit +qit pit −1(3.8) Por otro lado, particularizando para ciertos valores, se obtienen otras divergencias como la Alfa (para α+β=1 ), o la divergencia Beta (para α=0 ). Por lo tanto, podemos decir que la divergencia AB es una medida de similitud general, a partir de la cual es posible obtener muchas de las divergencias más utilizadas. 3.3.2 Propiedades Considerando el operador P[r] como la transformación que eleva todos los elementos de la matriz P al valor de r, es decir, pr it , se definen las siguientes propiedades de la divergencia AB: - Dualidad: esta propiedad implica que una permutación en los parámetros, provoca una permutación en las matrices. D(α,β) AB (PkQ) = D(α,β) AB (QkP)(3.9) - Inversión: un cambio de signo en los parámetros, se traduce en la inversión de los elementos de las matrices: D(−α,−β) AB (PkQ) = D(α,β) AB P[−1]kQ[−1](3.10) - Escalado de parámetros:las propiedades anteriores pueden ser consideradas como casos particulares de la propiedad de escalado de los parámetros α y β por un factor común ω∈R\{0} . La divergencia cuyos parámetros han sido re-escalados es proporcional a la divergencia original con ambos argumentos elevados al factor común, es decir: D(ωα,ωβ) AB (PkQ) = 1 ω2D(α,β) AB P·[ω]kQ·[ω].(3.11) Esta propiedad puede verse como un "zoom-in" a los argumentos de P y Q cuando ω<1 . Dicho "zoom" da más relevancia a los valores pequeños frente a los mayores. Al contrario, cuando ω>1 se
32 Capítulo 3. Modelo de Destino Recurrente (Common Fate Model, CFM) produce un efecto "zoom-out" donde los valores pequeños pierden relevancia en detrimento de los valores grandes (ver Figura 3.4). Figura 3.4 Ilustración gráfica de las propiedades de inversión y dualidad en la divergencia-AB. En el plano α−β están indicados como casos importantes divergencias particulares con puntos y líneas, especialmente la divergencia Kullback-Leibler DKL , Distancia de Hellinger DH , distancia Euclídea DE, distancia de Itakura-Saito DIS, Alfa-divergencia Dα A, y Beta-divergencia D(β) B[12]. Todas estas propiedades permiten reformular la definición de la divergencia AB y expresarla de distintas formas, por ejemplo en términos de otras divergencias combinadas con zooms de los parámetros. 3.3.3 Justifcación del estudio El hecho expuesto anteriormente y en [ 12 ], que implica que la divergencia AB incluye varias de las divergencias más populares en función del valor de los parámetros α y β , hace que esta divergencia sea muy interesante para examinar la eficiencia de nuestro algoritmo ya que da una gran versatilidad a la solución. También resulta muy interesante para resolver problemas de NMF ya que la divergencia AB cumple la restricción de no negatividad, al cumplirla las divergencias que abarca. Por otro lado, la divergencia AB es notablemente más robusta que otras, frente a los errores y al ruido, gracias al uso de los hiperparámetros α y β . El modelo factorizado puede verse como una función vectorial de una serie de parámetros θ , donde cada uno de sus elementos qit (θ)>0 es no negativo para un rango de parámetros determinado. Entonces, el estimador ˆ θ entre dos medidas discretas positivas P y Q , para la
3.3 Estudio de las Alfa-Beta divergencias 33 divergencia AB, es una solución de la ecuación: ∂D(α,β) AB (PkQ) ∂θ =−∑ it ∂qit ∂θ qα+β−1 it ln1−α(pit /qit ) = 0,(3.12) donde la función ln1−αes el logaritmo deformado, definido como: ln1−α(z) = (zα−1 α,si α6=0 ln(z),si α=0 Atendiendo a la ecuación 3.12, puede verse la influencia que tiene cada parámetro sobre la estimación. En el caso de α , se controlan los valores individuales de los cocientes pit /qit , lo que puede ser interpretado como un "zoom", donde se dará más importancia a valores altos, en el caso de α>1 , o por el contrario a valores pequeños, si α<1 . De igual forma, el parámetro β puede controlar la influencia del cociente pit /qit , normalmente se buscan valores que permitan un buen compromiso entre la robustez, para valores β>1 , y la eficiencia, para β=0 . En definitiva, el parámetro α , controla la influencia de los cocientes en el estimador, mientras que β , controla la ponderación de dichos cocientes dependiendo de los valores que mejor se ajusten al modelo. Figura 3.5 Ilustración gráfica de cómo los parámetros establecidos α y β pueden controlar la influencia de los ratios individuales pit /qit . La línea de puntos y rayas ( α+β=1 ) muestra la región donde el factor de ponderación multiplicativo qα+β−1 it en la ecuación de estimación es constante y unitario. La línea de rayas ( α=1 ) muestra la región donde el orden del logaritmo deformado de pit /qit es constante e igual al de la divergencia Kullback-Leibler estándar [12]. Resumiendo, la elección de la divergencia AB para este estudio se justifica gracias a su no negatividad, lo que permite que pueda usarse para resolver el problema de NMF, su versatilidad y su eficiencia, así como su robustez frente a ruidos y errores.
4 Simulaciones En este capítulo, se presenta todo el trabajo práctico realizado en Matlab ® para implementar el modelo expuesto en el Capítulo 3 y los resultados de las diversas simulaciones realizadas, con el fin de exponer el trabajo elaborado. 4.1 Datos de entrada Para poder comparar nuestras simulaciones con las de [ 53 ], se ha usado el mismo dataset, que consta de 10 mezclas diferentes de 5 instrumentos: viola, flauta travesera, violonchelo, saxo tenor y corno inglés. Todas las mezclas siguen la misma estructura: Instrumento A ( 3s ) | Instrumento B ( 3s ) | Mezcla de Instrumento A + Instrumento B (3s). Dicha estructura se aprecia en la Figura 4.1. Figura 4.1 Estructura de las pistas de audio de entrada del algoritmo. Todas las pistas de audio se han generado mediante un software de edición de audio 1 , renderizando la nota C4 (Do central) correspondiente a los 261.63 Hz . Los datos se han codificado en archivos WAV a 44.1kHz y 16 bits. Todas las pistas son mono canal. 1VIENNA SYMPHONIC LIBRARY (https:// vsl.co.at) 35
36 Capítulo 4. Simulaciones 4.2 Algoritmo paso a paso 1. Para empezar, es esencial pasar las pistas de audio a matrices de datos para poder trabajar con ellos. De esto se encarga la función audioread. 2. Una vez obtenida la matriz de datos de entrada, se calcula su STFT mediante la función spectrogram. La STFT se calcula con una ventana hamming de 1024 muestras de longitud, un solape del 50% y una NFFT de 1024 puntos, siendo estos los puntos que se usan para calcular la Transformada Discreta de Fourier (DFT), necesaria para el cálculo de la STFT, como se explica en la Sección 3.2.1. 3. Con la STFT calculada y guardada en una variable a la que hemos llamado XSTFT, llamamos a la función CFM, que ha de ser creada anteriormente. Esta función recibe la STFT de los datos de entrada, la frecuencia de muestreo que devuelve la función audioread y los valores de α y β , y devuelve la SAR (Relación Señal a Artefacto), SDR (Relación Señal a Distorsión) y SIR (Relación Señal a Interferencia). A partir de este punto, se detalla paso a paso lo que realiza la función CFM. 4. Lo primero que hace la función, es calcular el número de filas y columnas que tendrá la matriz X , que corresponden a las filas y columnas de XSTFT. También se inicializan las variables que marcan el número de filas y columnas que tendrá cada parche, en nuestro caso Na=4 y Nb=64 respectivamente. Por último, se calcula el número de parches que tendrá cada fila y cada columna y se guarda en las variables NfyNt. 5. Se calcula el tensor G , para ello se anidan dos bucles for para recorrer la matriz XSTFT e ir guardando en G parches de 4×64 , teniendo en cuenta que los parches tienen un 50% de solape, se obtiene un tensor de dimensión 4×64 ×256 ×25. 6. Una vez obtenido el tensor G , se calcula la matriz X realizando la 2D-DFT a cada parche de G . A continuación y siguiendo el Algoritmo 1, se inicializa V como el valor absoluto de X elevado a α , esto se ha hecho para simplificar los siguientes bucles iterativos. 7. Siguiendo lo dictado por el Algoritmo 1, tenemos que inicializar el tensor Paj (en el algoritmo equivale a ˆ Pα ). Este tensor es de tamaño Na×Nb×Nf×Nt×J donde J es una variable que almacena el número de fuentes, en nuestro caso J=2 , la función de este tensor se explica en la Sección 3.2.4 y no es más que una variable de un modelo de factorización para calcular las densidades de modulación de las fuentes. Para poder dar valores al tensor, primero tenemos que inicializar de forma aleatoria, con la función randn, el tensor A (en el algoritmo Aj(a,b,f) ), de tamaño Na×Nb×Nf×J y la matriz H (en el algoritmo Hj(t) ). Una vez hecho esto, anidamos dos bucles for, el primero se recorre tantas veces como fuentes tengamos y el segundo se recorre Nt veces. En cada iteración se realiza el producto de A por una entrada de Hy se almacena en una entrada de Paj. 8. En este paso, se calcula la densidad de modulación, que aparece en la Ecuación (3.1) como Pα , en nuestro algoritmo se almacenan en la variable Pa . Para este cálculo se ejecuta el Algoritmo 1, lo que significa que estamos ajustando los parámetros NMF del Modelo de Destino Recurrente, esto puede hacerse en este paso ya que previamente hemos inicializado todas las variables. El algoritmo se ejecuta dentro de un bucle iterativo, el cual se itera 100 veces, esta condición de parada es la recomendada por [ 53 ]. Una vez se llegue a la condición de parada, se toman los valores de la última actualización de Paj y se guardan en el tensor Pa como la suma de las Paj para ambas fuentes, teniendo así Pa una dimensión menos que Paj. Llegados a este punto, podemos decir que hemos acabado con el modelo de factorización y a partir de ahora empezaremos con el proceso propio de separación de fuentes. 9. Para la separación de fuentes lo que tenemos que hacer es implementar la Ecuación (3.2) en Matlab ®. Los resultados los guardaremos en una variable llamada S. 10. Lo siguiente que queremos, es obtener las formas de ondas correspondientes a los datos que hemos obtenido tras la separación. Lo primero que tendremos que hacer, será calcular la 2D-DFT inversa a cada parche de Sque guardaremos en el tensor iS.
4.3 Evaluación de los resultados 37 Una vez hecho eso, en un bucle tendremos que pasar el tensor iS de cuatro dimensiones, cuyos parches tienen un solape del 50% , a una matriz de dos dimensiones sin solape. Esto se hace simplemente recorriendo el tensor iS con los valores adecuados y guardando los datos en el tensor s , que es un tensor porque contiene las matrices de las dos fuentes. Por último, tendremos que calcular la STFT inversa de los datos correspondientes a cada fuente, con los mismos parámetros que usamos en el punto 2 para el cálculo de la STFT y guardarlos en las variables x1 y x2 respectivamente. Estas variables, x1 y x2 , podrían pasarse a pistas de audio con la función audiowrite. En la Figura 4.2 podemos ver un ejemplo para dos señales recuperadas correspondientes a una viola (Instrumento A) y a un saxo tenor (Instrumento B). Figura 4.2 Viola y Saxo recuperados tras el proceso de separación de fuentes. 4.3 Evaluación de los resultados Una vez hemos obtenido los resultados con el algoritmo expuesto en la sección anterior, es necesario medir la calidad de los resultados. Existen diferentes medidas para evaluar la calidad de los resultados obtenidos, como la distorsión o la cantidad de señal original que se ha conseguido separar. En nuestro caso nos vamos a centrar en los métodos objetivos [ 57 ], en concreto en una serie de medidas denominadas como medidas orientadas a la calidad de audio (Audio Quality Oriented, AQO), por ser estas las más usadas de entre las medidas objetivas. En estos métodos, se supone que cada fuente estimada produce un modelo, en el que el error total cometido, se divide en tres términos relacionados con tres tipos de error, el modelo se expresa de la siguiente forma [56]: ˆs(t) = sob j(t)+ einter f (t)+earte f (t)(4.1) donde sob j(t) es una deformación permitida de la fuente objetivo si(t) , einter f representa la interferencia que ejercen las fuentes no deseadas y earte f es el error generado en la propia separación. En otros casos, también habría que contar con otro error denominado eruido , que considera el ruido acústico, nosotros no lo tendremos en cuenta, debido a que las pistas de audio se han generado de forma sintética. A partir de este modelo de distorsión, se definen las siguientes medidas para evaluar la separación de las fuentes: •SDR (Signal to Distortion Ratio): compara las fuentes estimadas con las originales (error total), por lo que es la medida más usada para determinar la calidad de la separación de forma global. SDR :=10log10 ksob j k2 keinter f +earte f k2(4.2)
38 Capítulo 4. Simulaciones •SIR (Signal to Interference Ratio): mide la distorsión relativa causada por la interferencia de otras fuentes sobre la fuente objetivo. SIR :=10log10 ksob j k2 keinter f k2(4.3) •SAR (Signal to Artifacts Ratio): mide la distorsión relativa generada por el algoritmo al realizar la separación. SAR :=10log10 ksob j +einter f k2 kearte f k2(4.4) 4.4 Simulación 1 Esta primera simulación, se ha hecho usando las 10 mezclas de audio comentadas al principio del capítulo y ejecutando el algoritmo 5 veces para cada mezcla, como se hace en [ 53 ]. De estas ejecuciones se han obtenido 10 valores de SAR, SDR Y SIR para cada mezcla, formándose así una matriz de tamaño 100×3 , donde cada columna corresponde a una medida y cada fila a un valor. En esta primera simulación, todas las ejecuciones se han hecho fijando los parámetros αyβal valor 1, al igual que se hace en [53]. Para obtener estos valores hemos usado el toolbox de Matlab ® llamado BSS Eval, que se ha convertido en el estándar para medir la eficiencia de los algoritmos de BSS. Este toolbox fue presentado en [ 56 ] y nos hemos servido de [ 22 ] 2 para poder ejecutar de forma correcta sus funciones. Entre las diversas funciones que tiene este toolbox, en este trabajo se ha usado la función bss_eval_sources, que según [ 22 ], sirve para evaluar las señales estimadas de fuentes monocanal. Esta función recibe dos matrices, una con las fuentes estimadas y otra con las fuentes originales y devuelve 4 vectores de tamaño № de fuentes ×1 correspondientes a la SAR, SDR, SIR, y por último un vector que indica a qué fuente j original corresponde la fuente j estimada. En la Figura 4.3 podemos apreciar un diagrama de cajas orientativo de los resultados obtenidos tras la ejecución de nuestro algoritmo. Este diagrama de cajas se obtiene ejecutando la función boxplot en Matlab ® , función que recibe la matriz que contiene todos los valores calculados de SAR, SDR y SIR. Figura 4.3 Diagrama de cajas que contiene 100 valores de las SAR, SDR y SIR de las 10 mezclas de audio para α=1yβ=1. 2http:// bass-db.gforge.inria.fr/ bss_eval/
4.4 Simulación 1 39 Tabla 4.1 Valores correspondientes al diagrama de cajas de la Figura 4.3. SAR SDR SIR Máximo 14,95 13,16 18,05 Media 10,84 8,42 12,05 Mínimo 6,2 3,74 8,18 Para poder evaluar con una referencia los resultados obtenidos, hemos representado los expuestos en [ 53 ] en la Figura 4.4. Gracias a que el señor Fabian-Robert Stöter nos ofreció las pistas de audio resultante de sus simulaciones en [ 53 ], se ha podido comparar estas pistas con las de los instrumentos usando de nuevo BSS Eval. Figura 4.4 Diagrama de cajas que contiene 100 valores de las SAR, SDR y SIR de las 10 mezclas de audio para α=1yβ=1. Tabla 4.2 Valores correspondientes al diagrama de cajas de la Figura 4.4. SAR SDR SIR Máximo 15,54 14 19,31 Media 11,16 8,16 12,45 Mínimo 6,78 4,5 9,24 Los resultados son satisfactorios ya que se asemejan considerablemente, e incluso en nuestra simulación, la media de la SDR es un poco mejor. Podemos afirmar, que el modelo presentado en [ 53 ] se ha implementado de forma correcta en Matlab ® . En ambos resultados aparece una notable dispersión de los valores tanto en SAR, SDR como en SIR, hecho que se razonará en las próximas simulaciones ya que no se comenta en [ 53 ].
46 Capítulo 4. Simulaciones 4.8 Simulación 5 En esta ocasión, lo que se ha hecho es cambiar la estructura de las pistas de audio que se le han pasado al algoritmo. Estas pistas, tienen una duración de 3 segundos y contienen solo la mezcla de dos instrumentos tocando la misma nota y ejecutando un vibrato, ver Figura 4.12. Figura 4.12 Estructura de las pistas de audio de 3 segundos. Se han vuelto a usar las 10 mezclas anteriores, ahora recortadas a 3 segundos y los resultados de la separación se aprecian en la Figura 4.13. Para una óptima comparación, se han usado los mismos valores de los parámetros αyβque en la primera simulación. Figura 4.13 Diagrama de cajas que contiene 100 valores de las SAR, SDR y SIR de las 10 mezclas de audio de 3 segundos para α=1yβ=1. Tabla 4.10 Valores correspondientes al diagrama de cajas de la Figura 4.13. SAR SDR SIR Máximo 5,01 -0,79 7,41 Media 1,51 -2,23 1,53 Mínimo -2,66 -8,88 -3,49
4.9 Simulación 6 47 Se aprecia en la Figura 4.13 y se corrobora con los valores de la Tabla 4.10, que la separación es muy mala. El 100% de los 100 valores que hemos obtenido para la SDR son negativos. Con esta simulación, se ha demostrado que el algoritmo realiza un proceso de aprendizaje cuando recibe los sonidos originales, como en la primera simulación, donde detecta claramente las diversas componentes y eso le sirve para separar. En este caso, el algoritmo es incapaz de detectar dichas componentes y por tanto, no ofrece una solución válida para la separación. 4.9 Simulación 6 Se ha probado el algoritmo con instrumentos diferentes a los usados en [ 53 ], también incluidos en el dataset aportado por Fabian-Robert Stöter. En este caso, se han elegido 4 instrumentos: contrabajo, clarinete, guitarra eléctrica y órgano. Las mezclas tienen la misma estructura que las usadas en la primera simulación y todos los instrumentos se encuentran tocando una nota C4 mientras ejecutan un vibrato. Figura 4.14 Diagrama de cajas que contiene 60 valores de las SAR, SDR y SIR de las 6 mezclas de audio para α=1yβ=1. Tabla 4.11 Valores correspondientes al diagrama de cajas de la Figura 4.14. SAR SDR SIR Máximo 19,52 16,48 20,16 Media 11,17 9,12 12,79 Mínimo 6,1 3,31 7,51 Comparando con los valores de la primera simulación, se aprecia que hemos obtenido unos resultados ligeramente mejores. Centrándonos en los valores medios de la SDR, la media ha mejorado en 0,7 puntos. Se aprecia en este caso una mayor dispersión de los resultados frente a la Simulación 1 (Sección 4.4).
48 Capítulo 4. Simulaciones También se ha hecho una simulación para los valores óptimos de alfa y beta, calculados en la segunda simulación para la mezcla de violonchelo y saxo tenor, como se hiciera en la Simulación 4 (Sección 4.7). En este caso, si que se aprecia una mayor mejora de la SDR, con un aumento de 1,44 puntos frente a los valores expuestos en la Tabla 4.7 de la tercera simulación. Figura 4.15 Diagrama de cajas que contiene 60 valores de las SAR, SDR y SIR de las 6 mezclas de audio para α=1.8yβ=0.8. Tabla 4.12 Valores correspondientes al diagrama de cajas de la Figura 4.15. SAR SDR SIR Máximo 19,68 16,86 20,2 Media 10,44 9,28 14,73 Mínimo 6,25 4,84 11,31 Estos resultados, nos reafirman en el hecho de que el algoritmo es dependiente no solo de los instrumentos, si no que también lo es de los parámetros αyβ.
4.10 Simulación 7 49 4.10 Simulación 7 La última simulación realizada, se ha basado en ejecutar un algoritmo de separación NMF con las 10 mezclas utilizadas en la Simulación 1 (Sección 4.4). Así, podremos comparar el funcionamiento de CFM frente a NMF. El algoritmo usado se ha tomado de [ 32 ] 3 . Se ha ejecutado siguiendo la recomendación de los autores y una vez hecha la separación, se ha utilizado el toolbox BSS Eval para calcular la SAR, SDR y SIR. Se han hecho dos simulaciones para dos divergencias diferentes, la primera para Itakura-Saito y la segunda para Kullback-Leibler. Figura 4.16 Diagrama de cajas que contiene 20 valores de las SAR, SDR y SIR de las 10 mezclas de audio para la divergencia de Itakura-Saito. Tabla 4.13 Valores correspondientes al diagrama de cajas de la Figura 4.16. SAR SDR SIR Máximo 21,85 8,03 14 Media 7,62 3,39 6,84 Mínimo 4,66 -2,5 -1,53 Frente a la primera simulación, la media de la SDR ha bajado 4,97 puntos para la divergencia Itakura-Saito y 5,06 puntos para Kullback-Leibler, mientras que para la simulación donde se han optimizado los parámetros α y β el descenso de la SDR ha sido de 5,23 y 5,32 respectivamente. Con estos resultados, queda demostrado que NMF no es un modelo eficaz para fuentes unísonas mono-canal moduladas en tiempo y frecuencia, tal como se ha explicado teóricamente en la introducción del Capítulo 3. 3https:// github.com/ EliasKokkinis/ audio-source-separation
50 Capítulo 4. Simulaciones Figura 4.17 Diagrama de cajas que contiene 20 valores de las SAR, SDR y SIR de las 10 mezclas de audio para la divergencia de Kullback-Leibler. Tabla 4.14 Valores correspondientes al diagrama de cajas de la Figura 4.17. SAR SDR SIR Máximo 25,95 7,48 12,32 Media 7,76 3,3 5,67 Mínimo 4,39 -2,15 -1,84 Se aprecia en estas simulaciones una fuerte dispersión, que nos indica que NMF también es dependiente de los instrumentos. En cuanto a la dependencia a los parámetros α y β , los resultados varían considerablemente menos que en las simulaciones realizadas con CFM.
5 Conclusiones y Líneas Futuras 5.1 Trabajo realizado y conclusiones En este proyecto hemos trabajado sobre un método para explotar texturas de modulación recurrentes para el problema de la Separación Ciega de Fuentes, basándonos durante todo el trabajo en el presentado por F.-R. Stöter et al. en [53]. Para dar sentido a todo nuestro trabajo, se ha partido de un estudio teórico general de BSS y se han ido desarrollando de manera más extensa algunas técnicas y soluciones existentes, como NMF y sus propiedades, o las divergencias más utilizadas en la estimación de los parámetros. Tras la introducción teórica, se ha hecho un estudio teórico-práctico para obtener las fórmulas del algoritmo usado para resolver el problema se separación de fuentes. Para ello se ha introducido el concepto de divergencia AB, sus propiedades y sus ventajas respecto a otras divergencias conocidas. Finalmente en este estudio, se ha usado un algoritmo multiplicativo, basado en el estudio teórico anterior. Finalmente, se ha implementado el algoritmo en Matlab ® y se han realizado diversas simulaciones, presentando los resultados obtenidos. En función de los resultados obtenidos y su análisis, se han ido extrayendo varias conclusiones que presentamos a continuación. En primer lugar y como conclusión más importante, los valores de nuestros resultados para mezclas de dos instrumentos tocando la misma nota mientras ejecutan un vibrato, indican que este método funciona bien en este desafiante escenario. Lo que implica directamente, que la implementación del modelo en Matlab ® se ha realizado de forma correcta. En segundo lugar, del estudio de las divergencias AB, concluimos que su uso ha aportado mucha versatilidad al algoritmo, ya que al tratarse de una familia de divergencias que engloban a otras conocidas, pueden modificarse de manera sencilla las fórmulas para la estimación de las fuentes. También nos brinda la posibilidad de realizar un screening para determinar las zonas o puntos del plano alfa/beta donde los resultados son más favorables, como se ha hecho en la Sección 4.5. Se ha demostrado que para diferentes mezclas los valores óptimos de alfa y beta son diferentes, por lo que el estudio de las divergencias AB es de gran ayuda. Por último, se puede afirmar que con el estudio de las divergencias AB se han encontrado combinaciones de alfa/beta que mejoran los resultados de [ 53 ], donde solo se presentan resultados para α=1 y β=1 , por lo que podemos considerar la realización del trabajo y sus resultados como altamente satisfactorios. 51
52 Capítulo 5. Conclusiones y Líneas Futuras 5.2 Líneas futuras Tras la realización de este trabajo, consideramos que una de las mayores restricciones del mismo, es el tipo de pistas de audio con el que hemos trabajado. Sin duda, una de las líneas futuras se basa en probar el modelo con otras pistas de audio, sintéticas con diferente estructura y también con grabaciones reales, para ver cómo actúa el modelo frente al ruido. Por otra parte, sería de gran interés mejorar la eficiencia de nuestro código en Matlab ® , especialmente para el estudio de las AB divergencias ya que ha sido uno de los puntos negativos de este trabajo. Debido a los resultados obtenidos en la Simulación 5 (Sección 4.8), es necesario en un futuro estudiar a fondo el aprendizaje realizado por el algoritmo y cómo optimizar la separación cuando no se tienen las pistas originales previas a la mezcla. También podría estudiarse la aplicación del modelo al problema de la separación del habla. Ya se han usado algoritmos basados en NMF para este problema [ 11 ]. En la separación del habla, también se suele encontrar audio modulado en frecuencia y en amplitud, como características propias del habla, ya que el hablante no suele usar el mismo tono ni el mismo volumen durante una conversación, algo que puede ayudar a un correcto funcionamiento del modelo. Dichas modulaciones serían seguramente más variables que las estudiadas en este trabajo (que son siempre de 5 Hz), hecho que puede resultar una dificultad añadida para el modelo, además, las pistas de audio que recibiría no tendrían la misma estructura que en este estudio, algo que la Simulación 5 (Sección 4.8) ha demostrado que influye de forma muy negativa en la separación. Es común en los algoritmos de separación del habla, que haya una primera fase de entrenamiento del algoritmo con audios del conjunto de datos del problema a resolver, esto podría añadirse o ampliarse al actual modelo.
Índice de Figuras 2.1 Modelo BSS lineal instantáneo [60] 4 2.2 Modelo de mezcla convolutiva [14] 6 2.3 Respuesta impulsiva de una sala 10 2.4 Espectrograma de una melodía tocada en un xilófono [55] 11 2.5 Descomposición NMF multinivel del espectrograma de la Figura 2.4 13 2.6 Tensor de N=3 [13] 14 2.7 Modelo NMF bilineal 16 2.8 Esquema NMF con offset [13] 17 2.9 Esquema NMF multicapa [13] 17 2.10 Esquema NMF Proyectiva [13] 18 2.11 Esquema NMF Convexa [13] 18 2.12 Esquema NMF Convolutiva [13] 19 2.13 Esquema NMF Superpuesta [13] 19 2.14 Aproximación Large-Scale NMF 24 2.15 Esquema básico de la separación de fuentes de audio mediante NMF [19] 25 3.1 Transformada de Destino Recurrente, CFT [53] 28 3.2 Transformada de Destino Recurrente,CFT 28 3.3 Modelo de Destino Recurrente, CFM [53] 30 3.4 Ilustración gráfica de las propiedades de inversión y dualidad en la divergencia-AB 32 3.5 Ilustración gráfica de cómo los parámetros establecidos αyβpueden controlar la influencia de los ratios individuales pit /qit 33 4.1 Estructura de las pistas de audio de entrada del algoritmo 35 4.2 Viola y Saxo recuperados tras el proceso de separación de fuentes 37 4.3 Diagrama de cajas que contiene 100 valores de las SAR, SDR y SIR de las 10 mezclas de audio para α=1yβ=138 4.4 Diagrama de cajas que contiene 100 valores de las SAR, SDR y SIR de las 10 mezclas de audio para α=1yβ=139 4.5 Representación gráfica de los valores de la SAR, SDR y SIR en función de los parámetros αyβ para la mezcla de viola y saxo tenor 40 4.6 Representación gráfica de los valores de la SDR en función de los parámetros αyβpara las 10 mezclas 41 4.7 Representación gráfica de los valores de la SDR en función de los parámetros αyβpara la mezcla de viola y saxo tenor 41 4.8 Diagrama de cajas que contiene 100 valores de las SAR, SDR y SIR de las 10 mezclas de audio para αyβóptimos de cada mezcla 42 4.9 Diagrama de cajas que contiene 10 valores de SDR de cada mezcla 43 4.10 Diagrama de cajas que contiene 100 valores, 10 de cada mezcla para α=1.8yβ=0.844 4.11 Diagrama de cajas donde cada caja representa 10 valores SDR de cada una de las 10 mezclas α=1.8yβ=0.845 4.12 Estructura de las pistas de audio de 3 segundos 46 53
54 Índice de Figuras 4.13 Diagrama de cajas que contiene 100 valores de las SAR, SDR y SIR de las 10 mezclas de audio de 3 segundos para α=1yβ=146 4.14 Diagrama de cajas que contiene 60 valores de las SAR, SDR y SIR de las 6 mezclas de audio para α=1yβ=147 4.15 Diagrama de cajas que contiene 60 valores de las SAR, SDR y SIR de las 6 mezclas de audio para α=1.8yβ=0.848 4.16 Diagrama de cajas que contiene 20 valores de las SAR, SDR y SIR de las 10 mezclas de audio para la divergencia de Itakura-Saito 49 4.17 Diagrama de cajas que contiene 20 valores de las SAR, SDR y SIR de las 10 mezclas de audio para la divergencia de Kullback-Leibler 50
Índice de Tablas 4.1 Valores correspondientes al diagrama de cajas de la Figura 4.3 39 4.2 Valores correspondientes al diagrama de cajas de la Figura 4.4 39 4.3 Valores óptimos de los parámetros αyβpara cada mezcla y valor SDR máximo 41 4.4 Valores correspondientes al diagrama de cajas de la Figura 4.8 42 4.5 Valores correspondientes al diagrama de cajas de la Figura 4.9 43 4.6 Desviación típica de los valores de la Tabla 4.5 43 4.7 Valores correspondientes al diagrama de cajas de la Figura 4.10 44 4.8 Valores correspondientes al diagrama de cajas de la Figura 4.11 45 4.9 Desviación típica de los valores de la Tabla 4.8 45 4.10 Valores correspondientes al diagrama de cajas de la Figura 4.13 46 4.11 Valores correspondientes al diagrama de cajas de la Figura 4.14 47 4.12 Valores correspondientes al diagrama de cajas de la Figura 4.15 48 4.13 Valores correspondientes al diagrama de cajas de la Figura 4.16 49 4.14 Valores correspondientes al diagrama de cajas de la Figura 4.17 50 55
62 Bibliografía [52] P. Smaragdis, Convolutive speech bases and their application to supervised speech separation, IEEE Transactions on Audio, Speech, and Language Processing 15 (2007), no. 1, 1–12. [53] Fabian-Robert Stöter, Antoine Liutkus, Roland Badeau, Bernd Edler, and Paul Magron, Common Fate Model for Unison Source Separation, 41st International Conference on Acoustics, Speech and Signal Processing (ICASSP) (Shanghai, China), Proceedings of the 41st International Conference on Acoustics, Speech and Signal Processing (ICASSP), IEEE, 2016. [54] F. R. Stöter, S. Bayer, B. Edler, and P. Magron, Unison source separation, 17th International Conference on Digital Audio Effects, September 2014, pp. 235–241. [55] E. Vincent, N. Bertin, R. Gribonval, and F. Bimbot, From blind to guided audio source separation: How models and side information can improve the separation of sound, IEEE Signal Processing Magazine 31 (2014), no. 3, 107–115. [56] E. Vincent, R. Gribonval, and C. Fevotte, Performance measurement in blind audio source separation, IEEE Transactions on Audio, Speech, and Language Processing 14 (2006), no. 4, 1462–1469. [57] Emmanuel Vincent, Improved perceptual metrics for the evaluation of audio source separation, Latent Variable Analysis and Signal Separation (Berlin, Heidelberg) (Fabian Theis, Andrzej Cichocki, Arie Yeredor, and Michael Zibulevsky, eds.), Springer Berlin Heidelberg, 2012, pp. 430–437. [58] Emmanuel Vincent, Shoko Araki, Fabian Theis, Guido Nolte, Pau Bofill, Hiroshi Sawada, Alexey Ozerov, Vikrham Gowreesunker, Dominik Lutter, and Ngoc Q. K. Duong, The signal separation evaluation campaign (2007-2010): Achievements and remaining challenges, Signal Process. 92 (2012), no. 8, 1928–1936. [59] T. Virtanen, Monaural sound source separation by nonnegative matrix factorization with temporal continuity and sparseness criteria, IEEE Transactions on Audio, Speech, and Language Processing 15 (2007), no. 3, 1066–1074. [60] Xinling Wen, Research and simulation of linear instantaneous blind signal separation algorithm, Advances in Computer Science, Environment, Ecoinformatics, and Education (Berlin, Heidelberg) (Song Lin and Xiong Huang, eds.), Springer Berlin Heidelberg, 2011, pp. 119–124. [61] O. Yilmaz and S. Rickard, Blind separation of speech mixtures via time-frequency masking, IEEE Transactions on Signal Processing 52 (2004), no. 7, 1830–1847. [62] G. Zhou, Q. Zhao, Y. Zhang, T. Adalı, S. Xie, and A. Cichocki, Linked component analysis from matrices to high-order tensors: Applications to biomedical data, Proceedings of the IEEE 104 (2016), no. 2, 310–331.