I am trying to use discrete 2d-convolution to estimate continuous double convolution. The convolution integral is
$$g(x,y)=(f\ast h)(x,y) = \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} f(u,v) h(x-u, y-u) \mathrm{d}u \mathrm{d}v$$
Following the approximation techniques for one-dimensional convolution, I assume that $f$ and $h$ are smooth enough so that there exists small enough $t$ such that the surface around $f(it, jt) h((m-i)t, (n-j)t)$ for $i,j,m,n\in\mathbb{Z}$ is "flat". Then we can estimate the convolution integral as follows:
$$\hat{g}(mt,nt) = t^2 \sum_{i=-\infty}^{\infty} \sum_{j=-\infty}^{\infty} f(it, jt) h((m-i)t, (n-j)t)\tag{1}$$
The double-sum is a discrete double convolution (for example, it is implemented by conv2 in MATLAB). However, I ran into trouble when I implemented a toy example. Let
$$\begin{align}f(x,y)&=\exp\left[-\pi(x^2+y^2)\right]\\ h_r(x,y)&=\frac{J_1(2\pi r \sqrt{x^2+y^2})}{\pi\sqrt{x^2+y^2}} \end{align}$$
Thus, $f(x,y)$ is a "standard" 2-d Gaussian function while $h_r(x,y)$ is the 'jinc' function (a.k.a. the "sombrero function"), the 2-d Fourier transform of the circle function parametrized by the radius $r$: $\text{circle}_r(x,y)=\{1 \text{ if }\sqrt{x^2+y^2}\leq r, 0 \text{ otherwise}\}$. The 2-d Fourier transform of $f(x,y)$ is defined as $F(u,v)=\int\int f(x,y)e^{-2\pi i(ux+vy)}\mathrm{d}x\mathrm{d}y$. For this toy example I let radius $r=0.1$.
These particular $f$ and $h$ were selected on purpose. The two-dimensional Fourier transform of $f$ is $f$, and the two-dimensional Fourier transform of $h$ is circle function. Thus, the convolution between the two is the two-dimensional Fourier transform of a two-dimensional standard Gaussian truncated on a circle of radius 0.1, which makes for an easy sanity check. Here is how the "heat map" of $f$ looks like:
MATLAB code to generate this:
[X, Y]=meshgrid(-1:0.025:1,-1:0.025:1);
f=exp(-(X.^2+Y.^2)*pi);
colormap('hot')
imagesc(f)
colorbar
axis square
Some simplifications yield the following formula for the 2-d Fourier transform of the tructated standard Gaussian:
$$g(x,y)=\int_0^r \exp(-\pi u^2) J_0(2\pi u\sqrt{x^2+y^2})u \mathrm{d}u$$
where in this example $r=0.1$ and $J_\nu$ is the $\nu$-order Bessel function of the first kind (also used to define 'jinc' above). Here is the heat map:
This clearly illustrates that the truncation in the (2-d) time domain results in spreading in the (2-d) frequency domain. MATLAB code to generate the picture is:
[X, Y]=meshgrid(-1:0.025:1,-1:0.025:1);
for ii=1:81
for jj=1:81
g(ii,jj)=integral(@(u) exp(-u.^2.*pi).*besselj(0,2*sqrt(X(ii,jj).^2+Y(ii,jj).^2).*pi.*u).*u,0,0.1);
end
end
colormap('hot')
imagesc(g)
colorbar
axis square
Now, I am interested in studying the center of the above heatmap (that is, I am not concerned about the edge effects), so the convolution in (1) should approximate it well. However, when I perform it using MATLAB's built-in function I get nonsense. Here is the code:
[X, Y]=meshgrid(-1:0.025:1,-1:0.025:1);
f=exp(-(X.^2+Y.^2)*pi);
h=besselj(1,2.*pi.*0.1.*sqrt(X.^2+Y.^2))./sqrt(X.^2+Y.^2);
ghat=conv2(f,h);
colormap('hot')
imagesc(ghat)
colorbar
axis square
The black square inside contains NaNs.
Is my methodology incorrect? Am I using discrete convolution incorrectly? Is there a problem with MATLAB (unlikely)? Any help is appreciated. I would like to learn how to do this on a simple example above so that I could try different (and more difficult) $f$ and $h$ later.
UPDATE
I found the reason for the black box of NaNs: I didn't manually remove the discontinuity in the `jinc' function $h_r(x,y)$ at the origin $(x=0, y=0)$ (I also forgot some scaling factors, including $r$ in the numerator of $h_r(x,y)$ and $2\pi$ in $g(x,y)$. When all of these are corrected the figure obtained via approximation to the convolution matches that obtained by numerical integration. So this seems to work.