How to implement the fast fourier transform to correlate two 2d arrays?

Question

I want to do this with larger arrays, but so that I can understand it, i'll use a smaller example.

Let's say I have the following array:

    A = [[0, 0, 100],
         [0, 0, 0],
         [0, 0, 0]] 

If I want to calculate the correlation of this array with another, by multiplying the corresponding entries. For example, A*A would be equal to A (1*1 = 1, zeros everywhere else). I read that fast fourier transforms can be used to speed this up with large arrays. From what I read, I got the impression if I wanted to multiply two arrays A and B, as in A*B, that I could do it quicker with the following (using numpy in python):

a = np.conj(np.fft.fftn(A))
b = np.fft.fftn(B)
c = np.fft.ifft(a*b)

So in effect, take the fft of A, take the fft of B, multiply the two results, and then get the inverse of that result. However, I tried it with the case of A given above, multiplying itself. I had hoped the inverse multiplication would give me

[[0, 0, 10000],
 [0, 0, 0    ],
 [0, 0, 0    ]]

However I got something a bit different, closer to

[[10000, 0, 0],
 [10000, 0, 0],
 [10000, 0, 0]]

does anyone know what's going on? Sorry, I'm guessing there's something I'm misunderstanding about the fft.

Answer 1

If you must implement your own, you should be aware that multiplication in the frequency domain corresponds to circular convolution in the time domain. To obtain the desired linear convolution you need to pad both arrays with zeros, to a length at least twice the size of your original matrix¹:

s = [2*x for x in np.shape(A)]
a = np.conj(np.fft.fftn(A,s))
b = np.fft.fftn(B,s)
c = np.fft.ifftn(a*b)

_{¹ strictly speaking a size of 2n-1 (instead of 2n) will do, but FFTs tend to perform better when operating on sizes that are multiples of small prime factors.}

Answer 2

You should use scipy.signal.fftconvolve instead.

It is already implemented and has been extensively tested, particularly regarding the handling the boundaries. The only additional step necessary to go from the convolution to the correlation operator in 2D is to rotate the filter array by 180° (see this answer).