Convolution computations in Numpy/Scipy

So I tested this out and can now confirm a few things:

1) numpy.convolve is not circular, which is what the fft code is giving you:

2) FFT does not internally pad to a power of 2. Compare the vastly different speeds of the following operations:

x1 = np.random.uniform(size=2**17-1)x2 = np.random.uniform(size=2**17)np.fft.fft(x1)np.fft.fft(x2)

3) Normalization is not a difference -- if you do a naive circular convolution by adding up a(k)*b(i-k), you will get the result of the FFT code.

The thing is padding to a power of 2 is going to change the answer. I've heard tales that there are ways to deal with this by cleverly using prime factors of the length (mentioned but not coded in Numerical Recipes) but I've never seen people actually do that.

scipy.signal.fftconvolve does convolve by FFT, it's python code. You can study the source code, and correct you mix1 function.

As mentioned before, the scipy.signal.convolve function does not perform a circular convolution. If you want a circular convolution performed in realspace (in contrast to using fft's) I suggest using the scipy.ndimage.convolve function. It has a mode parameter which can be set to 'wrap' making it a circular convolution.

for idx, row in enumerate(outputs):    outputs[idx] = sp.ndimage.convolve(signal1[idx], signal2[idx], mode='wrap')