Convolución circular rápida (aplicación en reverbs)

Si consideramos un recinto o un determinado lugar como un sistema LTI[1] y de alguna forma obtenemos su respuesta impulsiva, podemos convolucionar esta última con cualquier grabación[2] y obtener como resultado como se "escucharía" en ese lugar.

La fórmula para una convolución discreta es la siguiente: donde es la entrada, la respuesta impulsiva e la salida del sistema.

En Matlab por ejemplo podriamos hacer:

y = conv(



x,ri); % x y ri vectores de señal de entrada y respuesta impulsiva 

Pero a medida que los vectores son más largos (más muestras, más segundos de duración por ejemplo) esta operación cada vez tarda más.

Una forma de hacerlo mucho más rápido (la diferencia es muy grande), es pasar todo esto al dominio de la frecuencia, donde la convolución se "transforma" en una multiplicación, operación computacionalmente mucho más sencilla[3].

Queda algo asi:

(cada letra en mayúscula representa la transformada de fourier discreta de cada vector antes mencionado)

Y en Matlab por ejemplo se reduce a:

X = fft(



 x, N ); % transformada de Fourier de la señal de entrada
H = fft( h, N 



); % transformada de Fourier del impulso
Y = H .* X; % multiplicación de espectros (elemento a elemento)
y = ifft( Y 



); % se vuelve al dominio del tiempo 

(fft es un algoritmo rápido para calcular la transformada discreta)

Donde (cantidad de puntos de la transformada) debe ser mayor que la suma de las longitudes de y de - 1, ya que la IDFT (la transformada inversa) de la multiplicacion es la convolución circular y para que esta coincida con la tradicional se debe pedir .

Por otra parte mientras más grande sea N el proceso será más lento. También debe ser potencia de 2, y el vector si hace falta es rellenado con ceros hasta llegar a N puntos.
Este método por ejemplo se usa mucho para realizar reverbs por convolución. Las respuestas impulsivas de varios lugares pintorescos como catedrales y bosques o recintos con determinadas características se pueden conseguir fácilmente "en internet".

Por ejemplo, para hacer "Mayo, Los sonidos de la Plaza", se que primero grabaron la respuesta impulsiva de la plaza tirando petardos (por su sonido fuerte y corto, que se asemeja a las carterísticas de un impulso) para asi obtener las características del sistema formado por la plaza y sus alrededores, y con esos datos generaron grabaciones para ser reproducidas en ese lugar adecuadamente, sin ecos, ni reverberaciones excesivas.

[1] que se comporta de forma lineal e invariante en el tiempo.

[2] preferiblemente grabada en una cámara anecoica o lo más cercano a esto posible.

[3] además, en programas como Matlab/ Octave, la multiplicación de vectores esta superoptimizada.

Código en Matlab, procesa archivos stereo y mono:

function fast_conv( archivo, impulso, salida )
% Convolución rápida de una señal de audio con la respuesta impulsiva de algún sistema modelado como LTI
% USO: fast_conv( 'archivo_a_procesar.wav', 'respuesta_al_impulso.wav', 'archivo_de_salida.wav' );
 
h_stereo = true;
x_stereo = true;
 
% respuesta impulsiva:
[ h, Fs1 ] = 



wavread( impulso ); % matriz de 2 columnas, una por cada canal
h1 = h(:,1); % columna 1 -> Canal izquierdo
h_channels = size



(h); h_channels = h_channels(2);
if( h_channels == 2 )
h2 = h(:,2); % columna 2 -> Canal derecho
else
h2 = h(:,1); % columna 1 -> Canal izquierdo
h_stereo = false;
end
 
% archivo a procesar:
[ x, Fs2 ] = 



wavread( archivo );
x1 = x(:,1); % columna 1 -> Canal izquierdo
x_channels = size



(x); x_channels = x_channels(2);
if( x_channels == 2 )
x2 = x(:,2); % columna 2 -> Canal derecho
else
x2 = x(:,1); % columna 1 -> Canal izquierdo
x_stereo = false;
end
 
% si las frecuencias de muestreo coinciden
if( Fs1 == Fs2 )
 
% busca la cantidad de puntos que debo utilizar para realizar la FFT
L = length(



h1) + length



(x1) - 1; % el tamaño de la convolución lineal
 
disp('



Procesando...');
 
% el siguiente código se puede reemplazar por la función nextpow2 (si esta disponible)
% N = 2^nextpow(L);
N = 2;
while( N
N = N*2;
end
 
H1 = fft( h1, N 



); % transformada de Fourier del impulso
X1 = fft( x1, N 



); % transformada de Fourier de la señal de entrada
 
% FFT(X,N) es la FFT de N puntos, rellenada con ceros si X tiene más de N puntos y truncada si tiene demás.
F_y1 = H1 .* X1; % multiplicación de espectros
 
y1 = ifft(



 F_y1 ); % se vuelve al dominio del tiempo
 
% normalización del audio: si "y = y/max(y)" -> "los valores fuera del rango [-1,+1] son recortados"
y1 = y1/( max



(y1)*1.01 ); % para que no "clipee"
 
% si alguno de los 2 archivos es stereo también procesa el otro canal
if( h_stereo || x_stereo )
H2 = fft( h2, N 



); % transformada de Fourier del impulso
X2 = fft( x2, N 



); % transformada de Fourier de la señal de entrada
 
F_y2 = H2 .* X2; % multiplicación de espectros
y2 = ifft(



 F_y2 ); % se vuelve al dominio del tiempo
 
y2 = y2/( max



(y2)*1.01 ); % para que no "clipee"
 
y = [ y1 y2 ]; % se define la salida stereo
else
y = y1; % se define la salida mono
end
 
wavwrite(



 y, Fs1, salida );
 




disp( 'Convolución realizada con éxito. Fin.' );
 
else
 
disp('



Error: los archivos no tienen la misma frecuencia de sampleo. No se puede realizar la convolución.');
 
end

(bajar archivo fast_conv.m)

Ejemplo1: Original - Resultado

Ejemplo2: Original - Resultado

Nota: los samples tienen licencias libres (ver links), la respuesta impulsiva (es de una conocida catedral) no la puedo distribuir, si consigo alguna libre, hago los ejemplos con esa. Escuchen y comparen.