The Discrete Hankel Transform

Abstract

The Hankel transform is an integral transform and is also known as the Fourier-Bessel transform. Until recently, there was no established discrete version of the transform that observed the same sort of relationship to its continuous counterpart as the discrete Fourier transform does to the continuous Fourier transform. Previous definitions of a discrete Hankel transform (DHT) only focused on methods to approximate the integrals of the continuous Hankel integral transform. Recently published work has remedied this and this chapter presents this theory. Specifically, this chapter presents a theory of a DHT that is shown to arise from a discretization scheme based on the theory of Fourier-Bessel expansions. The standard set of shift, modulation, multiplication, and convolution rules are shown. In addition to being a discrete transform in its own right, this DHT can approximate the continuous forward and inverse Hankel transform.

1. Introduction

The Hankel transform has seen applications in many areas of science and engineering. For example, there are applications in propagation of beams and waves, the generation of diffusion profiles and diffraction patterns, imaging and tomographic reconstructions, designs of beams, boundary value problems, etc. The Hankel transform also has a natural relationship to the Fourier transform since the Hankel transform of zeroth order is a 2D Fourier transform of a rotationally symmetric function. Furthermore, the Hankel transform also appears naturally in defining the 2D Fourier transform in polar coordinates and the spherical Hankel transform also appears in the definition of the 3D Fourier transform in spherical polar coordinates [1, 2].

As useful as the Hankel transform may be, there is no discrete Hankel transform (DHT) that exists that has the same relationship to the continuous Hankel transform in the same way that the discrete Fourier transform (DFT) exists alongside the continuous Fourier transform. By this, we mean that the discrete transform exists as a transform in its own right, has its own mathematical theory of the manipulated quantities, and finally as an added bonus, can be used to approximate the continuous version of the transform, with relevant sampling and interpolation theories. Until recently, a discrete Hankel transform merely implied an attempt to discretize the integral(s) of the continuous Hankel transform, rather than an independent discrete transform in its own right.

Such a theory of a DHT was recently proposed [3]. Thus, goal of this chapter is to outline the mathematical theory for the DHT. The mathematical standard set of “DFT-like” rules of shift, modulation, multiplication and convolution rules are derived and presented. A Parseval-like theorem is presented, as are sampling and interpolation theorems. The manner in which this DHT can be used to approximate the continuous Hankel transform is also explained.

2. Hankel transforms and Bessel series

To start, we define the Hankel transform and Fourier-Bessel series as used in this work.

2.1 Hankel transform

The nth-order Hankel transform $F\left(\rho \right)$ of the function $f\left(r\right)$ of a real variable, $r\ge 0$ , is defined by the integral [4]

$$F\left(\rho \right)={\mathbb{H}}_n\left(f\left(r\right)\right)=\underset{0}{\overset{\infty }{\int }}f\left(r\right)J_n\left(\mathit{\rho r}\right)\mathit{rdr},$$ E1

where $J_n\left(z\right)$ is the nth-order Bessel function of the first kind. If n is real and $n>-1/2$ , the transform is self-reciprocating and the inversion formula is given by

$$f\left(r\right)=\underset{0}{\overset{\infty }{\int }}F\left(\rho \right)J_n\left(\mathit{\rho r}\right)\mathit{\rho d\rho }.$$ E2

Thus, Hankel transforms take functions in the spatial r domain and transform them to functions in the spatial frequency $\rho$ domain $f\left(r\right)\iff F\left(\rho \right)$ . The notation $\iff$ is used to indicate a Hankel transform pair.

2.2 Fourier Bessel series

It is known that functions defined on a finite portion of the real line $\left[0,R\right]$ , can be expanded in terms of a Fourier Bessel series [5] given by

$$f\left(r\right)=\sum _{k=1}^{\infty }{f}_k{J}_n\left(\frac{j_{\mathit{nk}}r}{R}\right),$$ E3

where the order, n, of the Bessel function is arbitrary and $j_{\mathit{nk}}$ denotes the kth root of the nth Bessel function $J_n\left(z\right)$ . The Fourier Bessel coefficients $f_k$ of the function $f\left(r\right)$ are given by

$$f_k=\frac{2}{R^2{J}_{n+1}^2\left(j_{\mathit{nk}}\right)}\underset{0}{\overset{R}{\int }}f\left(r\right)J_n\left(\frac{j_{\mathit{nk}}r}{R}\right)r\mathit{dr}.$$ E4

Eqs. (3) and (4) can be considered to be a transform pair where the continuous function $f\left(r\right)$ is forward-transformed to the discrete vector $f_k$ given in (4). The inverse transform is then the operation which returns $f\left(r\right)$ if given $f_k$ , and is given by the summation in Eq. (3). The Fourier Bessel series has the same relationship to the Hankel transform as the Fourier series has to the Fourier transform.

3. Sampling and interpolation theorems for band-limited and space-limited functions

Sampling and interpolation theorems supply the answers to some important questions. For example, given a bandlimited function in frequency space, a sampling theorem answers the question of which samples of the original function are required in order to determine the function completely. The interpolation theorem shows how to interpolate those samples to recover the original function completely. Here, a band-limit means boundedness in frequency. In many applications such as tomography, the notion of a band-limit is not necessarily a property of a function, but rather a limitation of the measurement equipment used to acquire measurements. These measurements are then used to reconstruct some desired function. Thus, the sampling theorem can also answer the question of how band-limits (frequency sensitivities) of measurement equipment determine the resolution of those measurements.

Given a space-limited function, the sampling theorem answers the question of which samples in frequency space determine the function completely, i.e., those that are required to reconstruct the original function. In other words, the sampling theorem dictates which frequency measurements need to be made. As before, the interpolation theorem will give a formula for interpolating those samples to reconstruct the continuous function completely.

3.1 Sampling theorem for a band-limited function

We state here the sampling theorem in the same way that Shannon stated it for functions in time and frequency: if a spatial function $f\left(r\right)$ contains no frequencies higher than $W$ cycles per meter, then it is completely determined by a series of sampling points given by evaluating $f\left(r\right)$ at $r=\frac{j_{\mathit{nk}}}{W_{\rho }}$ where $W_{\rho }=2\mathit{\pi W}$ .

Proof: suppose that a function $f\left(r\right)$ is band-limited in the frequency Hankel domain so that its spectrum $F\left(\rho \right)$ is zero outside of an interval $\left[0,2\mathit{\pi W}\right]$ . The interval is written in this form since $W$ would typically be quoted in units of Hz (cycles per second) if using temporal units or cycles per meter if using spatial units. Therefore, the multiplication by $2\pi$ ensures that the final units are in ${\mathrm{s}}^{-1}$ or ${\mathrm{m}}^{-1}$ . Let us define $W_{\rho }=2\mathit{\pi W}$ . Since the Hankel transform $F\left(\rho \right)$ is defined on a finite portion of the real line $\left[0,W_{\rho }\right]$ , it can be expanded in terms of a Fourier Bessel series as

$$F\left(\rho \right)=\sum _{k=1}^{\infty }{F}_k{J}_n\left(\frac{j_{\mathit{nk}}\rho }{W_{\rho }}\right).$$ E5

where the Fourier Bessel coefficients can be found from Eq. (4) as

$$\begin{array}{l}{F}_k=\frac{2}{W_{\rho }^2{J}_{n+1}^2\left(j_{\mathit{nk}}\right)}{\int }_0^{W_{\rho }}F\left(\rho \right)J_n\left(\frac{j_{\mathit{nk}}\rho }{W_{\rho }}\right)\rho \mathit{d\rho }\\ =\frac{2}{W_{\rho }^2{J}_{n+1}^2\left(j_{\mathit{nk}}\right)}f\left(\frac{j_{\mathit{nk}}}{W_{\rho }}\right)\end{array}$$ E6

In (6), we have used the fact that $f\left(r\right)$ can be written in terms of its inverse Hankel transform, Eq. (2), in combination with the fact that the function is assumed band-limited.

Hence, a function bandlimited to $\left[0,W_{\rho }\right]$ can be written as

$$F\left(\rho \right)=\left\{\begin{array}{ll}\sum _{k=1}^{\infty }\frac{2}{W_{\rho }^2{J}_{n+1}^2\left(j_{\mathit{nk}}\right)}f\left(\frac{j_{\mathit{nk}}}{W_{\rho }}\right)J_n\left(\frac{j_{\mathit{nk}}\rho }{W_{\rho }}\right)& \rho <W_{\rho }\\ 0& \rho \ge W_{\rho }\end{array}\right.$$ E7

Eq. (7) states that the samples $f\left(\frac{j_{\mathit{nk}}}{W_{\rho }}\right)$ determine the function $f\left(r\right)$ completely since (i) $F\left(\rho \right)$ is determined by Eq. (7), and (ii) $f\left(r\right)$ is known if $F\left(\rho \right)$ is known. Another way of looking at this is that band-limiting a function to $\left[0,W_{\rho }\right]$ results in information about the original function at samples $r_{\mathit{nk}}=\frac{j_{\mathit{nk}}}{W_{\rho }}$ . So, Eq. (7) is the statement of the sampling theorem.

To verify that this sampling theorem is consistent with expectations, we recognize that the zeros of $J_n\left(z\right)$ are almost evenly spaced at intervals of $\pi$ and that the spacing becomes exactly $\pi$ in the limit as $z\to \infty$ . To determine the (bandlimited) function $f\left(r\right)$ completely, we need to sample the function at $f\left(\frac{j_{\mathit{nk}}}{W_{\rho }}\right)=f\left(\frac{j_{\mathit{nk}}}{2\mathit{\pi W}}\right)$ and these samples are (eventually) multiples of $\raisebox{1ex}{$\pi $}\!\left/ \!\raisebox{-1ex}{$(2\pi \mathit{W})$}\right.=\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$\left(2W\right)$}\right.$ apart, which is consistent with the standard Shannon sampling theorem which requires samples at multiples of $\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$\left(2W\right)$}\right.$ [6].

