DSP: Convolución eficiente que no introduce latencia

Introducción

Este artículo está basado en Efficient Convolution without Input-Output Delay.

Conceptos básicos.

Qué es lo que hace BruteFIR

Los conceptos básicos del funcionamiento de BruteFIR están explicados aquí:

A modo de resumen, BruteFIR hace la convolución usando "convolución FFT", en donde:

conv(x,h) = iFFT( FFT(x) * FFT(h))

Pero BruteFIR trabaja por bloques (buffer de memoria). Obviamente, antes de realizar la operación hay que esperar a que se llene el buffer. Si la música está ya en el disco duro, el tiempo puede ser despreciable. Pero si hay alguien está tocando la guitarra en directo, mientras no se llene el buffer, tardará un tiempo en reproducir el sonido. Por ejemplo, si el buffer es de 8000 muestras y el muestreo es de 48000hz, la latencia mínima introducida será de: 8000/48000 = 0.16sg

Cómo conseguir que el delay sea nulo

La respuesta es obvia si pensamos en lo que se llama en la DSPguide: Output Side Algorithm. En esta sección, se interpreta la matemática de la convolución como una máquina que se desplaza hacia la derecha. Esta "máquina" es capaz de generar un sample a medida que llegan los samples de la señal. Éste es el algoritmo normal de la convolución.

En la siguiente imagen, vemos la "máquina". Llamamos x[n] a la señal original, h[n] es el filtro e y[n] es la salida. Para generar la muestra de color verde en la salida, se ven involucradas las muestras azules de x[n] y las muestras amarillas de h[n] (observar que h[n] está colocado en orden inverso: la muestra 0 a la derecha y la M a la izquierda).

convolution_machine_01.svg

¿Por qué no se usa usa este algoritmo? Porque con filtros largos es muy poco eficiente, como se ve en la gráfica de la DSPguide (tened en cuenta que nuestros filtros pueden tener del orden de las 32000 muestras y la gráfica sólo llega a 1024).

¿Cómo conseguir tener "cero delay" y buena eficiencia?

Obviamente esto nos lleva a pensar en un método mixto que mezcle la convolución directa y la convolución por bloques. Esto es posible gracias a que el sistema es lineal (ver la superposición).

Para ello, lo que haremos será combinar la señal x[n], con bloques de h[n]: unos usando convolución directa y otros usando convolución FFT.

En la imagen, vemos que hemos dividido h en varios bloques: h1, h2, h3. Hacemos la convolución de h1 con x obteniendo y1, la convolución de h2 con x obteniendo y2, … El resultado final es la suma de y1, y2 e y3. Para el bloque amarillo (h1), usamos convolución directa. Para el resto usamos convolución FFT por bloques.

convolution_machine_02.svg

Observaciones:

  • La convolución de h1 con x, genera y1, que vemos que tiene un delay de 0 segundos. (Según entra x, se genera y1)
  • ¿Qué ocurre cuando "S>M1"?: podríamos pensar que cuando se gasta h1 se introduce delay, pero no es así. Visualmente, vemos que cuando se ha usado por completo h1, necesitaremos h2. Pero vemos visualmente que al convolucionar h2 se hace con muestras "antiguas" de "x", no con las nuevas que van entrando.
  • ¿Qué valor debería tener M1? Según la DSPguide, vemos que hasta 60 muestras la convolución directa es más rápida que la convolución FFT.
  • Con la convolución FFT, interesa usar bloques grandes, ya que procesar el doble de muestras no supone el doble de tiempo de procesado, sino mucho menos.
  • ¿Hay algún límite en el tamaño que puede tener h2? Sí lo hay. Vimos que h1 (que tiene M1 muestras) se combina como máximo con las primeras M1 muestras de "x". Eso significa que en la siguiente muestra de "y" se necesitará información de y2 (además de y1). Para generar "y2" se usa "h2" y las primeras muestras de "x" (tendremos "M1" muestras de "x" disponibles). Intuitivamente (pensando en la máquina), vemos que como sólo disponemos de M1 muestras de "x" el tamaño de h2 será como máximo de M1. Para el siguiente bloque (h2) será la suma de los bloques previos: M1+M1 = 2* M1. Para el siguiente (h3), la suma de los bloques previos: M1 + M1 + 2M1 = 4 M1. Resumiendo:

