# Analytical Expression for Convolution of Two 2D DFTs

I am trying to do an image analysis problem on some images, and I want to calculate some things in the frequency domain which are giving me problems. Say I have an (discrete) image $I(x,y)$ of size $N_x\times N_y$ (in practice this size is about $150\times40$ pixels). I want to calculate

$$DFT\left[\frac{\partial I}{\partial x}\cdot\frac{\partial I}{\partial y}\cdot\frac{\partial^2 I}{\partial x \partial y}\right]$$

and express this in terms of the DFT of $I$, which is an $N_x\times N_y$ matrix with entries $\hat{I}_{k,\ell}$. The product is to be interpreted pointwise. I discretize the derivatives with central differences, e.g.

$$\frac{\partial I}{\partial x} = \frac{I(x+1,y) - I(x-1,y)}{2\Delta x}$$ and so I can compute the DFT of various derivatives using the shift properties of DFT. I'm doing this in the frequency domain because other parts of my problem involve 4th derivatives which I think will be way too sensitive to compute in real-space, and I hope the frequency-space calculations will be more stable/robust.

I've attempted to calculate a slightly simpler case, involving only one derivative. My attempted solution is to use the convolution duality properties of the DFT and express multiplication in real space as convolution in frequency space. My attempted expression is below:

$$DFT\left[\left(\frac{\partial I}{\partial x}\right)^2\right]_{k,\ell} = DFT\left[\frac{\partial I}{\partial x}\right]*DFT\left[\frac{\partial I}{\partial x}\right] = \frac{-1}{N_x(\Delta x)^2}\times\sum_{m=0}^{N_x-1}\sin\left(\frac{2\pi m}{N_y}\right)\sin\left(\frac{2\pi (k-m)}{N_y}\right)\hat{I}_{m,\ell}\hat{I}_{(k-m),\ell}$$

However, I'm not sure if I'm properly generalizing circular convolution to more than one dimension, and even if I am, I don't know how to interpret the convolution property with a mix of derivatives (e.g. $DFT[\frac{\partial I}{\partial x}\frac{\partial I}{\partial y}]$). What am I summing to, $N_x-1$ or $N_y-1$? I suspect I might have something wrong in generalizing circular convolution from 1 to 2 dimensions.

In general, how should I interpret circular convolution in 2D? And how can that shed light on calculating this particular DFT?