Fast convolution

Fast convolution of an audio signal with the impulse response of a system modeled as LTI (Linear Time Invariant) one.

Uses FFT to perform the convolution by multiplication in the frequency domain. To reduce in this way, the computation time.

Works with mono or stereo input signals.

%
%   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( file, impulse, output )
% Fast convolution of an audio signal with the impulse response of a system modeled as LTI (Linear Time Invariant) one.
% USE: fast_conv( 'source.wav', 'impulse_response.wav', 'output.wav' );

h_stereo = true;
x_stereo = true;

clip_factor = 1.01; % clip factor default value

% impulse response:
[ h, Fs1 ] = wavread( impulse ); % 2 columns matrix, one per channel
h1 = h(:,1); % col 1 -> left channel
h_channels = size(h); h_channels = h_channels(2);
if( h_channels == 2 )
	h2 = h(:,2); % col 2 -> right channel
else
	h2 = h(:,1); % col 1 -> left channel
	h_stereo = false;
end

% file to process:
[ x, Fs2 ] = wavread( file );
x1 = x(:,1); % col 1 -> left channel
x_channels = size(x); x_channels = x_channels(2);
if( x_channels == 2 )
	x2 = x(:,2); % col 2 -> right channel
else
	x2 = x(:,1); % col 1 -> left channel
	x_stereo = false;
end

% if sample rates are the same
if( Fs1 == Fs2 )

	% searches for the amount of points required to perform the FFT
	L = length(h1) + length(x1) - 1; % linear convolution length

	disp('Processing...');

	
	% Use nextpow if you have it to replace below code -> N = 2^nextpow(L);
	N = 2;
	while( N<L ) % first 2 pow greater than L
		N = N*2;
	end
	

	% Note: N>=L is needed because the IDFT of the multiplication is the circular convolution and to match it to the
	% common one, N>=L is required (where L=N1+N2-1;N1=length(x);N2=length(h))
	

	% FFT(X,N) is the N points FFT, with zero padding if X has less than N points and truncated if has more.
	H1 = fft( h1, N ); % Fourier transform of the impulse
	X1 = fft( x1, N ); % Fourier transform of the input signal
	
	F_y1 = H1 .* X1; % spectral multiplication

	y1 = ifft( F_y1 ); % time domain again

	% audio normalization: if "y = y/max(y)" -> "values out of [-1,+1] range are clipped"
	y1 = y1/( max(y1)*clip_factor ); % to avoid clipping

	% if one of the files is stereo, process the other channel
	if( h_stereo || x_stereo )
		H2 = fft( h2, N ); % Fourier transform of the impulse
		X2 = fft( x2, N ); % Fourier transform of the input signal

		F_y2 = H2 .* X2; % spectral multiplication
		y2 = ifft( F_y2 ); % time domain again

		y2 = y2/( max(y2)*clip_factor ); % to avoid clipping

		y = [ y1 y2 ]; % stereo output (two channels)
	else
		y = y1; % mono output
	end

	wavwrite( y, Fs1, output );

	disp( 'Convolution success.' );

else

	disp('Error: files has different sample rate.');

end