[Numpy-discussion] reformulating a series of correlations as one fft, ifft pair

Paul Northug pnorthug@gmail....
Mon Apr 19 23:51:09 CDT 2010


 <josef.pktd <at> gmail.com> writes:

>
> On Mon, Apr 19, 2010 at 2:36 AM, Paul Northug <pnorthug <at> gmail.com> wrote:
<snip>
> >
> > #
> > import numpy as np
> > from scipy.signal import fftn, ifftn, correlate
> >
> > M, N, P, T = 2, 3, 4, 5
> >
> > phi = np.random.randn(M, N, P)
> > a = np.random.randn(N, P+T-1)
> > x = np.zeros((M, T))
> >
> > # 1. correlate
> >
> > x = np.array([correlate(phi[i], a, 'valid').squeeze() for i in range(M)])
> >
> > # 2. dot products
> >
> > y = np.zeros_like(x)
> > for b in range(P): y += np.dot(phi[:,:,b], a[:,b:b+T])
> >
> > # 3. a series of ffts (code adapted from scipy convolve)
> >
> > s1 = np.array(phi.shape[1:])
> > s2 = np.array(a.shape)
> > size = s1+s2-1
> > z = np.empty_like(y)
> > fsize = (2**np.ceil(np.log2(size))).astype(np.int) # power of 2 pad
> > af = fftn(a, fsize)
> > fslice = tuple([slice(0, sz) for sz in size])
> > for i in range(M):
> >    phif = fftn(phi[i], fsize).conj()
> >    z[i] = np.real(ifftn(phif * af)[fslice])[0,:T]
> >
> > # 4. <BROKEN> single fft, ifft pair - some rows are correct
> >
> > s1 = np.array(phi.shape)
> > s2 = np.array(a[np.newaxis].shape)
> > size = s1+s2-1
> > fsize = (2**np.ceil(np.log2(size))).astype(np.int)
> > af = fftn(a[np.newaxis], fsize)
> > phif = fftn(phi, fsize).conj()
> > fslice = tuple([slice(0, sz) for sz in size])
> > w = np.real(ifftn(phif * af)[fslice])[:,:,:T]
>
> I only checked this in your example
>
> >>> np.squeeze(signal.fftconvolve(phi[:,::-1,::-1], a[None,...],mode='valid'))
> array([[ 0.42695078,  4.3433924 , -2.53084395, -4.22672281, -5.16705981],
>        [-2.35108765, -2.15060513, -1.23528122,  1.40386599,  1.38189896]])
>
> >>> y
> array([[ 0.42695078,  4.3433924 , -2.53084395, -4.22672281, -5.16705981],
>        [-2.35108765, -2.15060513, -1.23528122,  1.40386599,  1.38189896]])
> >>> x
> array([[ 0.42695078,  4.3433924 , -2.53084395, -4.22672281, -5.16705981],
>        [-2.35108765, -2.15060513, -1.23528122,  1.40386599,  1.38189896]])
>
> Josef
>

Thanks Josef. You saved me a lot of work. I was missing the middle ::-1.

The fft approach is orders of magnitude slower than method #2 in this
case. Using CUDA fft3, it is about 15 times slower for the desired
problem size. It is not a fair comparison as the convolution is
degenerate along 2 dimensions but it is useful to know.


More information about the NumPy-Discussion mailing list