2D FFT what to do after converting both matrix into FFT-ed form?

Question

Assume that I have 2 matrix: image, filter; with size MxM and NxN.

My regular convolution looks like this and produces matrix output size (M-N+1)x(M-N+1). Basically it places the top-left corner of a filter on a pixel, convolute, then assign the sum onto that pixel:

for (int i=0; i<M-N; i++)
for (int j=0; j<M-N; j++)
{
    float sum = 0;
    for (int u=0; u<N; u++)
    for (int v=0; v<N; v++) 
        sum += image[i+u][j+v] * filter[u][v];
    output[i][j] = sum;
}

Next, to perform FFT:

1. Apply zero-padding to both image, filter to the right and bottom (that is, adding more zero columns to the right, zero rows to the bottom). Now both have size (M+N)x(M+N); the original image is at image[0->M-1][0-M-1].

2. (Do the same for both matrix) Calculate the FFT of each row into a new matrix, then calculate the FFT of each column of that new matrix.

Now, I have 2 matrices imageFreq and filterFreq, both size (M+N)x(M+N), which is the FFT-ed form of the image and the filter.

But how can I get the convolution values that I need (as described in the sample code) from them?

Answer 1

convolution between A,B using FFT is done by per element multiplication in the frequency domain so in 1D something like this:

1. convert A,B by FFT

   assuming the sizes are N,M of A[N],B[M] first zero pad to common size Q which is a power of 2 and at least M+N in size and then apply FFT:

   Q = exp2(ceil(log2(M+N)));
   zeropad(A,Q);
   zeropad(B,Q);
   a = FFT(A);
   b = FFT(B);
   

2. convolute

   in frequency domain use just element wise multiplication:

   for (i=0;i<Q;i++) a[i]*=b[i];
   

3. reconstruct result

   simply apply IFFT (inverse of FFT)...

   AB = IFFT(a); // crop to first N (real) elements
   

   and use only the first N element (unless algorithm used need more depends on what you are doing...)

For 2D you can either convolute directly in 2D (using 2 nested for loops) or convolve each axis separately. Beware that separating axises need also to normalize the result by some constant (which depends on dimensionality, resolution and kernel used)

So when put together (also assuming the same resolution NxN and MxM) first zero pad to (QxQ) and then:

Q = exp2(ceil(log2(M+N)));
zeropad(A,Q,Q);
zeropad(B,Q,Q);
a = FFT(A);
b = FFT(B);
for (i=0;i<Q;i++)
 for (j=0;j<Q;j++) a[i][j]*=b[i][j];
AB = IFFT(a); // crop to first NxN (real) elements

And again crop to AB to NxN size (unless ...) for more info see:

• How to compute Discrete Fourier Transform?

and all sublinks there... Also here at the end is 1D convolution example using NTT (its a special form of FFT) to compute bignum multiplication:

• Modular arithmetics and NTT (finite field DFT) optimizations

Also if you want real result then just use only the real parts of the result (ignore imaginary part).