2

I have a set of of 2D arrays that I have to compute the 2D correlation of. I have been trying many different things (even programming it in Fortran), but I think the fastest way will be calculating it using FFT.

Based on my tests and on this answer I can use scipy.signal.fftconvolve and it works fine if I'm trying to reproduce the output of scipy.signal.correlate2d with boundary='fill'. So basically this

scipy.signal.fftconvolve(a, a[::-1, ::-1], mode='same')

is equal to this (with the exception of a slight shift)

scipy.signal.correlate2d(a, a, boundary='fill', mode='same')

The thing is that the arrays should be computed in wrapped mode, since they are 2D periodic arrays (i.e., boundary='wrap'). So if I'm trying to reproduce the output of

scipy.signal.correlate2d(a, a, boundary='wrap', mode='same')

I can't, or at least I don't see how to do it. (And I want to use the FFT method, since it's way faster.)

Apparently Scipy used to have something like that that might have done the trick, but apparently it got left behind and I can't find it, so I think Scipy might have dropped support for it.

Anyway, is there a way to use scipy's or numpy's FFT routines to calculate this correlation of period arrays?

Community
  • 1
  • 1
TomCho
  • 3,204
  • 6
  • 32
  • 83
  • how does your `a` look like? – kmario23 Apr 19 '17 at 19:56
  • @kmario23 It's private experimental data so I can't really share here. But it's 200x200 arrays that are periodic in x and y. – TomCho Apr 19 '17 at 19:58
  • 1
    I don't have a good answer at the moment, but that "left behind" code you mentioned can be found [here](https://github.com/scipy/scipy/blob/maintenance/0.7.x/scipy/stsci/convolve/lib/Convolve.py#L144) – jakevdp Apr 19 '17 at 20:44

1 Answers1

4

The wrapped correlation can be implemented using the FFT. Here's some code to demonstrate how:

In [276]: import numpy as np

In [277]: from scipy.signal import correlate2d

Create a random array a to work with:

In [278]: a = np.random.randn(200, 200)

Compute the 2D correlation using scipy.signal.correlate2d:

In [279]: c = correlate2d(a, a, boundary='wrap', mode='same')

Now compute the same result, using the 2D FFT functions from numpy.fft. (This code assumes a is square.)

In [280]: from numpy.fft import fft2, ifft2

In [281]: fc = np.roll(ifft2(fft2(a).conj()*fft2(a)).real, (a.shape[0] - 1)//2, axis=(0,1))

Verify that both methods give the same result:

In [282]: np.allclose(c, fc)
Out[282]: True

And as you point out, using the FFT is much faster. For this example, it is about 1000 times faster:

In [283]: %timeit c = correlate2d(a, a, boundary='wrap', mode='same')
1 loop, best of 3: 3.2 s per loop

In [284]: %timeit fc = np.roll(ifft2(fft2(a).conj()*fft2(a)).real, (a.shape[0] - 1)//2, axis=(0,1))
100 loops, best of 3: 3.19 ms per loop

And that includes the duplicated computation of fft2(a). Of course, fft2(a) should only be computed once:

In [285]: fta = fft2(a)

In [286]: fc = np.roll(ifft2(fta.conj()*fta).real, (a.shape[0] - 1)//2, axis=(0,1))
Warren Weckesser
  • 110,654
  • 19
  • 194
  • 214