Consider a 2D function $f(x,y)$ on $\mathbb{R}^2$, which is finite and decays on some finite interval. I don't have a nice analytical/closed-form expression for $f(x,y)$, but can evaluate it at any $(x,y)$, and have thus computed it numerically in a discrete lattice of size $N_x$ and $N_y$, with finite limits $x \in [-a, a]$ and $y \in [-b, b]$. Here, the limits are chosen large enough such that $f(x,y) \to 0$ at the borders.
Question: How can I numerically evaluate the Fourier integral $$F(k_x,k_y) = \int_{-\infty}^{\infty}dx\int_{-\infty}^{\infty}dy e^{-i\mathbf{k}\cdot\mathbf{R}}f(x,y)$$ on a discrete lattice for $k_x$ and $k_y$?
Although these are continuous functions, I guess I could potentially express the answer in terms of discrete FT or FFT, which have many convenient implementations in most computer languages. I think I have to include an additional phase factor, but I don't know how to do this (at least not in 2D). Whatever I try to search for online is hidden by tons and tons of search results which seem irrelevant to me.
My approach so far has been to simply just discretize the integrals (using e.g. trapezoidal method or Simpson method). This is on the right level of complexity that I'm after, but feels crude and potentially error-prone.
The given two-dimensional integration problem can be solved by multiple using of one-dimensional integral Fourier transforms (FT).
To calculate the FT, Filon's integration formula (M.Abramowitz - I.Stegun, p.890) is intended. Taking in account zero boundary conditions and the complexity of the third derivatives estimation, looks convenient using of detalized lattice and the simplified integration expression in the form of
$$\int\limits_{x_0}^{x_{2n}} f(x)\, e^{-itx}\,\text dx \approx h\,\big(\beta(t\,h)F_{2n}(t)+\gamma(t\,h)F_{2n-1}(t)\big),\tag1$$ where $$\beta(\theta) = 2\left(\dfrac{1+\cos^2\theta}{\theta^2} - \dfrac{\sin 2\theta}{\theta^3}\right),\quad \gamma(\theta) = 4\left(\dfrac{\sin \theta}{\theta^3} -\dfrac{\cos\theta}{\theta^2}\right),\tag2$$ $$F_{2n}(t) = \sum\limits_{k=0}^{n-1} f_{2k}\,e^{-2ikt},\quad F_{2n-1}(t) = e^{-it}\sum\limits_{k=0}^{n-1} f_{2k+1}\, e^{-2ikt}.\tag3$$ Since formulas $(3)$ can be calculated via the FFT algorithm, then resulting calculation complexity of the proposed algorithm looks not big.