Fast discrete convolution in \(\mathbb{R}^2\) with radial kernels using non-uniform fast Fourier transform with nonequispaced frequencies (Q2287853): Difference between revisions

Let \(z_k\in {\mathbb R}^2\) be given points of the unit disk and \(f_k\in \mathbb C\), \(k=1,\ldots,N\). Further, \(\|x\|\) is the Euclidean norm of \(x\in {\mathbb R}^2\). Let \(g:(0,\infty) \to \mathbb R\) be given. Later \(g(t)\) is chosen as \(\log t\), \(t^2 \log t\), and \(t^{-2}\), \(t>0\). In this paper, the author presents a new algorithm, the so-called \textit{efficient Bessel decomposition} (EBD), for the fast evaluation of discrete convolutions with the radial kernel \(g(\|x\|)\) of the form \[ q_k = \sum_{\ell=1}^N g(\|z_k -z_{\ell}\|)\,f_{\ell}\,, \quad k = 1,\ldots,N\,. \] Using a trigonometric approximation of \(g(\|x\|)\), the discrete convolution can be rapidly computed by non-uniform fast Fourier transforms, see [\textit{D. Potts} et al., Numer. Math. 98, No. 2, 329--351 (2004; Zbl 1056.65146)]. The new EBD method approximates \(g(\|x\|)\) by a Bessel series outside the origin by a small number \(P\) of terms \[ g(\|x\|) \approx g(1) + \sum_{n=1}^P \alpha_n\,J_0(\rho_n \|x\|)\,, \quad \delta \le \|x\|\le 1\,, \] where \(J_0\) is the Bessel function of first kind, \(\rho_n\) is the \(n\)th positive zero of \(J_0\), and \(\delta >0\) is a cutoff parameter. Applying the trapezoidal rule to the integral \[ J_0(\rho_n \|x\|) = \frac{1}{2\pi} \int_{\|y\|=1} {\mathrm e}^{{\mathrm i}\rho_n\,x\cdot y}\,{\mathrm d}y\,, \] the author obtains a trigonometric approximation of \(g(\|x\|)\). An important consequence of EBD is the fact that much less terms last for a given accuracy. This allows for a faster evaluation of the discrete convolution at the cost of longer precomputation. A Matlab code of this method is available online. The computational cost and the error of the EBD method are estimated. Numerical results are presented too.