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 % <firstname.lastname@example.org> % % 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