3.2 Interpolation theorem for a band-limited function

It follows from Eq. (7) that $f\left(r\right)$ can be found by inverse Hankel transformation to give

$$\begin{array}{ll}f\left(r\right)& =\underset{0}{\overset{W_{\rho }}{\int }}\left\{\sum _{k=1}^{\infty }\frac{2}{W_{\rho }^2{J}_{n+1}^2\left(j_{\mathit{nk}}\right)}f\left(\frac{j_{\mathit{nk}}}{W_{\rho }}\right)J_n\left(\frac{j_{\mathit{nk}}\rho }{W_{\rho }}\right)\right\}{J}_n\left(\mathit{\rho r}\right)\mathit{\rho d\rho }\\ & =\sum _{k=1}^{\infty }\frac{2}{W_{\rho }^2{J}_{n+1}^2\left(j_{\mathit{nk}}\right)}f\left(\frac{j_{\mathit{nk}}}{W_{\rho }}\right)\underset{0}{\overset{W_{\rho }}{\int }}{J}_n\left(\frac{j_{\mathit{nk}}\rho }{W_{\rho }}\right)J_n\left(\mathit{\rho r}\right)\mathit{\rho d\rho }.\end{array}$$ E8

From Watson ([7], p. 134), we have the following result

$$\int J_n\left(\mathit{\alpha z}\right){\overline{J}}_n\left(\mathit{\beta z}\right)\mathit{zdz}=\frac{z\left\{\alpha J_{n+1}\left(\mathit{\alpha z}\right){\overline{J}}_n\left(\mathit{\beta z}\right)-\beta J_n\left(\mathit{\alpha z}\right){\overline{J}}_{n+1}\left(\mathit{\beta z}\right)\right\}}{{\alpha }^2-{\beta }^2}$$ E9

Eq. (9) can be used to simplify (8) to give

$$f\left(r\right)=\sum _{k=1}^{\infty }f\left(\frac{j_{\mathit{nk}}}{W_{\rho }}\right)\frac{2j_{\mathit{nk}}}{J_{n+1}\left(j_{\mathit{nk}}\right)}\frac{J_n\left({\mathit{rW}}_{\rho }\right)}{j_{\mathit{nk}}^2-r^2{W}_{\rho }^2}$$ E10

Eq. (10) gives the formula for interpolating the samples $f\left(\frac{j_{\mathit{nk}}}{W_{\rho }}\right)$ to reconstruct the continuous band-limited function $f\left(r\right)$ . Each term used to reconstruct the original function $f\left(r\right)$ consists of the samples of the function $f\left(r\right)$ at $r=\left(\frac{j_{\mathit{nk}}}{W_{\rho }}\right)$ multiplied by a reconstructing function of the form

$$\frac{2j_{\mathit{nk}}}{J_{n+1}\left(j_{\mathit{nk}}\right)}\frac{J_n\left({\mathit{rW}}_{\rho }\right)}{j_{\mathit{nk}}^2-r^2{W}_{\rho }^2}.$$ E11

3.3 Interpretation in terms of a jinc

Eq. (8) states

$$f\left(r\right)=\sum _{k=1}^{\infty }\frac{2}{W_{\rho }^2{J}_{n+1}^2\left(j_{\mathit{nk}}\right)}f\left(\frac{j_{\mathit{nk}}}{W_{\rho }}\right)\underset{0}{\overset{W_{\rho }}{\int }}{J}_n\left(\frac{j_{\mathit{nk}}\rho }{W_{\rho }}\right)J_n\left(\mathit{\rho r}\right)\mathit{\rho d\rho }$$ E12

In other research work [8], the generalized shift operator $R^{r_0}$ indicating a shift of $r_0$ acting on the function $f\left(r\right)$ has been defined by the formula

$$R^{r_0}f\left(r\right)=\underset{0}{\overset{\infty }{\int }}F\left(\rho \right)J_n\left(\rho r_0\right)J_n\left(\mathit{\rho r}\right)\mathit{\rho d\rho }.$$ E13

With this definition of a generalized shift operator, we recognize the integral in Eq. (12) as the inverse Hankel transform of the Boxcar function shifted by $\left(\frac{j_{\mathit{nk}}}{W_{\rho }}\right)$ . More explicitly,

$$\begin{array}{ll}\underset{0}{\overset{W_{\rho }}{\int }}{J}_n\left(\frac{j_{\mathit{nk}}\rho }{W_{\rho }}\right)J_n\left(\mathit{\rho r}\right)\mathit{\rho d\rho }& =\underset{\begin{array}{l}\text{generalized}\\ \text{shift of}\frac{j_{\mathit{nk}}}{W_{\rho }}\end{array}}{\underbrace{R^{\frac{j_{\mathit{nk}}}{W_{\rho }}}}} \underset{\begin{array}{l}\text{inverse Hankel transform of}\\ {\mathrm{\Pi }}_{W_{\rho }}\left(\rho \right)\end{array}}{\underbrace{\left\{\underset{0}{\overset{\infty }{\int }}{\mathrm{\Pi }}_{W_{\rho }}\left(\rho \right)J_n\left(\mathit{\rho r}\right)\mathit{\rho d\rho }\right\}}}\\ & =\frac{j_{\mathit{nk}}{W}_{\rho }^2}{j_{\mathit{nk}}^2-{\left({\mathit{rW}}_{\rho }\right)}^2}{J}_{n+1}\left(j_{\mathit{nk}}\right)J_n\left({\mathit{rW}}_{\rho }\right)\end{array}$$ E14

where

