The Discrete Hankel Transform

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.


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.

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

Hankel transform
The nth-order Hankel transform F ρ ð Þ of the function f r ð Þ of a real variable, r ≥ 0, is defined by the integral [4] where J n z ð Þ 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 Thus, Hankel transforms take functions in the spatial r domain and transform them to functions in the spatial frequency ρ domain f r ð Þ ⇔ F ρ ð Þ. The notation ⇔ is used to indicate a Hankel transform pair.

Fourier Bessel series
It is known that functions defined on a finite portion of the real line 0; R ½ , can be expanded in terms of a Fourier Bessel series [5] given by where the order, n, of the Bessel function is arbitrary and j nk denotes the kth root of the nth Bessel function J n z ð Þ. The Fourier Bessel coefficients f k of the function f r ð Þ are given by f r ð ÞJ n j nk r R r dr: Eqs. (3) and (4) can be considered to be a transform pair where the continuous function f r ð Þ is forward-transformed to the discrete vector f k given in (4). The inverse transform is then the operation which returns f r ð Þ 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.

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.

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 r ð Þ contains no frequencies higher than W cycles per meter, then it is completely determined by a series of sampling points given by evaluating f r ð Þ at r ¼ suppose that a function f r ð Þ is band-limited in the frequency Hankel domain so that its spectrum F ρ ð Þ is zero outside of an interval 0; 2πW ½ . 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π ensures that the final units are in s À1 or m À1 . Let us define W ρ ¼ 2πW. Since the Hankel transform F ρ ð Þ is defined on a finite portion of the real line 0; W ρ Â Ã , it can be expanded in terms of a Fourier Bessel series as where the Fourier Bessel coefficients can be found from Eq. (4) as In (6), we have used the fact that f r ð Þ 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 0; W ρ Â Ã can be written as

Interpretation in terms of a jinc
Eq. (8) states In other research work [8], the generalized shift operator R r 0 indicating a shift of r 0 acting on the function f r ð Þ has been defined by the formula 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 More explicitly, where 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 À1=2; À1=2 ½ 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 The function 2J 1 r ð Þ = r 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 Hence, for a zeroth-order Fourier Bessel transform, Eq. (12), the expansion for f r ð Þ reads 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.

Sampling theorem for a space-limited function
Now consider a space-limited function f r ð Þ so that f r ð Þ is zero outside of an interval 0; R ½ . It then follows that it can be expanded on 0; R ½ in terms of a Fourier Bessel series so that where the Fourier Bessel coefficients can be found from Here, we have used the definition of the Hankel transform F ρ ð Þ, Eq. (1), in the right hand side of Eq. (20). Hence, the function can be written as From Eq. (21), it is evident that the samples F j nk R determine the function f r ð Þ and hence its transform F ρ ð Þ completely. Another way of looking at this is that space limiting a function to 0; R ½ implies discretization in spatial frequency space, at

Interpolation theorem for a space-limited function
The Hankel transform of the function can then be found from a forward Hankel transformation as Using Eq. (9), Eq. (22) can be simplified to give Eq. (23) gives the formula for interpolating the samples F j nk R to give the continu-

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.

Forward transform
We demonstrated above that a band-limited function, with ρ , W ρ ¼ 2πW can be written as Evaluating the previous Eq. (24) at the sampling points ρ nm ¼ For m , N, then ρ nm ¼ , W ρ , and Eq. (25), summing over infinite k, is exact.
≥ W ρ and by the assumption of the bandlimited nature of the function, F ρ nm ð Þ ¼ 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 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 lim k!∞ f j nk W ρ ¼ 0, which means that for any arbitrarily small ε, there must exist an integer N for which f 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 In this case, the truncated sum in Eq. (26) does not represent F ρ nm ð Þ 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 j nk W ρ terms do not contribute significantly.

Inverse transform
Concomitantly, we know that for any space-limited function then for r , R, we can write More specifically, suppose that we follow the logic from the previous section that the function f r ð Þ that was bandlimited but also "effectively space-limited" due the truncation of the series in Eq. (25) at N. In that case then R ¼ j nN W ρ and the band-limit implies that F ρ ð Þ ¼ 0 for ρ . W ρ . Following this logic and using R ¼ where we truncated the series in Eq. (28) at N by using the fact that F ρ ð Þ ¼ 0 for Compare Eq. (29) to the "forward" transform from Eq. (26), repeated here for convenience, where we found that for m ¼ 1: Eqs. (29) and (30) are the fundamental relations used for the discrete Hankel transform proposed in the following sections.

