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.
#!/usr/bin/python # -*- coding: utf-8 -*- # # Copyright (C) 2006 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. # # Fast convolution of an audio signal with the impulse response of a system modeled as LTI (Linear Time Invariant) one. # Use: fast_conv.py source.wav impulse_response.wav output.wav from wav_array import * from FFT import fft, inverse_fft from pylab import size import sys def nextpow2( L ): N = 2 while N < L: N = N * 2 return N def fast_conv_vect( x, h ): # searches for the amount of points required to perform the FFT L = size(h) + size(x) - 1 # linear convolution length N = nextpow2( L ) # 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. H = fft( h, N ) # Fourier transform of the impulse X = fft( x, N ) # Fourier transform of the input signal Y = H * X # spectral multiplication y = inverse_fft( Y ) # time domain again return y source = sys.argv impulse = sys.argv conv_output = sys.argv clip_factor = 1.01 # clip factor default value [ h1, Fs1, h_bits ] = wavread( impulse ) # impulse response [ x1, Fs2, x_bits ] = wavread( source ) # file to process if Fs1 == Fs2 : # if sample rates are the same print "Processing..." y1 = fast_conv_vect( x1, h1 ).real # takes the real part to avoid a too small complex part (around e-18) # audio normalization: if "y = y/max(y)" -> "values out of [-1,+1] range are clipped" y1 = y1/( max(y1)*clip_factor ) # to avoid clipping wavwrite( y1, Fs1, conv_output ) print "Output file:", conv_output print "Convolution success." else: print "Error: files has different sample rate."