Simulador de recintos tridimensionales para el análisis del campo acústico en algoritmos de procesado de la señal

UNIVERSIDAD POLITECNICA DE VALENCIA ESCUELA POLITECNICA SUPERIOR DE GANDIA Master en Ingeniería Acústica “Simulador de recintos tridimensionales para el análisis del campo acústico en algoritmos de procesado de la señal” TESIS DE MASTER Autor: Jorge Machín Mindán Director: D. Miguel Ferrer Contreras GANDIA, 2010 Simulador de recintos tridimensionales para el análisis del campo acústico en algoritmos de procesado de la señal. 2 Autor: Jorge Machín Mindán, email: [email protected] Fecha de entrega: 16-10-2010 ÍNDICE RESUMEN .......................................................................................................................................... 3 I. OBJETIVOS ..................................................................................................................................... 4 II. INTRODUCCIÓN .......................................................................................................................... 4 III. DECISIONES ................................................................................................................................ 5 III.1. SIMULACION FDTD .......................................................................................................................................5 III.2. CONDICIONES DE CONTORNO ...............................................................................................................6 III.3. PASO TEMPORAL ...........................................................................................................................................6 III.4. CREACIÓN DE LAS VARIABLES Y MATRICES ..................................................................................7 III.5. SEÑALES DE EXCITACIÓN .........................................................................................................................7 III.6. FUENTES DE RUIDO .....................................................................................................................................8 III.7. ALGORITMO PARA EL CAR .......................................................................................................................8 III.8. RESPUESTA AL IMPULSO ..........................................................................................................................9 IV. CREACIÓN DEL INTERFAZ.....................................................................................................10 IV.1. INTERFAZ 2 D ............................................................................................................................................... 10 IV.2. INTERFAZ 3 D ............................................................................................................................................... 11 IV.3. RESTRICCIONES .......................................................................................................................................... 11 V. SIMUALCIONES EN DOS DIMENSIONES ...............................................................................13 V.1. ESTUDIO DEL CORRECTO FUNCIONAMIENTO DEL SISTEMA .............................................. 13 V.2. ESTUDIO DE LA INFLUENCIA DEL NÚMERO DE FUENTES Y SENSORES ........................ 15 V.3. EXPRESIÓN DE LOS RESULTADOS ...................................................................................................... 19 V.4. VARIACIÓN DEL ALFA (ABSORCIÓN DEL PML)............................................................................ 20 V.5. IDEA: ENCERRAR LA FUENTE DE RUIDO ........................................................................................ 22 V.6. SIMULACIÓN MULTIFUENTE ................................................................................................................. 24 VI. SIMULACIONES EN TRES DIMENSIONES ...........................................................................26 VI.1. COMPARACIÓN CON 2 D ......................................................................................................................... 26 VI.2. CAR EN TRES DIMENSIONES ................................................................................................................ 27 VII. CONCLUSIONES ......................................................................................................................35 VIII. LÍNEAS FUTURAS..................................................................................................................36 IX. BIBLIOGRAFÍA .........................................................................................................................37 ANEXO I. DEFINICIÓN DE VARIABLES .....................................................................................38 ANEXO II. EJEMPLO ......................................................................................................................40 ANEXO III. GUÍA PARA USUARIO. ..............................................................................................44 Simulador de recintos tridimensionales para el análisis del campo acústico en algoritmos de procesado de la señal. 3 Autor: Jorge Machín Mindán, email: [email protected] Fecha de entrega: 16-10-2010 Simulador de recintos tridimensionales para el análisis del campo acústico en algoritmos de procesado de la señal. Autor: Jorge Machín Mindán Director: Miguel Ferrer Contreras RESUMEN Por todos es conocido que el uso de Finite-Difference Time-Domain (FDTD) es cada vez más usado para las simulaciones acústicas, tanto en el caso de dos como en el tres dimensiones. Por lo que, en este proyecto se pretende diseñar una herramienta para pode utilizar esta técnica, y comprobar el correcto funcionamiento de diferentes algoritmos de procesado de la señal, como pueden ser el control activo de ruido (CAR) y/o Wave Field Synthesis (WFS). La técnica FDTD, al tratarse de una simulación en tiempo, nos da información del campo acústico resultante en cada instante, por lo que nos proporciona una idea del comportamiento de los algoritmos y de la sala tanto en el transitorio como en el estacionario (del filtro adaptativo). Las comprobaciones que se va a llevar a cabo en este trabajo van a ser: - Robustez y facilidad de manejo de las herramientas diseñadas. - El Estudio del comportamiento de los algoritmos en el estacionario, el cual nos da información del campo acústico resultante, y así poder ser comparado con la simulación sin el algoritmo. ABSTRACT Everyone knows that the use of Finite-Difference Time-Domain (FDTD) is increasingly used for acoustic simulations both in the case of two as in three dimensions. Herefore, the objectives of this project is the realization of an interface for the use of this technique and verify the correct operation of different signal processing algorithms, such as active noise control (ANC) and Wave Field Synyhesis (WFS). FDTD technique give us information of thr acoustics field at any moment, by giving us information about the bahavior of the algotihms and the room in the transient and in the steady-state (filter). Checks to be carried out in this work will be: - Robustness and ease of management tools designed. - The Study of the behavior of algorithms in the steady-state, which gives us the sound field information, so they can be compared to the simulation without the algorithm. Simulador de recintos tridimensionales para el análisis del campo acústico en algoritmos de procesado de la señal. 4 Autor: Jorge Machín Mindán, email: [email protected] Fecha de entrega: 16-10-2010 I. OBJETIVOS El objetivo de esta tesina es doble: por un lado se pretende diseñar una herramienta que permita simular la propagación del campo acústico en un recinto tridimensional, una vez definida la geometría del recinto y las características de las discontinuidades (paredes), y por el otro lado incluir el procesado discreto de diferentes algoritmos de Control activo de ruido y/o Wave Field Synthesis que gobiernen la generación de señales acústicas en cualquier punto del recinto simulado. De esta forma, se podrá analizar el campo acústico existente en cualquier plano del recinto (o incluso en cualquier región tridimensional), no sólo ante determinadas señales que se propagan de forma libre en el recinto sino también cuando estos algoritmos están en funcionamiento. II. INTRODUCCIÓN Esta tesina tenía como objetivo la creación de un interfaz para la simulación de la propagación acústica, y fue por ello por lo que la técnica elegida para esta tipo de simulación fue Finite Diference Time Domine (FDTD), una técnica cada vez más utilizada en acústica. Elegimos esta técnica porque nos permite visionar como actúan los diferentes algoritmos aplicados con el paso del tiempo. Es interesante ver como varía la presión sonora en cada instante de tiempo para comprender mejor los algoritmos de control activo de ruido. Sin embargo esta técnica también tiene un gran inconveniente, y es el gran coste computacional que conlleva este tipo de simulación, sobretodo en el caso de tres dimensiones. Es por esto, por lo que en este proyecto se ha decidido crear un interfaz para la simulación en dos dimensiones como una primer aproximación al problema. Por otro lado tenemos el control activo de ruido (CAR), el cual consiste en reducir el ruido, presente en un punto o zona, a partir de una señal (en contrafase al ruido) generada electrónicamente. Se basa en el principio de superposición de campos sonoros, más concretamente en el de interferencia destructiva. Simulador de recintos tridimensionales para el análisis del campo acústico en algoritmos de procesado de la señal. 5 Autor: Jorge Machín Mindán, email: [email protected] Fecha de entrega: 16-10-2010 III. DECISIONES A lo largo de todo el desarrollo del proyecto se ha visto la necesidad de tomar varias decisiones, a continuación se explican alguna de las más importantes. III.1. SIMULACION FDTD Al tratarse de una simulación acústica debemos utilizar la ecuación de ondas para obtener el nivel de presión en cada punto. A continuación vamos a explicar brevemente las ecuaciones y simplificaciones elegidas para llevar a cabo nuestras simulaciones. Ecuación de ondas en dos dimensiones[1]: (1) Donde: p, presión (escalar) x, coordenada espacial x y, coordenada espacial y ux, velocidad de la partícula en el eje x (vector) uy, velocidad de la partícula en el eje y (vector) t, es el tiempo c, velocidad de propagación del sonido en el aire ρ, es la densidad del aire Para la resolución de las ecuaciones anteriores aproximamos las primeras derivadas tanto en tiempo como en espacio, por el cociente de la diferencia entre muestras espaciales y el incremento espacial (derivada primera en espacio). La derivada en tiempo se aproxima de la misma manera, pero el incremento del denominador es el temporal y la diferencia es entre 2 instantes de tiempo consecutivos en el mismo punto. 0 0 0 2                             y u x u c t p t u y p t u x p y x y x    Simulador de recintos tridimensionales para el análisis del campo acústico en algoritmos de procesado de la señal. 6 Autor: Jorge Machín Mindán, email: [email protected] Fecha de entrega: 16-10-2010 Expresión de la aproximación utilizada para la solución de las derivadas parciales: Derivada en espacio (2) Derivada en tiempo (3) El método elegido, FDTD, calcula los valores en instantes alternados; calculamos ux y uy en el mismo instante mientras que “p” la calculamos en otro instante, con la finalidad de conseguir un error de segundo orden debido a las aproximaciones utilizadas en las derivadas temporales. Este cálculo se realiza dentro de un bucle, en cada interacción aumentamos el paso correspondiente a un incremento temporal. Expresiones a resolver en cada iteración [2]: Presión: (4) Velocidad: (5) Este vector velocidad se resuelve dos veces, una para cada eje (x e y). III.2. CONDICIONES DE CONTORNO Se decidió introducir al final de nuestro dominio de simulación unas condiciones de contorno absorbentes, para ver las diferencias entre distintas absorciones en las paredes. Es por lo que se eligió que la caracterización de las discontinuidades (paredes) fueran definidas mediante Perfectly Matched Layer (PML) [3] [4]. A su vez también se decidió tomar la condición de contorno, en la paredes, en la que la velocidad de las partículas en la superficie sea nula (ui = 0), donde ī es perpendicular a la dirección del contorno [5]. III.3. PASO TEMPORAL En nuestras simulaciones también se decidió utilizar un número de Courant (CFL) igual a 1 (si se desea utilizar un incremento temporal distinto, sólo hay que rellenar la casilla de Incremento Temporal con el valor deseado). h tiptip t tiuttiu tt      ),(),1( ),(),( 22 2 1 2 1  x tiutiu cttiptip tt     ),(),( ),(),( 2 1 2 1 2 22  t ttfttf t f x xxfxxf x f tt xx             )()( )()( 2 1 0 2 1 0 2 1 0 2 1 0 0 0 Simulador de recintos tridimensionales para el análisis del campo acústico en algoritmos de procesado de la señal. 7 Autor: Jorge Machín Mindán, email: [email protected] Fecha de entrega: 16-10-2010 El número de Courant (CFL) debe ser menor o igual a la unidad para que la simulación sea estable, ya que se trata de una condición de convergencia de ecuaciones diferenciales en derivadas parciales. Como consecuencia de esta condición, el paso temporal debe ser inferior a un cierto valor, de lo contrario la simulación producirá resultados incorrectos. III.4. CREACIÓN DE LAS VARIABLES Y MATRICES Una vez tomadas las decisiones correspondientes al PML y a las condiciones de contorno pasamos a definir las variables y las matrices. El primer paso es definir el contorno con los puntos cargados en la GUI, dicha interface será diseñada en Matlab (más adelante se explicará como esta diseñada). Con estos puntos definiremos una matriz que guarde los contornos. Será esta matriz la que aplicaremos como máscara a la hora de la simulación, la cual nos marcará las condiciones de contorno (velocidad de las partículas nula). A la vez que definimos esta máscara también haremos lo propio con las discontinuidades (las paredes), las cuales serán definidas mediante PMLs. Otra decisión importante es la de utilizar, por defecto, un incremento temporal tal que el número de Courant sea 1. No obstante, como se ha comentado antes, si se desea utilizar un incremento temporal menor, este se puede utilizar rellenando el hueco correspondiente a “Incremento Temporal” en el interfaz. La explicación de cómo se realizan los pasos anteriores esta explicado en el anexo I. III.5. SEÑALES DE EXCITACIÓN Con el fin de poder realizar diferentes simulaciones se ha decidido elegir varios tipos de señales de excitación, estas señales son aplicables a cada una de las fuentes de ruido: - Señal senoidal. Tono puro (frecuencia fundamental definida por el usuario). - Señal ricker. Excitación Ricker, frecuencia de corte definida por el usuario. - Ruido filtrado. El ruido generado tiene una función de probabilidad Gaussiana con una densidad de potencia constante hasta la frecuencias de corte (Wn) definida por el usuario, la media de esta señal es cero y su desviación típica uno. Todas las fuentes de ruido blanco que emitan a la vez tienen la misma señal. Esto nos permite aproximar las superficies que generen sonido mediante fuentes puntuales emitiendo en fase. - Archivo wav. Para este paso hay que realizar un remuestreo de la señal, para que la frecuencia de muestreo de dicha señal coincida con la inversa del incremento temporal de la simulación. Simulador de recintos tridimensionales para el análisis del campo acústico en algoritmos de procesado de la señal. 8 Autor: Jorge Machín Mindán, email: [email protected] Fecha de entrega: 16-10-2010 III.6. FUENTES DE RUIDO En primer lugar debemos comentar que en nuestras simulaciones usamos fuentes puntuales. Por todos es conocido que estas fuentes se caracterizan por comportarse como fuentes omnidireccionales. Sabemos que esto no es real, pero lo podemos considerar como una aproximación. Este “problema” se podría solucionar definiendo las fuentes como pistones planos (formados por un conjunto de fuentes puntuales). III.7. ALGORITMO PARA EL CAR En esta tesina se decidió implementar el sistema de Control Activo de Ruido basado en el algoritmo LMS con la 󳮁estructura de filtrado-x, ya que se trata de uno de los más usados para este fin [6] [7] [8]. Este algoritmo nos permite utilizar filtros adaptativos, es decir, que se recalculen sus coeficientes en cada interacción del algoritmo. Dichos coeficientes cambian conforme varía la señal de error (los sensores en nuestro caso) intentando que la potencia de la señal de error sea mínima, en nuestro caso la señal de error es la señal en los sensores. Una vez decidido ya el algoritmo a utilizar se consideró la posibilidad de aplicarlo con varias fuentes secundarías y varios sensores a la vez, pero sólo con una señal de referencia (una única fuente de ruido). Después de diseñar este algoritmo se vieron otras posibles aplicaciones, como se verá más adelante, por lo que se decidió realizar también el algoritmo con varias señales de referencia. Para el correcto funcionamiento de la simulación deben calcularse tantos filtros adaptativos como pares de señal de referencia – fuentes secundarias haya en la simulación. Hay varios métodos de calcular los nuevos coeficientes del filtro adaptativo, pero en nuestro caso se ha elegido el cálculo de los nuevos coeficientes como la media de todas las señales de error (ek) por su correspondiente señal de referencia filtrada (xfijk) [6]. Los nuevos coeficientes del filtro se calculan conforme a la siguiente expresión: wij=wij-2*mu*(xfij1*e1+xfij2*e2+...xfijK*ek) (6) A su vez se generarán tantas señales anti-ruido como fuentes secundarias haya (j): yj=x1*w1j+x2*w2j+...+xI*wIJ (7) Tendremos tantas señales de referencia como fuentes de ruido haya (i): xi Simulador de recintos tridimensionales para el análisis del campo acústico en algoritmos de procesado de la señal. 9 Autor: Jorge Machín Mindán, email: [email protected] Fecha de entrega: 16-10-2010 A su vez estas señales de referencia serán filtradas por las plantas acústicas B_i(j,:): xfijk= xi*Bij (8) Y se monitorizará la señal de error (señal a minimizar) en tantos puntos como sensores haya (k): ek Otra decisión, casi obligatoria, es que el usuario pueda elegir el paso de convergencia (mu), que como hemos visto anteriormente es un parámetro que regula la velocidad de convergencia del filtro adaptativo. La elección de este parámetro tiene especial interés ya que es el que puede hacer que el sistema converja lentamente o que incluso se vuelva inestable. La señal (o señales) de referencia que obtenemos para el control activo de ruido suponemos que se toma de la fuente de ruido de una forma mecánica (acelerómetro o similar), por lo que no se trata de un sensor acústico (micrófono). Esto se hace así para evitar retroalimentaciones acústicas en nuestro sistema. III.8. RESPUESTA AL IMPULSO Una vez ya elegido el algoritmo para el CAR, ya conocemos los cálculos que debemos realizar, y uno de estos cálculos es el de las plantas, o caminos acústicos. Como ya se ha comentado, el algoritmo LMS con filtrado x necesita, para su correcto funcionamiento, de tantas plantas acústicas como pares fuente-sensor haya: Las plantas relacionadas con cada fuente secundaria son guardadas como una variable nueva (B_i), mientras que cada sensor se guarda como una fila diferente B_i(j,:). La señal elegida para obtener estas plantas acústicas es un impulso. Una vez emitido este impulso y guardado las señales recibidas en los sensores durante la simulación, ya tenemos la información necesaria para describir el camino acústico. El número de muestras de la planta es definido por el usuario para buscar una buena “relación” entre fiabilidad y tiempo de cálculo. El tamaño de la planta también determina cuando el CAR empieza a funcionar, es decir, el CAR tiene que almacenar tantas muestras como tamaño tengan los filtros para poder realizar sus cálculos. Todas las decisiones anteriores son extrapolables para el caso de simulaciones en tres dimensiones. Simulador de recintos tridimensionales para el análisis del campo acústico en algoritmos de procesado de la señal. 16 Autor: Jorge Machín Mindán, email: [email protected] Fecha de entrega: 16-10-2010 4º caso. La fuente anterior emitiendo. Y un sistema de CAR con tres fuentes antiruido y con tres sensores. 5º caso. La fuente anterior emitiendo. Y un sistema de CAR con tres fuentes antiruido y con cuatro sensores. 6º caso. La fuente anterior emitiendo. Y un sistema de CAR con cuatro fuentes antiruido y con cuatro sensores. En la siguiente tabla se muestra la posición de las fuentes y los sensores en la simulación: Fig.11. Tabla que representa la posición de los sensores y fuentes secundarias en el sistema de CAR Caso 1 (posición) Caso 2 (posición) Caso 3 (posición) Caso 4 (posición) Caso 5 (posición) Caso 6 (posición) Fuente 1 (0.5 0.5) 1 (0.5 0.5) 1 (0.5 0.5) 1 (0.5 0.5) 1 (0.5 0.5) 1 (0.5 0.5) Fuente Antiruido 0 1 (1.5 1.8) 2 (1.5 1.8) (1 1.4) 3 (1.5 1.8) (1 1.4) (1.7 0.8) 3 (1.5 1.8) (1 1.4) (1.7 0.8) 4 (1.5 1.8) (1 1.4) (1.7 0.8) (1 1) Sensor 0 1 (2 1.5) 2 (2 1.5) (2.2 1.2) 3 (2 1.5) (2.2 1.2) (2.3 1.8) 4 (2 1.5) (2.2 1.2) (2.3 1.8) (2.5 1.5) 4 (2 1.5) (2.2 1.2) (2.3 1.8) (2.5 1.5) Simulador de recintos tridimensionales para el análisis del campo acústico en algoritmos de procesado de la señal. 17 Autor: Jorge Machín Mindán, email: [email protected] Fecha de entrega: 16-10-2010 La siguiente imagen muestra el esquema de la simulación: Fig.12. Esquema que representa el montaje de la simulación. El primer paso para comprobar el comportamiento del algoritmo LMS con filtrado x multicanal es la comparación entre las señales grabadas en el mismo punto en distintas condiciones. Con esta primera comparación queremos mostrar que sucede en el sensor uno (2 1.5) conforme aumentamos el número de estos y de fuentes de ruido. Fig.13. Representación de la señal grabada por un sensor en los diferentes casos. En el eje x tiempo de simulación (s), eje y presión instantánea (arbitraría). Simulador de recintos tridimensionales para el análisis del campo acústico en algoritmos de procesado de la señal. 18 Autor: Jorge Machín Mindán, email: [email protected] Fecha de entrega: 16-10-2010 Como vemos en la gráfica la presión en el punto disminuye conforme aumentamos el número de sensores y de fuentes, esto no siempre es así, hay casos en los que el sistema CAR da mayores atenuaciones (en los sensores) en los casos de una fuente secundaria y un único sensor de error. Pero el campo acústico resultante en los casos con más fuentes y/o sensores se observa una mayor superficie de actuación (atenuación). Un ejemplo de lo anteriormente comentado se puede ver en el anexo II. Una vez que hemos visto las señales en los sensores de error queremos observar el campo acústico resultante de la simulación para conocer la atenuación en cada caso, ya que como hemos visto antes una mayor atenuación en un punto no tiene porque significar una mayor zona de actuación. Para expresar los resultados utilizamos la potencia media de la presión acústica. Usamos esta potencia para hacernos una idea del campo acústico resultante (no queremos ver el campo del estacionario del CAR, es simplemente una forma de comparar los diferentes casos). Dado que todas las simulaciones tienen la misma duración y la fuente es la misma en todos los casos, las diferencias se deben únicamente al CAR. Fig.14. Representación gráfica de las potencias medias obtenidas durante las simulaciones: de izquierda a derecha y de arriba abajo. Solo fuente, CAR 1 f y 1s, CAR 2f y 2s, CAR 3f y 3s, CAR 3 f y 4s, CAR 4f y 4s. Tras ver la gráfica podemos observar que cuantos más pares de sensores – fuentes secundarias tenga nuestro sistema mayor será el área de influencia del CAR. Como hemos comentado antes, debemos encontrar la mejor relación entre atenuación y superficie afectada dependiendo de cada caso. Así como la mejor posición posible de los transductores. Simulador de recintos tridimensionales para el análisis del campo acústico en algoritmos de procesado de la señal. 19 Autor: Jorge Machín Mindán, email: [email protected] Fecha de entrega: 16-10-2010 V.3. EXPRESIÓN DE LOS RESULTADOS Una vez que vemos que el sistema funciona correctamente, intentamos averiguar la mejor forma de expresar los resultados para que se pueda mostrar la diferencia entre una sala con CAR y sin él. Pensamos que la mejor forma de representar esta diferencia es la de mostrar los campos acústicos en las dos situaciones. En el caso sin CAR se realizará la potencia media de la presión en un intervalo concreto de la simulación, mientras que el caso del CAR el resultado será la potencia media de la presión del mismo intervalo de tiempo, pero con la condición de que se haya alcanzado el estacionario. La siguiente gráfica muestra la comparación de potencias medias con y sin CAR, del mismo intervalo de tiempo una vez ya en el estacionario (llega a las 350 ms) durante un intervalo de 150 ms. La representación se realiza en escala logarítmica arbitraría. Fig.15. Comparación de potencias medias entre una sala con y sin CAR. Se muestran las dos salas juntas para tener la misma referencia a la hora de calcular la escala logarítmica. La sala sin CAR esta desde 2 hasta 4 en el eje y, mientras que la sala con CAR esta entre 0 y 2. Para poder visualizar mejor la diferencia se ha calculado la atenuación entre las dos simulaciones como el logaritmo del cociente de la potencia media sin CAR y la potencia media con CAR. Simulador de recintos tridimensionales para el análisis del campo acústico en algoritmos de procesado de la señal. 20 Autor: Jorge Machín Mindán, email: [email protected] Fecha de entrega: 16-10-2010 Fig.16. Diferencia entre las potencias medias con y sin CAR. Escala logarítmica referencia la potencia media de la sala sin el CAR. Como podemos observar la zona de actuación producida por el CAR es bastante extensa y cubre la zona en la que queríamos atenuar el ruido. En la gráfica se observa que tenemos un máximo de atenuación de 48 dB, y la disminución del ruido es superior a los 25 dB en la zona de interés, zona en la que están situados los sensores. V.4. VARIACIÓN DEL ALFA (ABSORCIÓN DEL PML) En este apartado se pretende comprobar los resultados de las simulaciones con diferentes absorciones del PML. Dado que los ejemplos anteriores se han llevado a cabo con un alfa (absorción del PML) de 0.5, se ha querido comprobar como se comportaba el mismo algoritmo en condiciones más adversas, con un alfa menor. Como nuevo coeficiente de absorción del PML se ha elegido 0.1. A continuación podemos ver la comparación de potencias medias (una vez que el CAR está en el estacionario) entre la sala con y sin CAR. Simulador de recintos tridimensionales para el análisis del campo acústico en algoritmos de procesado de la señal. 21 Autor: Jorge Machín Mindán, email: [email protected] Fecha de entrega: 16-10-2010 Fig. 17. Representación de las potencias medias con y sin CAR escala logarítmica arbitraría. La representación sigue la misma estructura que la Fig.11. Cálculo de la atenuación en cada punto como la división, en dB, de la potencia media sin CAR entre la potencia media con CAR. Fig.18. Diferencia entre las potencias medias con y sin CAR. Escala logarítmica referencia la potencia media de la sala sin el CAR. Simulador de recintos tridimensionales para el análisis del campo acústico en algoritmos de procesado de la señal. 22 Autor: Jorge Machín Mindán, email: [email protected] Fecha de entrega: 16-10-2010 En esta gráfica se observa que no sólo hay puntos en los que la señal es atenuada al activar el CAR, sino que también existen regiones de la sala en los que la presión aumenta (en torno a 10 dB), esto es mucho más real que el caso anterior en el que el CAR conseguía atenuaciones en toda la sala. Todo ello se debe a que el anterior alfa del PML era muy elevado, es decir, absorbía mucha energía, consiguiendo condiciones casi anecoicas en la sala. Por el contrario el nuevo alfa no tiene este efecto, por lo que se asemeja mejor a una sala real. En lo referente a la atenuación en este nuevo caso, vemos como el CAR sigue atenuando nuestra zona de interés, llegando a atenuaciones de 48 dB y siendo superior a 10 dB en la zona que nos interesa. De todas estas observaciones podemos sacar la conclusión de que hay que tener especial cuidado al elegir el alfa del PML, para simular lo mejor posible nuestro caso a estudio. También hay que prestar atención en no confundir el alfa del PML con el coeficiente de absorción tal y como lo entendemos en acústica. V.5. IDEA: ENCERRAR LA FUENTE DE RUIDO Una vez que ya hemos visto que la herramienta diseñada funcionaba correctamente, sabemos como expresar los resultados y tras observar las simulaciones anteriores, pasamos al siguiente paso, buscar alguna utilidad a este sistema. Viendo el ejemplo anterior se nos ocurrió intentar “encerrar” la fuente de ruido acústicamente, es decir, diseñar un sistema de control activo de ruido rodeando la fuente de ruido (senoidal de 100 Hz) con fuentes secundarias. Con la posición de los sensores de error no se pretende que el ruido sea mínimo en esos puntos, sino que el ruido sea también mínimo en todos los puntos de la sala. A continuación se muestra un esquema de cómo se monto la simulación: Fig.19. Esquema de la simulación. Simulador de recintos tridimensionales para el análisis del campo acústico en algoritmos de procesado de la señal. 23 Autor: Jorge Machín Mindán, email: [email protected] Fecha de entrega: 16-10-2010 En esta simulación se eligió un alfa (absorción PML) de 0.1, lo cual quiere decir que la absorción por parte del PML es baja, se diseño así para intentar simular mejor las condiciones de un sala cualquiera. En la siguiente gráfica se muestra la comparación de potencias medias durante 500 ms, ya llegado el estacionario del CAR. Con una señal de excitación de 100 Hz. Fig. 20. Representación de las potencias medias con y sin CAR escala logarítmica arbitraría. La representación sigue la misma estructura que la Fig.11. En la imagen anterior observamos como este sistema consigue elevadas atenuaciones en toda la sala. También se puede observar que las atenuaciones más altas no se encuentran en los sensores de error, pero esto no es importante, ya que el objetivo de este montaje era atenuar el ruido en toda la sala. Como hemos visto esta configuración es válida para baja frecuencia (100 Hz), pero que sucedería con frecuencias más altas, a continuación se muestra una imagen con el mismo diseño anterior, pero cambiando la frecuencia de emisión a 500 Hz, en ambos casos se eligen señales senoidales porque es más fácil ver su atenuación. Simulador de recintos tridimensionales para el análisis del campo acústico en algoritmos de procesado de la señal. 24 Autor: Jorge Machín Mindán, email: [email protected] Fecha de entrega: 16-10-2010 Fig. 21. Representación de las potencias medias con y sin CAR escala logarítmica arbitraría. Señal de excitación senoidal de 500 Hz. La representación sigue la misma estructura que la Fig.11. Como podemos ver en la gráfica anterior, el ruido también se ha atenuado en la sala, sin llegar a conseguir atenuaciones tan elevadas como en el caso anterior. Como se puede ver este diseño no esta preparado para frecuencias tan elevadas, pero este podría ser modificado para que actuase con mayor eficiencia: añadiendo más fuentes secundarias y colocándolas más cerca entre sí. V.6. SIMULACIÓN MULTIFUENTE De esta última aplicación surgió otra idea, la de intentar “encerrar” varias fuentes, para simular el ruido proveniente de una superficie vibrante, ventana o similar. Esta superficie emisora se simula con varias fuentes de ruido puntuales emitiendo en fase. Para realizar esta simulación se decidió posicionar cinco fuentes de ruido cerca de la pared, y seis fuentes antiruido, todas ellas emitiendo una señal senoidal de 100 Hz en fase. También se decidió posicionar 5 sensores de error, estos sensores se posicionaron a medio metro de las fuentes secundarias formando una línea vertical, con la idea de que si el ruido es mínimo en estos puntos, aquellos puntos que estén detrás de ellos también se verán atenuados. En las siguientes líneas se va a mostrar la diferencia de atenuaciones entre las dos simulaciones, con y sin CAR. A continuación observamos la comparación de las señales grabadas en un sensor en las dos simulaciones, la gráfica se encuentra en unidades lineales arbitrarias. Simulador de recintos tridimensionales para el análisis del campo acústico en algoritmos de procesado de la señal. 25 Autor: Jorge Machín Mindán, email: [email protected] Fecha de entrega: 16-10-2010 Fig.22. Comparación entre las señales grabadas con el CAR y sin él, en los sensores de error. Para poder comprobar los efectos del CAR en toda la sala mostramos, en la siguiente imagen, los campos acústicos resultantes de las dos simulaciones. Fig. 23. Representación de las potencias medias con y sin CAR escala logarítmica arbitraría. La representación sigue la misma estructura que la Fig.11. Como vemos se han conseguido atenuaciones en toda la sala, lo cual era nuestro objetivo. Simulador de recintos tridimensionales para el análisis del campo acústico en algoritmos de procesado de la señal. 32 Autor: Jorge Machín Mindán, email: [email protected] Fecha de entrega: 16-10-2010 Matriz: 100x100x100 100 ms Con CAR (4 F y 4 S) 73 min. Fig.34. Tabla que expresa los tiempos de cálculo de diferentes simulaciones. Viendo el problema que supone el coste computacional en tres dimensiones, sobretodo si aplicamos el CAR, se decidió intentar minimizar un poco este coste probando diferentes métodos de cálculo de los nuevos coeficientes del filtro adaptativo, sin dejar de lado el algortimo LMS con filtrado x. Se propusieron dos posibles soluciones, la técnica mininax y la técnica scanning, que a continuación se explican. De forma similar a la importancia de la elección del tipo de algoritmo a utilizar, la técnica de cálculo de los coeficientes es igual de importante, es decir, si conseguimos llegar al estacionario del CAR (que converja el filtro) tenemos que llegar al mismo resultado, independientemente del camino escogido. La misma idea puede ser aplicada a la estrategia para calcular los coeficientes del filtro, puede tardar más o menos, pero el resultado debe ser el mismo con las todas posibles técnicas [9]. - La estrategia “minimax” consiste en calcular los coeficientes del filtro adaptativo con el valor máximo de los sensores de error en valor absoluto. e_kmax= max(abs(e1),abs(e2), …., abs(ek)) Wij= Wij – 2*mu*(xf_ijk*e_kmax) (9) - La otra técnica, la de “scanning”, consiste en ir alternando la captura de datos de un sensor a otro en cada iteración. Wij=Wij-2*mu*(xf_ij1*e1)󳮁 (10) En la siguiente iteración: wij(n+1)=wij(n)-2*mu*(xf_ij2*e2) Así sucesivamente. Una vez definidas y programadas las técnicas a utilizar pasamos a las simulaciones. Utilizaremos la herramienta de dos dimensiones ya que el tiempo de cálculo se ve disminuido notablemente. A continuación mostraremos la evolución temporal de los coeficientes de los filtros adaptativos para ver si las tres técnicas de cálculo nos dan los mismos resultados. Y cuanto tiempo tardan en converger cada una de ellas. Simulador de recintos tridimensionales para el análisis del campo acústico en algoritmos de procesado de la señal. 33 Autor: Jorge Machín Mindán, email: [email protected] Fecha de entrega: 16-10-2010 Fig.35. Comparación de la evolución temporal de los coeficientes del filtro adaptativo con las distintas técnicas. Como podemos observar en la gráfica anterior, los dos primeros algoritmos convergen y llegan pronto a su solución, mientras que la tercera imagen parece que no llega a converger, siempre mantiene una oscilación alrededor de la solución. A continuación mostraremos los campos acústicos resultantes de la misma simulación con las diferentes técnicas. De arriba abajo técnica de media, scanning y minimax. Simulador de recintos tridimensionales para el análisis del campo acústico en algoritmos de procesado de la señal. 34 Autor: Jorge Machín Mindán, email: [email protected] Fecha de entrega: 16-10-2010 Fig.36. Comparación de potencias medias con distintas técnicas de cálculo. De arriba abajo, técnica media, sacanning y minimax. Tiempo de cálculo: Tamaño matriz Sistema CAR Media Scanning Minimax 150 x 100 500 ms 4 F y 4 S 430’’ 415’’ 420’’ 100 x 100 x 100 50 ms 4 F y 4 S 36’ 28’ 35’ Fig.37. Tabla que expresa el tiempo de cálculo en las distintas simulaciones. Tras ver los resultados de las tres técnicas y ver también el tiempo que lleva su procesado, recomendamos el uso de la técnica scanning para el cálculo de los nuevos coeficientes del filtro adaptativo, ya que las tres técnicas ofrecen resultados muy parecidos y es en el tiempo de cálculo en lo que se diferencian, siendo la técnica de scanning la que menos tiempo consume. Simulador de recintos tridimensionales para el análisis del campo acústico en algoritmos de procesado de la señal. 35 Autor: Jorge Machín Mindán, email: [email protected] Fecha de entrega: 16-10-2010 VII. CONCLUSIONES Creemos que hemos conseguido cumplir con los objetivos marcados, desarrollar una herramienta para la simulación tridimensional de diferentes algoritmos de procesado de la señal, y con ella poder observar los campos acústicos resultantes. Y tras ver los resultados de las simulaciones anteriores, con diferentes configuraciones, creemos que se trata de una herramienta robusta y fácil de utilizar. Podemos concluir que FDTD es una herramienta apropiada para ver como se comporta el CAR con el paso del tiempo, nos permite ver el campo acústico en cada momento, la evolución de los filtros del CAR, la señales en los sensores, el campo acústico en el transitorio, etc. Tras ver la comparación podemos llegar a la conclusión de que el sistema CAR puede ser utilizado para intentar aislar las fuentes de ruido de una forma acústica. Si bien es cuierto que debemos tener en cuenta la propagación del sonido para elegir la mejor colocación posible de nuestras fuentes secundarias y sensores de error. De todos modos un objetivo de este interfaz es la posibilidad de realizar distintas simulaciones para descubrir estas posiciones y así poder minimizar los efectos del ruido. Hay que prestar atención al alfa del PML, ya que no se puede comparar con el coeficiente de absorción tal y como lo entendemos en acústica, sino que es la absorción del PML. Mucho del tiempo invertido en la simulación es debido a las operaciones que tiene que realizar el CAR para funcionar, es por esto por lo que se propone el uso de la técnica scanning para ganar tiempo de cálculo. Que como hemos podido ver minimizan las operaciones a realizar en el algoritmo del CAR, si bien es cierto que pueden tardar más en converger los filtros adaptativos con estas técnicas, si el parámetro de convergencia “mu” no esta optimizado. Simulador de recintos tridimensionales para el análisis del campo acústico en algoritmos de procesado de la señal. 36 Autor: Jorge Machín Mindán, email: [email protected] Fecha de entrega: 16-10-2010 VIII. LÍNEAS FUTURAS Es cierto que estas interfaces se han diseñado con objeto de simular diferentes algoritmos, pero también se le pueden dar otras aplicaciones, como puede ser la caracterización de salas mediante la obtención de la respuesta al impulso, sobretodo a baja frecuencia debido a la relación entre tamaño de la malla y frecuencia máxima a utilizar. Otra aplicación puede ser didáctica con la explicación de fenómenos acústicos, difracción, superposición de ondas, efecto Doppler, etc. Una de las líneas futuras más interesantes puede ser la de realizar distintas simulaciones con diferentes algoritmos de procesado de la señal, ya que al final sólo se puedo realizar un algoritmo de Control Activo de Ruido en este trabajo, debido al poco tiempo disponible para la realización de este proyecto. Otra línea interesante sería la de bajar el tiempo de cálculo de las simulaciones, para ello debería depurarse el código, optimizar el cálculo numérico, usar hardware basado en GPU’s, etc. Por otro lado creo que una buena forma de bajar el tiempo de simulación es la de crear líneas de simetría para reducir el espacio de cálculo. Simulador de recintos tridimensionales para el análisis del campo acústico en algoritmos de procesado de la señal. 37 Autor: Jorge Machín Mindán, email: [email protected] Fecha de entrega: 16-10-2010 IX. BIBLIOGRAFÍA [1] Acoustics. L. L. Beranek. Acoustics Society of Ameria, American Institute of Physics. 1993. [2] Adaptation of FDTD tcheniques to acoustics modeling. J.G. Maloney y K.E. Cummings. 11 th Annual REview of Progress in Applied Computational Electromagnetics. 1981. [3] Finite Difference Time Domain Tutorial. Dr Ian Drumm, The University of Salford. [4] Evaluation of the Perfectly Matched Layer for Computational Acoustics. Thomas L. Geers. 2002. [5] Wave Field Synthesis simulation by means of finite-difference time-domain technique. J. Escolano. [6] Multichannel adaptive filtering for signal enhacement. B. Widrow. 1981. [7] Analysis of the filtered-X LMS algorithm. E. Bjarson. 1995. [8] New block filtered-X LMS algorithms for active noise control Systems. S.M. Kuo. 2007. [9] Filtrado multicanal para control local sonoro basado en algoritmos de proyección afín. M. Ferrer. 2008 Simulador de recintos tridimensionales para el análisis del campo acústico en algoritmos de procesado de la señal. 38 Autor: Jorge Machín Mindán, email: [email protected] Fecha de entrega: 16-10-2010 ANEXO I. Definición de Variables DEFINICIÓN DE LA MATRIZ DE PUNTOS DE CONTORNO El primer paso para la creación de las matrices que vamos a utilizar es definir una matriz donde guardaremos todos los puntos del contorno. Estos puntos son guardados en orden, cada punto, con sus respectivas coordenadas, en una fila. Ejemplo de definición del primer punto del contorno: datos.pos_puntos(1,:) = str2num(get(hObject,'String')); Se elige esta forma de definir los puntos para poder guardar todos los puntos en la misma variable y así poder definir el número de puntos que se necesite, inclusive un número mayor de puntos de los que nos permite definir la GUI. Otra razón que nos lleva a esta decisión es la de evitar la utilización de bucles a la hora de leer esta información. DEFINICIÓN DE LAS MATRICES DE CONTORNO, PMLs Y MATRICES DE SIMULACIÓN 2D Como ya hemos comentado antes los puntos del contorno se definen en orden, por lo que ya conocemos las dimensiones máximas de nuestra matriz, la cual representa la superficie de simulación. No sólo se crea la matriz en la que se desarrollará la simulación, también se crean otras matrices que utilizaremos como máscara a la hora del cálculo. Estas matrices son las encargadas de definir las condiciones de contorno (Pared rígida. Velocidad = 0). De una forma similar serán definidas las matrices para definir las paredes mediante PML, en la creación de estas matrices podemos elegir el alfa máximo de la matriz así como la longitud del PML. El primer paso para definir las matrices es calcular las dimensiones máximas de la matriz de simulación. Se calculan como los máximos en las dos dimensiones x e y. Con estos máximos se puede definir ya las matrices que guardarán los datos. puntos= datos.pos_puntos / datos.malla; % calculo de numero de nodos= tamaño/malla dimensiones= max(puntos); % es solo para el máximo %% Matrices para el PML nx= dimensiones(1); % En nx guardamos el máximo en el eje X ny= dimensiones(2); % En ny guardamos el máximo en el eje Y alfaux=zeros(nx+1,ny); % Matriz para guardar información del contorno en el eje X alfauy=zeros(nx,ny+1); matrizX= ones(dimensiones(1)+1,dimensiones(2)); % Matriz con las dimensiones máximas para el eje X matrizY= ones(dimensiones(1),dimensiones(2)+1); Simulador de recintos tridimensionales para el análisis del campo acústico en algoritmos de procesado de la señal. 39 Autor: Jorge Machín Mindán, email: [email protected] Fecha de entrega: 16-10-2010 Para realizar la matriz de contorno y de los PMLs, se cogen los puntos de contorno en parejas, de dos en dos, para posteriormente comprobar sí los puntos están en la misma recta x o y, y sabiendo en que recta está el contorno podemos definir ya el PML de esa pared. La matriz de contorno la definimos como una matriz de unos, menos los contornos que son 0, para poder multiplicarlos luego por la matriz de la velocidad (condición de contorno de velocidad nula en las paredes). Para definir los PMLs, creamos los PMLs hacia el interior, si la pared esta cerca del mínimo o del máximo, ya sabemos hacia que lado se define el PML. En el caso en que la pared sea interior definimos el PML de una forma triangular, es decir, creamos un PML con el máximo en la pared y decreciente hacia ambos extremos. Una vez realizado todo lo anterior, ya tenemos nuestras matrices definidas y preparadas para la simulación. Simulador de recintos tridimensionales para el análisis del campo acústico en algoritmos de procesado de la señal. 40 Autor: Jorge Machín Mindán, email: [email protected] Fecha de entrega: 16-10-2010 ANEXO II. EJEMPLO En estas líneas se va a intentar demostrar un fenómeno que se da en algunos sistemas de control activo de ruido. Este fenómeno consiste en que en algunas ocasiones se obtienen mayores atenuaciones en los sistemas con una fuente antiruido y un sensor, pero dicha atenuación se da en una pequeña zona alrededor del sensor, mientras que si añadimos más sensores al sistema, la atenuación es menor pero se da en una mayor superficie. Para llevar a cabo esta demostración se han realizado dos simulaciones, en la primera de ellas hay una fuente secundaría y un sensor de error, mientras que en la segunda hay una fuente y dos sensores de error. Esquema de las simulaciones: Fig.38y 39. Diseño de las simulaciones sutilizadas para el ejemplo. A continuación vamos a mostrar la señal grabada durante las dos simulaciones en el mismo punto. Fig.40. Comparación señal grabada en el sensor 1 durante las dos simulaciones. Simulador de recintos tridimensionales para el análisis del campo acústico en algoritmos de procesado de la señal. 41 Autor: Jorge Machín Mindán, email: [email protected] Fecha de entrega: 16-10-2010 La señal azul se corresponde con el sensor uno del caso 2 y la señal roja con el sensor del caso 1. Como se observa la señal de error es menor en el primer caso, tal y como habíamos comentado, mientras que como se observa en la siguiente imagen, el CAR tiene una mayor superficie de actuación (consigue menores atenuaciones, pero en un espacio mayor). Fig.41. Comparación de la potencia media en los dos casos. Imagen superior caso 1, imagen inferior caso 2. A continuación se procedió a calcular la atenuación conseguida al añadir un sensor, se calculó como el logaritmo de la división de la potencia media en el caso 1 entre la potencia media del caso 2: