A Python Program for convolution

Here is an example of fast convolution in a Python script. It's nice that the language allows us to show the process quite well. It requires pylab, scipy and numpy. For playback I am using the sndfile lib sndfile-play program.

# fast convolution example
# (c) V Lazzarini, 2010
# GNU Public License
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

import sys
from scipy import *
from pylab import *
from scipy.io import wavfile


# read impulse and signal input
(sr,impulse) = wavfile.read(sys.argv[1])
(sr,signalin) = wavfile.read(sys.argv[2])

S = len(impulse) # impulse length
L = len(signalin) # signal length

# find fftsize as the next power of 2
# beyond S*2-1
N = 2
while(N <= S*2-1): N *= 2

# signal blocks for FFT
imp = zeros(N)
sig = zeros(N)

# output signal & amplitude
sigout = zeros(L+S-1)
amp = max(signalin)

# copy impulse for FFT
imp[0:S] = impulse
spec_imp = fft(imp)

for i in range(0, L/S):

p = S*i
# get block from input
sig[0:S] = signalin[p:p+S]
# perform convolution and overlap-add
sigout[p:p+2*S] += ifft(spec_imp*fft(sig))[0:2*S]

# write file to output, scaling it to original amp
wavfile.write(sys.argv[3], sr,
array(amp*sigout/max(sigout), dtype='int16'))

# play it using a libsndfile utility
import os
try:

os.spawnlp(os.P_WAIT, 'sndfile-play',
'sndfile-play', sys.argv[3])

except:

pass