Convolution in the frequency domain

Let \(\mathbf{u} \in \mathbb{R}^{N_r \times N_c}\) and \(\mathbf{h} \in \mathbb{R}^{L_r \times L_c}\) represent variables with zero-pad boundary conditions, i.e.

\[\begin{split}\begin{array}{rcl} u_{n_1,n_2} = \left\{ \begin{array}{lc} 0 & n_1,\,n_2 < 0 \\ u_{n_1,n_2} & \\ 0 & n_1 \geq N_r, \, n_2 \geq N_c \end{array}\right. & \; \mbox{ and } \; & h_{n_1,n_2} = \left\{ \begin{array}{lc} 0 & n_1,\,n_2 < 0 \\ h_{n_1,n_2} & \\ 0 & n_1 \geq L_r, \, n_2 \geq L_c \end{array}\right. \end{array} \end{split}\]

Then

\(\begin{array}{l} \mbox{Linear Convolution} \\ \hline \begin{array}{rcl}\mathbf{b} & = & \mathbf{u} * \mathbf{h} \\ b_{k_1,k_2}& = & \displaystyle \sum_{n_1=0}^{L_r} \sum_{n_2=0}^{L_c} h_{n_1 , n_2} \cdot u_{n_1 - k_1, n_2 - k_2} \end{array} \\ \; \\ \mbox{Circular Convolution} \\ \hline \begin{array}{rcl}\mathbf{b} & = & \mathbf{u} \circledast \mathbf{h} \\ \tilde{b}_{k_1,k_2}& = & \displaystyle \sum_{n_1=0}^{L_r} \sum_{n_2=0}^{L_c} h_{n_1 , n_2} \cdot \tilde{u}_{n_1 - k_1, n_2 - k_2} \end{array} \\ \mbox{where} \\ \begin{array}{rcl} \tilde{u}_{n_1,n_2} & = & u_{(n_1)_{N_r},(n_2)_{N_c}} \qquad \color{darkgray} \mbox{ circular extension of } \mathbf{u}\\ (a)_{b} & = & \mbox{remainder}(a,b) \end{array} \end{array}\)

Convolution theorem

Recall that

\[\begin{split} \begin{array}{rcl} % X & = & {\cal{F}}\{ \mathbf{x} \} \\ % X_{k_1,k_2} & = & \displaystyle \sum_{k_1=0}^{N_r-1}\sum_{k_2=0}^{N_c-1} x_{m_1,m_2} \cdot W^{m_1\cdot k_1}_{N_r} \cdot W^{m_2\cdot k_2}_{N_c} \qquad \color{darkgray} W^a_b = \mbox{e}^{-j\cdot 2\pi \cdot \frac{a}{b}} % \end{array} \end{split}\]

and that, by definition,

\({\cal{F}}\{ \mathbf{x} \} = {\cal{F}}\{ \tilde{\mathbf{x}} \} \qquad \color{darkgray}(\mbox{equivalent to } \tilde{\mathbf{x}} = {\cal{F}}^{-1}\{ X \})\)

then

\[ {\cal{F}}\{ \mathbf{u} \circledast \mathbf{h} \} = U \odot H \]

where

\({\cal{F}}\{ \cdot \}\) represents the Fourier transform;
\(H = {\cal{F}}\{ \mathbf{h}, \mbox{shape}= \mathbf{x}\}. \qquad \color{darkgray} \mathbf{h} \mbox{ is zero-padded to match } \mathbf{x}\mbox{'s size}\)

1. Proof

Summary

Let \(\mathbf{x} \in \mathbb{R}^{N_r \times N_c}\) and \(\mathbf{h} \in \mathbb{R}^{L_r \times L_c}\)

\(U = {\cal{F}}\{ \mathbf{u} \}\)
\(H = {\cal{F}}\{ \mathbf{h}, \mbox{shape}= \mathbf{x}\}\)
\(B = X \odot H\).

The proof is based on [DM84].

Step 1

\[\begin{split} % \begin{array}{rcl} \mathbf{b} & = & \displaystyle{\cal{F}}^{-1}\{ U \odot H \} \\ b_{n_1,n_2}& = & \displaystyle\frac{1}{N_r\cdot N_c} \sum_{k_1=0}^{N_r-1}\sum_{k_2=0}^{N_c-1} H_{k_1,k_2} \cdot \bbox[white,4px,border:2px solid white]{U_{k_1,k_2}} \cdot W^{\scriptscriptstyle -n_1\cdot k_1}_{N_r} \cdot W^{-n_2\cdot k_2}_{N_c} \qquad \color{darkgray} W^a_b = \mbox{e}^{-j\cdot 2\pi \cdot \frac{a}{b}} \\ % \end{array} \end{split}\]

\(\color{white}\mbox{Considering } U_{k_1,k_2} = \bbox[white,6px,border:2px solid white]{ \displaystyle \bbox[white,4px,border:2px solid white]{ \sum_{k_1=0}^{N_r-1}\sum_{k_2=0}^{N_c-1} \tilde{u}_{m_1,m_2} } \cdot \bbox[white,4px,border:2px solid white]{ W^{m_1\cdot k_1}_{N_r} \cdot W^{m_2\cdot k_2}_{N_c} } \; }\)

\(\color{white}\therefore\)

\[\begin{split}\color{white}\begin{array}{rcl} b_{n_1,n_2} & = & \displaystyle\sum_{k_1=0}^{N_r-1}\sum_{k_2=0}^{N_c-1} \tilde{u}_{m_1,m_2} \cdot \bbox[white,4px,border:2px solid white]{ \displaystyle\frac{1}{N_r\cdot N_c} \sum_{k_1=0}^{N_r-1}\sum_{k_2=0}^{N_c-1} H_{k_1,k_2} \cdot W^{(m_1-n_1)\cdot k_1}_{N_r} \cdot W^{(m_2-n_2)\cdot k_2}_{N_c} }\\ % & = & \displaystyle\sum_{k_1=0}^{N_r-1}\sum_{k_2=0}^{N_c-1} \tilde{u}_{m_1,m_2} \cdot \tilde{h}_{n_1 - m_1, n_2 - m_2} \\ % & = & h_{n_1,n_2} \circledast u_{n_1,n_2} \end{array} \end{split}\]

Step 2

\[\begin{split} % \begin{array}{rcl} \mathbf{b} & = & \displaystyle{\cal{F}}^{-1}\{ U \odot H \} \\ b_{n_1,n_2}& = & \displaystyle\frac{1}{N_r\cdot N_c} \sum_{k_1=0}^{N_r-1}\sum_{k_2=0}^{N_c-1} H_{k_1,k_2} \cdot \bbox[white,4px,border:2px solid green]{U_{k_1,k_2}} \cdot W^{\scriptscriptstyle -n_1\cdot k_1}_{N_r} \cdot W^{-n_2\cdot k_2}_{N_c} \qquad \color{darkgray} W^a_b = \mbox{e}^{-j\cdot 2\pi \cdot \frac{a}{b}} \\ % \end{array} \end{split}\]

\(\mbox{Considering } U_{k_1,k_2} = \bbox[white,6px,border:2px solid green]{ \displaystyle \bbox[white,4px,border:2px solid lightblue]{ \sum_{k_1=0}^{N_r-1}\sum_{k_2=0}^{N_c-1} \tilde{u}_{m_1,m_2} } \cdot \bbox[white,4px,border:2px solid red]{ W^{m_1\cdot k_1}_{N_r} \cdot W^{m_2\cdot k_2}_{N_c} } \; }\)

\(\color{white}\therefore\)

\[\begin{split}\color{white}\begin{array}{rcl} b_{n_1,n_2} & = & \displaystyle\sum_{k_1=0}^{N_r-1}\sum_{k_2=0}^{N_c-1} \tilde{u}_{m_1,m_2} \cdot \bbox[white,4px,border:2px solid white]{ \displaystyle\frac{1}{N_r\cdot N_c} \sum_{k_1=0}^{N_r-1}\sum_{k_2=0}^{N_c-1} H_{k_1,k_2} \cdot W^{(m_1-n_1)\cdot k_1}_{N_r} \cdot W^{(m_2-n_2)\cdot k_2}_{N_c} }\\ % & = & \displaystyle\sum_{k_1=0}^{N_r-1}\sum_{k_2=0}^{N_c-1} \tilde{u}_{m_1,m_2} \cdot \tilde{h}_{n_1 - m_1, n_2 - m_2} \\ % & = & h_{n_1,n_2} \circledast x_{n_1,n_2} \end{array} \end{split}\]

Step 3

\[\begin{split} % \begin{array}{rcl} \mathbf{b} & = & \displaystyle{\cal{F}}^{-1}\{ U \odot H \} \\ b_{n_1,n_2}& = & \displaystyle\kern-0.5em\bbox[lightblue,4px]{^{\;}}\kern-0.22em\frac{1}{N_r\cdot N_c} \sum_{k_1=0}^{N_r-1}\sum_{k_2=0}^{N_c-1} H_{k_1,k_2} \cdot \bbox[white,4px,border:2px solid white]{\phantom{U_{k_1,k_2}}\kern-1.625em \bbox[red,4px]{^{\;\;}}\;\;} \cdot W^{\scriptscriptstyle -n_1\cdot k_1}_{N_r} \cdot W^{-n_2\cdot k_2}_{N_c} \qquad \color{darkgray} W^a_b = \mbox{e}^{-j\cdot 2\pi \cdot \frac{a}{b}} \\ % \end{array} \end{split}\]

\(\mbox{Considering } U_{k_1,k_2} = \bbox[white,6px,border:2px solid green]{ \displaystyle \bbox[white,4px,border:2px solid lightblue]{ \sum_{k_1=0}^{N_r-1}\sum_{k_2=0}^{N_c-1} \tilde{u}_{m_1,m_2} } \cdot \bbox[white,4px,border:2px solid red]{ W^{m_1\cdot k_1}_{N_r} \cdot W^{m_2\cdot k_2}_{N_c} } \; }\)

\(\therefore\)

\[\begin{split}\begin{array}{rcl} b_{n_1,n_2} & = & \displaystyle\sum_{k_1=0}^{N_r-1}\sum_{k_2=0}^{N_c-1} \tilde{u}_{m_1,m_2} \cdot \bbox[white,4px,border:2px solid white]{ \displaystyle\frac{1}{N_r\cdot N_c} \sum_{k_1=0}^{N_r-1}\sum_{k_2=0}^{N_c-1} H_{k_1,k_2} \cdot W^{(m_1-n_1)\cdot k_1}_{N_r} \cdot W^{(m_2-n_2)\cdot k_2}_{N_c} }\\ % & \color{white}= & \color{white}\displaystyle\sum_{k_1=0}^{N_r-1}\sum_{k_2=0}^{N_c-1} \tilde{u}_{m_1,m_2} \cdot \tilde{h}_{n_1 - m_1, n_2 - m_2} \\ % & \color{white}= & \color{white}h_{n_1,n_2} \circledast x_{n_1,n_2} \end{array} \end{split}\]

Step 4

\[\begin{split} % \begin{array}{rcl} \mathbf{b} & = & \displaystyle{\cal{F}}^{-1}\{ U \odot H \} \\ b_{n_1,n_2}& = & \displaystyle\frac{1}{N_r\cdot N_c} \sum_{k_1=0}^{N_r-1}\sum_{k_2=0}^{N_c-1} H_{k_1,k_2} \cdot \bbox[white,4px,border:2px solid white]{ U_{k_1,k_2} } \cdot W^{\scriptscriptstyle -n_1\cdot k_1}_{N_r} \cdot W^{-n_2\cdot k_2}_{N_c} \qquad \color{darkgray} W^a_b = \mbox{e}^{-j\cdot 2\pi \cdot \frac{a}{b}} \\ % \end{array} \end{split}\]

\(\mbox{Considering } U_{k_1,k_2} = \bbox[white,6px,border:2px solid white]{ \displaystyle \bbox[white,4px,border:2px solid white]{ \sum_{k_1=0}^{N_r-1}\sum_{k_2=0}^{N_c-1} \tilde{u}_{m_1,m_2} } \cdot \bbox[white,4px,border:2px solid white]{ W^{m_1\cdot k_1}_{N_r} \cdot W^{m_2\cdot k_2}_{N_c} } \; }\)

\(\therefore\)

\[\begin{split}\begin{array}{rcl} b_{n_1,n_2} & = & \displaystyle\sum_{k_1=0}^{N_r-1}\sum_{k_2=0}^{N_c-1} \tilde{u}_{m_1,m_2} \cdot \bbox[white,4px,border:2px solid green]{ \displaystyle\frac{1}{N_r\cdot N_c} \sum_{k_1=0}^{N_r-1}\sum_{k_2=0}^{N_c-1} H_{k_1,k_2} \cdot W^{(m_1-n_1)\cdot k_1}_{N_r} \cdot W^{(m_2-n_2)\cdot k_2}_{N_c} }\\ % & = & \displaystyle\sum_{k_1=0}^{N_r-1}\sum_{k_2=0}^{N_c-1} \tilde{u}_{m_1,m_2} \cdot \tilde{h}_{n_1 - m_1, n_2 - m_2} \\ % & = & h_{n_1,n_2} \circledast x_{n_1,n_2} \end{array} \end{split}\]

2. Additional comments

\(\tilde{\mathbf{h}}\) is the circular extension of the zero-padded version of \(\mathbf{h}\)
To be consintant, \(b_{n_1,n_2}\) should be replaced by \(\tilde{b}_{n_1,n_2}\) to emphasize that it represents a circular extended vector.

Matrix extensions: zero-pad, circurlar and symmetric

The result of convolving image (matrix) \(\mathbf{u}\) with a filter \(\mathbf{h}\) is influenced by the boundary conditions of the former; this is specially important for deconvolution (for instance, see [SerraCapizzano04]).

The collapsible items below show the visual effect of three different boundary conditions:

---------------------------------------------------------------------------
ModuleNotFoundError                       Traceback (most recent call last)
/tmp/ipykernel_356/785914097.py in <module>
----> 1 get_ipython().run_line_magic('load_ext', 'tikzmagic')
      2 
      3 preamble = '''
      4     \definecolor{OliveGreen}{RGB}{60,128,49}
      5     \definecolor{BrickRed}{RGB}{182,50,28}

~/checkouts/readthedocs.org/user_builds/f2o/envs/latest/lib/python3.7/site-packages/IPython/core/interactiveshell.py in run_line_magic(self, magic_name, line, _stack_depth)
   2405                 kwargs['local_ns'] = self.get_local_scope(stack_depth)
   2406             with self.builtin_trap:
-> 2407                 result = fn(*args, **kwargs)
   2408             return result
   2409 

~/checkouts/readthedocs.org/user_builds/f2o/envs/latest/lib/python3.7/site-packages/decorator.py in fun(*args, **kw)
    230             if not kwsyntax:
    231                 args, kw = fix(args, kw, sig)
--> 232             return caller(func, *(extras + args), **kw)
    233     fun.__name__ = func.__name__
    234     fun.__doc__ = func.__doc__

~/checkouts/readthedocs.org/user_builds/f2o/envs/latest/lib/python3.7/site-packages/IPython/core/magic.py in <lambda>(f, *a, **k)
    185     # but it's overkill for just that one bit of state.
    186     def magic_deco(arg):
--> 187         call = lambda f, *a, **k: f(*a, **k)
    188 
    189         if callable(arg):

~/checkouts/readthedocs.org/user_builds/f2o/envs/latest/lib/python3.7/site-packages/IPython/core/magics/extension.py in load_ext(self, module_str)
     31         if not module_str:
     32             raise UsageError('Missing module name.')
---> 33         res = self.shell.extension_manager.load_extension(module_str)
     34 
     35         if res == 'already loaded':

~/checkouts/readthedocs.org/user_builds/f2o/envs/latest/lib/python3.7/site-packages/IPython/core/extensions.py in load_extension(self, module_str)
     78             if module_str not in sys.modules:
     79                 with prepended_to_syspath(self.ipython_extension_dir):
---> 80                     mod = import_module(module_str)
     81                     if mod.__file__.startswith(self.ipython_extension_dir):
     82                         print(("Loading extensions from {dir} is deprecated. "

~/checkouts/readthedocs.org/user_builds/f2o/envs/latest/lib/python3.7/importlib/__init__.py in import_module(name, package)
    125                 break
    126             level += 1
--> 127     return _bootstrap._gcd_import(name[level:], package, level)
    128 
    129 

~/checkouts/readthedocs.org/user_builds/f2o/envs/latest/lib/python3.7/importlib/_bootstrap.py in _gcd_import(name, package, level)

~/checkouts/readthedocs.org/user_builds/f2o/envs/latest/lib/python3.7/importlib/_bootstrap.py in _find_and_load(name, import_)

~/checkouts/readthedocs.org/user_builds/f2o/envs/latest/lib/python3.7/importlib/_bootstrap.py in _find_and_load_unlocked(name, import_)

ModuleNotFoundError: No module named 'tikzmagic'

1. Zero-pad

2. Circular extension

3. Symmetric extension

Spatial and Fourier based convolution: Matching their results

In order to match the computation of the convolution of image \(\mathbf{u} \in \mathbb{R}^{N_r \times N_c}\) and filter \(\mathbf{h} \in \mathbb{R}^{L_r \times L_c}\) in the spatial domain

\[ \mathbf{b} = \mathbf{u} * \mathbf{h}, \]

and in the Fourier domain, the convolution theorem must be taken into account; the Fourier-based computation of the convolution can be broken down into three steps:

1. Pre-processing

The input is extended / zero-padded depending on the boundary conditions; the filter is zero-pad to match the size of the pre-processed image.

Figure 1.a Preprocessing

2. Product in the Fourier domain

This step includes the computation of the forward Fourier transform of the pre-processed versions of the input and filter, and the inverse Fourier tranform of the Hadamard product between them.

\( U_F = {\cal{F}}\{ \mathbf{u}_{ZP} \}\)
\( H_F = {\cal{F}}\{ \mathbf{h}_{ZP} \}\)

\[ \begin{eqnarray*} \mathbf{b}_{BC} & = & {\cal{F}}^{-1}\{ U_F \odot H_F \} \end{eqnarray*} \]

3. Post-processing

In order to match the output of the spatial convolution, the previous result must be cropped out.

Figure 1.b Postprocessing

The collapsible items below include details and simulations associated with three particular boundary conditions usually found when using the convolution operator:

Linear convolution

Figure 2. Linear convolution.

Code snippet

# Pre-processing  -- (np: numpy)
Lr0 = np.int( (Lr - np.remainder(Lr,2))/2 )
Lr1 = Lr - Lr0 - 1

padShp = (u.shape[0]+Lr-1,u.shape[1]+Lc-1)

# Product in the Fourier domain (see additional comments)
UF = rfft2( u, s=padShp, axes=(0,1) )
HF = rfft2( h, s=padShp, axes=(0,1) )

bBC = irfft2( UF*HF, s=padShp, axes=(0,1) )

# Post-processing
b = bBC[Lr0:u.shape[0]+Lr0,Lc0:u.shape[1]+Lc0] 

Additional comments

If pyfftw is to be used, then

import pyfftw

rfft2 = pyfftw.interfaces.numpy_fft.rfft2
irfft2 = pyfftw.interfaces.numpy_fft.irfft2

if JAX is to be used, then

import jax
from jax import numpy as jnp

rfft2 = jnp.fft.rfft2
irfft2 = jnp.fft.irfft2

Stand-alone simulation

from PIL import Image 
import numpy as np
from scipy import signal

import matplotlib.pyplot as PLT
%matplotlib inline

import pyfftw
        
# Imports needed to read images from an URL address
import requests
from io import BytesIO

        
# Test images from the SIPI image database
fname = {0: requests.get('http://sipi.usc.edu/database/misc/5.2.10.tiff'),   # bridge (grayscale)
         1: requests.get('http://sipi.usc.edu/database/misc/boat.512.tiff'), # boat (grayscale)
         2: requests.get('http://sipi.usc.edu/database/misc/4.2.03.tiff'),   # mandrill (color)
        }     

# Read input image      
u = np.asarray( Image.open(BytesIO(fname.get(0).content)) ).astype(float) / 256.0
     
# Generate a random filter     
Lr = 7
Lc = 8
h = np.random.randn(Lr,Lc)
      
# Linear convolution (sptail domain)
bSpt = signal.convolve2d(u, h, boundary='fill', mode='same')

# --- Linear convolution (frequency domain) ---
# Freq. domain: pre-processing
Lr0 = np.int( (Lr - np.remainder(Lr,2))/2 )
Lr1 = Lr - Lr0 - 1
Lc0 = np.int( (Lc - np.remainder(Lc,2))/2 )
Lc1 = Lc - Lc0 - 1

# Freq. domain: product in Fourier
padShp = (u.shape[0]+Lc-1,u.shape[1]+Lr-1)
UF = pyfftw.interfaces.numpy_fft.rfft2( u, s=padShp, axes=(0,1) )
HF = pyfftw.interfaces.numpy_fft.rfft2( h, s=padShp, axes=(0,1) )

bBC = pyfftw.interfaces.numpy_fft.irfft2( UF*HF, s=padShp, axes=(0,1) )

# Freq. domain: post-processing
bFrq = bBC[Lc1:u.shape[0]+Lc1,Lr1:u.shape[1]+Lr1] 

# ---  ---

# Plot results

diff = np.linalg.norm(bSpt-bFrq) 

figure, f = PLT.subplots(ncols=3, figsize=(12, 18)) 
f[0].imshow( u , cmap='gray') 
f[0].set_title('original') 
f[0].set_axis_off()
f[1].imshow( bSpt, cmap='gray') 
f[1].set_title(r'$\mathbf{b}_{spt} = conv(\mathbf{u},\mathbf{h})$') 
f[1].set_axis_off()
f[2].imshow( bSpt-bFrq , cmap='gray') 
txt = r'$\mathbf{b}_{spt} - \mathbf{b}_{frq}$, ($\|\| \cdot \|\|_2 =$'+'{:01.2g})'.format(diff)
f[2].set_title(txt) 
f[2].set_axis_off()

PLT.show()

F2O-based simulation

Convolution with circular B.C.

Figure 3. Convolution with circular B.C.

Code snippet

# Pre-processing  -- (np: numpy)
Lr0 = np.int( (Lr - np.remainder(Lr,2))/2 )
Lr1 = Lr - Lr0 - np.remainder(Lr,2)
Lc0 = np.int( (Lc - np.remainder(Lc,2))/2 )
Lc1 = Lc - Lc0 - np.remainder(Lc,2)

padShp = u.shape

# Product in the Fourier domain (see additional comments)
UF = rfft2( u, s=padShp, axes=(0,1) )
HF = rfft2( h, s=padShp, axes=(0,1) )

bBC = irfft2( UF*HF, s=padShp, axes=(0,1) )

# Post-processing
bFrq = np.roll(bBC, (-Lr0 + 1 - np.remainder(Lr,2),-Lc0 + 1 - np.remainder(Lc,2)), axis=(0,1) )

Additional comments

If pyfftw is to be used, then

import pyfftw

rfft2 = pyfftw.interfaces.numpy_fft.rfft2
irfft2 = pyfftw.interfaces.numpy_fft.irfft2

if JAX is to be used, then

import jax
from jax import numpy as jnp

rfft2 = jnp.fft.rfft2
irfft2 = jnp.fft.irfft2

Stand-alone simulation

from PIL import Image, ImageOps
import numpy as np
from scipy import signal

import matplotlib.pyplot as PLT
%matplotlib inline

import pyfftw
        
# Imports needed to read images from an URL address
import requests
from io import BytesIO

        
# Test images from the SIPI image database
fname = {0: requests.get('http://sipi.usc.edu/database/misc/5.2.10.tiff'),   # bridge (grayscale)
         1: requests.get('http://sipi.usc.edu/database/misc/boat.512.tiff'), # boat (grayscale)
         2: requests.get('http://sipi.usc.edu/database/misc/4.2.03.tiff'),   # mandrill (color)
        }     

# Read input image      
u = np.asarray( ImageOps.grayscale( Image.open(BytesIO(fname.get(2).content)) ) ).astype(float) / 256.0
     
# Generate a random filter     
Lr = 7
Lc = 8
h = np.random.randn(Lr,Lc)
      
# Convolution with circular BC (sptail domain)
bSpt = signal.convolve2d(u, h, boundary='wrap', mode='same')

# --- Convolution with circular BC (frequency domain) ---
# Freq. domain: pre-processing
Lr0 = np.int( (Lr - np.remainder(Lr,2))/2 )
Lr1 = Lr - Lr0 - np.remainder(Lr,2)
Lc0 = np.int( (Lc - np.remainder(Lc,2))/2 )
Lc1 = Lc - Lc0 - np.remainder(Lc,2)


padShp = u.shape

# Freq. domain: product in Fourier
UF = pyfftw.interfaces.numpy_fft.rfft2( u, s=padShp, axes=(0,1) )
HF = pyfftw.interfaces.numpy_fft.rfft2( h, s=padShp, axes=(0,1) )

bBC = pyfftw.interfaces.numpy_fft.irfft2( UF*HF, s=padShp, axes=(0,1) )

# Freq. domain: post-processing
bFrq = np.roll(bBC, (-Lr0 + 1 - np.remainder(Lr,2),-Lc0 + 1 - np.remainder(Lc,2)), axis=(0,1) )

# ---  ---

# Plot results

diff = np.linalg.norm(bSpt-bFrq) 

figure, f = PLT.subplots(ncols=3, figsize=(12, 18)) 
f[0].imshow( u , cmap='gray') 
f[0].set_title('original') 
f[0].set_axis_off()
f[1].imshow( bSpt, cmap='gray') 
f[1].set_title(r'$\mathbf{b}_{spt} = conv(\mathbf{u},\mathbf{h})$') 
f[1].set_axis_off()
f[2].imshow( bSpt-bFrq , cmap='gray') 
txt = r'$\mathbf{b}_{spt} - \mathbf{b}_{frq}$, ($\|\| \cdot \|\|_2 =$'+'{:01.2g})'.format(diff)
f[2].set_title(txt) 
f[2].set_axis_off()

PLT.show()

F2O-based simulation

Convolution with symmetric B.C.

Figure 4. Convolution with symmetric B.C.

Code snippet

# Pre-processing  -- (np: numpy)
Lr0 = np.int( (Lr - np.remainder(Lr,2))/2 )
Lr1 = Lr - Lr0 - np.remainder(Lr,2)
Lc0 = np.int( (Lc - np.remainder(Lc,2))/2 )
Lc1 = Lc - Lc0 - np.remainder(Lc,2)

uBC = np.pad(u, ((Lr0,Lc0),(Lr1,Lc1)), mode='symmetric' )

padShp = (uBC.shape[0]+Lr-np.remainder(Lr,2), uBC.shape[1]+Lc-np.remainder(Lc,2))

# Product in the Fourier domain (see additional comments)
UF = rfft2( u, s=padShp, axes=(0,1) )
HF = rfft2( h, s=padShp, axes=(0,1) )

bBC = irfft2( UF*HF, s=padShp, axes=(0,1) )

# Post-processing
b = bBC[Lr-1:u.shape[0]+Lr-1,Lc-1:u.shape[1]+Lc-1] 

Additional comments

If pyfftw is to be used, then

import pyfftw

rfft2 = pyfftw.interfaces.numpy_fft.rfft2
irfft2 = pyfftw.interfaces.numpy_fft.irfft2

if JAX is to be used, then

import jax
from jax import numpy as jnp

rfft2 = jnp.fft.rfft2
irfft2 = jnp.fft.irfft2

Stand-alone simulation

from PIL import Image 
import numpy as np
from scipy import signal

import matplotlib.pyplot as PLT
%matplotlib inline

import pyfftw
        
# Imports needed to read images from an URL address
import requests
from io import BytesIO

        
# Test images from the SIPI image database
fname = {0: requests.get('http://sipi.usc.edu/database/misc/5.2.10.tiff'),   # bridge (grayscale)
         1: requests.get('http://sipi.usc.edu/database/misc/boat.512.tiff'), # boat (grayscale)
         2: requests.get('http://sipi.usc.edu/database/misc/4.2.03.tiff'),   # mandrill (color)
        }     

# Read input image      
u = np.asarray( Image.open(BytesIO(fname.get(1).content)) ).astype(float) / 256.0
     
# Generate a random filter     
Lr = 7
Lc = 8
h = np.random.randn(Lr,Lc)
      
# Convolution with symmetric BC (sptail domain)
bSpt = signal.convolve2d(u, h, boundary='symm', mode='same')

# --- Convolution with symmetric BC (frequency domain) ---
# Freq. domain: pre-processing
Lr0 = np.int( (Lr - np.remainder(Lr,2))/2 )
Lr1 = Lr - Lr0 - np.remainder(Lr,2)
Lc0 = np.int( (Lc - np.remainder(Lc,2))/2 )
Lc1 = Lc - Lc0 - np.remainder(Lc,2)

uBC = np.pad(u, ((Lr0,Lr0),(Lc1,Lc1)), mode='symmetric' )

padShp = (uBC.shape[0]+Lr-np.remainder(Lr,2), uBC.shape[1]+Lc-np.remainder(Lc,2))

# Freq. domain: product in Fourier
UF = pyfftw.interfaces.numpy_fft.rfft2( uBC, s=padShp, axes=(0,1) )
HF = pyfftw.interfaces.numpy_fft.rfft2( h, s=padShp, axes=(0,1) )

bBC = pyfftw.interfaces.numpy_fft.irfft2( UF*HF, s=padShp, axes=(0,1) )

# Freq. domain: post-processing
bFrq = bBC[Lr-1:u.shape[0]+Lr-1,Lc-1:u.shape[1]+Lc-1]

# ---  ---

# Plot results

diff = np.linalg.norm(bSpt-bFrq) 

figure, f = PLT.subplots(ncols=3, figsize=(12, 18)) 
f[0].imshow( u , cmap='gray') 
f[0].set_title('original') 
f[0].set_axis_off()
f[1].imshow( bSpt, cmap='gray') 
f[1].set_title(r'$\mathbf{b}_{spt} = conv(\mathbf{u},\mathbf{h})$') 
f[1].set_axis_off()
f[2].imshow( bSpt-bFrq , cmap='gray') 
txt = r'$\mathbf{b}_{spt} - \mathbf{b}_{frq}$, ($\|\| \cdot \|\|_2 =$'+'{:01.2g})'.format(diff)
f[2].set_title(txt) 
f[2].set_axis_off()

PLT.show()

F2O-based simulation