Intuitive discretization scheme and kernel
The preceding development shows that a "natural," N-dimensional discretization scheme in finite space 0; R ½ and finite frequency space 0; W ρ Â Ã is given by The relationship W ρ R ¼ j 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 It is pointed out in [10] that the zeros of J n z ð Þ are almost evenly spaced at intervals of π and that the spacing becomes exactly π in the limit as z ! ∞. In fact, it is shown in [10] that a simple asymptotic form for the Bessel function is given by Eq. (32) becomes a better approximation to J n z ð Þ as z ! ∞. The zeros of the cosine function are at odd multiples of π=2. Therefore, an approximation to the Bessel zero, j nk is given by Using this approximation, then For larger values of N as would typically be used in a discretization scheme, then we can write approximately 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 ¼ 2WT 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 ¼ 2WR À n=2 ð Þnumbers and it will certainly be true that specifying N ¼ 2WR numbers will specify the function, in exact analogy to Shannon's result.

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 2 J 2 nþ1 j nk À ÁJ n j nk j nm j nN : 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 where j nm represents the mth zero of the nth-order Bessel function J n x ð Þ, and δ mi is the Kronecker delta function, defined as

Transformation matrix
With inspiration from the notation in [11], and an additional scaling factor of 1=j nN , we define an N À 1 ð ÞÂ N À 1 ð Þtransformation matrix with the (m,k)th entry given by 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 J n j ni j nk = j nN J n j nk j nm = j nN Eq. (40) states that the rows and columns of the matrix Y nN m, k are orthonormal and can be written in matrix form as where I is the N À 1 dimensional identity matrix and we have written the N À 1 square matrix Y nN m, k as Y nN . The forward and inverse truncated and discretized transforms given in Eqs. (26) and (29) can be expressed in terms of Y nN . The forward transform, Eq. (26), can be written as Similarly, the inverse transform, Eq. (29), can be written as

Another choice of transformation matrix
Following the notation in [12], we can also define a different N À 1 ð ÞÂ N À 1 ð Þ transformation matrix with the (m,k)th entry given by 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 nN m, k ¼ T nN k, m so that T nN is a real, symmetric matrix. The relationship between the T nN m, k and Y nN m, k matrices is given by The orthogonality relationship, Eq. (37), can be written as Eq. (40) states that the rows and columns of the matrix T nN are orthonormal so that T nN is an orthogonal matrix. Furthermore, T nN is also symmetric. Eq. (46) can be written in matrix form as Therefore, the T nN matrix is unitary and furthermore orthogonal since the entries are real. Using the symmetric, orthogonal transformation matrix T nN , the forward transform from Eq. (26) can be written in as Similarly, the inverse discrete transform of Eq. (29) can be written as

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 nN m, k or T nN m, k . Y nN m, k is closer to the discretized version of the continuous Hankel transform that we hope the discrete version emulates. However, T nN m, k 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: 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 ≤ m, k , N À 1. This can be written in matrix form as where F is any N À 1 dimensional column vector and f is also any column vector, defined in the same manner. The inverse discrete Hankel transform (IDHT) is then given by This can also be written in matrix form as We note that the forward and inverse transforms are the same.

Proof
We show the proof for the Y nN formulation, but it proceeds similarly if Y nN is replaced with T nN . Substituting Eq. (52) into the right hand side of (50) gives Switching the order of the summation in Eq. (54) gives The inside summations as indicated in Eq. (55) are recognized as yielding the Diracdelta function, the orthogonality property of Eq. (40) (or Eq. (46) if using T nN ), which in turn yields the desired result. This proves that the DHT given by (50) can be inverted by (52).

Generalized Parseval theorem
Inner products are preserved and thus energies are preserved under the T nN matrix formulation. To see this, consider any two vectors given by the transform The Y nN matrix formulation does not directly preserve inner products: However, under the Y nN formulation, the inner product between g k J nþ1 j nk ð Þ and h k J nþ1 j nk ð Þ is preserved. To see this, we calculate the inner product between them as Making use of the now-present Dirac-delta function, Eq. (58) simplifies to give a modified Parseval relationship In other words, under a DHT using the Y 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 nN , the T nN based DHT is energy preserving, meaning that where the overbar indicates a conjugate transpose and the superscript T indicates a transpose.
For the formulation with Y nN as the transformation kernel, the equivalent expression is It is obvious from Eq. (59) that if we define the "scaled" vector then by straighforward substitution of scaled vectors and their conjugates, it follows that

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 nN is used but all expressions apply equally if Y nN is replaced with T nN .

Transform of Kronecker-Delta function
The discrete counterpart of the Dirac-delta function is the Kronecker-delta function, δ kk 0 . We recall that the DHT as defined above is a matrix transform from a N À 1 dimensional vector to another. The vector δ kk o is interpreted as the vector as having zero entries everywhere except at position k ¼ k 0 (k 0 fixed so δ kk 0 is a vector), or in other words, the k 0 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 The symbol H Á ð Þ 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 Here, f k ⇔ F m denotes a transform pair and Y nN m, k 0 is k 0 th column of the matrix Y nN .

