# How to manually implement convolution with FFTs?

I'm trying to manually implement a convolution using FFTs and it isn't working as expected. I know I'm probably missing some subtlety with padding, shifting, or conjugation, (all of which I've tried playing with) but can't find it. Here's what I've got:

import numpy as np
import scipy.signal as sig
import matplotlib.pyplot as plt

# Both of these work as expected, and give the same result
def myconvolve(f, g) :
#return sig.fftconvolve(f, g, mode='same')
return sig.convolve2d(f, g, 'same', boundary='fill', fillvalue=0)

# Here's my version, which doesn't work
def myconvolve2(f, g) :
p1, p2 = (np.r_[f.shape]-g.shape).astype(int)//2
F = np.fft.fft2(f)
G = np.fft.fft2(g)
FG = F*np.conj(G)
res = np.real(np.fft.ifft2(FG))
return res

# Generate some test functions
X,Y = np.mgrid[-1:1:20j,-1:1:20j]
g = np.exp( -X**2 -Y**2 )**2
np.random.seed(1)
f = np.random.rand(100,70)
f[f<.998] = 0

# Convolve
fc = myconvolve2(f,g)

# Plot
plt.figure(1,figsize = (10,15), dpi=80); plt.clf(); plt.imshow(f.T)
plt.figure(2,figsize = (10,15), dpi=80); plt.clf(); plt.imshow(g.T)
plt.figure(3,figsize = (10,15), dpi=80); plt.clf(); plt.imshow(fc.T)


Output:

left to right: f, g, f*g (correct), f*g (my version)

• only had a cursory glance but I see inside your function definition you never actually do anything with F1 and F2? – DamienBradley Jan 11 '20 at 2:52
• You're right. In trimming things down for SE, I accidentally removed that line! Fixed. – argentum2f Jan 12 '20 at 16:56
• Also, does dsp not support syntax highlighting like other SE sites do? – argentum2f Jan 13 '20 at 14:54

I figured out my problem. The kernel needs to be shifted so the 'center' is on the corner of the image (which acts as the origin in an FFT). The built-in ifftshift function works just fine for this. (Note, there are some subtleties here depending on whether you have even or odd shapes or differences in shape that can result in one row or column shifts. I picked ifftshift in this case because it matches the result of the comparison function.) There's also a very faint but visible difference between my implementation and the library function (a faint streaking) that could just be numerical errors or something.

import numpy as np
import scipy.signal as sig
import matplotlib.pyplot as plt

def myconvolve(f, g) :
return sig.fftconvolve(f, g, mode='same')
return sig.convolve2d(f, g, 'same', boundary='fill', fillvalue=0)

def myconvolve2(f, g) :
# Pad g to equal size of f. This assumes f is smaller in both dimensions
# and that the difference between f and g dimensions are even
p1, p2 = (np.r_[f.shape]-g.shape).astype(int)//2

# Shift g to 'center' on top left corner (the 'origin' in an fft)

# Multiply spectra

# Generate some test functions
X,Y = np.mgrid[-1:1:21j,-1:1:21j]
g = np.exp( -X**2*2 -Y**2*2 )**2
np.random.seed(1)
f = np.random.rand(101,71)
f[f<.998] = 0

# Convolve