Req1: Un bloque tendrá como máximo la suma de muestras de todos los bloques anteriores.

Además:

NOTA: Si a cada bloque se le asigna la suma de muestras de los bloques anteriores, se obtiene el algoritmo de mínimo coste. Tiene la ventaja de ser la forma en la que el consumo de procesador es mínima, y la desventaja de que hace que el uso del procesador no sea homogéneo.

Uso del procesador

Es obvio que podríamos recalcular la convolución "h2" a medida que nos llegan nuevos samples, sin embargo, lo más eficiente es que el número de convoluciones calculadas sea el mínimo posible. ¿Cuándo hay que calcular la convolución con "h2"?

Req2: el cálculo de la convolución de "h2" tendrá que terminar antes de que se consuma el sample [M1+1].

El Req2 es obvio. Para el sample [M1+1], necesitamos información de la convolución con h2. Lo más efectivo sería hacer la convolución entre las M1 muestras de "x" y las M2-M1 muestras de "h2". El tiempo disponible para realizar ese cálculo será el de la duración de un sample (entre el sample M1 y el sample M1+1).
Cada cuánto tiempo hay que recalcular esta convolución. Para realizar el mínimo número de operaciones posibles, se necesitaría otra conv(x,h) cuando se tengan otras "M1" muestras de "x".

Consideramos "L = M2-M1" el númeo de muestras que tiene "h2" y "J=M1" el número de muestras que tiene "h1". Para minimizar el número de convoluciones que realizamos con h2, tendremos que tomar el mismo número de muestras de "x" que tengamos en h2, es decir, "L" muestras de "x". El cálculo tendrá que estar realizado cuando se tengan "J" muestras de "x". El requerimiento equivalente a Req2 será:

Req2: L<=J.

En la gráfica vemos este requerimiento. El cálculo de la convolución entre "h2" y "L" samples de "x" puede comenzar cuando se tienen "J" muestras de "x". A partir de ese momento, se dispone hasta la muestra "L" para tener terminado el cálculo. Si queremos maximizar ese tiempo disponible haremos "L = 2*J". Con esa condición, el tiempo disponible para realizar el cálculo será máximo (y en principio permite un uso del procesador más eficiente, dado que elimina los tiempos en los que el procesador no puede ser usado. Vemos también que si L fuese mayor que J, no habría tiempo para realizar los cálculos, por los que se introduciría latencia.

uso_procesador.svg

Especificación del algoritmo.

  • REQ1: El filtro "h" irá dividido en bloques.
  • REQ2: El primer bloque se calculará usando convolución directa.
  • REQ3: El resto de bloques se calculará mediante convolución FFT.
  • REQ4: El tamaño del primer bloque deberá ser pequeño (por ejemplo 64 o 32 muestras).
  • REQ5: El segundo bloque tendrá un número de muestras menor o igual al primer bloque.
  • REQ6: Del tercer bloque en adelante tendrán un número de muestras menor o igual a la suma de todas las muestras de los bloques anteriores.
  • REQ7: El número de muestras del segundo bloque no debería ser menor que la mitad del número de muestras del primer bloque. (No estrictamente un requerimiento, pero evita recalcular información ya disponible.

Observaciones

El documento mencionado al comienzo del artículo, considera el caso de mínimo coste (L = J) y el de uso del procesador más homogéneo (2L = J). Yo prefiero permitir el margen: L <= J <= 2L, y permitir que la disponibilidad de la información sea quien haga que nos acerquemos más a un criterio u a otro.

Salvo que se diga otra cosa, el contenido de esta obra está bajo la licencia: Creative Commons Reconocimiento NoComercial CompartirIgual 2.5 España.