[SciPy-user] Re: Help on performance of signal.convolve

[Travis Oliphant]

Yichun Wei wrote:

> Hello Experts,
>
> Could you give some more detailed explanation about 
> fftpack.convolve.init_convolution_kernel ? It is somewhat hard for me 
> to guess its behavior... It seems not designed for general purpose 
> convolution, but for constructing filters.


I'm not sure about fftpack.convolve but I think it is limited to 1-d 
convolution.

Your large convolutions are usually done using the Fourier Transform (as 
the direct method implemented by convolveND will be slow for large data 
-- though it currently could use some optimizations).

Basically, using the Fourier transform takes advantage of the fact that 
multiplication in the (discrete) Fourier domain is the same as 
*periodic* convolution in the originating domain.   To make this the 
same as linear convolution, you have to zero pad to N1 + N2 - 1 where N1 
and N2 are the lengths of the same dimensions of the two arrays. 

So, something like the following untested code should get you started:

def convolvefft(arr1,arr2):
       s1 = array(arr1.shape())
       s2 = array(arr2.shape())
       fftsize = s1 + s2 - 1
       # finds the closest power of 2 in each dimension (you may comment 
this out and compare speeds)
       fftsize= pow(2,ceil(log2(fftsize))) 

       ARR1 = fftn(arr1,fftsize)
       ARR2 = fftn(arr2,fftsize)
      
       RES = ifftn(ARR1*ARR2) 
       #RES=RES[validpart]  # I'm not sure how to get the correct part 
--- first try would be to just truncate to the shape you wanted
       return RES