%
%   Copyright (C) 2004 by Hernán Ordiales
%   <audiocode@uint8.com.ar>
%
%   This program is free software; you can redistribute it and/or modify
%   it under the terms of the GNU General Public License as published by
%   the Free Software Foundation; either version 2 of the License, or
%   (at your option) any later version.
%

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;

clip_factor = 1.01; % valor por defecto


% 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<L ) % se busca la primer potencia de 2 mayor a L
		N = N*2;
	end

	% Nota: se pide N>=L ya que la IDFT de la multiplicacion es la convolución circular y para que esta
	% equivalga a la tradicional se debe pedir N>=L (donde L=N1+N2-1;N1=longitud(x);N2=longitud(h))

	% 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.
	H1 = fft( h1, N ); % transformada de Fourier del impulso
	X1 = fft( x1, N ); % transformada de Fourier de la señal de entrada
	
	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)*clip_factor ); % 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)*clip_factor ); % 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