Potenciales Ventajas del Método Pseudoespectral
en Auralizaciones
Potential Advantages of the Pseudospectral Method in Auralizations
Jorge Petrosino
#1
, Lucas Landini
#
, Andrés Bonino Reta
#
, Georgina Lizaso
#
#
Departamento de Humanidades y Artes. Universidad Nacional de Lanús
29 de Septiembre 3901, Remedios de Escalada, Buenos Aires, Argentina
1
jorgepetrosino@gmail.com
Abstract— Auralization is a term introduced to describe the
recreation of the experience of acoustic phenomena a listener
would perceive in a specific soundfield. Sound propagation in a
soundfield can be simulated with geometric based models or
wave based models. Each one offers particular advantages and
disadvantages. For wave based models, the finite element
method, the boundary element method or the finite difference
method are widely mentioned. They are characterized for
achieving very precise results for individual frequencies
applied to small and moderately sized rooms. Geometric
methods lead to the ray tracing method or the image source
method. These methods achieve good results for high
frequencies and are efficient in large rooms and complex
structures, but are not able to represent in a simple manner
specific wave phenomena such as diffraction. Commercial
software used to produce auralizations is usually based on a
hybrid model combining ray tracing and image sources. This
paper proposes an exploration on possible advantages and
challenges on the use of the k-space pseudospectral method for
wave based auralizations.
Keywords: auralization; pseudospectral method; numerical
methods
Resumen— El término auralización se refiere a la
recreación en los oídos de un oyente de la sensación que
percibiría en un espacio acústico determinado. La propagación
del sonido en un espacio acústico puede simularse mediante
procesos basados en el modelo geométrico o en el modelo
ondulatorio. Cada uno ofrece ventajas y dificultades
particulares. Entre los modelos basados en ondas pueden
mencionarse principalmente el método de elementos finitos, el
de contornos finitos o el de diferencias finitas. Se caracterizan
por lograr resultados muy precisos para frecuencias únicas
aplicadas a recintos de tamaño pequeño o medio. Por otro
lado, los modelos geométricos dan lugar al método del trazado
de rayos o al método de las fuentes imagen. Estos métodos
logran buenos resultados en frecuencias altas y resultan
eficientes en salas de gran tamaño con estructuras complejas,
pero no pueden dar cuenta en forma sencilla de fenómenos
específicamente ondulatorios como la difracción. Los
programas comerciales utilizados para obtener auralizaciones
suelen utilizar un modelo híbrido combinando el trazado de
rayos y las fuentes imagen. El presente trabajo tiene la
intención de iniciar una exploración sobre las posibles ventajas
y desafíos del uso del método pseudoespectral del espacio k
para obtener auralizaciones basadas en el modelo ondulatorio.
Palabras clave: auralización; método pseudoespectral;
métodos numéricos
I. INTRODUCTION
Vorländer define auralización como la técnica de crear
señales en el rango audible mediante métodos numéricos, a
partir de información simulada o medida que describa el
espacio acústico, las fuentes y posición del oyente [1]. Los
primeros intentos de auralización incluyeron modelos
físicos a escala [2].
El proceso se basa en la determinación de la respuesta al
impulso del espacio acústico a ser representado. Para ello se
requiere seleccionar un punto en el que estará localizada la
fuente en el espacio acústico. Dependiendo de la
complejidad de la simulación, el modelo podría incluir
patrones de directividad de la fuente sonora y de algunas
otras características específicas, como la respuesta al
impulso del sistema acústico o electroacústico que generará
el sonido. Es necesario, además, disponer de un modelo
matemático que represente el espacio acústico sobre el que
se simulará la propagación del sonido, así como
información sobre la ubicación del receptor y sus
características.
La respuesta al impulso del modelo suele obtenerse
mediante procesos basados en el modelo geométrico o en el
modelo ondulatorio, cada uno de los cuales presenta
ventajas y dificultades particulares.
Los modelos que se basan en acústica geométrica se
implementan mediante el método del trazado de rayos (ray
tracing method: RTM) o el método de las fuentes imagen
(image source method: ISM). Algunos paquetes comerciales
utilizan un método híbrido que combina la técnica de
trazado de rayos con la de fuentes imagen [3]. Estos
métodos logran buenos resultados en frecuencias altas y
resultan eficientes en salas de gran tamaño con estructuras
complicadas. Simulan la emisión de un impulso y sus
reflexiones en las distintas superficies hasta alcanzar el
punto de registro. La respuesta al impulso se obtiene en
forma de un registro de esas múltiples reflexiones que
forman lo que se conoce como ecograma. Estrictamente
hablando no incluye todas las características esperables de
una respuesta al impulso ya que el modelo geométrico no
puede dar cuenta en forma sencilla de fenómenos
específicamente ondulatorios como la difracción.
Entre los modelos basados en ondas pueden mencionarse
principalmente el método de elementos finitos (finite
element method: FEM) y el de contornos finitos (boundary
Recibido: 08/12/18; Aceptado: 18/06/19
Revista elektron, Vol. 3, No. 2, pp. 67-74 (2019)
ISSN 2525-0159
67
element method: BEM). Estos se caracterizan por lograr
resultados muy precisos para frecuencias únicas aplicadas a
recintos de tamaño pequeño o medio [1]. Estos métodos
resuelven las ecuaciones diferenciales convirtiendo las
variables continuas en variables discretas, lo que provoca un
efecto de dispersión numérica. Esto vuelve impráctico el
estudio de casos en los que la fuente emite una onda de
audio compleja (como la palabra hablada) o un impulso, y
por ello quedan limitados a simulaciones para frecuencias
únicas. Para poder obtener la respuesta al impulso se
requiere un complejo trabajo de solución independiente para
frecuencias puntuales o bandas estrechas de frecuencia
cuyos resultados luego deben combinarse para conformar
una respuesta al impulso global, lo que conlleva un aumento
considerable del tiempo de cómputo.
El presente trabajo propone el uso de un método diferente
de resolución de las ecuaciones diferenciales de
propagación de ondas para generar auralizaciones. El
método pseudoespectral del espacio k (k-space
pseudospectral method) podría permitir realizar
simulaciones tanto para ondas de audio complejas como
para un impulso en una única sesión de cómputo,
ubicándolo como una interesante alternativa a los métodos
ondulatorios tradicionales
A. Zonas de Validez de los Distintos Modelos de
Comportamiento Acústico de Recintos
A mediados del siglo XX, Manfred Schroeder analizó la
región de transición de la respuesta en frecuencia de
recintos que separa la región de baja frecuencia dominada
por modos separados y la región de frecuencias dominada
por una gran superposición de modos, en donde sería
posible considerar ciertas variables de tipo estadístico.
Esta frecuencia de transición depende del tiempo de
reverberación y del volumen del recinto. En 1954 [4],
Schroeder propuso obtener el valor de dicha frecuencia
como se muestra en la Ec. (1)
(1)
En la publicación original, Schoeder utiliza un valor de
K
s
= 4000. En 1962 [5], Schroeder y Kuttruff revisan dicha
propuesta, reduciendo el factor a la mitad K
s
= 2000. Fuchs
indica que es difícil considerar un valor inequívoco para
este factor, por lo que propone considerar una franja
indeterminada de frecuencias entre ambos valores de Ks
citando diversos autores que argumentan respecto de la
utilización de uno u otro valor [6].
Debido a la limitada capacidad de cálculo disponible en
la época de la publicación original, la única manera de tratar
con recintos de estructuras complejas era utilizando
métodos estadísticos, lo que volvía imprescindible conocer
los límites de validez de dichos métodos. La propuesta de
Schroeder y Kuttruff permitía determinar una frecuencia
límite. Los métodos estadísticos podían ofrecer resultados
adecuados para frecuencias superiores a fs. Por debajo de
esa frecuencia era necesario utilizar la teoría ondulatoria,
que sólo podía resolverse para casos limitados al uso de
formas geométricas sencillas. Por otra parte, el rango de
frecuencias por sí solo no es garantía de validez de un
modelo estadístico ya que se requiere incorporar otros
criterios, como por ejemplo el de la distribución uniforme
de materiales acústicos en el recinto.
La disponibilidad de computadoras hizo posible el
trabajo con modelos basados en acústica geométrica, que
asumen un comportamiento regido por reflexiones
especulares para analizar el comportamiento de recintos.
Para que los resultados obtenidos mediante estos modelos
puedan considerarse adecuados es necesario que las
longitudes de onda sean menores que el tamaño de las
superficies sobre las que impactan. Esta condición impone
restricciones importantes al rango de frecuencias ya que
para tener en cuenta la influencia de irregularidades
pequeñas en las superficies, el estudio queda limitado a
valores muy elevados de frecuencias. La transición entre el
comportamiento como onda y como rayo es gradual y
resulta difícil de determinar con un único valor de
frecuencia.
Es importante mencionar que tanto el modelo geométrico
como el estadístico tienen una frecuencia mínima de validez,
sin embargo los modelos ondulatorios no se encuentran
limitados en frecuencia. Su alcance sólo podría verse
comprometido por la capacidad de cómputo o por
dificultades en la definición de modelos precisos de los
elementos físicos incorporados al mismo.
La Fig. 1 muestra valores de fs para distintos tiempos de
reverberación y volúmenes entre 100 m
2
y 1200 m
3
.
Fig. 1. Frecuencia de Schroeder (fs) para distintos tiempos de
reverberación, considerando Ks = 2000.
B. Modelos Basados en Acústica Geométrica
Existen principalmente dos métodos que se basan en el
modelo geométrico: el de trazado de rayos y el de las
fuentes imagen.
El método del trazado de rayos genera un haz de rayos
que parte de la fuente. La densidad de rayos utilizada en el
análisis depende de la implementación particular, lo que
puede dar lugar en casos poco densos a que terminen
omitiéndose las reflexiones en alguna superficie de poca
extensión. Cuando se detecta que un rayo interseca un
elemento se genera una reflexión especular a partir de dicho
punto. Los métodos de rayos no tienen en cuenta la longitud
de onda específica, ni los valores instantáneos de presión o
velocidad de partículas, sino que toman en consideración la
energía transmitida. La reflexión en una determinada
superficie puede ser afectada por un coeficiente de
absorción. En casos más detallados es posible simular
superficies difusas, considerando que en el punto de
intersección existe una nueva fuente de energía sonora
Revista elektron, Vol. 3, No. 2, pp. 67-74 (2019)
ISSN 2525-0159
68
http://elektron.fi.uba.ar
emitiendo en todas direcciones según la ley del coseno de
Lambert. Este método no puede dar cuenta de los
fenómenos puramente ondulatorios, como interferencia o
difracción, aunque existen modos de incorporar procesos
complementarios para suplir esta falencia. La cantidad de
variantes propuestas es demasiado grande para dar cuenta
aquí de todas las alternativas. Como ejemplo de algunas de
estas variantes, podemos mencionar que existen versiones
que utilizan procesos estocásticos para emular la absorción,
haciendo que algunos rayos o partículas incidentes se
reflejen y otros sean eliminados.
El método de las fuentes imagen reemplaza cada
superficie reflectante por una fuente virtual a distancia,
como si se tratase del reflejo en un espejo plano. En
principio esto permite concentrar la capacidad de cómputo
solamente en aquellos rayos que unen a las fuentes imagen
con el punto de registro. Las fuentes creadas por la primera
reflexión constituyen las fuentes imagen de primer orden,
resultando necesario utilizar un orden tan alto como
cantidad de reflexiones deseen simularse. Los coeficientes
de absorción pueden tenerse en cuenta atenuando la
intensidad de las fuentes imagen. Cuanto más tiempo de
simulación se requiera, más fuentes imagen de orden
superior será necesario incorporar. La principal diferencia
en eficiencia con el trazado de rayos es que requiere menos
capacidad de cómputo para tiempos cortos, como los
correspondientes a fuentes imagen de órdenes bajos, pero
crece en forma exponencial superando al requerido por el
método de rayos al simular tiempos mayores.
El método de fuentes imagen no resulta adecuado para
trabajar con superficies curvas, para incluir reflexión difusa,
ni trabajar con espacios de formas complejas que den lugar
a fuentes imagen ocultas, como por ejemplo al considerar
recintos acoplados. La difracción no puede ser representada
mediante este método. Podría sin embargo incorporarse
información de fase para simular interferencia. Esta
aparente ventaja no termina de ser suficientemente valiosa,
ya que los efectos de interferencia tienen mayor influencia
en la zona de bajas frecuencias, donde no se cumplen los
postulados del modelo geométrico.
El método de las fuentes imagen tiene mejor resolución
temporal, sin embargo el de trazado de rayos posibilita la
inclusión del fenómeno de difusión. Es por esto que un
método híbrido, como adecuada combinación de ambos,
puede proporcionar mejores aproximaciones a la respuesta
al impulso de un recinto real. En general los programas de
simulación acústica comerciales actuales utilizan las fuentes
imagen para determinar las reflexiones tempranas y el
trazado de rayos para las tardías, siendo un tema de
investigación científica actual el determinar el tiempo
adecuado de transición entre ambos [1].
Para poder tener en cuenta características que dependan
de la frecuencia (como la absorción o la difusión) el trazado
de rayos debe repetirse para cada una de las bandas de
frecuencia determinadas en el análisis. Cada banda de
frecuencia da lugar a un reflectograma (o ecograma) que
posteriormente debe combinarse con el resto para obtener la
respuesta al impulso del recinto. Esta respuesta al impulso
luego se convoluciona con un archivo de audio para generar
una señal representativa de aquello que escucharía un
oyente ubicado en el lugar seleccionado.
C. Modelos Basados en Acústica Ondulatoria
Entre los modelos ondulatorios comúnmente utilizados
pueden citarse el método de elementos finitos (FEM) y el de
contornos finitos (BEM). Estos modelos permiten en
principio tener en cuenta todos los fenómenos ondulatorios,
por lo que no tienen limitación en cuanto a la utilización de
superficies curvas, ni a la simulación de efectos de
difracción, difusión o interferencia. Los elementos presentes
en el recinto no se representan mediante coeficientes como
los de absorción o difusión, sino en términos de sus
propiedades físicas como la densidad, la velocidad de
propagación del sonido interna, o el módulo de
compresibilidad.
Con estos métodos se obtienen resultados muy precisos
pero a un alto costo computacional. Uno de los
inconvenientes de los métodos numéricos de solución de las
ecuaciones diferenciales es que el tipo de procesamiento
realizado genera dispersión numérica en las soluciones [7].
Esto provoca que la velocidad de propagación varíe con la
frecuencia, lo que obliga a realizar simulaciones para
frecuencias puntuales y luego superponer los resultados para
obtener la respuesta al impulso final del recinto.
Tanto en el método de elementos finitos, como en el de
contornos finitos el espacio acústico se divide en pequeños
elementos de tamaño adecuado a las necesidades de la
simulación. Para cada elemento se establecen relaciones
entre las presiones y los desplazamientos (a través de
ecuaciones diferenciales en derivadas parciales), que luego
dan lugar en conjunto a la solución global aproximada.
Ambos reemplazan la tarea de resolver ecuaciones
complejas sobre volúmenes o superficies complejas, por la
de resolver una gran cantidad de ecuaciones simples sobre
una gran cantidad de elementos simples.
La diferencia relativa entre ambos es que el método de
elementos finitos resuelve la ecuación de onda en todo el
espacio, mientras que el de contornos finitos resuelve las
ecuaciones en el contorno que delimita el espacio de interés
lo que significa menos elementos de cálculo. Sin embargo,
esta aparente ventaja relativa es opacada en cuanto a su
eficiencia computacional porque las matrices necesarias
para elementos finitos son poco densas y simétricas,
mientras que las de contornos finitos son densas y
asimétricas [8].
Por otra parte, al resolver la ecuación de onda en todo el
espacio, el tiempo de procesamiento requerido en elementos
finitos no varía con la cantidad de puntos de registro, por lo
que pueden obtenerse resultados para diferentes lugares del
recinto en forma conjunta.
Estos métodos se utilizan actualmente para el cálculo en
recintos de tamaño pequeño o mediano y para frecuencias
bajas [9]. La referencia a recintos pequeños o medianos se
hace porque es allí donde las frecuencias de Schroeder son
relativamente elevadas como muestra la figura 1 y en ellos
los métodos geométricos no resultan satisfactorios. Para
recintos de mayor tamaño la frecuencia de Schroeder es
razonablemente baja y el problema podría resolverse con los
métodos geométricos que son más eficientes.
Es importante prestar atención a que tanto los métodos
geométricos como los ondulatorios mencionados
anteriormente requieren, por motivos diferentes, realizar
procesos de simulación considerando frecuencias únicas o
bandas estrechas de frecuencias. Que luego componen los
Revista elektron, Vol. 3, No. 2, pp. 67-74 (2019)
ISSN 2525-0159
69
http://elektron.fi.uba.ar
resultados que dan lugar a la respuesta al impulso. Desde un
punto de vista computacional, un posible método que
permitiese obtener la respuesta al impulso en forma
completa durante un único proceso de simulación podría
presentarse como una interesante alternativa, aún cuando el
método en sí mismo no ofreciera particulares ventajas desde
los requerimientos generales de cálculo.
II. MÉTODO PSEUDOESPECTRAL DEL ESPACIO K
Los métodos espectrales son una clase de técnicas
numéricas que resuelven ecuaciones diferenciales. La idea
central es la de obtener la solución como suma de ciertas
funciones básicas (sinusoidales, por ejemplo) y luego
calcular los coeficientes que mejor satisfagan la ecuación
diferencial [10].
Tienen relación cercana con los métodos de elementos
finitos. La principal diferencia es que los métodos
espectrales utilizan funciones básicas que son no nulas en
todo el dominio, mientras que los elementos finitos utilizan
funciones que son no nulas en pequeños subdominios.
Podría decirse que los métodos espectrales adoptan un
enfoque global mientras que los de elementos finitos
utilizan un enfoque local [11] [12].
Los métodos pseudoespectrales son una subclase de los
métodos espectrales. La diferencia principal es que han sido
adaptados para representar las funciones en una grilla
rectangular, lo que simplifica la evaluación de ciertos
operadores y puede incrementar la velocidad de cálculo
permitiendo el uso de la transformada rápida de Fourier. El
nivel de precisión de los métodos espectrales y los
pseudoespectrales es muy similar [13] [14].
Cuando una onda acústica pasa a través de un medio se
producen variaciones en la presión, la densidad, la
temperatura, la velocidad de las partículas y otras variables.
Estos cambios pueden describirse mediante ecuaciones
diferenciales de primer orden acopladas que tienen en
cuenta la conservación del momento, de la masa y de la
energía en el medio. Es común en acústica combinar estas
ecuaciones en relación con una variable (que puede ser la
presión sonora o la velocidad de las partículas) formando
una ecuación diferencial de segundo orden.
El método pseudoespectral del espacio k aplicado a la
propagación de ondas acústicas resuelve las ecuaciones
diferenciales de primer orden, en lugar de utilizar la
ecuación combinada de segundo orden [15]. Esto se hace
por dos motivos. En primer lugar, porque permite incluir de
manera sencilla fuentes de masa y de fuerza en las
ecuaciones discretas. En segundo lugar, porque habilita a
utilizar una capa especial anisotrópica conocida como capa
perfectamente acoplada (PML: perfectly matched layer) con
el fin de absorber las ondas acústicas que alcanzan los
bordes del dominio utilizado para el cálculo [16].
El método pseudoespectral del espacio k ha sido utilizado
profusamente en los últimos años en aplicaciones
biomédicas relacionadas con la propagación de ondas
ultrasónicas en tejidos [17] [18]. También ha sido utilizado
para el estudio de la alteración producida por variaciones de
la temperatura atmosférica en la propagación de ondas
mecánicas [19]. Hemos hallado muy pocos trabajos que
utilicen métodos espectrales en acústica de recintos [20].
Hace poco menos de una década Treeby, Cox (University
College London) y Jaros (Brno University of Technology)
desarrollaron un paquete complementario (toolbox) para
MATLAB u Octave que implementa de modo eficiente el
método subespectral del espacio k. El paquete fue diseñado
para permitir simulaciones en el dominio del tiempo de
ondas acústicas propagándose en medios complejos. Si bien
su uso puede ser más general, fue diseñado especialmente
para el estudio de la propagación de ondas ultrasónicas en
tejido biológico en relación con imaginería médica [21].
Nuestro equipo viene trabajando con el paquete k-wave
en diversas aplicaciones acústicas. Resulta muy sencilla la
incorporación de un número ilimitado de fuentes sonoras
con comportamientos temporales que pueden ser
arbitrariamente definidos por el usuario utilizando cualquier
tipo de función de MATLAB/Octave. Permite realizar
simulaciones en una, dos o tres dimensiones. Pueden
combinarse fuentes de presión y de velocidad. También
pueden incorporarse una cantidad arbitraria de puntos de
registro de presión y/o velocidad, sin límite y sin que esto
implique mayor tiempo computacional, ya que el método
resuelve siempre las ecuaciones en todos los puntos de su
grilla rectangular. Por otra parte, pueden especificarse
valores de densidad, velocidad de propagación y coeficiente
de atenuación en cada punto de la grilla en forma
independiente.
Los métodos numéricos tradicionales basados en el
modelo ondulatorio producen errores en la estimación del
gradiente, lo que da lugar a una dependencia de la velocidad
del sonido con la frecuencia, generando una clase de
dispersión numérica sin sentido físico. Los métodos
pseudoespectrales reemplazan el cálculo de los gradientes
espaciales mediante diferencias finitas por el cálculo
espectral del gradiente por medio de series de Fourier. Esto
evita la dispersión debida a la discretización espacial pero
no aquella que surge de la aproximación de gradientes
temporales por diferencias finitas. Sin embargo, el error
cometido por el salto discreto temporal puede corregirse
para un medio homogéneo hasta el límite de Nyquist [22].
Estas características resultan suficientemente interesantes
como para iniciar un análisis exploratorio sobre la
potencialidad de su utilización en auralizaciones intentando
obtener la respuesta al impulso en un solo ciclo de
simulación.
III. PROPAGACIÓN DE UN IMPULSO EN 3D
Los métodos clásicos que implementan el modelo
ondulatorio trabajan simulando el comportamiento ante
frecuencias puras o bandas estrechas para luego componer
los resultados en una única respuesta al impulso general.
Esto implica generalmente repetir el proceso innumerables
veces, como queda claramente expresado por Aretz y cols.
[23] al indicar que: "los resultados fueron calculados para el
rango de frecuencias entre 200 Hz y 2500 Hz con saltos de
1 Hz". Asumiendo un tiempo de cómputo de 5 minutos por
frecuencia [1] esto tomaría más de 190 horas continuas de
cómputo.
En esta sección se explora la posibilidad de obtener la
respuesta al impulso en forma directa simulando la
propagación en un espacio de tres dimensiones de una
fuente impulsiva en una sola sesión de cómputo. La
hipótesis de trabajo es que la simulación de la propagación
de un impulso en un espacio de tres dimensiones desde la
fuente hasta un sensor ubicado a cierta distancia provocará
Revista elektron, Vol. 3, No. 2, pp. 67-74 (2019)
ISSN 2525-0159
70
http://elektron.fi.uba.ar
un registro que mantendrá las características del impulso
original, con poca o nula deformación.
Tanto en los métodos pseudoespectrales como en los de
elementos finitos se analiza la estabilidad de las soluciones
mediante el coeficiente descripto en la Ec. (2) conocido
como CFL (Courant - Friedrichs - Lewy). Corresponde a la
razón entre la distancia que una onda avanzará en un salto
temporal elemental y el salto de grilla.
(2)
Los valores particulares de CFL que corresponden a
condiciones de estabilidad cambian cuando se comparan
modelos numéricos diferentes, como los pseudoespectrales
y los de elementos finitos. Por ejemplo, un CFL = 0.3 con 2
puntos por longitud de onda resultará en un salto temporal
Δt cinco veces mayor al que corresponde a un modelo de
diferencias finitas con 10 puntos por longitud de onda y el
mismo valor de CFL [24] [25]. En las simulaciones se
utilizó un valor CFL = 0.3 por presentar un adecuado
balance entre precisión y velocidad de cómputo para medios
homogéneos o ligeramente heterogéneos. El valor de salto
temporal que utiliza internamente el k-wave se obtiene
mediante la Ec. (3)
(3)
Considerando valores de CFL = 0.3, Δx = 1 cm, cmax =
343 m/s, obtenemos un Δt = 8.75 µs. Este Δt no se
corresponde con un período de muestreo en el sentido
utilizado normalmente en procesamiento de audio. La
frecuencia máxima de validez para las funciones temporales
que corresponden a las fuentes o a los registros se obtiene
mediante la Ec. (4), donde cmin es la velocidad mínima de
propagación en el medio.
(4)
Se llevó a cabo la simulación de la propagación de una
señal impulsiva en un espacio tridimensional de
256x256x256 puntos de grilla, lo que da un total de más de
16 millones de nodos, representando un espacio cúbico de
2,56 m de lado con un salto discreto de 1 cm para cada una
de las dimensiones espaciales. La frecuencia máxima de
validez obtenida mediante la Ec. (4) fue de 17150 Hz. La
fuente se ubicó en las coordenadas (0 m; -1 m; 0 m), con el
sensor ubicado en (0.5 m; 0.5 m; 0.5 m). Esto da una
distancia entre la fuente y el sensor de 1.66 m. El tiempo de
cómputo en una computadora personal (procesador i7 de 3.4
GHz con 16 GHz de RAM) fue de 32 min. La Fig. 2
muestra el espectro del registro del sensor. El caso ideal
implicaría que el registro tuviese un espectro de módulo
constante para todas las frecuencias. Los cálculos previos
predicen una frecuencia máxima válida para la simulación
de 17150 Hz. Puede observarse en nuestros resultados que
si se considera aceptable una variación menor a 1 dB en la
respuesta al impulso, el rango de frecuencia debería
limitarse en este caso a 14 kHz.
Estos resultados son perfectamente compatibles con la
hipótesis de que es posible obtener la respuesta al impulso
completa en un sólo proceso de simulación.
Fig. 2. Espectro de la respuesta al impulso (módulo de la función de
transferencia del sistema)
El tiempo de cómputo reportado en esta sección sólo tuvo
en cuenta la propagación desde la fuente hasta el sensor, con
el fin de poner a prueba la hipótesis de que el impulso no se
altera al propagarse, a diferencia de lo que sucede en los
métodos de elementos finitos. Dicho tiempo de cómputo
aumentará si se intenta obtener auralizaciones en un recinto
para poder registrar las sucesivas reflexiones. De todas
maneras, estos tiempos resultan inferiores que los
esperables por la repetición de simulaciones mediante
elementos finitos para cada frecuencia por separado.
Es importante insistir aquí que dado que se resuelven las
ecuaciones diferenciales en todo el espacio, el tiempo de
cómputo no depende en forma directa de la cantidad de
sensores incorporados a la simulación, por lo cual pueden
utilizarse gran cantidad de puntos de registro, obteniendo en
la misma sesión de cómputo múltiples auralizaciones.
IV. PROPAGACIÓN DE ONDAS COMPLEJAS EN 2D
Como se mencionó anteriormente, el proceso de
auralización consiste en obtener la respuesta al impulso en
un recinto para convolucionar luego esta respuesta al
impulso con una señal de audio. Este proceso genera un
audio final con las características de lo que sería registrado
por un oyente ubicado en donde se encuentra el sensor en la
simulación.
Claramente simular solamente la propagación de un
impulso supone una ventaja en tiempo de cómputo respecto
de lo que implicaría simular la propagación de todo el audio
ya que el registro de la respuesta al impuso sólo requiere
simular el tiempo necesario para que las reflexiones
disminuyan hasta un nivel considerado suficiente para la
simulación. Este tiempo de simulación es independiente de
la duración del audio. Se asume que la convolución de la
respuesta al impuso con el archivo de audio sin procesar
genera un resultado similar a la que provocaría el audio si
fuese emitido desde la fuente. Para que esto sea válido es
necesario que el sistema sea lineal e independiente en el
tiempo (sistema LTI, por linear time invariant). La entrada
al sistema es la señal emitida por la fuente, la salida es el
registro y el sistema es todo el espacio por el cual las ondas
se propagan.
Con el método pseudoespectral del espacio k fue posible
realizar una prueba de concepto para evaluar la linealidad
del sistema al considerar la propagación de ondas. Esta
prueba consistió en obtener una auralización de la
Revista elektron, Vol. 3, No. 2, pp. 67-74 (2019)
ISSN 2525-0159
71
http://elektron.fi.uba.ar
propagación de todo el audio completo correspondiente a un
registro de la palabra "hola" y comparar este resultado con
la obtención de una respuesta al impulso y su posterior
convolución con el registro de la palabra original.
La simulación de la propagación de un audio complejo
no resulta posible con los métodos tradicionales, pero puede
en cambio desarrollarse con el método pseudoespectral. El
tiempo de cómputo para simular la propagación de un audio
completo es muy elevado, por lo cual se decidió para esta
prueba de concepto utilizar una versión limitada a dos
dimensiones, lo que disminuye enormemente la cantidad de
nodos a considerar.
A. Propagación Directa de un Archivo de Audio
En la simulación realizada se utilizó una única velocidad
de propagación en medio homogéneo (aire), con lo cual
cmin = cmax = 343 m/s. Con los valores seleccionados se
obtuvo una frecuencia fmax = 17150 Hz
Se definió un espacio de propagación de 2,56 m de lado.
Estos valores fueron elegidos para dar lugar a una grilla
cuadrada de 256 puntos de lado (65536 puntos en la grilla)
que resulta computacionalmente eficiente para el método
pseudoespectral. El toolbox k-wave dispone de una función
denominada checkFactors() que facilita tomar decisiones
respecto de valores de grilla eficientes para el cálculo.
Se colocó una fuente en x0 = 0 m; y0 = - 1 m y un sensor
de presión sonora en x0 = 0 m; y0 = 1 m. Se impuso a la
fuente de presión la condición de ser tipo Dirichlet
(obligando de este modo a que el valor de presión en el
punto correspondiente a la fuente sea exactamente el
especificado). En esta experiencia se pretendía analizar
solamente el comportamiento de propagación por lo que el
sonido se propaga en el aire sin obstáculos ni paredes
Fig. 3. Representación espacial de un instante en la propagación de una
señal de audio a lo largo del espacio definido en la simulación
El desarrollo temporal de la fuente fue especificado a
partir de la información registrada previamente en un
archivo de audio de palabra hablada. Se registró una palabra
breve (hola) cuya duración temporal fue de 317 ms. El
audio fue registrado con una frecuencia de muestreo de 48
kHz.
Dado que k-wave opera con un Δt muy inferior al del
período de muestreo, fue necesario remuestrear el audio
para que se adapte al Δt utilizado en la solución de las
ecuaciones diferenciales de propagación. Se aplicó la
función resample() para llevar la señal de audio original
muestreada a 48 kHz a una señal muestreada a 114 kHz.
Este audio remuestreado se constituyó en la información
de desarrollo temporal de la fuente de presión tipo Dirichlet.
La propagación de toda la información de audio a lo largo
de los 2 metros que separaban la fuente del sensor tomó un
tiempo de cálculo de 4 min en una computadora personal de
escritorio (procesador i7 de 3.4 GHz con 16 GHz de RAM).
El toolbox ofrece al comienzo de la simulación una
estimación del tiempo de cómputo total que requerirá la
simulación. La Fig. 3 muestra un instante de la simulación
de propagación del audio.
El audio registrado por el sensor fue reproducido
directamente en MATLAB con una frecuencia de muestreo
de 114 kHz. Resultó perfectamente claro, con un ligero
incremento de frecuencias bajas.
Fig. 4. Diagramas temporales del audio original emitido por la fuente y el
audio registrado por el sensor.
La Fig. 4 muestra los diagramas temporales del audio
emitido por la fuente y el audio captado por el sensor
ubicado a 2 m de distancia. Puede notarse un breve retardo
en el registro correspondiente al tiempo de propagación
entre la fuente y el sensor.
B. Propagación de un Impulso
Para llevar a cabo la segunda simulación se mantuvo la
configuración anterior utilizando una fuente impulsiva. Se
eligió un tiempo total de procesamiento ligeramente
superior al de propagación del impulso (t_final = 6.3 ms).
En este caso el tiempo total de cómputo fue de sólo 6
segundos. La Fig. 5 muestra el diagrama temporal y el
espectral de la respuesta al impulso. El incremento en
frecuencias bajas que se observa en la respuesta en
frecuencia es producto de la utilización de una fuente
puntual de tipo Dirichlet cuando se realiza una simulación
en dos dimensiones.
Revista elektron, Vol. 3, No. 2, pp. 67-74 (2019)
ISSN 2525-0159
72
http://elektron.fi.uba.ar
Fig. 5. Diagramas temporales de la señal registrada (respuesta al impulso)
y su espectro.
C. Convolución de la Respuesta al Impulso con el Audio
Se convolucionó la respuesta al impulso con el audio
original remuestreado a 114 kHz (para respetar el Δt de
procesamiento de k-wave). El resultado fue nuevamente un
audio claro y perfectamente inteligible, con un ligero
incremento en frecuencias bajas. Tal como era de esperar,
el audio recuperado posee un retardo de poco menos de 6
ms respecto del audio original, idéntico al de la simulación
de propagación del audio completo.
Fig. 6. Diagramas temporales del audio original emitido por la fuente y el
audio registrado por el sensor.
V. CONCLUSIONES Y TRABAJO FUTURO
Los resultados reportados resultan alentadores para
continuar la exploración en profundidad de los métodos
pseudoespectrales con miras a utilizarlos en auralizaciones.
Las pruebas realizadas muestran que es posible obtener una
respuesta al impluso en un solo ciclo de cómputo, y que el
resultado de la propagación de una onda de audio completa
coincide con la propagación de un impulso y su posterior
convolución con el audio original, en forma consistente con
lo que se espera de un sistema LTI. La función de
transferencia en 3D resulta prácticamente plana en un rango
de frecuencias adecuado para sus aplicaciones en audio. En
las simulaciones en 2D se verifica un incremento en
frecuencias bajas que podría ser compensado mediante
filtros para obtener una respuesta plana.
Hemos realizado pruebas preliminares en ambientes de
mayor tamaño incluyendo una variedad de configuraciones
de superficies reflectantes con interesantes resultados. Estos
trabajos necesitan ser verificados y sometidos a un análisis
más cuidadoso antes de ser reportados.
Parte de nuestro equipo de investigación ha comenzado a
explorar con diferentes materiales (densidad, velocidad de
propagación, absorción), con la intención de repetir las
experiencias de propagación de audio en situaciones más
complejas. Resulta necesario continuar con una mayor
exploración de los límites sobre los cambios que pueden ser
introducidos en las características del medio sin que se
pierda la posibilidad de simular la propagación de un
impulso en un único proceso de cómputo.
AGRADECIMIENTOS
Los autores de este trabajo quieren agradecer muy
especialmente al equipo conformado por Damián Fernández,
Ianina Canalis, Ian Kuri y Nicolás Casais Dassie por sus
aportes durante todo el proceso.
El presente trabajo es parte es parte de un proyecto de
investigación que se encuentra en desarrollo en la
Universidad Nacional de Lanús.
REFERENCIAS
[1] M. Vorländer, Auralization: fundamentals of acoustics, modelling,
simulation, algorithms and acoustic virtual reality, Berlin, Alemania:
Springer Science & Business Media, 2008
[2] D. J. Meares, “The Use of Scale Models in the Acoustic Design of
Studios,” SMPTE J., vol. 87, no. 10, pp. 684-687, Oct. 1978
[3] G. M. Naylor, “ODEON-Another hybrid room acoustical model,”
Appl. Acoust., vol. 38, no. 2-4, pp. 131-143, 1993
[4] M. Schroeder, “Die statistischen Parameter der Frequenzkurven von
groβen Räumen,” Acta Acust. united Ac., vol. 4, no. 5, pp. 594-600,
Jan. 1954
[5] M. R. Schroeder, y K. H. Kuttruff, “On frequency response curves
in rooms. Comparison of experimental, theoretical, and Monte Carlo
results for the average frequency spacing between maxima,” J.
Acoust. Soc. Am., vol. 34, no. 1, pp. 76-80, 1962
[6] H. V. Fuchs, Applied Acoustics: Concepts, Absorbers, and Silencers
for Acoustical Comfort and Noise Control: Alternative Solutions-
Innovative Tools-Practical Examples, Berlin, Germany: Springer
Science & Business Media, 2013
[7] J. O. Smith, Physical audio signal processing: For virtual musical
instruments and audio effects. San Francisco, CA, USA: W3K
Publishing, 2018
[8] L. Chen, L. Liu, W. Zhao, y H. Chen, “2D acoustic design
sensitivity analysis based on adjoint variable method using different
types of boundary elements,” Acoust. Aust., vol. 44, no. 2, pp. 343-
357, Jul. 2016
[9] M. Vorländer, “Computer simulations in room acoustics: Concepts
and uncertainties,” J. Acoust. Soc. Am., vol. 133, no. 3, pp. 1203-
1213, Mar. 2013
[10] C. Canuto, M.Y. Hussaini, A. Quarteroni, y T. A. Zang, Spectral
methods: Fundamentals in single domains, Berlin, Germany:
Springer-Verlag, Berlin, 2006
[11] A. Solomonoff, y E. Turkel, “Global properties of pseudospectral
methods,” J. Comput. Phys., vol. 81, no. 2, pp. 239-276, Apr. 1989
[12] G. T. Huntington, y A. V. Rao, “A comparison between global and
local orthogonal collocation methods for solving optimal control
problems,” en 2007 Amer. Contr. Conf. IEEE, pp.1950-1957.
[13] S. A. Orszag, “Comparison of pseudospectral and spectral
approximation,” Stud. Appl. Math., vol. 51, no. 3, pp. 253-259, 1972
[14] J. Virieux, H. Calandra, y R. E. Plessix, “A review of the spectral,
pseudo-spectral, finite-difference and finite-element modelling
Revista elektron, Vol. 3, No. 2, pp. 67-74 (2019)
ISSN 2525-0159
73
http://elektron.fi.uba.ar
techniques for geophysical imaging,” Geophys. Prospect., vol. 59,
no. 5, pp. 794-813, Sep. 2011
[15] K. Firouzi, B. Cox, B. Treeby, y N. Saffari, “A first-order k-space
model for elastic wave propagation in heterogeneous media,” J.
Acoust. Soc. Am., vol. 129, no. 4, pp. 2611-2611, Apr. 2011
[16] Q. H. Liu, “The PSTD algorithm: A time ‐ domain method
combining the pseudospectral technique and perfectly matched
layers,” J. Acoust. Soc. Am., vol. 101, no. 5, pp. 3182, May. 1997
[17] Y. Jing, F. C. Meral, y G. T. Clement, “Time-reversal transcranial
ultrasound beam focusing using a k-space method,” Phys. Med. Biol.,
vol. 57, no. 4, pp. 901, Jan. 2012
[18] B. Treeby, M. Tumen, y B. Cox, “Time domain simulation of
harmonic ultrasound images and beam patterns in 3D using the k-
space pseudospectral method,” en Int. Conf. Medical Image
Computing Computer-Assisted Intervention, Berlin, 2011, pp. 363-
370.
[19] J. Hargreaves, P. Kendrick, S. Von Hünerbein, “Simulating acoustic
scattering from atmospheric temperature fluctuations using a k-space
method,” J. Acoust. Soc. Am., vol. 135, no. 1, pp. 83-92, Jan. 2014
[20] F. Pind, M. S. Mejling, A. P. Engsig-Karup, C. H. Jeong, y J.
Strømann-Andersen, “Room Acoustic Simulations using High-Order
Spectral Element Methods,” en Euronoise 2018, pp. 2085-2092.
[21] B. E. Treeby, y B. T. Cox, “k-Wave: MATLAB toolbox for the
simulation and reconstruction of photoacoustic wave fields,” J.
Biomed. Opt., vol. 15, no. 2, Mar. 2010, Art. no. 021314.
[22] J. L. Robertson, B. T. Cox, y B. E. Treeby, “Quantifying numerical
errors in the simulation of transcranial ultrasound using
pseudospectral methods,” en IEEE Int. Ultrasonics Symposium
(IUS), 2014, pp. 2000-2003.
[23] M. Aretz, R. Nöthen, M. Vorländer, y D. Schröder, “Combined
broadband impulse responses using FEM and hybrid ray-based
methods,” en Proc. EAA Symp. Auralization, Jun. 2009, pp. 201-206.
[24] Q. H. Liu, “The pseudospectral time-domain (PSTD) algorithm for
acoustic waves in absorptive media,” IEEE T. Ultrason. Ferr., vol.
45, no. 4, pp. 1044-1055, Jul. 1998
[25] G. Zhao, y Q. H. Liu, “The unconditionally stable pseudospectral
time-domain (PSTD) method,” IEEE Microw. Wirel. Co., vol. 13, no.
11, pp. 475-477, Nov. 2003
Revista elektron, Vol. 3, No. 2, pp. 67-74 (2019)
ISSN 2525-0159
74
http://elektron.fi.uba.ar
Enlaces de Referencia
- Por el momento, no existen enlaces de referencia
Copyright (c) 2019 Jorge Petrosino
This work is licensed under a Creative Commons Attribution-NonCommercial-NoDerivatives 4.0 International License.
Revista elektron, ISSN-L 2525-0159
Facultad de Ingeniería. Universidad de Buenos Aires
Paseo Colón 850, 3er piso
C1063ACV - Buenos Aires - Argentina
revista.elektron@fi.uba.ar
+54 (11) 528-50889