using Fourier transforms to do convolution?

Question

According to the Convolution theorem, we can convert the Fourier transform operator to convolution.

Using Python and Scipy, my code is below but not correct. Can you help me and explain it?

import tensorflow as tf
import sys
from scipy import signal

from scipy import linalg
import numpy as np

x = [[1 , 2] , [7 , 8]]

y = [[4 , 5] , [3 , 4]]


print "conv:" ,  signal.convolve2d(x , y , 'full')

new_x = np.fft.fft2(x)
new_y = np.fft.fft2(y)


print "fft:" , np.fft.ifft2(np.dot(new_x , new_y))

The result of code:

conv: [[ 4 13 10]
 [31 77 48]
 [21 52 32]]

fft: [[  20.+0.j   26.+0.j]
 [ 104.+0.j  134.+0.j]]

I'm confused!

Answer 1

The problem may be in the discrepancy between the discrete and continuous convolutions. The convolution kernel (i.e. y) will extend beyond the boundaries of x, and these regions need accounting for in the convolution.

scipy.signal.convolve will by default pad the out of bounds regions with 0s, which will bias results: https://docs.scipy.org/doc/scipy-0.18.1/reference/generated/scipy.signal.convolve2d.html

The Fourier multiplication will not do this by default - you could test this by making padded x, y arrays and comparing the results.

The discrepancy between such techniques should diminish as the kernel size becomes much less than the image dimensions.

As a further note - you should not use the dot product between new_x, new_y. Instead, just multiply the arrays with the * operator.

Hope this helps.

Answer 2

I answer my question. The correct code.

import sys
from scipy import signal

from scipy import linalg
import numpy as np

x = [[1 , 0 , 0 , 0] , [0 , -1 , 0 , 0] , [0 , 0 , 3 , 0] , [0 , 0 , 0 , 1]]
x = np.array(x)
y = [[4 , 5] , [3 , 4]]
y = np.array(y)

print "conv:" ,  signal.convolve2d(x , y , 'full')

s1 = np.array(x.shape)
s2 = np.array(y.shape)

size = s1 + s2 - 1


fsize = 2 ** np.ceil(np.log2(size)).astype(int)
fslice = tuple([slice(0, int(sz)) for sz in size])


new_x = np.fft.fft2(x , fsize)


new_y = np.fft.fft2(y , fsize)
result = np.fft.ifft2(new_x*new_y)[fslice].copy()

print "fft for my method:" , np.array(result.real , np.int32)

print "fft:" , np.array(signal.fftconvolve(x ,y) , np.int32)