Fast audio 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.

#!/usr/bin/python
# -*- coding: utf-8 -*-

#
# Copyright (C) 2006 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.
#

# 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[1]
impulse = sys.argv[2]
conv_output = sys.argv[3]

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."