$${\mathrm{\Pi }}_{W_{\rho }}\left(\rho \right)=\left\{\begin{array}{ll}1& 0\le \rho \le W_{\rho }\\ 0& \text{otherwise}\end{array}\right.$$ E15

The boxcar function is a generalized version of the standard Rect function. The Rect function is usually defined as the function which has value 1 over the interval $\left[-1/2,-1/2\right]$ and is zero otherwise. Now this is interesting specifically because of the interpretation of Eq. (14). Had we been working in the standard Fourier domain instead of the Hankel domain, the Boxcar function would be replaced with the Rect function and the Hankel transform would be replaced with a standard Fourier transform. Proceeding with this line of thinking, the inverse Fourier transform of the Rect function would be a sinc function, which is the standard interpolation function of the classical Shannon interpolation formula. Hence, the Fourier equivalent interpretation of Eq. (14) is a shifted sinc function. Thus, the formulation above follows exactly the standard Shannon Interpolation formula (see the original publication [9], or the classic paper reprint [6]).

For the relatively simple case of the zeroth-order Hankel transform, the inverse Hankel transform of the Boxcar function is given by

$$\begin{array}{ll}\underset{0}{\overset{\infty }{\int }}{\mathrm{\Pi }}_{W_{\rho }}\left(\rho \right)J_0\left(\mathit{\rho r}\right)\mathit{\rho d\rho }& =\underset{0}{\overset{W_{\rho }}{\int }}{J}_0\left(\mathit{\rho r}\right)\mathit{\rho d\rho }\\ & =\frac{W_{\rho }}{r}{J}_1\left(W_{\rho }r\right)=W_{\rho }^2\frac{J_1\left(W_{\rho }r\right)}{W_{\rho }r}.\end{array}$$ E16

The function $\raisebox{1ex}{$2J_1\left(r\right)$}\!\left/ \!\raisebox{-1ex}{$r$}\right.$ is often termed the jinc or sombrero function and is the polar coordinate analog of the sinc function, so Eq. (16) is a scaled version of a jinc function.

In fact, from Eqs. (13), (14) and (16), it follows that the generalized shifted version of the jinc function is given by

$$R^{\frac{j_{0k}}{W_{\rho }}}\left\{\frac{2J_1\left(W_{\rho }r\right)}{W_{\rho }r}\right\}=\frac{2j_{0k}{J}_1\left(j_{0k}\right)}{j_{0k}^2-{\left({\mathit{rW}}_{\rho }\right)}^2}{J}_0\left({\mathit{rW}}_{\rho }\right).$$ E17

Hence, for a zeroth-order Fourier Bessel transform, Eq. (12), the expansion for $f\left(r\right)$ reads

$$f\left(r\right)=\sum _{k=1}^{\infty }f\left(\frac{j_{0k}}{W_{\rho }}\right)\frac{1}{J_1^2\left(j_{0k}\right)}\underset{=R^{\frac{j_{0k}}{W_{\rho }}}\left\{\frac{2J_1\left(W_{\rho }r\right)}{W_{\rho }r}\right\}}{\underbrace{\frac{2j_{0k}{J}_1\left(j_{0k}\right)J_0\left({\mathit{rW}}_{\rho }\right)}{j_{0k}^2-{\left({\mathit{rW}}_{\rho }\right)}^2}}}$$ E18

Eq. (18) says that the interpolating function is a shifted jinc function, in analogy with a scaled sinc being the interpolating function for the sampling theorem used for Fourier transforms.

3.4 Sampling theorem for a space-limited function

Now consider a space-limited function $f\left(r\right)$ so that $f\left(r\right)$ is zero outside of an interval $\left[0,R\right]$ . It then follows that it can be expanded on $\left[0,R\right]$ in terms of a Fourier Bessel series so that

$$f\left(r\right)=\sum _{k=1}^{\infty }{f}_k{J}_n\left(\frac{j_{\mathit{nk}}r}{R}\right),$$ E19

where the Fourier Bessel coefficients can be found from

$$f_k=\frac{2}{R^2{J}_{n+1}^2\left(j_{\mathit{nk}}\right)}\underset{0}{\overset{R}{\int }}f\left(r\right)J_n\left(\frac{j_{\mathit{nk}}r}{R}\right)r\mathit{dr}=\frac{2}{R^2{J}_{n+1}^2\left(j_{\mathit{nk}}\right)}F\left(\frac{j_{\mathit{nk}}}{R}\right).$$ E20

Here, we have used the definition of the Hankel transform $F\left(\rho \right)$ , Eq. (1), in the right hand side of Eq. (20). Hence, the function can be written as

$$f\left(r\right)=\left\{\begin{array}{ll}\sum _{k=1}^{\infty }\frac{2}{R^2{J}_{n+1}^2\left(j_{\mathit{nk}}\right)}F\left(\frac{j_{\mathit{nk}}}{R}\right)J_n\left(\frac{j_{\mathit{nk}}r}{R}\right)& r<R\\ 0& r\ge R\end{array}\right.$$ E21

From Eq. (21), it is evident that the samples $F\left(\frac{j_{\mathit{nk}}}{R}\right)$ determine the function $f\left(r\right)$ and hence its transform $F\left(\rho \right)$ completely. Another way of looking at this is that space limiting a function to $\left[0,R\right]$ implies discretization in spatial frequency space, at frequencies ${\rho }_{\mathit{nk}}=\frac{j_{\mathit{nk}}}{R}$ .

3.5 Interpolation theorem for a space-limited function

The Hankel transform of the function can then be found from a forward Hankel transformation as

$$\begin{array}{ll}F\left(\rho \right)& =\underset{0}{\overset{\infty }{\int }}f\left(r\right)J_n\left(\mathit{\rho r}\right)\mathit{rdr}\\ & =\sum _{k=1}^{\infty }\frac{2}{R^2{J}_{n+1}^2\left(j_{\mathit{nk}}\right)}F\left(\frac{j_{\mathit{nk}}}{R}\right)\underset{0}{\overset{R}{\int }}{J}_n\left(\frac{j_{\mathit{nk}}r}{R}\right)J_n\left(\mathit{\rho r}\right)\mathit{rdr}\end{array}$$ E22

Using Eq. (9), Eq. (22) can be simplified to give

$$F\left(\rho \right)=\sum _{k=1}^{\infty }F\left(\frac{j_{\mathit{nk}}}{R}\right)\frac{2j_{\mathit{nk}}}{J_{n+1}\left(j_{\mathit{nk}}\right)}\frac{J_n\left(\mathit{\rho R}\right)}{j_{\mathit{nk}}^2-{\left(\mathit{\rho R}\right)}^2}$$ E23

Eq. (23) gives the formula for interpolating the samples $F\left(\frac{j_{\mathit{nk}}}{R}\right)$ to give the continuous function $F\left(\rho \right)$ .

4. Intuitive discretization scheme for the Hankel transform

Based on the sampling theorems above, in this section we explore how assuming that a function can be simultaneously band-limited and space-limited will naturally lead to a discrete Hankel transform. Although it is known that it is not possible to fulfill both of these conditions exactly, it is possible to keep the spectrum within a given frequency band, and to have the space function very small outside some specified spatial interval (or vice-versa). Hence, it is possible for functions to be “effectively” space and band-limited.

4.1 Forward transform

We demonstrated above that a band-limited function, with $\rho <W_{\rho }=2\mathit{\pi W}$ can be written as

$$F\left(\rho \right)=\sum _{k=1}^{\infty }\frac{2}{W_{\rho }^2{J}_{n+1}^2\left(j_{\mathit{nk}}\right)}f\left(\frac{j_{\mathit{nk}}}{W_{\rho }}\right)J_n\left(\frac{j_{\mathit{nk}}\rho }{W_{\rho }}\right).$$ E24

Evaluating the previous Eq. (24) at the sampling points ${\rho }_{\mathit{nm}}=\frac{j_{\mathit{nm}}{W}_{\rho }}{j_{\mathit{nN}}}$ (for any integer N) gives for $m<N$

$$F\left(\frac{j_{\mathit{nm}}{W}_{\rho }}{j_{\mathit{nN}}}\right)=\sum _{k=1}^{\infty }\frac{2}{W_{\rho }^2{J}_{n+1}^2\left(j_{\mathit{nk}}\right)}f\left(\frac{j_{\mathit{nk}}}{W_{\rho }}\right)J_n\left(\frac{j_{\mathit{nk}}}{W_{\rho }}\frac{j_{\mathit{nm}}{W}_{\rho }}{j_{\mathit{nN}}}\right) m<N.$$ E25

For $m<N$ , then ${\rho }_{\mathit{nm}}=\frac{j_{\mathit{nm}}{W}_{\rho }}{j_{\mathit{nN}}}<W_{\rho }$ , and Eq. (25), summing over infinite k, is exact. For $m\ge N$ , then ${\rho }_{\mathit{nm}}=\frac{j_{\mathit{nm}}{W}_{\rho }}{j_{\mathit{nN}}}\ge W_{\rho }$ and by the assumption of the bandlimited nature of the function, $F\left({\rho }_{\mathit{nm}}\right)=0$ .

Now, suppose that in addition to being band-limited, the function is also effectively space limited. As mentioned above, it is known that a function cannot be finite in both space and spatial frequency (equivalently it cannot be finite in both time and frequency if using the Fourier transform). However, if a function is effectively space limited, this means that there exists an integer N for which $f\left(\frac{j_{\mathit{nk}}}{W_{\rho }}\right)\approx 0$ for $k>N$ . In other words, we can find an interval beyond which the space function is very small. In fact, since the Fourier-Bessel series in (24) is known to converge, it is known that $\underset{k\to \infty }{\text{lim}}f\left(\frac{j_{\mathit{nk}}}{W_{\rho }}\right)=0$ , which means that for any arbitrarily small $\epsilon$ , there must exist an integer N for which $f\left(\frac{j_{\mathit{nk}}}{W_{\rho }}\right)<\epsilon$ for $k>N$ .

Hence, using the argument of “effectively space limited” in the preceding paragraph, we can terminate the series in Eq. (25) at a suitably chosen $N$ that ensures the effective space limit. Terminating the series at $k=N$ is the same as assuming that $f\left(r\right)\approx 0$ for $r>R=\frac{j_{\mathit{nN}}}{W_{\rho }}$ . Noting that at $k=N$ , the last term in (25) is $J_n\left(\frac{j_{\mathit{nN}}{j}_{\mathit{nm}}}{j_{\mathit{nN}}}\right)=J_n\left(j_{\mathit{nm}}\right)=0$ , then after terminating at N, Eq. (25) becomes for $m=1..N-1$

$$F\left(\frac{j_{\mathit{nm}}{W}_{\rho }}{j_{\mathit{nN}}}\right)=\sum _{k=1}^{N-1}\frac{2}{W_{\rho }^2{J}_{n+1}^2\left(j_{\mathit{nk}}\right)}f\left(\frac{j_{\mathit{nk}}}{W_{\rho }}\right)J_n\left(\frac{j_{\mathit{nk}}{j}_{\mathit{nm}}}{j_{\mathit{nN}}}\right).$$ E26

In this case, the truncated sum in Eq. (26) does not represent $F\left({\rho }_{\mathit{nm}}\right)$ exactly due to the truncation at N terms, but should provide a reasonably good approximation since the Fourier-Bessel series is known to converge and we are assuming that we have terminated the series at the point where additional $f\left(\frac{j_{\mathit{nk}}}{W_{\rho }}\right)$ terms do not contribute significantly.

4.2 Inverse transform

Concomitantly, we know that for any space-limited function then for $r<R$ , we can write

$$f\left(r\right)=\sum _{m=1}^{\infty }\frac{2}{R^2{J}_{n+1}^2\left(j_{\mathit{nm}}\right)}F\left(\frac{j_{\mathit{nm}}}{R}\right)J_n\left(\frac{j_{\mathit{nm}}r}{R}\right).$$ E27

More specifically, suppose that we follow the logic from the previous section that the function $f\left(r\right)$ that was bandlimited but also “effectively space-limited” due the truncation of the series in Eq. (25) at N. In that case then $R=\frac{j_{\mathit{nN}}}{W_{\rho }}$ and the band-limit implies that $F\left(\rho \right)=0$ for $\rho >W_{\rho }$ . Following this logic and using $R=\frac{j_{\mathit{nN}}}{W_{\rho }}$ , then Eq. (27) becomes

$$f\left(r\right)=\sum _{m=1}^{N-1}\frac{2W_{\rho }^2}{j_{\mathit{nN}}^2{J}_{n+1}^2\left(j_{\mathit{nm}}\right)}F\left(\frac{j_{\mathit{nm}}{W}_{\rho }}{j_{\mathit{nN}}}\right)J_n\left(\frac{j_{\mathit{nm}}{W}_{\rho }}{j_{\mathit{nN}}}r\right)$$ E28

where we truncated the series in Eq. (28) at N by using the fact that $F\left(\rho \right)=0$ for $\rho \ge W_{\rho }$ to deduce that $F\left(\frac{j_{\mathit{nm}}{W}_{\rho }}{j_{\mathit{nN}}}\right)=0$ for $m\ge N$ . Evaluating (28) at $r_{\mathit{nk}}=\frac{j_{\mathit{nk}}R}{j_{\mathit{nN}}}=\frac{j_{\mathit{nk}}}{W_{\rho }}$ gives for $k=1..N-1$

$$f\left(\frac{j_{\mathit{nk}}}{W_{\rho }}\right)=\sum _{m=1}^{N-1} \frac{2W_{\rho }^2}{j_{\mathit{nN}}^2{J}_{n+1}^2\left(j_{\mathit{nm}}\right)}F\left(\frac{j_{\mathit{nm}}{W}_{\rho }}{j_{\mathit{nN}}}\right)J_n\left(\frac{j_{\mathit{nm}}{j}_{\mathit{nk}}}{j_{\mathit{nN}}}\right).$$ E29

Compare Eq. (29) to the “forward” transform from Eq. (26), repeated here for convenience, where we found that for $m=1..N-1$

$$F\left(\frac{j_{\mathit{nm}}{W}_{\rho }}{j_{\mathit{nN}}}\right)=\sum _{k=1}^{N-1}\frac{2}{W_{\rho }^2{J}_{n+1}^2\left(j_{\mathit{nk}}\right)}f\left(\frac{j_{\mathit{nk}}}{W_{\rho }}\right)J_n\left(\frac{j_{\mathit{nk}}{j}_{\mathit{nm}}}{j_{\mathit{nN}}}\right).$$ E30

Eqs. (29) and (30) are the fundamental relations used for the discrete Hankel transform proposed in the following sections.

4.3 Intuitive discretization scheme and kernel

The preceding development shows that a “natural,” N-dimensional discretization scheme in finite space $\left[0,R\right]$ and finite frequency space $\left[0,W_{\rho }\right]$ is given by

$$r_{\mathit{nk}}=\frac{j_{\mathit{nk}}}{W_{\rho }}=\frac{j_{\mathit{nk}}R}{j_{\mathit{nN}}}{\rho }_{\mathit{nk}}=\frac{j_{\mathit{nk}}}{R}=\frac{j_{\mathit{nk}}{W}_{\rho }}{j_{\mathit{nN}}}k=1\dots N-1.$$ E31

The relationship $W_{\rho }R=j_{\mathit{nN}}$ can be used to change from finite frequency domain to a finite space domain and vice-versa. The size of the transform N, can be determined from $W_{\rho }R=j_{\mathit{nN}}$ .

It is pointed out in [10] that the zeros of $J_n\left(z\right)$ are almost evenly spaced at intervals of $\pi$ and that the spacing becomes exactly $\pi$ in the limit as $z\to \infty$ . In fact, it is shown in [10] that a simple asymptotic form for the Bessel function is given by

$$J_n\left(z\right)\approx \sqrt{\frac{2}{\mathit{\pi z}}}\text{cos}\left[z-\left(n+\frac{1}{2}\right)\frac{\pi }{2}\right]$$ E32

Eq. (32) becomes a better approximation to $J_n\left(z\right)$ as $z\to \infty$ . The zeros of the cosine function are at odd multiples of $\pi /2$ . Therefore, an approximation to the Bessel zero, $j_{\mathit{nk}}$ is given by

$$j_{\mathit{nk}}\approx \left(2k+n-\frac{1}{2}\right)\frac{\pi }{2}.$$ E33

Using this approximation, then $W_{\rho }R=j_{\mathit{nN}}$ becomes

$$2\mathit{\pi WR}=j_{\mathit{nN}}\approx \left(2N+n-\frac{1}{2}\right)\frac{\pi }{2}$$ E34

For larger values of N as would typically be used in a discretization scheme, then we can write approximately

$$2\mathit{WR}\approx \left(N+\frac{n}{2}\right)$$ E35

This is exactly analogous to the corresponding expression for Fourier transforms. Specifically, for temporal Fourier transforms Shannon [6] argued that “If the function is limited to the time interval T and the samples are spaced 1/(2 W) seconds apart (where W is the bandwidth), there will be a total of 2WT samples in the interval. All samples outside will be substantially zero. To be more precise, we can define a function to be limited to the time interval T if, and only if, all the samples outside this interval are exactly zero. Then we can say that any function limited to the bandwidth W and the time interval T can be specified by giving $N=2\mathit{WT}$ numbers”. Following this line of thinking, Eq. (35) states that for an nth-order Hankel transform, any function limited to the bandwidth W and the space interval R can be specified by giving $N=\left(2\mathit{WR}-n/2\right)$ numbers and it will certainly be true that specifying $N=2\mathit{WR}$ numbers will specify the function, in exact analogy to Shannon’s result.

4.4 Proposed kernel for the discrete transform

The preceding sections show that both forward and inverse discrete versions of the transforms contain an expression of the form

$$\frac{2}{J_{n+1}^2\left(j_{\mathit{nk}}\right)}{J}_n\left(\frac{j_{\mathit{nk}}{j}_{\mathit{nm}}}{j_{\mathit{nN}}}\right).$$ E36

This leads to a natural choice of kernel for the discrete transform, as shall be outlined in the next section. To aid in in the choice of kernel for the discrete transform, we present a useful discrete orthogonality relationship shown in [11] that for $1\le m,i\le N-1$

$$\sum _{k=1}^{N-1}\frac{4J_n\left(\frac{j_{\mathit{nm}}{j}_{\mathit{nk}}}{j_{\mathit{nN}}}\right)J_n\left(\frac{j_{\mathit{nk}}{j}_{\mathit{ni}}}{j_{\mathit{nN}}}\right)}{J_{n+1}^2\left(j_{\mathit{nk}}\right)}=j_{\mathit{nN}}^2{J}_{n+1}^2\left(j_{\mathit{nm}}\right){\delta }_{\mathit{mi}}$$ E37

where $j_{\mathit{nm}}$ represents the mth zero of the nth-order Bessel function $J_n\left(x\right)$ , and ${\delta }_{\mathit{mi}}$ is the Kronecker delta function, defined as

$${\delta }_{\mathit{mn}}=\left\{\begin{array}{ll}1& \text{if}m=n\\ 0& \text{otherwise}\end{array}\right..$$ E38

If written in matrix notation, then the Kronecker delta of Eq. (38) is the identity matrix.

Fisk-Johnson discusses the analytical derivation of Eq. (37) in the appendix of [11]. Eq. (37) is exactly true in the limit as $N\to \infty$ and is true for $N>30$ within the limits of computational error $\left(\pm {10}^{-7}\right)$ . For smaller values of $N$ , Eq. (37) holds with the worst case for the smallest value of N giving $\pm {10}^{-3}$ .

5. Transformation matrices

5.1 Transformation matrix

With inspiration from the notation in [11], and an additional scaling factor of $1/j_{\mathit{nN}}$ , we define an $\left(N-1\right)\times \left(N-1\right)$ transformation matrix with the (m,k)th entry given by

$$Y_{m,k}^{\mathit{nN}}=\frac{2}{j_{\mathit{nN}}{J}_{n+1}^2\left(j_{\mathit{nk}}\right)}{J}_n\left(\frac{j_{\mathit{nm}}{j}_{\mathit{nk}}}{j_{\mathit{nN}}}\right)1\le m,k\le N-1.$$ E39

In Eq. (39), the superscripts n and N refer to the order of the Bessel function and the dimension of the space that are being considered, respectively. The subscripts m and k refer to the (m,k)th entry of the transformation matrix.

Furthermore, the orthogonality relationship, Eq. (37), states that

$${\sum }_{k=1}^{N-1}{Y}_{i,k}^{\mathit{nN}}{Y}_{k,m}^{\mathit{nN}}={\sum }_{k=1}^{N-1}4\frac{J_n\left(\raisebox{1ex}{$j_{\mathit{ni}}{j}_{\mathit{nk}}$}\!\left/ \!\raisebox{-1ex}{$j_{\mathit{nN}}$}\right.\right)J_n\left(\raisebox{1ex}{$j_{\mathit{nk}}{j}_{\mathit{nm}}$}\!\left/ \!\raisebox{-1ex}{$j_{\mathit{nN}}$}\right.\right)}{j_{\mathit{nN}}^2{J}_{n+1}^2\left(j_{\mathit{nk}}\right)J_{n+1}^2\left(j_{\mathit{nm}}\right)}={\delta }_{\mathit{im}}.$$ E40

Eq. (40) states that the rows and columns of the matrix $Y_{m,k}^{\mathit{nN}}$ are orthonormal and can be written in matrix form as

$$Y^{\mathit{nN}}{Y}^{\mathit{nN}}=\mathbf{{\rm I}},$$ E41

where $\mathbf{I}$ is the $N-1$ dimensional identity matrix and we have written the $N-1$ square matrix $Y_{m,k}^{\mathit{nN}}$ as $Y^{\mathit{nN}}$ . The forward and inverse truncated and discretized transforms given in Eqs. (26) and (29) can be expressed in terms of $Y^{\mathit{nN}}$ . The forward transform, Eq. (26), can be written as

$$F\left({\rho }_{\mathit{nm}}\right)=\frac{j_{\mathit{nN}}}{W_{\rho }^2}\sum _{k=1}^{N-1}{Y}_{m,k}^{\mathit{nN}}f\left(r_{\mathit{nk}}\right).$$ E42

Similarly, the inverse transform, Eq. (29), can be written as

$$f\left(r_{\mathit{nk}}\right)=\frac{W_{\rho }^2}{j_{\mathit{nN}}}\sum _{m=1}^{N-1}{Y}_{k,m}^{\mathit{nN}}F\left({\rho }_{\mathit{nm}}\right).$$ E43

5.2 Another choice of transformation matrix

Following the notation in [12], we can also define a different $\left(N-1\right)\times \left(N-1\right)$ transformation matrix with the (m,k)th entry given by

$$T_{m,k}^{\mathit{nN}}=2\frac{J_n\left(\raisebox{1ex}{$j_{\mathit{nm}}{j}_{\mathit{nk}}$}\!\left/ \!\raisebox{-1ex}{$j_{\mathit{nN}}$}\right.\right)}{J_{n+1}\left(j_{\mathit{nm}}\right)J_{n+1}\left(j_{\mathit{nk}}\right)j_{\mathit{nN}}}1\le m,k\le N-1.$$ E44

In Eq. (44), the superscripts n and N refer to the order of the Bessel function and the dimension of the space that are being considered, respectively. The subscripts m and k refer to the (m,k)th entry of the matrix. From (39), it can be seen that $T_{m,k}^{\mathit{nN}}=T_{k,m}^{\mathit{nN}}$ so that $T^{\mathit{nN}}$ is a real, symmetric matrix. The relationship between the $T_{m,k}^{\mathit{nN}}$ and $Y_{m,k}^{\mathit{nN}}$ matrices is given by

$$T_{m,k}^{\mathit{nN}}\frac{J_{n+1}\left(j_{\mathit{nm}}\right)}{J_{n+1}\left(j_{\mathit{nk}}\right)}=Y_{m,k}^{\mathit{nN}}.$$ E45

The orthogonality relationship, Eq. (37), can be written as

$$\sum _{k=1}^{N-1}4\frac{J_n\left(\raisebox{1ex}{$j_{\mathit{nm}}{j}_{\mathit{nk}}$}\!\left/ \!\raisebox{-1ex}{$j_{\mathit{nN}}$}\right.\right)J_n\left(\raisebox{1ex}{$j_{\mathit{nk}}{j}_{\mathit{ni}}$}\!\left/ \!\raisebox{-1ex}{$j_{\mathit{nN}}$}\right.\right)}{J_{n+1}^2\left(j_{\mathit{nm}}\right)J_{n+1}^2\left(j_{\mathit{nk}}\right)j_{\mathit{nN}}^2}=\sum _{k=1}^{N-1}{T}_{m,k}^{\mathit{nN}}{T}_{k,i}^{\mathit{nN}}={\delta }_{\mathit{mi}}.$$ E46

Eq. (40) states that the rows and columns of the matrix $T^{\mathit{nN}}$ are orthonormal so that $T^{\mathit{nN}}$ is an orthogonal matrix. Furthermore, $T^{\mathit{nN}}$ is also symmetric. Eq. (46) can be written in matrix form as

$$T^{\mathit{nN}}{T}^{\mathit{nN}}=T^{\mathit{nN}}{\left(T^{\mathit{nN}}\right)}^T=\mathbf{I}.$$ E47

Therefore, the $T^{\mathit{nN}}$ matrix is unitary and furthermore orthogonal since the entries are real.

Using the symmetric, orthogonal transformation matrix $T^{\mathit{nN}}$ , the forward transform from Eq. (26) can be written in as

$$\begin{array}{ll}F\left({\rho }_{\mathit{nm}}\right)& =\frac{R^2}{j_{\mathit{nN}}}\sum _{k=1}^{N-1}{T}_{m,k}^{\mathit{nN}}\frac{J_{n+1}\left(j_{\mathit{nm}}\right)}{J_{n+1}\left(j_{\mathit{nk}}\right)}f\left(r_{\mathit{nk}}\right)\\ & =\frac{j_{\mathit{nN}}}{W_{\rho }^2}\sum _{k=1}^{N-1}{T}_{m,k}^{\mathit{nN}}\frac{J_{n+1}\left(j_{\mathit{nm}}\right)}{J_{n+1}\left(j_{\mathit{nk}}\right)}f\left(r_{\mathit{nk}}\right)\end{array}$$ E48

Similarly, the inverse discrete transform of Eq. (29) can be written as

$$\begin{array}{ll}f\left(r_{\mathit{nk}}\right)& =\frac{j_{\mathit{nN}}}{R^2}\sum _{m=1}^{N-1}{T}_{k,m}^{\mathit{nN}}\frac{J_{n+1}\left(j_{\mathit{nk}}\right)}{J_{n+1}\left(j_{\mathit{nm}}\right)}F\left({\rho }_{\mathit{nm}}\right)\\ & =\frac{W_{\rho }^2}{j_{\mathit{nN}}}\sum _{m=1}^{N-1}{T}_{k,m}^{\mathit{nN}}\frac{J_{n+1}\left(j_{\mathit{nk}}\right)}{J_{n+1}\left(j_{\mathit{nm}}\right)}F\left({\rho }_{\mathit{nm}}\right).\end{array}$$ E49

6. Discrete forward and inverse Hankel transform

From the previous section is it clear that the two natural choices of kernel for a DHT are either $Y_{m,k}^{\mathit{nN}}$ or $T_{m,k}^{\mathit{nN}}$ . $Y_{m,k}^{\mathit{nN}}$ is closer to the discretized version of the continuous Hankel transform that we hope the discrete version emulates. However, $T_{m,k}^{\mathit{nN}}$ is an orthogonal and symmetric matrix, therefore it is energy preserving and will be shown to lead to a Parseval-type relationship if chosen as the kernel for the DHT. Thus, to define a discrete Hankel transform (DHT), we can use either formulation:

$$F_m=\sum _{k=1}^{N-1}{Y}_{m,k}^{\mathit{nN}}{f}_k\text{or}{F}_m=\sum _{k=1}^{N-1}{T}_{m,k}^{\mathit{nN}}{f}_k.$$ E50

Here, the transform is of any $N-1$ dimensional vector $f_k$ to any $N-1$ dimensional vector $F_m$ for the integers $m,k$ where $1\le m,k<N-1$ . This can be written in matrix form as

$$\mathbf{F}=Y^{\mathit{nN}}\mathbf{f}\text{or}\mathbf{F}=T^{\mathit{nN}}\mathbf{f}$$ E51

where $\mathbf{F}$ is any $N-1$ dimensional column vector and $\mathbf{f}$ is also any column vector, defined in the same manner.

The inverse discrete Hankel transform (IDHT) is then given by

$$f_k=\sum _{m=1}^{N-1}{Y}_{k,m}^{\mathit{nN}}{F}_m\text{or}{f}_k=\sum _{m=1}^{N-1}{T}_{k,m}^{\mathit{nN}}{F}_m.$$ E52

This can also be written in matrix form as

$$\mathbf{f}=Y^{\mathit{nN}}\mathbf{F}\text{or}\mathbf{f}=T^{\mathit{nN}}\mathbf{F}.$$ E53

We note that the forward and inverse transforms are the same.

Proof

We show the proof for the $Y^{\mathit{nN}}$ formulation, but it proceeds similarly if $Y^{\mathit{nN}}$ is replaced with $T^{\mathit{nN}}$ . Substituting Eq. (52) into the right hand side of (50) gives

$$\sum _{k=1}^{N-1}{Y}_{m,k}^{\mathit{nN}}{f}_k=\sum _{k=1}^{N-1}{Y}_{m,k}^{\mathit{nN}}\left[\sum _{p=1}^{N-1}{Y}_{k,p}^{\mathit{nN}}{F}_p\right].$$ E54

Switching the order of the summation in Eq. (54) gives

$$\sum _{p=1}^{N-1}\underset{{\delta }_{\mathit{mp}}}{\underbrace{\sum _{k=1}^{N-1}{Y}_{m,k}^{\mathit{nN}}{Y}_{k,p}^{\mathit{nN}}}}{F}_p=\sum _{p=1}^{N-1}{\delta }_{\mathit{mp}}{F}_p=F_m$$ E55

The inside summations as indicated in Eq. (55) are recognized as yielding the Dirac-delta function, the orthogonality property of Eq. (40) (or Eq. (46) if using $T^{\mathit{nN}}$ ), which in turn yields the desired result. This proves that the DHT given by (50) can be inverted by (52).

7. Generalized Parseval theorem

Inner products are preserved and thus energies are preserved under the $T^{\mathit{nN}}$ matrix formulation. To see this, consider any two vectors given by the transform $\mathbf{g}=T^{\mathit{nN}}\mathbf{G}$ , $\mathbf{h}=T^{\mathit{nN}}\mathbf{H}$ then

$${\mathbf{g}}^T\mathbf{h}={\left(T^{\mathit{nN}}\mathbf{G}\right)}^T{T}^{\mathit{nN}}\mathbf{H}={\mathbf{G}}^T\underset{=I}{\underbrace{{\left(T^{\mathit{nN}}\right)}^T{T}^{\mathit{nN}}}}\mathbf{H}={\mathbf{G}}^T\mathbf{H}.$$ E56

The $Y^{\mathit{nN}}$ matrix formulation does not directly preserve inner products:

$${\mathbf{g}}^T\mathbf{h}={\left(Y^{\mathit{nN}}\mathbf{G}\right)}^T{Y}^{\mathit{nN}}\mathbf{H}={\mathbf{G}}^T{\left(Y^{\mathit{nN}}\right)}^T{Y}^{\mathit{nN}}\mathbf{H}.$$ E57

However, under the $Y^{\mathit{nN}}$ formulation, the inner product between $\frac{g_k}{J_{n+1}\left(j_{\mathit{nk}}\right)}$ and $\frac{h_k}{J_{n+1}\left(j_{\mathit{nk}}\right)}$ is preserved. To see this, we calculate the inner product between them as

$$\begin{array}{ll}\sum _{k=1}^{N-1}\frac{g_k}{J_{n+1}\left(j_{\mathit{nk}}\right)}\frac{h_k}{J_{n+1}\left(j_{\mathit{nk}}\right)}& =\sum _{k=1}^{N-1}\frac{1}{J_{n+1}^2\left(j_{\mathit{nk}}\right)}\sum _{p=1}^{N-1}{Y}_{k,p}^{\mathit{nN}}{G}_p\sum _{q=1}^{N-1}{Y}_{k,q}^{\mathit{nN}}{H}_q\\ & =\sum _{p=1}^{N-1}\sum _{q=1}^{N-1}\frac{1}{J_{n+1}^2\left(j_{\mathit{np}}\right)}\underset{{\delta }_{\mathit{pq}}}{\underbrace{\sum _{k=1}^{N-1}\frac{4\left(\raisebox{1ex}{$j_{\mathit{nk}}{j}_{\mathit{np}}$}\!\left/ \!\raisebox{-1ex}{$j_{\mathit{nN}}$}\right.\right)J_n\left(\raisebox{1ex}{$j_{\mathit{nk}}{j}_{\mathit{nq}}$}\!\left/ \!\raisebox{-1ex}{$j_{\mathit{nN}}$}\right.\right)}{j_{\mathit{nN}}^2{J}_{n+1}^2\left(j_{\mathit{nk}}\right)J_{n+1}^2\left(j_{\mathit{nq}}\right)}}}{H}_q{G}_p\end{array}$$ E58

Making use of the now-present Dirac-delta function, Eq. (58) simplifies to give a modified Parseval relationship

$$\sum _{k=1}^{N-1}\left(\frac{g_k}{J_{n+1}\left(j_{\mathit{nk}}\right)}\right)\left(\frac{h_k}{J_{n+1}\left(j_{\mathit{nk}}\right)}\right)=\sum _{p=1}^{N-1}\left(\frac{H_p}{J_{n+1}\left(j_{\mathit{np}}\right)}\right)\left(\frac{G_p}{J_{n+1}\left(j_{\mathit{np}}\right)}\right).$$ E59

In other words, under a DHT using the $Y^{\mathit{nN}}$ matrix, inner products of the scaled functions are preserved but not the inner products of the functions themselves.

As a consequence of the orthogonality property of $T^{\mathit{nN}}$ , the $T^{\mathit{nN}}$ based DHT is energy preserving, meaning that

$${\overline{\mathbf{F}}}^T\mathbf{F}={\overline{\mathbf{f}}}^T\mathbf{f}.$$ E60

where the overbar indicates a conjugate transpose and the superscript T indicates a transpose.

For the formulation with $Y^{\mathit{nN}}$ as the transformation kernel, the equivalent expression is

$${\overline{\mathbf{F}}}^T\mathbf{F}={\left(Y^{\mathit{nN}}\overline{\mathbf{f}}\right)}^T{Y}^{\mathit{nN}}\mathbf{f}={\overline{\mathbf{f}}}^T{\left(Y^{\mathit{nN}}\right)}^T{Y}^{\mathit{nN}}\mathbf{f}.$$ E61

It is obvious from Eq. (59) that if we define the “scaled” vector

$$f_k^{\mathit{\text{Scaled}}}=\frac{f_k}{J_{n+1}\left(j_{\mathit{nk}}\right)}\text{and}{F}_p^{\mathit{\text{Scaled}}}=\frac{F_p}{J_{n+1}\left(j_{\mathit{np}}\right)},$$ E62

then by straighforward substitution of scaled vectors and their conjugates, it follows that

$${\overline{\left({\mathbf{F}}^{\mathit{\text{Scaled}}}\right)}}^T{\mathbf{F}}^{\mathit{\text{Scaled}}}={\overline{\left({\mathbf{f}}^{\mathit{\text{Scaled}}}\right)}}^T{\mathbf{f}}^{\mathit{\text{Scaled}}}.$$ E63

8. Transform rules

In keeping with the development of the well-known discrete Fourier transform, we develop the standard toolkit of rules for the discrete Hankel transform. In the following, $Y^{\mathit{nN}}$ is used but all expressions apply equally if $Y^{\mathit{nN}}$ is replaced with $T^{\mathit{nN}}$ .

8.1 Transform of Kronecker-Delta function

The discrete counterpart of the Dirac-delta function is the Kronecker-delta function, ${\delta }_{{\mathit{kk}}_0}$ . We recall that the DHT as defined above is a matrix transform from a $N-1$ dimensional vector to another. The vector ${\delta }_{{\mathit{kk}}_o}$ is interpreted as the vector as having zero entries everywhere except at position $k=k_0$ ( $k_0$ fixed so ${\delta }_{{\mathit{kk}}_0}$ is a vector), or in other words, the k₀th column of the $N-1$ sized identity matrix. The DHT of the Kronecker-delta can be found from the definition of the forward transform via

$$\mathbb{H}\left({\delta }_{{\mathit{kk}}_o}\right)=\sum _{k=1}^{N-1}{Y}_{m,k}{\delta }_{{\mathit{kk}}_o}=Y_{m,k_0}^{\mathit{nN}}$$ E64

The symbol $\mathbb{H}\left(\cdot \right)$ is used to denote the operation of taking the discrete Hankel transform. This gives us our first DHT transform pair of order n dimension $N-1$ , and we denote this relationship as

$${\delta }_{{\mathit{kk}}_o}\iff Y_{m,k_0}^{n,N}$$ E65

Here, $f_k\iff F_m$ denotes a transform pair and $Y_{m,k_0}^{\mathit{nN}}$ is k₀th column of the matrix $Y^{\mathit{nN}}$ .

8.2 Inverse transform of the Kronecker Delta function

From Eq. (65), we can deduce the vector $f_k$ that transforms to the Kronecker-delta vector ${\delta }_{{\mathit{mm}}_o}$ function. Namely, we take the forward transform of

$$f_k=Y_{k,m_0}^{n,N}.$$ E66

As before, $Y_{k,m_0}^{\mathit{nN}}$ represents the m₀th column of the transformation matrix $Y^{\mathit{nN}}$ . From the forward definition of the transform, Eq. (50), the transform of $Y_{k,m_0}^{n,N}$ is given by

$$F_m=\sum _{k=1}^{N-1}{Y}_{m,k}^{\mathit{nN}}{f}_k=\sum _{k=1}^{N-1}{Y}_{m,k}^{\mathit{nN}}{Y}_{k,m_0}^{\mathit{nN}}={\delta }_{{\mathit{mm}}_0},$$ E67

where we have used the orthogonality relationship (40). This gives us another DHT pair:

$$Y_{k,m_0}^{n,N}\iff {\delta }_{{\mathit{mm}}_o}.$$ E68

8.3 The generalized shift operator

For a one-dimensional Fourier transform, one of the known transform rules is the shift rule, which states that

$$f\left(x-a\right)={\mathbb{F}}^{-1}\left\{e^{-\mathit{ia\omega }}\widehat{f}\left(\omega \right)\right\}=\frac{1}{2\pi }\underset{-\infty }{\overset{\infty }{\int }}\left\{e^{-\mathit{ia\omega }}\widehat{f}\left(\omega \right)\right\}{e}^{\mathit{i\omega t}}\mathit{d\omega }.$$ E69

In Eq. (69), $\widehat{f}\left(\omega \right)$ is the Fourier transform of $f\left(x\right)$ , ${\mathbb{F}}^{-1}$ denotes an inverse Fourier transform and $e^{-\mathit{ia\omega }}$ is the kernel of the Fourier transform operator. Motivated by this result, we define a generalized-shift operator by finding the inverse DHT of the DHT of the function multiplied by the DHT kernel. This is a discretized version of the definition of a generalized shift operator as proposed by Levitan [8] (he suggested the complex conjugate of the Fourier operator, which for Fourier transforms is the inverse transform operator). We thus propose the definition of the generalized-shifted function to be given by

$$f_{k,k_o}^{\mathit{\text{shift}}}=\sum _{p=1}^{N-1}{Y}_{k,p}^{\mathit{nN}}\underset{\begin{array}{l}\text{shift in Hankel}\\ \text{domain}\end{array}}{\underbrace{\left\{Y_{p,k_o}^{\mathit{nN}}{F}_p\right\}}},$$ E70

where $1\le k,k_o\le N-1$ . For a single, fixed value of $k_o$ , then $f_{k,k_o}^{\mathit{\text{shift}}}$ is another $N-1$ vector, with the notation $f_{k,k_o}^{\mathit{\text{shift}}}$ implying a $k_0$ -shifted version of $f_k$ . This generalizes the notion of the shift, usually denoted $f_{k-k_o}$ , which inevitably encounters difficulty when the subscript $k-k_o$ falls outside of the range $\left[1,N-1\right]$ . We note that if all possible shifts $k_o$ are considered, then $f_{k,k_o}^{\mathit{\text{shift}}}$ is a $N-1$ square matrix (in other words, a two dimensional structure), whereas the original un-shifted $f_k$ is an $N-1$ vector. For the discrete Fourier transform, when the shifted subscript $k-k_o$ falls outside the range of the indices, is it usually interpreted modulo the size of the DFT. However, the kernel of the Fourier transform is periodic so this does not create difficulties for the DFT. The Bessel functions are not periodic so the same trick cannot be used with the Hankel transform. In fact, this lack of periodicity and lack of simple relationship between $J_n\left(x-y\right)$ and $J_n\left(x\right)$ is the reason that the continuous Hankel transform does not have a convolution-multiplication rule [13]. Thus, the notation $f_{k-k_o}$ would not make mathematical sense when used with the DHT. With the definition given by Eq. (70), no such confusion arises since the definition is unambiguous for all allowable values of $k$ and $k_o$ .

The shifted function $f_{k,k_o}^{\mathit{\text{shift}}}$ can also be expressed in terms of the original un-shifted function $f_k$ . Using the definition of $F_m$ from Eq. (50) and a dummy change of variable, then Eq. (70) becomes

$$f_{k,k_o}^{\mathit{\text{shift}}}=\sum _{p=1}^{N-1}{Y}_{k,p}^{\mathit{nN}}{Y}_{p,k_o}^{\mathit{nN}}{F}_p=\sum _{p=1}^{N-1}{Y}_{k,p}^{\mathit{nN}}{Y}_{p,k_o}^{\mathit{nN}}\sum _{m=1}^{N-1}{Y}_{p,m}^{\mathit{nN}}{f}_m.$$ E71

Changing the order of summation gives

$$f_{k,k_o}^{\mathit{\text{shift}}}=\sum _{p=1}^{N-1}{Y}_{k,p}^{\mathit{nN}}{Y}_{p,k_o}^{\mathit{nN}}{F}_p=\sum _{m=1}^{N-1}\underset{\text{shift operator}}{\underbrace{\sum _{p=1}^{N-1}{Y}_{k,p}^{\mathit{nN}}{Y}_{p,k_o}^{\mathit{nN}}{Y}_{p,m}^{\mathit{nN}}}}{f}_m.$$ E72

As indicated in Eq. (72), the quantity in brackets can be considered to be a type of shift operator acting on the original unshifted function. We can define this as

$$S_{k,k_o,m}^{\mathit{nN}}=\sum _{p=1}^{N-1}{Y}_{k,p}^{\mathit{nN}}{Y}_{p,k_o}^{\mathit{nN}}{Y}_{p,m}^{\mathit{nN}}.$$ E73

It then follows that Eq. (72) can be written as

$$f_{k,k_o}^{\mathit{\text{shift}}}=\sum _{m=1}^{N-1}{S}_{k,k_o,m}^{\mathit{nN}}{f}_m.$$ E74

This triple-product shift operator is similar to previous definitions of shift operators for multidimensional Fourier transforms that rely on Hankel transforms [1, 2] and of generalized Hankel convolutions [14, 15, 16].

8.4 Transform of the generalized shift operator

We now consider the forward DHT transform of the shifted function $f_{k,k_o}^{\mathit{\text{shift}}}$ . From the definition, the DHT of the shifted function can be found from

$$\sum _{k=1}^{N-1}{Y}_{m,k}^{\mathit{nN}}{f}_{k,k_o}^{\mathit{\text{shift}}}=\sum _{k=1}^{N-1}{Y}_{m,k}^{\mathit{nN}}\sum _{p=1}^{N-1}{Y}_{k,p}^{\mathit{nN}}{Y}_{p,k_o}^{\mathit{nN}}{F}_p.$$ E75

Changing the order of summation gives

$$\sum _{p=1}^{N-1}\underset{={\delta }_{\mathit{mp}}}{\underbrace{\sum _{k=1}^{N-1}{Y}_{m,k}^{\mathit{nN}}{Y}_{k,p}^{\mathit{nN}}}}{Y}_{p,k_o}^{\mathit{nN}}{F}_p=\sum _{p=1}^{N-1}{\delta }_{\mathit{mp}}{Y}_{p,k_o}^{\mathit{nN}}{F}_p=Y_{m,k_o}^{\mathit{nN}}{F}_m.$$ E76

This yields another transform pair and is the shift-modulation rule. This rule analogous to the shift-modulation rule for regular Fourier transforms whereby a shift in the spatial domain is equivalent to modulation in the frequency domain:

$$f_{k,k_o}^{\mathit{\text{shift}}}\iff Y_{m,k_o}^{\mathit{nN}}{F}_m.$$ E77

Note that Eq. (77) does not imply a summation over the m index. For a fixed value of $k_o$ on the left hand side, the corresponding transformed value of $F_m$ is multiplied by the $\left(m,k_o\right)$ th entry of the $Y^{\mathit{nN}}$ matrix.

8.5 Modulation

We consider the forward DHT of a function “modulated” in the space domain $f_k=Y_{k,k_o}^{\mathit{nN}}{g}_k$ . Here, the interpretation of $f_k=Y_{k,k_o}^{\mathit{nN}}{g}_k$ is that the kth entry of the vector g is multiplied by the $\left(k,k_o\right)$ th entry of $Y^{\mathit{nN}}$ for a fixed value of $k_o$ . No summation is implied so this is not a dot product; both $f_k$ and $Y_{k,k_o}^{\mathit{nN}}{g}_k$ are $N-1$ vectors. Again, we implement the definition of the forward transform

$$\sum _{k=1}^{N-1}{Y}_{m,k}^{\mathit{nN}}{f}_k=\sum _{k=1}^{N-1}{Y}_{m,k}^{\mathit{nN}}{Y}_{k,k_o}^{\mathit{nN}}{g}_k,$$ E78

and write $g_k$ in terms of its inverse transform

$$g_k=\sum _{p=1}^{N-1}{Y}_{k,p}^{\mathit{nN}}{G}_p.$$ E79

Then Eq. (78) becomes

$$\sum _{k=1}^{N-1}{Y}_{m,k}^{\mathit{nN}}{f}_k=\sum _{k=1}^{N-1}{Y}_{m,k}^{\mathit{nN}}{Y}_{k,k_o}^{\mathit{nN}}{g}_k=\sum _{k=1}^{N-1}{Y}_{m,k}^{\mathit{nN}}{Y}_{k,k_o}^{\mathit{nN}}\sum _{p=1}^{N-1}{Y}_{k,p}^{\mathit{nN}}{G}_p.$$ E80

Interchanging the order of summation gives

$$\sum _{p=1}^{N-1}\underset{\text{shift operator}}{\underbrace{\sum _{k=1}^{N-1}{Y}_{m,k}^{\mathit{nN}}{Y}_{k,k_o}^{\mathit{nN}}{Y}_{k,p}^{\mathit{nN}}}}{G}_p=G_{m,k_o}^{\mathit{\text{shift}}}.$$ E81

By comparing Eq. (81) with Eqs. (72) and (73), we recognize the shift operator as indicated in (81). This produces a modulation-shift rule as would be expected so that the forward DHT of a modulated function is equivalent to a generalized shift in the frequency domain. This yields another transform pair:

$$Y_{k,k_o}^{\mathit{nN}}{g}_k\iff G_{m,k_o}^{\mathit{\text{shift}}}.$$ E82

In other words, Eq. (82) says that modulation in the space domain is equivalent to shift in the frequency domain, as would be expected for a (generalized) Fourier transform.

8.6 Convolution

We consider the convolution using the generalized shifted function previously defined. The convolution of two functions is defined as

$$f_k={\left(g^{\ast }h\right)}_k=\sum _{k_0=1}^{N-1}{g}_{k_o}{h}_{k,k_o}^{\mathit{\text{shift}}}.$$ E83

The meaning of Eq. (83) follows from the traditional definition of a convolution: multiply one of the functions by a shifted version of a second function and then sum over all possible shifts.

Subsequently, from the definition of the inverse transforms, we obtain

$$\begin{array}{ll}{f}_k=\sum _{k_0=1}^{N-1}{g}_{k_o}{h}_{k,k_o}^{\mathit{\text{shift}}}& =\sum _{k_0=1}^{N-1}\underset{g_{k_o}}{\underbrace{\sum _{q=1}^{N-1}{Y}_{k_o,q}^{\mathit{nN}}{G}_q}}\underset{h_{k,k_o}^{\mathit{\text{shift}}}}{\underbrace{\sum _{p=1}^{N-1}{Y}_{k,p}^{\mathit{nN}}{Y}_{p,k_o}^{\mathit{nN}}{H}_p}}\\ & =\sum _{q=1}^{N-1}\sum _{p=1}^{N-1}\underset{={\delta }_{\mathit{pq}}}{\underbrace{\sum _{k_0=1}^{N-1}{Y}_{p,k_o}^{\mathit{nN}}{Y}_{k_o,q}^{\mathit{nN}}}}{Y}_{k,p}^{\mathit{nN}}{H}_p{G}_q.\end{array}$$ E84

But from the orthogonality relationship (40), the summation over $k_0$ gives the Kronecker delta function, so that Eq. (84) becomes

$$\begin{array}{ll}{\left(g^{\ast }h\right)}_k& =\sum _{k_0=1}^{N-1}{g}_{k_o}{h}_{k,k_o}^{\mathit{\text{shift}}}=\sum _{q=1}^{N-1}\sum _{p=1}^{N-1}{\delta }_{\mathit{pq}}{Y}_{k,p}^{\mathit{nN}}{H}_p{G}_q\\ & =\sum _{p=1}^{N-1}{Y}_{k,p}^{\mathit{nN}}\left(H_p{G}_p\right)\end{array}$$ E85

The right hand side of Eq. (85) is clearly the inverse transform of the product of the transforms $H_p{F}_p$ . This gives us another transform pair

$${\left(g^{\ast }h\right)}_k=\sum _{k_0=1}^{N-1}{g}_{k_o}{h}_{k,k_o}^{\mathit{\text{shift}}}\iff H_m{G}_m.$$ E86

It follows from Eq. (85) that interchanging the roles of $g$ and $h$ will yield the same result, meaning

$$\sum _{k_0=1}^{N-1}{g}_{k,k_o}^{\mathit{\text{shift}}}{h}_{k_o}=\sum _{p=1}^{N-1}{Y}_{k,p}^{\mathit{nN}}{G}_p{H}_p.$$ E87

Therefore, it follows that

$${\left(h^{\ast }g\right)}_k=\sum _{k_0=1}^{N-1}{g}_{k,k_o}^{\mathit{\text{shift}}}{h}_{k_o}=\sum _{k_0=1}^{N-1}{g}_{k_o}{h}_{k,k_o}^{\mathit{\text{shift}}}={\left(g^{\ast }h\right)}_k.$$ E88

8.7 Multiplication

We now consider the forward transform of a product in the space domain $f_k=g_k{h}_k$ so that

$$\sum _{k=1}^{N-1}{Y}_{m,k}^{\mathit{nN}}{g}_k{h}_k=\sum _{k=1}^{N-1}{Y}_{m,k}^{\mathit{nN}}\underset{g_k}{\underbrace{\sum _{q=1}^{N-1}{Y}_{k,q}^{\mathit{nN}}{G}_q}} \underset{h_k}{\underbrace{\sum _{p=1}^{N-1}{Y}_{k,p}^{\mathit{nN}}{H}_p}}.$$ E89

Rearranging gives

$$\begin{array}{ll}\sum _{k=1}^{N-1}{Y}_{m,k}^{\mathit{nN}}{g}_k{h}_k& =\sum _{q=1}^{N-1}{G}_q\sum _{p=1}^{N-1}\underset{\text{shift operator}}{\underbrace{\sum _{k=1}^{N-1}{Y}_{m,k}^{\mathit{nN}}{Y}_{k,q}^{\mathit{nN}}{Y}_{k,p}^{\mathit{nN}}}}{H}_p\\ & =\sum _{q=1}^{N-1}{G}_q{H}_{m,q}^{\mathit{\text{shift}}}={\left(G^{\ast }H\right)}_m.\end{array}$$ E90

This gives us yet another transform pair that says that multiplication in the spatial domain is equivalent to convolution in the transform domain:

$$g_k{h}_k\iff \sum _{q=1}^{N-1}{G}_q{H}_{m,q}^{\mathit{\text{shift}}}={\left(G^{\ast }H\right)}_m.$$ E91

Interchanging the roles of $G$ and $H$ in Eq. (91) demonstrates that convolution in the transform domain also commutes:

$${\left(G^{\ast }H\right)}_m=\sum _{q=1}^{N-1}{G}_q{H}_{m,q}^{\mathit{\text{shift}}}=\sum _{q=1}^{N-1}{G}_{m,q}^{\mathit{\text{shift}}}{H}_q={\left(H^{\ast }G\right)}_m.$$ E92

9. Using the DHT to approximate the continuous transform

9.1 Approximation to the continuous transform

Eqs. (26) and (29) show how the DHT can be used to calculate the continuous Hankel transform at finite points. From Eqs. (26) and (29), it is clear that given a continuous function $f\left(r\right)$ evaluated at the discrete points $r_{\mathit{nk}}$ (given by Eq. (31)) in the space domain ( $1\le k\le N-1$ ), its nth-order Hankel-transform function $F\left(\rho \right)$ evaluated at the discrete points ${\rho }_{\mathit{nm}}$ (given in Eq. (31)) in the frequency domain ( $1\le m\le N-1$ ), can be approximately given by

$$F\left[m\right]=\alpha \sum _{k=1}^{N-1}{Y}_{m,k}^{\mathit{nN}}f\left[k\right]\Rightarrow \mathbf{F}=\alpha Y^{\mathit{nN}}\mathbf{f}$$ E93

where $\alpha$ is a scaling factor to be discussed below, and $F\left[m\right]=F\left({\rho }_{\mathit{nm}}\right)$ , $f\left[k\right]=f\left(r_{\mathit{nk}}\right)$ .

Similarly, given a continuous function $F\left(\rho \right)$ evaluated at the discrete points ${\rho }_{\mathit{nm}}$ in the frequency domain ( $1\le m\le N-1$ ), its nth-order inverse Hankel transform $f\left(r\right)$ evaluated at the discrete points $r_{\mathit{nk}}$ ( $1\le k\le N-1$ ), can be approximately given by

$$f\left[k\right]=\frac{1}{\alpha }\sum _{m=1}^{N-1}{Y}_{m,k}^{\mathit{nN}}F\left[m\right]\Rightarrow \mathbf{f}=\frac{1}{\alpha }{Y}^{\mathit{nN}}\mathbf{F}$$ E94

For both the forward and inverse transforms, $\alpha$ is a scaling factor and $\alpha =\frac{R^2}{j_{\mathit{nN}}}$ or equivalently $\alpha =\frac{j_{\mathit{nN}}}{W_{\rho }^2}$ , where R is the effective space limit and $W_{\rho }$ is the effective band limit (in m⁻¹). The scaling factor $\alpha$ chosen for using the DHT to approximate the CHT depends on whether information is known about the band-limit or space-limit. We already introduced the idea of an effective limit in the previous sections, where a function was defined as being “effectively limited in space by R” means that if $r>R$ , then $f\left(r\right)\approx 0$ for all $r>R$ . In other words, the function can be made as close to zero as desired by selecting an R that is large enough. The same idea can be applied to the spatial frequency domain, where the effective domain is denoted by $W_{\rho }$ .

The relationship $W_{\rho }R=j_{\mathit{nN}}$ , derived in the previous sections, holds between the ranges in space and frequency. Choosing $N$ determines the dimension (size) of the DHT and determines $j_{\mathit{nN}}$ . The determination of $j_{\mathit{nN}}$ (via choosing N) determines the range in one domain once the range in the other domain is chosen. In fact, any two of $R,W_{\rho },j_{\mathit{nN}}$ can be chosen but the third must follow from $W_{\rho }R=j_{\mathit{nN}}$ . A similar relationship applies when using the discrete Fourier transform, any two of the range in each domain and the size of the DFT can be chosen independently. In previous sections, we showed that the size of the DHT required can be quickly approximated from $2\mathit{WR}=\frac{W_{\rho }R}{\pi }\approx \left(N+\frac{n}{2}\right)$ .

9.2 Sampling points

In order to properly use the discrete transform to approximate the continuous transform, a function has to be sampled at specific discretization points. For a finite spatial range $\left[0,R\right]$ and a Hankel transform of order $n$ , these sampling points are given in the space domain as $r_{\mathit{nk}}$ and frequency domain by ${\rho }_{\mathit{nm}}$ , given in Eq. (31) and repeated here for convenience

$$r_{\mathit{nk}}=\frac{j_{\mathit{nk}}}{W_{\rho }}=\frac{j_{\mathit{nk}}R}{j_{\mathit{nN}}}{\rho }_{\mathit{nm}}=\frac{j_{\mathit{nm}}}{R}=\frac{j_{\mathit{nm}}{W}_{\rho }}{j_{\mathit{nN}}}k,m=1\dots N-1$$ E95

It is important to note that as in the case of the computation of the transformation matrix $Y^{\mathit{nN}}$ , the first Bessel zero $j_{n1}$ used in computing the discretization points is the first non-zero value. Eq. (95) demonstrates that some of the ideas known for the DFT also apply to the DHT. That is, making the spatial domain larger (larger R) implies making the sampling density tighter in frequency (the ${\rho }_{\mathit{nm}}$ get closer together). Similarly, making the frequency domain larger (larger $W_{\rho }$ ) implies a tighter sampling density (smaller step size) in the spatial domain. Although $j_{\mathit{nm}}$ are not equispaced, they are nearly so for higher values of m and for purposes of developing quick intuitions on ideas such as sampling density, if is convenient to approximately think of $j_{\mathit{nk}}\approx \left(k+\frac{n}{2}\right)\pi$ .

9.3 Implementation and availability of the software

The software used to calculate the DHT is based on the MATLAB programming language. The software can be downloaded from

• http://dx.doi.org/10.6084/m9.figshare.1453205

• https://github.com/uchouinard/DiscreteHankelTransform

The implementation of the discrete Hankel transform is decomposed into distinct functions. These functions consist of various steps that have to be performed in order to properly execute the transform. These steps are as follows:

• Calculate $N$ Bessel zeros of the Bessel function of order $n$

• Generate of $N$ sample points (if using the DHT to approximate the continuous transform)

• Sample the function (if needed)

• Create the $Y^{\mathit{nN}}$ transformation matrix

• Perform the matrix-function multiplication

The steps are the same regardless if the function is in the space or frequency domain. Furthermore, the $Y^{\mathit{nN}}$ transformation matrix is used for both the forward and inverse transform. The second and third steps in the list above are only needed if the function (vector) to be transformed is not already given as a set of discrete points. In the case of a continuous function, it is important to evaluate the function at the sampling points in Eq. (95). Failing to do so results in the function not being properly transformed since there is a necessary relationship between the sampling points and the transformation matrix $Y^{\mathit{nN}}$ . In order to perform the steps listed above, several Matlab functions have been developed. These functions are shown in Table 1.

Function name	Calling sequence	Description
besselzero	besselzero(n,k,kind)	Calculation of k Bessel zeros of the nth-order Bessel function of kind—developed in [17]
freqSampler	freqSampler(R,zeros)	Creation of sample points in the frequency domain (Eq. (95))
spaceSampler	spaceSampler(R,zeros)	Creation of sample points in the space domain (Eq. (95))
YmatrixAssembly	YmatrixAssembly(n,N,zeros)	Creation of $Y^{\mathit{nN}}$ matrix from the zeros

Table 1.

Set of available functions.

Additionally, the matlab script GuidetoDHT.m is included to illustrate the execution of the necessary computational steps.

9.4 Verification of the software

The software was tested by using the DHT to approximate the computation of both the continuous Hankel forward and inverse transforms and comparing the results with known (continuous) forward and inverse Hankel transform pairs. Different transform orders $n$ were evaluated.

For the purpose of testing the accuracy of the DHT and IDHT, the dynamic error was used, defined as [12]

$$e\left(v\right)=20{\text{log}}_{10}\left[\frac{\left|f\left(v\right)-f^{\ast }\left(\nu \right)\right|}{\text{max}\left|f^{\ast }\left(v\right)\right|}\right]$$ E96

This error function compares the difference between the exact function values $f\left(v\right)$ (evaluated from the continuous function) and the function values estimated via the discrete transform, $f^{\ast }\left(\nu \right)$ , scaled with the maximum value of the discretely estimated samples. The dynamic error uses the ratio of the absolute error to the maximum amplitude of the function on a log scale. Therefore, negative decibel errors imply an accurate discrete estimation of the true transform value. The transform was also tested for accuracy on itself by performing consecutive forward and then inverse transformation. This is done to verify that the transforms themselves do not add errors. For this evaluation, the average absolute error $\frac{1}{N}\sum _{i=1}^N\left|f_i-{f_i}^{\ast }\right|$ was used. The methodology of the testing is given in further detail in [18] and also in the theory paper [3].

10. Summary and conclusions

In this chapter, the theory of the discrete Hankel transform as a “standalone” transform was motivated and presented. The standard operating rules for multiplication, modulation, shift and convolution were also demonstrated. Sampling and interpolation theorems were shown. The theory and numerical steps to use the presented discrete theory for the purpose of approximating the continuous Hankel transform was also shown. Links to the publicly available, open-source numerical code were also included.

Acknowledgments

The author acknowledges the contributions of Mr. Ugo Chouinard, who developed and tested the numerical code to which links are provided in this chapter. This work was financially supported by the Natural Sciences and Engineering Research Council of Canada.

Conflict of interest

The author declares that there are no conflicting interests.