Inverse transform of the Kronecker Delta function
From Eq. (65), we can deduce the vector f k that transforms to the Kroneckerdelta vector δ mm o function. Namely, we take the forward transform of As before, Y nN k, m 0 represents the m 0 th column of the transformation matrix Y nN . From the forward definition of the transform, Eq. (50), the transform of Y n, N k, m 0 is given by where we have used the orthogonality relationship (40). This gives us another DHT pair:

The generalized shift operator
For a one-dimensional Fourier transform, one of the known transform rules is the shift rule, which states that In Eq. (69),f ω ð Þ is the Fourier transform of f x ð Þ, F À1 denotes an inverse Fourier transform and e Àiaω 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 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 x À y ð Þand J n x ð Þ 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 shift k, k o can also be expressed in terms of the original unshifted function f k . Using the definition of F m from Eq. (50) and a dummy change of variable, then Eq. (70) becomes Changing the order of summation gives 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 It then follows that Eq. (72) can be written as 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].

Transform of the generalized shift operator
We now consider the forward DHT transform of the shifted function f shift k, k o . From the definition, the DHT of the shifted function can be found from Changing the order of summation gives 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: 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 m; k o ð Þth entry of the Y nN matrix.

Modulation
We consider the forward DHT of a function "modulated" in the space domain f k ¼ Y nN k, k o g k . Here, the interpretation of f k ¼ Y nN k, k o g k is that the kth entry of the vector g is multiplied by the k; k o ð Þth entry of Y nN for a fixed value of k o . No summation is implied so this is not a dot product; both f k and Y nN k, k o g k are N À 1 vectors. Again, we implement the definition of the forward transform and write g k in terms of its inverse transform Then Eq. (78) becomes Interchanging the order of summation gives 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: 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.

Convolution
We consider the convolution using the generalized shifted function previously defined. The convolution of two functions is defined as 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 But from the orthogonality relationship (40), the summation over k 0 gives the Kronecker delta function, so that Eq. (84) becomes 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 It follows from Eq. (85) that interchanging the roles of g and h will yield the same result, meaning Therefore, it follows that

Multiplication
We now consider the forward transform of a product in the space domain Rearranging gives This gives us yet another transform pair that says that multiplication in the spatial domain is equivalent to convolution in the transform domain: Interchanging the roles of G and H in Eq. (91) demonstrates that convolution in the transform domain also commutes:  (29), it is clear that given a continuous function f r ð Þ evaluated at the discrete points r nk (given by Eq. (31)) in the space domain (1 ≤ k ≤ N À 1), its nth-order Hankel-transform function F ρ ð Þ evaluated at the discrete points ρ nm (given in Eq. (31)) in the frequency domain (1 ≤ m ≤ N À 1), can be approximately given by where α is a scaling factor to be discussed below, and F m ½ ¼ F ρ nm ð Þ, f k ½ ¼ f r nk ð Þ. Similarly, given a continuous function F ρ ð Þ evaluated at the discrete points ρ nm in the frequency domain (1 ≤ m ≤ N À 1), its nth-order inverse Hankel transform f r ð Þ evaluated at the discrete points r nk (1 ≤ k ≤ N À 1), can be approximately given by For both the forward and inverse transforms, α is a scaling factor and α ¼ R 2 where R is the effective space limit and W ρ is the effective band limit (in m À1 ). The scaling factor α 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 r ð Þ ≈ 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 ρ .
The relationship W ρ R ¼ j 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 nN . The determination of j 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 ρ , j nN can be chosen but the third must follow from W ρ R ¼ j 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 2WR ¼ Á .

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 0; R ½ and a Hankel transform of order n, these sampling points are given in the space domain as r nk and frequency domain by ρ nm , given in Eq. (31) and repeated here for convenience It is important to note that as in the case of the computation of the transformation matrix Y 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 ρ nm get closer together). Similarly, making the frequency domain larger (larger W ρ ) implies a tighter sampling density (smaller step size) in the spatial domain. Although j 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 nk ≈ k þ n 2 À Á π.

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 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 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 nN . In order to perform the steps listed above, several Matlab functions have been developed. These functions are shown in Table 1.
Additionally, the matlab script GuidetoDHT.m is included to illustrate the execution of the necessary computational steps.

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 v ð Þ ¼ 20 log 10 f This error function compares the difference between the exact function values f v ð Þ (evaluated from the continuous function) and the function values estimated via the discrete transform, f * ν ð Þ, 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 1 N ∑ N i¼1 f i À f i * was used. The methodology of the testing is given in further detail in [18] and also in the theory paper [3].

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 Creation of Y nN matrix from the zeros Table 1.
Set of available functions.