Procesamiento Digital de Señales
Diciembre 2014
APLICACIÓN DE SEÑALES ESS EN LA CARACTERIZACIÓN DE SISTEMAS IGNACIO CALDERÓN DE PALMA 1, MARIANO SERATTIN 2 Y MARIELA ALBA3 Universidad Nacional de Tres de Febrero, Ingeniería de Sonido, Procesamiento Digital de Señales, Buenos Aires, Argentina.
[email protected] 2
[email protected] [email protected]
RESUMEN – Este trabajo presenta fundamentos teóricos de señales conocidas como sweeps senoidales exponenciales (ESS) y su implementación en la caracterización de sistemas. Su aplicación no se limita a sistemas lineales e invariantes en el tiempo, sino que se extiende a aquellos que presentan alinealidades. Esto último demuestra la potencia e importancia de las técnicas aquí expuestas. La validez del método se demostró de manera práctica con aplicación a dos mediciones realizadas en dos sistemas físicos reales: un recinto y un filtro analógico. El procesamiento de señales se llevó a cabo en MATLAB®, donde se desarrollaron algoritmos destinados a las mediciones y a la comprensión teórica de conceptos relacionados con las ESS.
ABSTRACT – This work presents both theory and implementation of exponential sine sweep signals (ESS) in system characterization. Their application is not limited only to linear systems, but also to non-linear ones, this demonstrates the power behind the techniques developed throughout this paper. This method’s validity was expressed by making measurements over two physical systems: a room and an analogue filter. Signal processing was carried out in MATLAB®, where algorithms destined to measurements and theoretical concepts were developed.
1. INTRODUCCIÓN El análisis de sistemas no lineales es una tarea laboriosa a nivel teórico y numérico, por lo que, en general, a fines prácticos se los suele aproximar a modelos lineales o simplemente se asume su linealidad. Sin embargo, en este trabajo se presenta un método, ampliamente utilizado, basado en sweeps senoidales, comúnmente llamados chirp. El objetivo de este trabajo es demostrar la efectividad y aplicación de las señales del sweep en el análisis de sistemas reales que no necesariamente cumplen con las condiciones de linealidad e invariancia en el tiempo. La primera parte del trabajo desarrolla algunos conceptos teóricos claves para comprender el funcionamiento del método de ESS (Exponential Sine Sweep). La segunda parte describe ciertas ventajas y posibles problemas que se pueden encontrar a la hora de implementar el método. La tercera parte muestra la implementación real de las ESS en dos mediciones distintas, una midiendo la respuesta del filtro de una consola de audio y la otra describiendo los parámetros de sala de un aula de la Universidad Nacional de Tres de Febrero. Además se describen los algoritmos desarrollados en MATLAB® para realizar lo expuesto con anterioridad. Finalmente se discuten los resultados y se ofrece un Anexo con los códigos completos.
Calderón De Palma, Serattin, Alba
2. TRABAJOS PREVIOS El estudio de las ESS se remonta a trabajos de, por ejemplo, Griesinger [1], quienes establecieron que las mismas podrían ser más eficientes a la hora de analizar sistemas, que no cumplieran las condiciones LTI (Linear Time Invariant), frente a las señales más comúnmente utilizadas hasta el momento, como lo eran las MLS (Maximum Length Sequence). Sin embargo, debido al tiempo de procesamiento digital que requerían, su implementación y popularidad no se logró hasta la publicación de Farina [2] en el año 2000 y el posterior lanzamiento del software Aurora®, también del mismo autor. Desde entonces, se encuentran innumerables referencias al método de ESS. Algunos trabajos a destacar son los de Stan et. al [3], Muller y Massarini [4], Holters et. al [5] y por supuesto Farina [6]. Es recurrente la idea de comparar al método de ESS con otros, los resultados han sido ampliamente favorables para la implementación de las señales aquí tratadas.
3. CONCEPTOS TEÓRICOS E INTRODUCCIÓN A ESS 3.1. SISTEMAS LTI Y NO LINEALES
1
Procesamiento Digital de Señales
Diciembre 2014
Los sistemas denominados LTI son aquellos que cumplen con las condiciones de linealidad e invariancia en el tiempo. De la teoría de sistemas se considera que si un sistema es LTI, entonces, se lo puede caracterizar completamente mediante su respuesta al impulso ℎ(𝑡). Por lo tanto, la salida 𝑦(𝑡) del sistema frente a una señal de entrada 𝑥(𝑡) queda descripta por:
𝑦(𝑡) ∗ 𝑥 𝑖 (𝑡) = ℎ(𝑡) ∗ 𝑥(𝑡) ∗ 𝑥 𝑖 (𝑡) = ℎ(𝑡) ∗ 𝛿(𝑡) = ℎ(𝑡) (4)
+∞
𝑦(𝑡) = 𝑥(𝑡) ∗ ℎ(𝑡) = ∫
𝑥(𝜏)ℎ(𝑡 − 𝜏)𝑑𝜏
(1)
−∞
La operación en (1) muestra que la salida del sistema será la convolución entre la señal de entrada y la respuesta al impulso del sistema LTI. Los sistemas reales no suelen cumplir con las condiciones anteriores, pero se los suele aproximar o se asume que cumplen con ellas. En estas aproximaciones se basan, por ejemplo, los métodos de ruido impulsivo para la medición de parámetros acústicos de salas y auditorios. La modelización de los sistemas no lineales conlleva un desarrollo matemático más elaborado, dado que la convolución lineal presentada en (1) ya no ofrece la salida real del sistema. Para lograr la aproximación de 𝑦(𝑡), se utilizan los métodos de Volterra, Wiener y regresiones con kernel polinomiales [7]. La relación entrada-salida de un sistema invariante en el tiempo utilizando series de Volterra está dado por:
Figura 1: Esquema para la obtención de la respuesta al impulso de un sistema LTI. Es importante destacar que la elección del filtro inverso es un paso vital para el correcto procesamiento de los datos. Dado que la señal de sweep incrementa sus valores de frecuencia de un modo no constante en el tiempo, el espectro resultante en el dominio frecuencial no es plano, sino que tiene una caída de 3dB/octava. Por lo tanto para compensar esto, se propone un filtro inverso que sea la reversa temporal de la señal de entrada y con un término de ganancia, que asegure una pendiente de +3dB/octava en el dominio frecuencial. El filtro inverso luego de Novak et al. [9]:
𝑥 𝑖 (𝑡) =
∞ 𝑞
𝑦(𝑡) = 𝑘0 + ∑ 𝑘𝑞 ∗ 𝑥(𝑡)
(2)
ω ω1 ln (ω2 ) 1
2𝜋𝑇 [
(
𝑒
ω −𝑡 ln( 2 ) ω1 ) 𝑇
𝑥(−𝑡)
(6)
]
𝑞=1
Donde 𝑘𝑞 son los kernel o núcleos de Volterra de orden 𝑞. Además se impone la condición que la parte no lineal no tiene memoria (solo lo hace la parte lineal) y es puramente algebraica. En el caso en que 𝑘𝑞 = 0 ∀ 𝑞 ≠ 0 la respuesta del sistema se reduce a la de un sistema LTI.
3.2. MÉTODO ESS
La figura 2 muestra la respuesta en frecuencia de 𝑥(𝑡) y 𝑥 𝑖 (𝑡). Se puede apreciar la pendiente de -3dB/oct y +3dB/oct de la señal original y su correspondiente filtro inverso. Por otro lado la figura 3 son dos espectrogramas que corresponden a una señal de sweep, de duración 15 s. cubriendo el espectro audible de 20 Hz. a 20 kHz., y su correspondiente filtro inverso. Se puede apreciar el incremento y decremento en frecuencia de tipo exponencial para la señal y su respectivo filtro.
La definición de la señal de sweep exponencial logarítmica está dada por:
𝑥(𝑡) = sin [
𝜔1 𝑇
𝜔 ln( 2 ) 𝜔1
𝑡
(𝑒 𝑇
𝜔 ln( 2 ) 𝜔1
− 1)]
(3)
La señal 𝑥(𝑡) es una señal senoidal con una variación exponencial de su frecuencia angular, que va entre 𝜔1 y 𝜔2 , en un período de tiempo 𝑇. El método de ESS se apoya en un concepto clave del análisis de señales y sistemas, la búsqueda de un correcto filtro inverso 𝑥 𝑖 (𝑡), que garantice la obtención de una señal Delta de Dirac teórica, tras la convolución del mismo con la señal de entrada 𝑥(𝑡) en el dominio temporal. De este modo se puede recuperar la respuesta al impulso del sistema cuando realizamos la convolución entre la respuesta del sistema y el filtro inverso (Fig.1). Esto se puede demostrar de la siguiente forma:
Calderón De Palma, Serattin, Alba
Figura 2: Respuesta en frecuencia del sweep 𝑥(𝑡) (línea azul) y del filtro inverso 𝑥 𝑖 (𝑡) (línea verde).
2
Procesamiento Digital de Señales
Diciembre 2014
Figura 5: Acercamiento de la señal obtenida luego de la deconvolución de la señal de sweep con fade in/out. Figura 3: Espectrograma de una señal de sweep y su filtro inverso (rango de frecuencias de 20 Hz. a 20 kHz., duración: 15 s.)
Sin embargo la incorporación del fade in/out podría provocar otros problemas, como es el caso del Preringing.
Si se realiza la convolución entre la señal y su filtro inverso, el resultado esperado sería una Delta de Dirac perfecta. Sin embargo, al aplicar dicha operación entre las dos señales, el resultado no es el ideal. La señal es una función seno cardenal con una oscilación amortiguada en sus extremos. Una posible explicación de este fenómeno, es el hecho de que se trabaja en con un espectro acotado entre 𝜔1 y 𝜔2 . Para demostrar este fenómeno, se realizó la experiencia en MATLAB utilizando dos métodos de convolución, el primero una convolución circular utilizando zero-padding, para llevarla a una convolución lineal, y el segundo una convolución pasando por el dominio frecuencial, aprovechando el teorema de convolución de la teoría de sistemas. En la figura 4 podemos ver, en un acercamiento, los resultados de ambas operaciones y una pequeña oscilación que se extiende hacia los extremos.
4. VENTAJAS Y POSIBLES PROBLEMAS DE LAS ESS
Figura 4: Acercamiento de la señal obtenida luego de la deconvolución realizada con convolución circular (arriba) y aplicando el teorema de convolución (abajo). Al incluir un fade-in y fade-out en la señal de sweep, se observó como la respuesta temporal mejora en sus extremos y se disminuye la oscilación. La figura 5 muestra la mejora en el comportamiento de la señal temporal.
Calderón De Palma, Serattin, Alba
4.1. APLICACIÓN Y VENTAJAS Como se mencionó, las ESS son utilizadas en el análisis de sistemas que presentan alinealidades. La razón de esto es que cuando se recupera la respuesta al impulso, se separan los componentes de distorsión armónica de la respuesta puramente lineal, lo cual permite, analizar cada respuesta individualmente. Esto se logra también con las señales MLS, sin embargo, la gran ventaja de las señales de sweep es que colocan a las componentes de distorsión en forma de “picos” por detrás de la respuesta lineal a una distancia específica:
∆𝑡 = 𝑇
ln(𝑁) 𝜔 ln (𝜔2 ) 1
(7)
Donde ∆𝑡 representa el retardo entre la respuesta al impulso lineal y la componente de distorsión de orden 𝑁. Como consecuencia de la ubicación de las componentes de distorsión, la relación señal ruido, que representa la relación en dB entre el valor medio de la señal y las componentes de ruido y distorsión en la cola de la respuesta lineal luego de la deconvolución, es excelente, pudiendo superar en 20 dB a aquella de una señal MLS [3]. La diferencia se da principalmente debido a que al hacer la medición con señales MLS, luego de la deconvolución, las componentes de distorsión armónica aparecen en la cola de la respuesta lineal, mientras que en las ESS son separadas de la misma.
4.2. PROBLEMAS Y SOLUCIÓNES Existen algunos problemas a la hora de trabajar con las ESS. Farina propuso soluciones a los siguientes problemas [6]:
3
Procesamiento Digital de Señales
Diciembre 2014
Pre-ringing Ruidos impulsivos durante la medición Desajustes del clock digital Promediado de múltiples mediciones Ecualización generada (no deseada) por equipos de ensayo En el presente trabajo solo se discutirán dos de estos problemas particulares que podrían conllevar una medición errónea.
4.2.1. RUIDOS IMPULSIVOS DURANTE LA MEDICIÓN Los ruidos impulsivos que se producen durante la grabación de la señal de salida del sistema en estudio, contaminan la medición. Estos efectos no deseados pueden verse en un espectrograma como una línea vertical (Figura 6). Figura 7: IR filtrado por banda de octava (a 1KHz) contaminado con ruido impulsivo (arriba), y sin contaminación (abajo).
Figura 6: Medición de ESS contaminada por un ruido impulsivo Al convolucionar la salida contaminada con el filtro inverso, las impurezas se verán como un barrido en frecuencia descendente el cual se puede apreciar en la figura 7.
Otros parámetros como la Claridad, también serán afectados ante estos ruidos impulsivos espurios, por lo que los resultados también serán erróneos, aunque en menor medida. Existen variadas formas de eliminar este problema. Si bien lo primero que se podría pensar es en eliminar el defecto sonoro de la señal grabada, esto no resuelve el problema. Farina entonces, propone determinar a qué frecuencia del barrido senusoidal se produce el ruido impulsivo y aplicar a la grabación un filtro de banda estrecha en esa frecuencia exacta. De esta manera se logrará restaurar la forma de onda. En la figura 8 se muestra la reducción de la amplitud del efecto espuria y como la cola reverberante es prácticamente libre de distorsión.
Figura 7: Respuesta al impulso y Efecto no deseado La parte anterior a la respuesta al impulso lineal es irrelevante debido a que está dentro de la porción que se eliminará para obtener la respuesta al impulso lineal. El problema radica en la parte posterior de este barrido espuria. Uno de los problemas más importantes es que al modificar la última parte de la medición, se producen errores en el cálculo del Tiempo de Reverberación calculado según la norma ISO 3382. El error consiste en una sobreestimación del tiempo de reverberación. En la figura 7 se muestra esta sobreestimación del T30 debido a la impureza causada por el ruido impulsivo.
Calderón De Palma, Serattin, Alba
Figura 8: Resultado del filtro FFT (arriba), y la IR filtrada por banda de octava sin el evento impulsivo (abajo)
4.2.2. PROMEDIADO DE MÚLTIPLES MEDICIONES
4
Procesamiento Digital de Señales
En métodos anteriores como el MLS, el promediado de múltiples mediciones permitía mejorar la relación señal ruidos (S/R). Sin embargo, el promedio sólo es válido si se tiene un sistema LTI, lo cual es imposible de lograr en las mediciones acústicas de recintos. En el caso del ESS, no es necesario acudir a estos promedios para mejorar la relación S/R, solo basta con realizar una única medición de un barrido lo suficientemente largo (norma ISO 18233/2006). A pesar de ello, en ocasiones no es posible usar un barrido largo y se debe acudir a la utilización de un sweep corto. En estos casos se hace indispensable realizar el promedio de varias mediciones. El problema que se presentará al promediar múltiples respuestas al impulso obtenidas mediante el ESS, será la cancelación o subestimación de las frecuencias altas en la última parte de la cola reverberante. Esto se puede apreciar en la figura 9 donde se presenta un ejemplo de éste problema. Por arriba de los 350Hz la IR obtenida por promedio está subestimada con respecto a la obtenida a partir de un sweep largo. En este caso, entre las frecuencias de 5KHz a 6KHZ esa subestimación llega aproximadamente a los 10 dB.
Diciembre 2014
Junto con la subestimación de la energía en la caída, también se puede apreciar un empeoramiento en la relación S/R. Para compensar este problema, Farina utiliza un archivo de música que contenga en el canal izquierdo el sine sweep, y en el derecho la señal grabada del recinto. Luego esta señal estéreo la procesa con un plugin de Aurora llamado “Cross Functions” mediante el cual, se calcula la función de transferencia (𝐻1 ) del sistema en el dominio de las frecuencias. La función de transferencia está dada por la división entre el promedio cruzado (𝐺𝐿𝑅 ) y el promedio autoespectro (𝐺𝐿𝐿 ) en el dominio espectral (ecuación 8)
𝐻1 (𝑓) =
𝐺𝐿𝑅 𝐺𝐿𝐿
(8)
Luego la 𝐻1 se deconvoluciona y se obtiene la respuesta al impulso en el dominio temporal. En el ejemplo anterior, la respuesta al impulso recuperada se muestra en la figura 11.
Figura 11: IR filtrada por banda de octavas de 50 barridos de 1 seg. Utilizando “Cross Functions”
Figura 9: Espectro del barrido único de 50 seg. (azul), y promedio de 50 barridos de 1 seg. (magenta) Como consecuencia, la pendiente de la curva de caída se verá afectada y por ende el tiempo de reverberación calculado (figura 10).
Se puede apreciar que al procesar los datos en el dominio de la frecuencia, se obtiene una IR con una caída apenas menor a la obtenida con un solo sweep largo, y una mejora en la relación S/R. A pesar de ello, sin dudas, es mejor la utilización de un solo barrido largo.
5. DESARROLLO EXPERIMENTAL Se estudió la aplicación de las ESS realizando dos mediciones en la Universidad Nacional de Tres de Febrero. Se midió la respuesta en magnitud y fase de un filtro de una consola de audio analógica y se buscaron los parámetros acústicos de un aula de la Universidad. Todos los cálculos y procesamiento de las señales se llevaron a cabo en MATLAB®. Los códigos se detallan en el ANEXO 1.
5.1. MEDICIÓN DEL FILTRO ANALÓGICO
Figura 10: IR filtrada por banda de octava de un barrido único de 50 seg. (arriba), y promedio de 50 barridos de 1 seg. (abajo)
Calderón De Palma, Serattin, Alba
Para obtener la respuesta al impulso se generó una señal de sweep con las características de la tabla 1 y su correspondiente filtro inverso. Luego se procesó la señal de sweep a través de uno de los filtros de una consola de audio Mackie Onyx 1640.
5
Procesamiento Digital de Señales
Diciembre 2014
Características del sweep Frecuencia inferior 30 Hz Frecuencia superior 18 kHz Duración 10 s. Silencio 5 s. 0.1 s. Fade in/out (Ventana) (Hann) Frecuencia de muestreo 44100 Hz
Finalmente se graficaron los resultados, los cuales pueden apreciarse en las figuras 14 y 15.
Tabla 1: Datos de la señal de sweep utilizada en las mediciones Una vez obtenida la señal procesada se realizó la convolución entre dicha señal y el filtro inverso, a fin de recuperar la respuesta lineal del sistema. La figura 12 muestra el resultado obtenido, se puede ver, en el acercamiento, que la respuesta al impulso no es una delta de Dirac sino que posee perturbaciones.
Figura 12: Respuesta al impulso del filtro recuperada luego de la deconvolución Esta respuesta temporal se transformó al dominio frecuencial utilizando una FFT de N puntos, donde N es la longitud de la respuesta al impulso. Luego se obtuvo la magnitud de la respuesta en decibeles, utilizando:
|𝐻(𝜔)| = 10 log (√𝑅𝑒 2 (𝐻(𝜔)) + 𝐼𝑚2 (𝐻(𝜔)))
(9)
Además se recuperó la respuesta en fase utilizando la función de MATLAB Angle.m. La figura 13 muestra un diagrama básico del procesamiento.
Figura 14: Magnitud de la respuesta en frecuencia en dB del filtro (arriba) y respuesta de fase en radianes (abajo)
Figura 15: Acercamiento de la magnitud de la respuesta en frecuencia y fase, utilizando una escala lineal. Como se puede ver, los resultados demuestran que el filtro es un pasa banda cuya frecuencia central es 135 Hz. y su fase es de tipo lineal, las características se resumen en la tabla 2.
Características del filtro Frecuencia Central Frecuencia de Corte superior Frecuencia de Corte inferior Ancho de Banda Tipo de fase
135 Hz 185 Hz 100 Hz 85 Hz Lineal
Tabla 2: Características del filtro medido
5.2. MEDICIÓN DE PARÁMETROS ACÚSTICOS DE UN AULA
Figura 13: Diagrama de procesamiento para obtener la magnitud y fase del filtro a partir de la respuesta al impulso
Calderón De Palma, Serattin, Alba
Para medir los parámetros acústicos del recinto, primero se reprodujo la misma señal de sweep utilizada para el caso del filtro y se grabó la respuesta de la sala. La medición se llevó a cabo utilizando los siguientes equipos: Monitores Genelec 1037c Micrófono de medición Earthworks M50 Dodecaedro y Subwoofer Outline
6
Procesamiento Digital de Señales
Placa de audio Motu Traveler Mk3 Mixer Mackie Onyx 1640
Diciembre 2014
Es válido afirmar que la claridad es una propiedad que tiene una relación directa con el tiempo de reverberación. Se define el parámetro 𝐶50 en la siguiente ecuación [9]: 50𝑚𝑠
Figura 16 – Diagrama en bloque del proceso de obtención de la señal
ℎ2 (𝑡)𝑑𝑡 ∫ 𝐶50 = 10 log ( 0𝑚𝑠 ) ∞ ∫50𝑚𝑠 ℎ2 (𝑡) 𝑑𝑡
(10)
Mientras que el parámetro 𝐶80 se define como [9]: Luego se realizó la deconvolución de la señal medida de sala y se recuperó la respuesta al impulso. En la figura 17 se aprecia la respuesta al impulso lineal a la derecha y las componentes de distorsión armónica a la izquierda. Para el análisis se utilizó solo la componente lineal.
Figura 17: Respuesta al impulso de la sala y componentes de distorsión armónica (Amplitud normalizada) Una vez recuperada la respuesta puramente lineal, se la proceso para obtener los siguientes parámetros: Tiempo de Reverberación (TR20, TR30) Claridad (C50, C80) Initial Time Delay Gap (ITDG)
80𝑚𝑠
∫0𝑚𝑠 ℎ2 (𝑡)𝑑𝑡 𝐶80 = 10 log ( ∞ ) ∫80𝑚𝑠 ℎ2 (𝑡) 𝑑𝑡
(11)
5.2.3. Initial Time Delay Gap (ITDG) Se define al Initial Time Delay Gap como el tiempo de arribo entre el sonido directo en un punto determinado de la sala, y el arribo de la primer reflexión [10]. Es un parámetro que depende de la posición de la sala en donde se realiza la medición del mismo, ya que queda determinada por las reflexiones que inciden en ese punto. También se ve afectado por la posición de la fuente emisora. Valores pequeños de ITDG son deseables para una sala de dimensiones pequeñas. Por ejemplo, valores por debajo de los 15 ms obtenidos en el centro de la sala son aceptables.
5.2.4. Procesamiento y resultados A modo esquemático se muestra un diagrama de procesamiento de la señal ℎ[𝑛] en la figura 18.
5.2.1. TIEMPO DE REVERBERACIÓN Factores como la inteligibilidad de la palabra pueden a llegar verse perjudicadas por tiempos de reverberación elevados dentro de la sala. Sabiendo que el tiempo de reverberación juega un papel importante en las características acústicas de una sala, definimos al tiempo de reverberación, como el tiempo que transcurre en el que la intensidad del sonido decae 60 dB del nivel original. Este parámetro se define como 𝑇𝑅60 [8]. Análogamente, podemos definir a los parámetros 𝑇𝑅20 y 𝑇𝑅30 como el tiempo en el cual la intensidad de la señal decae 20 dB o 30 dB respectivamente.
5.2.2. Claridad (C50, C80) El parámetro claridad, se define como la relación de energía de la onda sonora (en este caso, de la energía de la respuesta al impulso de la sala) en un tiempo determinado, y luego de dicho tiempo. Representa una comparación entre la energía temprana y tardía.
Calderón De Palma, Serattin, Alba
Figura 18: Diagrama del procesamiento de la respuesta al impulso lineal del sistema La respuesta lineal se procesó por dos ramas distintas. La primera (rama superior en figura 18) filtra la señal con un filtro de banda de octava, utilizando la función FiltoctDSP.m. Luego, en base a esta señal filtrada se calculó el C50 y C80 por octava. La función a cargo de esta tarea fue CxDSP.m. Ésta función utiliza el método numérico de trapecios para calcular las siguientes integrales, representativas de los valores de C50 y C80. Con la misma señal filtrada por bandas de octava se calculó la integral de Schroeder con el fin de obtener
7
Procesamiento Digital de Señales
Diciembre 2014
una curva suave para el cálculo del TR. Por definición, la integral de Schroeder se calcula de la siguiente manera: ∞
𝑠(𝑡) = ∫ ℎ2 (𝑡)𝑑𝑡
(12)
𝑡
Donde ℎ(𝑡) es la respuesta al impulso y la integral se realiza en “reversa”. La función SchrointDSP.m obtiene la imagen de la integral en (12), utilizando la regla de trapecios y expresa el resultado normalizado en dB. Esta curva normalizada resulta tener una pendiente muy suave (principalmente en frecuencias medias y altas), lo cual permite realizar fácilmente una interpolación lineal que permita calcular los valores de TR. La función TRDSP.m busca los valores que correspondan a una caída de 20 dB y 30 dB, descartando los primeros 5 dB, pero sin interpolación lineal, dado que se consideró que la respuesta poseía una pendiente aproximadamente lineal luego de realizar la integral en (12). Esta última afirmación es válida para el rango de frecuencias de interés y una buena SNR. Dado que hablamos de un aula, se consideraron de mayor interés el rango de frecuencias correspondientes al habla. Teniendo en cuenta esto y sabiendo que la inteligibilidad de la palabra se corresponde principalmente con la comprensión de las consonantes, cuyas componentes frecuenciales son altas (entre 1 kHz y 4 kHz), la aproximación tiene su validez. El cálculo del TR20 y TR30 se deriva de las siguientes formulas:
𝑇𝑅30 = 2 (𝑡(−35 𝑑𝐵) − 𝑡(−5 𝑑𝐵))
(13)
𝑇𝑅20 = 3 (𝑡(−25 𝑑𝐵) − 𝑡(−5 𝑑𝐵))
(14)
Donde 𝑡(𝑖 𝑑𝐵) representa el instante en el que la función 𝑠(𝑡), normalizada y expresada en dB, registra la caída 𝑖 en decibeles. Finalmente, se evaluó el ITDG. En general el análisis de este parámetro se realiza de modo visual evaluando el reflectograma o la envolvente de la señal. Para éste trabajo se optó por recuperar la envolvente mediante la obtención de una señal analítica utilizando la transformada de Hilbert [11], de este modo:
𝐸𝑛𝑣𝑜𝑙𝑣𝑒𝑛𝑡𝑒 = √ℎ2 (𝑡) + (ℋ(ℎ(𝑡))2
(15)
Donde ℋ(ℎ(𝑡)) representa la transformada de Hilbert de la señal ℎ(𝑡). Luego, en base a la envolvente, la función ITDG.m calcula el tiempo que hay entre el sonido directo y la mayor reflexión que aparece en los primeros 30 ms, realizando una evaluación de la envolvente por pasos determinados por el usuario. El criterio de “mayor reflexión” es algo relativo, por lo que se le ofrece la posibilidad al usuario de confirmar la reflexión seleccionada o elegir una a mano. La figura 19 muestra la pantalla de selección de reflexiones para el cálculo de ITDG.
Calderón De Palma, Serattin, Alba
Figura 19: Pantalla de selección de reflexiones para calcular el ITDG, la flecha verde indica la reflexión seleccionada por ITDG.m Por último la tabla 3 muestra los resultados obtenidos.
Frecuencia TR20 TR30 Normalizada [s.] [s.] [Hz] 30 1,41 1,38 63 1,60 1,43 125 0,98 1,03 250 0,61 0,77 500 0,58 0,58 1000 0,46 0,47 2000 0,59 0,63 4000 0,69 0,74 8000 0,59 0,64
C80 [dB]
C50 [dB]
ITDG [ms.]
-6,30 -11,29 -1,34 -2,94 4,78 -0,72 8,38 4,68 9,81 7,17 1,56 12,09 8,27 10,00 6,51 9,94 6,44 10,73 7,51
Tabla 3: Resultados TR20, TR30, C80 y C50 por banda de octava Se puede observar que los valores de TR son bajos en frecuencias medias y altas, lo cual garantiza niveles altos de Claridad, en particular se le dio importancia a los valores de C50, el cual es un parámetro significativo en los espacios destinados al habla. Son valores esperados debido a la finalidad del recinto estudiado. Por otro lado el valor de ITDG es pequeño. Éste también era un valor esperado, principalmente por las dimensiones pequeñas y forma paralelepipédica de la sala. Es importante destacar que los valores de la banda de 16 kHz fueron desestimados y también podrían ser aquellos en la banda de 30 Hz, debido a que el recorrido en frecuencia del sweep no fue lo suficientemente grande como para caracterizar completamente a dichas bandas. Nuevamente, teniendo en cuenta esta medición en particular, no es necesario cubrir todo el espectro audible de 20 Hz a 20 kHz.
8
Procesamiento Digital de Señales
6. CONCLUSIONES En este trabajo se presentó el método de ESS y su implementación en la caracterización de sistemas LTI y no lineales. No solo se demostraron conceptos teóricos de las ESS, sino que además se realizaron dos mediciones para demostrar la validez y eficacia del método de medición de la respuesta al impulso. Aún quedan por discutir el resto de los problemas que se podrían encontrar y sus soluciones. Además se han hecho algunas simplificaciones en los códigos desarrollados para las mediciones, que se aplicaron a los casos estudiados. Sería ideal extender el alcance de los algoritmos para su posterior aplicación a casos más generales. Finalmente se puede concluir que el método expuesto en este trabajo representa grandes ventajas frente a otros, debido a su fácil implementación, velocidad de procesamiento y superioridad en la relación señal ruido.
7. REFERENCIAS [1] D. Griesinger. “Beyond MLS – Occupied Hall Measurements With FFT Techniques”. 101st AES Convention. 1996 [2] A. Farina. “Simultaneous Measurements of Impulse Response and Distortion with a Sine Sweep Technique”. 108th AES Convention. 2000 [3] G. Stan, J. Embrechts, D. Archambeau. “Comparison of different impulse response measurement techniques”. JAES Vol. 50, No. 4. 2002
Calderón De Palma, Serattin, Alba
Diciembre 2014
[4] S. Müller, P. Massarani. “Transfer-Function Measurement with Sweeps”, JAES Vol. 49, No 6. 2001 [5] M. Holters, T. Corbach, U. Zölzer. “Impulse Response Measurement Techniques and their applicability in the Real World”. 12th Int. Conference on Digital Audio Effects (DAFx-09). 2009 [6] A. Farina. “Advancements in impulse response measurements by sine sweeps”. 122nd AES Convention. 2007 [7] M. Franz, B. Schölkopf “A unifying view of Wiener and Volterra Theory and Polynomial Kernel Regression”. Neural Computation, December 2006, Vol. 18, No. 12 [8] T.D. Rossing. “Springer handbook of acoustics”. Springer Science+Business Media. 2007. ISBN 978-0387-30446-5. [9] ISO 3382. “Acoustics-Measurement of the reverberation time of rooms with reference to other acoustical parameters”. Second edition 1997-06-15. [10] L. Beranek. “Concert Halls and Opera Houses”. Springer Science & Business Media, 2004. ISBN 0387955240. [11] G. Lindgren. “Stationary Stochastic Processes: Theory and Application”. CRC Press. 2013. ISBN 9781-4665-5779-6.
9
Procesamiento Digital de Señales
Diciembre 2014
ANEXO 1: CÓDIGOS COMPLETOS DESARROLLADOS EN MATLAB 1. Función para la creación del sweep logarítmico sin fade in/out function [sweep,sweep_inverso] = logsinesweep(f1,f2,fs,T,silencio) % % % % % %
Los valores de entrada a la función son: f1 = Frecuencia inicial f2 = Frecuencia final fs = Frecuencia de muestreo de la señal T = Duración de la señal (sin incluir el silencio) Silencio = Silencio agregado al final de la señal
w1 = 2*pi*f1; %Freq angular inicial w2 = 2*pi*f2; %Freq angular final t = 0:1/fs:T; %Vector tiempo en s. L = T/log(w2/w1); %Tasa de incremento de frecuencia exponencial % Sweep Logaritmico (ESS - Exponential Sine Sweep) sweep = sin((w1*L)*((exp(t/L))-1)); %Sweep exponencial % El filtro inverso se realizó tomando en cuenta las consideraciones de % Antonín Novák, Laurent Simon, František Kadlec y Pierrick Lotton en % "Nonlinear System Identification Using Exponential Swept-Sine Signal" f1 = (w1/(2*pi*L))*exp(-t/L).*fliplr(sweep); %Filtro inverso sweep_inverso = f1./max(f1); %Filtro inverso normalizado % Agregar silencio al final sweep = [sweep zeros(1,silencio*fs)]; sweep_inverso = [sweep_inverso zeros(1,silencio*fs)]; % Escribir .wav wavwrite(sweep,fs,'ESS (Exponential Sine Sweep)'); wavwrite(sweep_inverso,fs,'IESS (Exponential Sine Sweep Inverso)'); end
2. Función para la creación de sweep logarítmico con fade in/out function [sweep,sweep_inverso] = logsweepfade(f1,f2,fs,T,silencio,Tfadein,Tfadeout) %Función para generar sweeps logarítmicos, comienza en una frecuencia f1 y %termina en una frecuencia f2. Luego se le agrega un silencio al final de %la señal cuya duración se especifica en segundos. Los controles de fade in %y fade out también se especifican en segundos. w1 = 2*pi*f1; %Freq angular inicial w2 = 2*pi*f2; %Freq angular final t = 0:1/fs:T; %Vector tiempo en s. L = T/log(w2/w1); %Tasa de incremento de frecuencia exponencial % Sweep exponencial sweep = sin((w1*L)*((exp(t/L))-1)); % Generar ventana para fade in/out Win1 = hann(fs*2*Tfadein,'symmetric'); Win2 = hann(fs*2*Tfadeout,'symmetric'); Win1 = Win1(1:1:end/2);
Calderón De Palma, Serattin, Alba
10
Procesamiento Digital de Señales
Diciembre 2014
Win2 = Win2(end/2:1:end); %Ventana para fade in+fade out WinFin = [Win1' ones(1,length(sweep)-length(Win1)-length(Win2)) Win2']; % Señal de sweep con fade in y fade out sweep = sweep.*WinFin; % El filtro inverso se realizó tomando en cuenta las consideraciones de % Antonín Novák, Laurent Simon, František Kadlec y Pierrick Lotton en % "Nonlinear System Identification Using Exponential Swept-Sine Signal" f1 = (w1/(2*pi*L))*exp(-t/L).*fliplr(sweep); %Filtro inverso sweep_inverso = f1./max(f1); %Filtro inverso normalizado % Agregar Silencio sweep = [sweep zeros(1,silencio*fs)]; sweep_inverso = [sweep_inverso zeros(1,silencio*fs)]; % Escribir .wav wavwrite(sweep,fs,'ESS (Exponential Sine Sweep con Fade IN-OUT)'); wavwrite(sweep_inverso,fs,'IESS (Exponential Sine Sweep Inverso con Fade IN-OUT)'); end
3. Verificación de señal de sweep sin fade in/out %% f1 f2 fs
I) Generar señales de sweep y filtro inverso = 20; %frec.inicial = 20e3; % frec. final = 44100; %Solo se tomo por ser un estándar, además ofrece la posibilidad %de exportar luego a Adobe Audition usando wavwrite T = 15; %Duración del sweep sil = 5; %Duración del silencio [s,invs] = logsinesweep(f1,f2,fs,T,sil); H1 = fft(s)/length(s); Hdb1 = 20*log10(abs(H1)); H2 = fft(invs)/length(invs); Hdb2 = 20*log10(abs(H2)); freq = (0:length(H1)-1)*fs/length(H1); a = round(f1*length(H1)/fs); %Indexa hasta f1 Hz aprox. b = round(f2*length(H1)/fs); %Indexa hasta f2 Hz aprox. figure(1) %Grafica de la respuesta de la señal de sweep y su filtro inverso semilogx(freq(a:b),Hdb1(a:b),'linewidth',2);hold on; semilogx(freq(a:b),Hdb2(a:b),'g','linewidth',2); xlabel('Frecuencia (Hz.)');ylabel('Amplitud (dB)'); legend('Señal de Sweep','Filtro Inverso'); grid on %% II) Espectrograma figure(2) subplot(2,1,1) [y,f,t,p]=spectrogram(s,(hann(2^12))',[],2^8,fs); surf(t,f,10*log10(abs(p)),'EdgeColor','none'); axis xy; axis([0 T 0 f2]); colormap(jet); view(0,90); xlabel('Tiempo (s.)'); ylabel('Frecuencia (Hz)'); title('ESS'); subplot(2,1,2); [y1,f1,t1,p1]=spectrogram(invs,(hann(2^12))',[],2^8,fs); surf(t1,f1,10*log10(abs(p1)),'EdgeColor','none'); axis xy; axis([0 T 0 f2]); colormap(jet); view(0,90);
Calderón De Palma, Serattin, Alba
11
Procesamiento Digital de Señales
Diciembre 2014
xlabel('Tiempo (s.)'); ylabel('Frecuencia (Hz)'); title('Filtro inverso'); axis([0 T 0 f2]);
%% III) Realizar convolución del filtro inverso y el sweep %Se usó la convolución circular por ser más rápida, primero se realiza zero %padding. También se realizo una convolución pasando por frecuencia y se %midió el tiempo de procesado de cada método. %CONVOLUCIÓN TEMPORAL tic sconv = [s zeros(1,length(s)+length(s)-1)]; % zero padding para que % cconv=conv invsconv = [invs zeros(1,length(invs)+length(invs)-1)]; h1 = cconv(sconv,invsconv); h1 = h1(1:1:length(s)); tconvl = toc; % ¿Cuanto tarda en realizar la convolución? %CONVOLUCIÓN PASANDO POR FRECUENCIA tic S = fft(s); INVS = fft(invs); H = S.*INVS; h2 = ifft(H); tconvf = toc; % ¿Cuanto tarda en realizar la convolución? %% IV) Visualizar deconvolución t = 0:1/fs:T+sil; %Vect. Tiempo mas silencio figure(3) subplot(2,1,1) plot(t,h1); xlabel('tiempo (s.)');ylabel('h1');axis([T-.01 T+.01... min(h1)+10 max(h1)+10]) subplot(2,1,2); plot(t,h2,'r'); xlabel('tiempo (s.)');ylabel('h2');axis([T-.01 T+.01... min(h2)+10 max(h2)+10]) %% V) Escribir .wav de la deconvolución para posible análisis en software comercial wavwrite(h1,fs,16,'RespImpt');wavwrite(h2,fs,16,'RespImpf');
4. Verificación de señal de sweep con fade in/out %% f1 f2 fs
I) Generar señales de sweep y filtro inverso = 20; %frec. inicial = 20e3; %frec. final = 44100; %Solo se tomo por ser un estándar, además ofrece la posibilidad %de exportar luego a Adobe Audition usando wavwrite T = 15; %Duracion sweep sil = 5; %Duración silencio fadein = 0.2; %Duración fade in fadeout = 0.2; % Duración fade out [s,invs] = logsweepfade(f1,f2,fs,T,sil,fadein,fadeout); %Visualizar magnitud de filtro inverso y señal de sweep H1 = fft(s)/length(s); Hdb1 = 20*log10(abs(H1)); H2 = fft(invs)/length(invs); Hdb2 = 20*log10(abs(H2)); freq = (0:length(H1)-1)*fs/length(H1); a = round(f1*length(H1)/fs); %Indexa hasta f1 Hz aprox. b = round(f2*length(H1)/fs); %Indexa hasta f2 Hz aprox.
Calderón De Palma, Serattin, Alba
12
Procesamiento Digital de Señales
Diciembre 2014
figure(1) %Grafica de la respuesta de la señal de sweep y su filtro inverso semilogx(freq(a:b),Hdb1(a:b),'linewidth',2);hold on; semilogx(freq(a:b),Hdb2(a:b),'g','linewidth',2); xlabel('Frecuencia (Hz.)');ylabel('Amplitud (dB)'); legend('Señal de Sweep','Filtro Inverso'); grid on %% II) Espectrograma figure(2) subplot(2,1,1) [y,f,t,p]=spectrogram(s,(hann(2^12))',[],2^8,fs); surf(t,f,10*log10(abs(p)),'EdgeColor','none'); axis xy; axis([0 T 0 f2]); colormap(jet); view(0,90); xlabel('Tiempo (s.)'); ylabel('Frecuencia (Hz)'); title('ESS'); subplot(2,1,2); [y1,f1,t1,p1]=spectrogram(invs,(hann(2^12))',[],2^8,fs); surf(t1,f1,10*log10(abs(p1)),'EdgeColor','none'); axis xy; axis([0 T 0 f2]); colormap(jet); view(0,90); xlabel('Tiempo (s.)'); ylabel('Frecuencia (Hz)'); title('Filtro inverso'); %% III) Realizar convolución del filtro inverso y el sweep %Se usó la convolución circular por ser más rápida, primero se realiza zero %padding. También se realizo una convolución pasando por frecuencia y se %midió el tiempo de procesado de cada método. %CONVOLUCIÓN TEMPORAL tic sconv = [s zeros(1,length(s)+length(s)-1)]; % zero padding para que % cconv=conv invsconv = [invs zeros(1,length(invs)+length(invs)-1)]; h1 = cconv(sconv,invsconv); h1 = h1(1:1:length(s)); tconvl = toc; %Duracion proceso de convolucion metodo 1 %CONVOLUCIÓN PASANDO POR FRECUENCIA tic S = fft(s); INVS = fft(invs); H = S.*INVS; h2 = ifft(H); tconvf = toc; %Duración proceso convolucion metodo 2 %% IV) Visualizar deconvolución t = 0:1/fs:T+sil; %Vect. Tiempo mas silencio figure(3) subplot(2,1,1) plot(t,h1); xlabel('tiempo (s.)');ylabel('h1');axis([T-.01 T+.01... min(h1)+10 max(h1)+10]) subplot(2,1,2); plot(t,h2,'r'); xlabel('tiempo (s.)');ylabel('h2');axis([T-.01 T+.01... min(h2)+10 max(h2)+10]) %% V) Escribir .wav de la deconvolución para posible análisis en software comercial wavwrite(h1,fs,16,'RespImptfade');wavwrite(h2,fs,16,'RespImpffade');
5. Análisis del filtro analógico %% I) Señal de Consola [s,fs] = wavread('consola_senal_grabada'); [invs,~] = wavread('filtro_inverso1'); s = [s' zeros(1,length(s)+length(invs)-1)];
Calderón De Palma, Serattin, Alba
13
Procesamiento Digital de Señales
Diciembre 2014
invs = [invs' zeros(1,length(s)- length(invs))... zeros(1,length(s)+length(invs)-1)]; h = cconv(s,invs); h = h(9*fs:1:11*fs)./max(h); %Seleccionamos la parte significativa de h t = 0:1/fs:(length(h)-1)/fs; %Vector tiempo %Gráfico de señal + acercamiento figure(1) subplot(1,2,1) plot(t,h,'k'); xlabel('Tiempo (s.)'); ylabel('Amplitud relativa normalizada') subplot(1,2,2) plot(t,h,'k'); xlabel('Tiempo (s.)'); axis([0.96 1.06 -0.14 0.14]) ylabel('Amplitud relativa normalizada') H = fft(h)/length(h); Hdb = 20*log10(abs(H)); %Magnitud logarítmica phase = angle(H); %Calcular fase freq = (0:length(H)-1)*fs/length(H); a = round(30*length(H)/fs); %Indexa hasta 30 Hz aprox. b = round(18e3*length(H)/fs); %Indexa hasta 18 kHz aprox. %Gráficos de Magnitud y fase figure(2) subplot(2,1,1); semilogx(freq(a:b),Hdb(a:b)); xlabel('Frecuencia (Hz.)'); ylabel('Magnitud (dB)');grid on subplot(2,1,2) semilogx(freq(a:b),phase(a:b)); xlabel('Frecuencia (Hz.)'); ylabel('Fase (rad.)')
6. Análisis de parámetros de sala %% I) Parámetros de sala clear all clc clf %Cargar datos de mediciones [s,fs] = wavread('aula_senal_grabada'); [invs,~] = wavread(filtro_inverso1'); %Recuperar respuesta al impulso s = [s' zeros(1,length(s)+length(invs)-1)]; invs = [invs' zeros(1,length(s)- length(invs))... zeros(1,length(s)+length(invs)-1)]; h = cconv(s,invs); h = h./max(abs(h)); %Convolución+normalización %En un programa más genérico se le puede dar la opción al usuario de %seleccionar que parte corresponde al sonido directo y componentes %armónicas: % t = 0:1/fs:(length(h)-1)/fs; % figure(1) % plot(t,h); % [X,~] = getpts(figure(1)); %Seleccionar los limites donde corto la señal % h1 = h(round(X(1)*fs):round(X(2)*fs)); % t1 = X(1):1/fs:X(2); % plot(t1,h1); % %Si además solo quiero el sonido directo % [X1,~] = getpts(figure(1)); % h2 = h1(round(X1(1)*fs):round(X1(2)*fs)); % t2 = X1(1):1/fs:X1(2); % figure(2) % plot(t2,h2);
Calderón De Palma, Serattin, Alba
14
Procesamiento Digital de Señales
Diciembre 2014
%Para este trabajo se realizó el análisis de esta señal en particular por %inspección h1 = h(8.5*fs:1:11*fs); % Respuesta al impulso lineal + comp. armonicas t = 0:1/fs:(length(h1)-1)/fs; %vector tiempo %Grafico 1, respuesta al impulso recuperada figure(1) plot(t,h1); xlabel('Tiempo (s.)'); grid on; ylabel('Amplitud relativa normalizada'); % Filtrar señal con un filtro de banda de octava que vaya de 30Hz a 18kHz % Primero hay que dejar la respuesta al impulso lineal (quitar componentes % de distorsión armonica) h2 = h(10*fs:1:11*fs); [~,nmax] = max(abs(h2)); h2 = h2(nmax:1:end); %Asumo que el maximo corresponde al sonido directo [y,f] = filtthirdoctDSP(h2,fs,1); % Realizar integral de Schroeder según ISO 3382 [~,N] = size(y); for i=1:N x = SchrointDSP(y(:,i)); z{1,i} = x; end % Calculo del TR, no se uso un polinomio de ajuste, solo la integral de % Schroeder for i=1:N [TR20{1,i},TR30{1,i}] = TRDSP(z{i},fs); end TR20 = cell2mat(TR20);TR30 = cell2mat(TR30); figure(2) semilogx(f,TR20,'-om',f,TR30,'-sg');xlabel('Frecuencia (Hz.)'); grid on ylabel('TR (s.)');legend('TR20','TR30'); axis([30 max(f) 0 max(max(TR20,TR30))+0.5]); %Calculo de ITDG %Recuperar envolvente y normalizar env = sqrt(hilbert(h2).^2); env_norm = abs(env./max(env)); [itdg,idx,~,~] = ITDG(env_norm,fs,0.005); %Método automatico t2 = 0:1/fs:(length(env_norm)-1)/fs; %vector tiempo gráfico envolvente figure(3) set(figure(3),'position',[1 1 1280 800]) %Que aparezca maximizada plot(t2,env_norm); axis([-0.01 0.2 0 1.05]); hold on xlabel('Tiempo (s.)');ylabel('Envolvente (Amplitud relativa normalizada)'); [pks,idx1] = findpeaks(env_norm,'MINPEAKHEIGHT',0.15,'NPEAKS',10,... 'MINPEAKDISTANCE',round(0.0005*fs)); %Buscar picos plot(t2(idx1),pks+0.008,'k+'); %Graficar picos %Graficar pico seleccionado por la funcion ITDGDSP plot(t2(idx),env_norm(idx)+0.03,'kv','markerfacecolor',[0 1 0]); [ITDG,~] = getpts(figure(3)); %Selección de reflexión por parte de usuario %Si está conforme con la elección del programa presiona "ENTER" y la salida %es el valor que se calculo con la función ITDGDSP, sino puede seleccionar %la reflexión que considera más adecuada (criterio personal) y presionar
Calderón De Palma, Serattin, Alba
15
Procesamiento Digital de Señales
Diciembre 2014
%"ENTER". En el caso que se haya equivocado con la selección, el usuario %puede presionar "Backspace" y elegir otra vez. %El valor seleccionado por ITDGDSP se indica con una flecha VERDE, los %picos relativos superiores a 0.15 (amplitud relativa) se marcan con un '+' if isempty(ITDG)==0 format shorteng; disp('------------------');disp(['ITDG = ' num2str(ITDG) 's.']); disp('------------------') else format shorteng; disp('------------------');disp(['ITDG = ' num2str(itdg) 's.']); disp('------------------') end %Calculo de C80 y C50 for i=1:N [C80(1,i),C50(1,i)] = CxDSP(y(:,i),fs); end % Gráficas de C50, C80 figure(4) semilogx(f,C50,'-om',f,C80,'-sg');xlabel('Frecuencia (Hz.)'); grid on ylabel('Cx (dB)');legend('C50','C80');
6.1 FiltthirdoctDSP.m function [Y,f] = filtthirdoctDSP(x,Fs,B) % Función que filtra la señal de entrada y devuelve un vector de columnas % cada una con la señal filtrada por tercio de octava o banda de octava. % B = tipo de filtro aplicado. Si B=1, filtro de banda de octava, si B=3 % filtro de tercio de octava. % Rango del filtro de tercio: 40 Hz - 15861 Hz (F. Centrales) % Rango del filtro de octava: 32 Hz - 8 kHz (F. Centrales) if B ~= 1 && B ~= 3 error('Error: B debe ser igual a 3 (Filtrado por tercio)o 1 (Filtrado por octava') end if B == 3 %Primero diseñar para 1kHz F0 = 1000; h = fdesign.octave(3,'Class 0', 'N,F0',8,F0,Fs); %Obtener frecuencias normalizadas F0 = validfrequencies(h); y = zeros(length(x),1); for i = 1:length(F0) h.F0 = F0(i); %Determinar propiedad F0 para h H(i) = design(h,'butter'); y(:,i) = filter(H(i),x); end Y = [y(:,3),y(:,4),y(:,5),y(:,6),y(:,7),y(:,8),y(:,9),y(:,10),y(:,11),... y(:,12),y(:,13),y(:,14),y(:,15),y(:,16),y(:,17),y(:,18),y(:,19),... y(:,20),y(:,21),y(:,22),y(:,23),y(:,24),y(:,25),y(:,26),... y(:,27),y(:,28)]; f = F0(3:1:28); else
Calderón De Palma, Serattin, Alba
16
Procesamiento Digital de Señales
Diciembre 2014
%Primero diseñar para 1kHz F0 = 1000; h = fdesign.octave(1,'Class 0', 'N,F0',8,F0,Fs); %Obtener frecuencias normalizadas F0 = validfrequencies(h); y = zeros(length(x),1); for i = 1:length(F0) h.F0 = F0(i); %Determinar propiedad F0 para h H(i) = design(h,'butter'); y(:,i) = filter(H(i),x); end Y = [y(:,1),y(:,2),y(:,3),y(:,4),y(:,5),y(:,6),y(:,7),... y(:,8),y(:,9)]; f = F0(1:1:end-1); end end
6.2 SchrointDSP.m function hdb = SchrointDSP(x) %Funcion que calcula la integral en reversa de Schroeder para una señal x h = cumtrapz(x(end:-1:1).^2); % Integración en reversa h = h(end:-1:1); % Dar vuelta la curva % Normalizar curva y dar valor en dB hdb = 10*log10(abs(h)+eps)-10*log10(abs(max(h))); end
6.3 TRDSP.m function [TR20,TR30] = TRDSP(x,fs) % Función que calcula TR20 y TR30 a partir de la integral de Schroeder, en % principio busca la caida de 5dB y luego las respectivas caidas de 30 db y % 20 db. Finalmente, utilizando la frecuencia de muestreo (fs), calcula el % tiempo de reverberación. % x = integral de Schroeder
r1 = find(x > -5,1,'last'); r2 = find(x > -25,1,'last'); r3 = find(x > -35,1,'last'); TR20 = 3*(r2-r1)/fs; TR30 = 2*(r3-r1)/fs; end
6.4 CxDSP.m function [C80,C50] = CxDSP(x,fs) %Función que calcula los parametros de claridad definidos en norma ISO3382 %Se utiliza regla de trapecios de matlab
Calderón De Palma, Serattin, Alba
17
Procesamiento Digital de Señales
Diciembre 2014
num80 = trapz(x(1:1:0.08*fs).^2); den80 = trapz(x(0.08*fs:1:end).^2); num50 = trapz(x(1:1:0.05*fs).^2); den50 = trapz(x(0.05*fs:1:end).^2); C80 = 10*log10(num80/den80); C50 = 10*log10(num50/den50); end
6.5 ITDG.m function [itdg,idx,maximo,picos] = ITDG(env,fs,delta) %Función a la que se ingresa la envolvente,fs y un parametro delta. %luego la analiza por pasos de tiempo (delta), buscando los maximos %locales y finalmente calcula el máximo local más grande y %lo utiliza para calcular ITDG %Tambien se añade el calculo de maximos locales (picos), evaluados cada %delta segundos. %IMPORTANTE: Se corta el primer milisegundo que se supone corresponde al %sonido directo y se analizan solo 30ms de la señal, debido a estudios de %Beranek que sugiere un ITDG
