My question is about quantifying the error that occurs by approximating the continuous Fourier transform on a finite domain by using a discretised version with resolution $N$ for a function of a given Sobolev order $k$. I'm interested in the $L^2$-norm on the full (but finite) discrete grid, which I can prove is at most $\mathcal{O}(N^{-(k-1)})$, but numerically, I see $\mathcal{O}(N^{-k})$. I'm looking for ways to prove this better bound. First let me introduce some notation to formulate the question more clearly.
Let
$$\Omega:=[0,L]^2, \qquad \hat\Omega=\mathbb{Z}^2,$$
as well as the finite discretisations (taking $N\in 2\mathbb{N}$ for convenience)
$$\Omega_\mathrm{fin}:=\left[0:\frac{L}{N}:L\right)^2 \quad \text{and} \quad \hat\Omega_\mathrm{fin}:=\left[-\frac{N}{2}:\frac{N}{2}\right)^2.$$
The notation $[a:b:c]=\{a,a+b,\ldots,c-b,c\}$ is Matlab-inspired (missing $b$ defaults to $1$), with round brackets indicating exclusion of the respective endpoint as for open/closed intervals.
Define
$$\mathcal{F}[f](\hat x, \hat y):= \frac{1}{L^2}\int_\Omega f(x,y) \mathrm{e}^{-2\pi\mathrm{i}(\hat x \frac{x}{L}+\hat y \frac{y}{L})} \mathrm{d} x\,\mathrm{d} y$$ and $$\mathcal{F}_\mathrm{fin}:=\frac{1}{N^2} \sum_{(x,y)\in\Omega_\mathrm{fin}} f(x,y) \mathrm{e}^{-2\pi\mathrm{i}(\hat x \frac{x}{L}+\hat y \frac{y}{L})},$$
as well as the usual Sobolev Space
$$H^k(\Omega):=\{f\in L^2(\Omega) : \int_\Omega \left\lvert\frac{\partial^{k_1+k_2} f}{\partial^{k_1} x \, \partial^{k_2}y}\right\rvert^2\mathrm{d} x\,\mathrm{d} y < \infty, \; k_1+k_2\le k\}.$$
This space can be characterised on the Fourier side as
$$H^k(\hat\Omega):=\{\hat f\in L^2(\hat\Omega) : \langle (\hat x,\hat y) \rangle^k \hat f(\hat x,\hat y)\in L^2(\hat\Omega)\},$$
where $\langle z \rangle:= \sqrt{1+\lvert z\rvert^2}$ is the regularised absolute value.
The norm I'm trying to bound is
$$\left\lVert \mathcal{F}[f]-\mathcal{F}_\mathrm{fin}[f\bigr|_{\Omega_\mathrm{fin}}] \right\rVert_{L^2(\hat\Omega_\mathrm{fin})},$$
which should be $\mathcal{O}(N^{-k})$ for $f\in H^k$, ideally.
I'd also be happy with a 1d-proof or a reference (all the references I found so far deal with fft-errors as opposed to the error of discretising the continuous transform).
This finishes the actual question, below are my attempts/reasoning for the rate I claim to prove, resp. what a proof for the better rate will likely have to show.
I started by analysing the pointwise difference
$$\left\lvert \mathcal{F}[f](\hat x_0, \hat y_0)-\mathcal{F}_\mathrm{fin}[f\bigr|_{\Omega_\mathrm{fin}}](\hat x_0, \hat y_0) \right\rvert.$$
Basically, the sum in $\mathcal{F}_\mathrm{fin}$ is the trapezoidal rule approximation to the integral in $\mathcal{F}$. An analysis in terms of Fourier coefficients seems necessary, because the trapezoidal rule in terms of differentiability only achieves $\mathcal{O}(N^{-2})$, regardless of higher $k$. The following idea is from [J. Waldvogel, Towards a general error theory of the trapezoidal rule; in Approximation and Computation, 2011].
By introducing $h(x,y):= f(x,y) \mathrm{e}^{-2\pi\mathrm{i}(\hat x_0 \frac{x}{L}+\hat y_0 \frac{y}{L})}$ and the trapezoidal quadrature operator
$$T[h]:= \frac{L^2}{N^2} \sum_{(x,y)\in\Omega_\mathrm{fin}} h(x,y),$$
we see that $\mathcal{F}[f](\hat x_0, \hat y_0)=\mathcal{F}[h](0,0)$. By inserting $h=\mathcal{F}^{-1}[\mathcal{F}[h]]$ into T, we see that
$$T[h]=\frac{L^2}{N^2} \sum_{(x,y)\in\Omega_\mathrm{fin}} \sum_{(\hat x,\hat y)\in\hat\Omega} \mathcal{F}[h](\hat x, \hat y)\mathrm{e}^{2\pi\mathrm{i}(\hat x \frac{x}{L}+\hat y \frac{y}{L})}\\
=L^2 \sum_{(\hat x,\hat y)\in\hat\Omega} \mathcal{F}[h](\hat x, \hat y) \delta(\hat x \mathrm{mod} N) \delta(\hat y \mathrm{mod} N)=L^2 \sum_{(\hat x,\hat y)\in\hat\Omega} \mathcal{F}[h](N\hat x, N\hat y),$$
by virtue of the summation property of the roots of unity (after interchanging the sums)
$$ \sum_{\ell=0}^{N-1} \mathrm{e}^{2\pi\mathrm{i} j \frac{\ell}{N}} = N\delta( j \mathrm{mod} N) \quad \forall j\in\mathbb{Z}.$$
In particular, we have that
$$\left\lvert \mathcal{F}[f](\hat x_0, \hat y_0)-\mathcal{F}_\mathrm{fin}[f\bigr|_{\Omega_\mathrm{fin}}(\hat x_0, \hat y_0)] \right\rvert = \biggl\lvert \mathcal{F}[h](0,0)-\sum_{(\hat x,\hat y)\in\hat\Omega} \mathcal{F}[h](N\hat x, N\hat y)\biggr\rvert \\
= \biggl\lvert \sum_{(\hat x,\hat y)\in\hat\Omega\setminus(0,0)} \mathcal{F}[f](N\hat x+\hat x_0, N\hat y+\hat y_0)\biggr\rvert.$$
By inserting the representation $\hat f(\hat x,\hat y)=\langle (\hat x, \hat y) \rangle^{-k}\hat g(\hat x, \hat y)$ for a $\hat g\in L^2(\hat\Omega)$, we can pull out $N^{-k}$ while the sum is still finite (both $\lvert \hat x_0 \rvert$ and $\lvert \hat y_0 \rvert$ are less than $\frac{N}{2}$). Summing this over the $N^2$ points in $\hat\Omega_\mathrm{fin}$ gives $N^{-(k-1)}$ for the $L^2$-norm.
Any further gain (using this avenue of proof) would have to come from bounding
$$\sum_{(\hat x_0,\hat y_0)\in\hat\Omega_\mathrm{fin}} \biggl\lvert \sum_{(\hat x,\hat y)\in\hat\Omega\setminus(0,0)} \langle (\hat x+\frac{\hat x_0}{N}, \hat y+\frac{\hat y_0}{N}) \rangle^{-k}\hat g(N\hat x+\hat x_0, N\hat y+\hat y_0)\biggr\rvert^2$$
with a lower power of $N$ than $2$. There is certainly some decay there, but harnessing it is not obvious. In particular, the knowledge that $\hat g \in L^2(\hat\Omega)$ is not easily applied.
Thanks for reading through this long question!
Best Answer
Also note a paper, but it is in Russian: Kurbatov A.V., Kurbatov V.G. Approximation of (integral) FT via DFT.
http://www.vestnik.vsu.ru/program/view/view.asp?sec=physmath&year=2012&num=02&f_name=2012-02-20