Set of available functions.
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.
- Hankel transform
- transform rules
- discrete transform
- polar coordinates
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 . 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
where is the
Thus, Hankel transforms take functions in the spatial
2.2 Fourier Bessel series
It is known that functions defined on a finite portion of the real line , can be expanded in terms of a Fourier Bessel series  given by
where the order,
Eqs. (3) and (4) can be considered to be a transform pair where the continuous function is forward-transformed to the discrete vector given in (4). The inverse transform is then the operation which returns if given , 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 contains no frequencies higher than cycles per meter, then it is completely determined by a series of sampling points given by evaluating at where .
where the Fourier Bessel coefficients can be found from Eq. (4) as
Hence, a function bandlimited to can be written as
Eq. (7) states that the samples determine the function completely since (i) is determined by Eq. (7), and (ii) is known if is known. Another way of looking at this is that band-limiting a function to results in information about the original function at samples . 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 are
3.2 Interpolation theorem for a band-limited function
It follows from Eq. (7) that can be found by inverse Hankel transformation to give
From Watson (, p. 134), we have the following result
Eq. (10) gives the formula for interpolating the samples to reconstruct the continuous band-limited function . Each term used to reconstruct the original function consists of the samples of the function at multiplied by a reconstructing function of the form
3.3 Interpretation in terms of a jinc
Eq. (8) states
In other research work , the generalized shift operator indicating a shift of acting on the function 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,
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 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 , or the classic paper reprint ).
For the relatively simple case of the zeroth-order Hankel transform, the inverse Hankel transform of the Boxcar function is given by
The function 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.
Hence, for a zeroth-order Fourier Bessel transform, Eq. (12), the expansion for 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.
3.4 Sampling theorem for a space-limited function
Now consider a space-limited function so that is zero outside of an interval . It then follows that it can be expanded on in terms of a Fourier Bessel series so that
where the Fourier Bessel coefficients can be found from
From Eq. (21), it is evident that the samples determine the function and hence its transform completely. Another way of looking at this is that space limiting a function to implies discretization in spatial frequency space, at frequencies .
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
4.1 Forward transform
We demonstrated above that a band-limited function, with can be written as
Evaluating the previous Eq. (24) at the sampling points (for any integer
For , then , and Eq. (25), summing over infinite
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
Hence, using the argument of “effectively space limited” in the preceding paragraph, we can terminate the series in Eq. (25) at a suitably chosen that ensures the effective space limit. Terminating the series at is the same as assuming that for . Noting that at , the last term in (25) is , then after terminating at
In this case, the truncated sum in Eq. (26) does not represent
4.2 Inverse transform
Concomitantly, we know that for any space-limited function then for , we can write
More specifically, suppose that we follow the logic from the previous section that the function that was bandlimited but also “effectively space-limited” due the truncation of the series in Eq. (25) at
4.3 Intuitive discretization scheme and kernel
The preceding development shows that a “natural,”
The relationship can be used to change from finite frequency domain to a finite space domain and vice-versa. The size of the transform
It is pointed out in  that the zeros of are
Eq. (32) becomes a better approximation to as . The zeros of the cosine function are at odd multiples of . Therefore, an approximation to the Bessel zero, is given by
Using this approximation, then becomes
For larger values of
This is exactly analogous to the corresponding expression for Fourier transforms. Specifically, for temporal Fourier transforms Shannon  argued that “If the function is limited to the time interval
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
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  that for
where represents the
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 . Eq. (37) is exactly true in the limit as and is true for within the limits of computational error . For smaller values of , Eq. (37) holds with the worst case for the smallest value of
5. Transformation matrices
5.1 Transformation matrix
With inspiration from the notation in , and an additional scaling factor of , we define an transformation matrix with the (
In Eq. (39), the superscripts
Furthermore, the orthogonality relationship, Eq. (37), states that
Eq. (40) states that the rows and columns of the matrix are orthonormal and can be written in matrix form as
where is the dimensional identity matrix and we have written the square matrix as . The forward and inverse truncated and discretized transforms given in Eqs. (26) and (29) can be expressed in terms of . The forward transform, Eq. (26), can be written as
Similarly, the inverse transform, Eq. (29), can be written as
5.2 Another choice of transformation matrix
Following the notation in , we can also define a different transformation matrix with the (
In Eq. (44), the superscripts
The orthogonality relationship, Eq. (37), can be written as
Therefore, the matrix is unitary and furthermore orthogonal since the entries are real.
Using the symmetric, orthogonal transformation matrix , the forward transform from Eq. (26) can be written in as
Similarly, the inverse discrete transform of Eq. (29) can be written as
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 or . is closer to the discretized version of the continuous Hankel transform that we hope the discrete version emulates. However, 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
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.
Switching the order of the summation in Eq. (54) gives
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
7. Generalized Parseval theorem
Inner products are preserved and thus energies are preserved under the matrix formulation. To see this, consider any two vectors given by the transform , then
The matrix formulation does not directly preserve inner products:
However, under the formulation, the inner product between and 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 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 , the based DHT is energy preserving, meaning that
where the overbar indicates a conjugate transpose and the superscript
For the formulation with 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
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, is used but all expressions apply equally if is replaced with .
8.1 Transform of Kronecker-Delta function
The discrete counterpart of the Dirac-delta function is the Kronecker-delta function, . We recall that the DHT as defined above is a matrix transform from a dimensional vector to another. The vector is interpreted as the vector as having zero entries everywhere except at position (fixed so is a vector), or in other words, the
The symbol is used to denote the operation of taking the discrete Hankel transform. This gives us our first DHT transform pair of order
Here, denotes a transform pair and is
8.2 Inverse transform of the Kronecker Delta function
From Eq. (65), we can deduce the vector that transforms
As before, represents the
where we have used the orthogonality relationship (40). This gives us another DHT pair:
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
In Eq. (69), is the Fourier transform of , denotes an inverse Fourier transform and 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  (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
where . For a single, fixed value of , then is another vector, with the notation implying a -shifted version of . This generalizes the notion of the shift, usually denoted , which inevitably encounters difficulty when the subscript falls outside of the range . We note that if
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].
8.4 Transform of the generalized shift operator
We now consider the forward DHT transform of the shifted function . 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
We consider the forward DHT of a function “modulated” in the space domain . Here, the interpretation of is that the
and write 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.
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
The right hand side of Eq. (85) is clearly the inverse transform of the product of the transforms . This gives us another transform pair
It follows from Eq. (85) that interchanging the roles of and will yield the same result, meaning
Therefore, it follows that
We now consider the forward transform of a product in the space domain so that
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 and in Eq. (91) demonstrates that convolution in the transform domain also commutes:
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 evaluated at the discrete points (given by Eq. (31)) in the space domain (), its
where is a scaling factor to be discussed below, and , .
Similarly, given a continuous function evaluated at the discrete points in the frequency domain (), its
For both the forward and inverse transforms, is a scaling factor and or equivalently , where
The relationship , derived in the previous sections, holds between the ranges in space and frequency. Choosing determines the dimension (size) of the DHT and determines . The determination of (via choosing
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 and a Hankel transform of order , these sampling points are given in the space domain as and frequency domain by , 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 , the first Bessel zero 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
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
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 Bessel zeros of the Bessel function of order
Generate of sample points (if using the DHT to approximate the continuous transform)
Sample the function (if needed)
Create the 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 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 . 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 |
|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 matrix from the zeros|
Additionally, the matlab script
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 were evaluated.
For the purpose of testing the accuracy of the DHT and IDHT, the dynamic error was used, defined as 
This error function compares the difference between the exact function values (evaluated from the continuous function) and the function values estimated via the discrete transform, , 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 was used. The methodology of the testing is given in further detail in  and also in the theory paper .
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.
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.