Matlab Code for the Discrete Hankel Transform

Natalie Baddour; Ugo Chouinard

(1) Overview

Introduction

There have been many attempts to define a Discrete Hankel Transform (DHT) in the literature, however prior work has focused on proposing methods to approximate the calculation of the continuous Hankel integral, for example as given in [, ]. This stands in stark contrast to the approach taken with the Fourier transform where the Discrete Fourier Transform (DFT) is a transform in its own right, with its own mathematical theory of the manipulated quantities. An additional feature of a carefully derived DFT is that it can be used to approximate the continuous Fourier transform, with relevant sampling and interpolation theories that can be used. Recently, a DHT was proposed as a complete and orthogonal transform that possesses its own mathematical theory, including the standard set of shift, modulation, multiplication and convolution rules []. In addition, this DHT can be used to approximate the continuous Hankel transform in the same manner that the Discrete Fourier transform is known to be able to approximate the continuous Fourier transform.

Overview of the Discrete Hankel Transform

The Continuous Hankel Transform

The forward Hankel transform of order n transforms a function f(r) in the spatial domain to a function F(ρ) in the spatial frequency domain and is given by [, p. 5.6]

(1)

ℍ n {f (r)} = F n (ρ) = ∫ 0 ∞ f (r) r J n (ρ r) d r

M1 \documentclass[10pt]{article} \usepackage{wasysym} \usepackage[substack]{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage[mathscr]{eucal} \usepackage{mathrsfs} \usepackage{pmc} \usepackage[Euler]{upgreek} \pagestyle{empty} \oddsidemargin -1.0in \begin{document} \[ {{\rm{H}}_n}\{ f(r)\} = {F_n}(\rho ) = \int_0^\infty {f(r)r{J_n}(\rho r){\rm{ }}dr} \] \end{document}

where J_n (z) is the nth order Bessel function of the first kind. The inverse transform is given by

(2)

ℍ n − 1 {F n (ρ)} = f (r) = ∫ 0 ∞ ρ F n (ρ) J n (ρ r) d ρ

M2 \documentclass[10pt]{article} \usepackage{wasysym} \usepackage[substack]{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage[mathscr]{eucal} \usepackage{mathrsfs} \usepackage{pmc} \usepackage[Euler]{upgreek} \pagestyle{empty} \oddsidemargin -1.0in \begin{document} \[ {{\rm{H}}_n}^{ - 1}\{ {F_n}(\rho )\} = f(r) = \int_0^\infty {\rho {F_n}(\rho ){J_n}(\rho r){\rm{ }}d\rho } \] \end{document}

More on the continuous transform can be found in [].

Discrete Hankel Transform

The nth order discrete Hankel transform (DHT) proposed in [] is defined as the transformation of the discrete vector f to vector F given by

(3)

F = Y n N f

M3 \documentclass[10pt]{article} \usepackage{wasysym} \usepackage[substack]{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage[mathscr]{eucal} \usepackage{mathrsfs} \usepackage{pmc} \usepackage[Euler]{upgreek} \pagestyle{empty} \oddsidemargin -1.0in \begin{document} \[ {\bf{F}} = {Y^{nN}}{\bf{f}} \] \end{document}

This discrete transform consists of taking an N – 1 vector f and a (N – 1) × (N – 1) square matrix of Hankel order n, Y^nN, to perform the matrix-vector multiplication and obtain the N – 1 DHT vector F. If the DHT as defined in (3) is used to approximate the CHT, then the vector f represents the sampled function to be transformed and the vector F represents the discrete function in the transformed (Hankel) domain. The Y^nN matrix in equation (3) is defined as having the m, k^th entry given by

(4)

Y m, k n N = 2 j n N J n + 1 2 (j n k) J n (j n m j n k j n N) 1 ≤ m, k ≤ N − 1

M4 \documentclass[10pt]{article} \usepackage{wasysym} \usepackage[substack]{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage[mathscr]{eucal} \usepackage{mathrsfs} \usepackage{pmc} \usepackage[Euler]{upgreek} \pagestyle{empty} \oddsidemargin -1.0in \begin{document} \[ Y_{m,k}^{nN} = \frac{2}{{{j_{nN}}J_{n + 1}^2({j_{nk}})}}{J_n}\left( {\frac{{{j_{nm}}{j_{nk}}}}{{{j_{nN}}}}} \right){\rm{ }}1 \le m,k \le N - 1{\rm{ }} \] \end{document}

where j_nk is the k^th zero of the Bessel function of the first kind of order n[]. Properties of the DHT as defined in equation (3) are shown in [].

Since the core of the tested discrete transform is the transformation matrix Y^nN, various properties have to be maintained. One of these properties is that the matrix Y^nN possesses orthogonality properties, where Y^nNY^nN = I. In order to preserve the requisite properties of Y^nN and therefore of the DHT itself, the first Bessel zero used in computing the entries of the Y^nN matrix is the first non-zero value of the Bessel zero of order n. If the Y^nN matrix is not assembled following this rule, the matrix loses its orthogonality property and thus performing the discrete transform leads to improper results. If the DHT is used to approximate a CHT, then this restriction also applies to the discretization of the continuous function, as shall be discussed further below.

The inverse discrete Hankel transform f of the vector F is then given by

(5)

f = Y n N F

M5 \documentclass[10pt]{article} \usepackage{wasysym} \usepackage[substack]{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage[mathscr]{eucal} \usepackage{mathrsfs} \usepackage{pmc} \usepackage[Euler]{upgreek} \pagestyle{empty} \oddsidemargin -1.0in \begin{document} \[ {\bf{f}} = {Y^{nN}}{\bf{F}} \] \end{document}

The discrete forward and inverse Hankel transforms as given in equations (3) and (5) have been shown to possess the standard set of shift, modulation, multiplication and convolution rules. In addition, this DHT can be used to approximate the continuous Hankel transform in the same manner that the Discrete Fourier transform is known to be able to approximate the continuous Fourier transform at certain discrete points.

Discrete Hankel Transform to approximate the Continuous Hankel Transform

Given a continuous function f(r) evaluated at the discrete points r_nk in the space domain (1 ≤ k ≤ N – 1), its nth order Hankel-transform function F(ρ) evaluated at the discrete points ρ_nm (1 ≤ m ≤ N – 1), can be approximately given by []

(6)

F [m] = α ∑ k = 1 N − 1 Y m, k n N f [k] ⇒ F = α Y n N f

M6 \documentclass[10pt]{article} \usepackage{wasysym} \usepackage[substack]{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage[mathscr]{eucal} \usepackage{mathrsfs} \usepackage{pmc} \usepackage[Euler]{upgreek} \pagestyle{empty} \oddsidemargin -1.0in \begin{document} \[ F\left[ m \right] = \alpha \sum\limits_{k = 1}^{N - 1} {Y_{m,k}^{nN}f\left[ k \right]} \;\;\;\; \Rightarrow \;\;\;{\bf{F}} = \alpha {Y^{nN}}{\bf{f}} \] \end{document}

where α is a scaling factor to be discussed below, and F[m] = F(ρ_nm), f[k] = f(r_nk).

Conversely, 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 function f(r) evaluated at the discrete points r_nk (1 ≤ k ≤ N – 1), can be approximately given by

(7)

f [k] = 1 α ∑ m = 1 N − 1 Y m, k n N F [m] ⇒ f = 1 α Y n N F

M7 \documentclass[10pt]{article} \usepackage{wasysym} \usepackage[substack]{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage[mathscr]{eucal} \usepackage{mathrsfs} \usepackage{pmc} \usepackage[Euler]{upgreek} \pagestyle{empty} \oddsidemargin -1.0in \begin{document} \[ f[k] = \frac{1}{\alpha }\sum\limits_{m = 1}^{N - 1} {Y_{m,k}^{nN}F[m]} \;\;\; \Rightarrow \;\;\;\;{\bf{f}} = \frac{1}{\alpha }{Y^{nN}}{\bf{F}} \] \end{document}

For both the forward and inverse transforms, α is a scaling factor which depends on the function properties and shall be discussed further below. The choice of discretization points r_nk and ρ_nm is also discussed below. The full theory of the discrete Hankel transform is given in [].

Discretization Points

In order to properly use the discrete transform to approximate the continuous transform, a function has to be discretized at specific sampling 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

(8)

r n k = j n k j n N R for 1 ≤ k ≤ N − 1

M8 \documentclass[10pt]{article} \usepackage{wasysym} \usepackage[substack]{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage[mathscr]{eucal} \usepackage{mathrsfs} \usepackage{pmc} \usepackage[Euler]{upgreek} \pagestyle{empty} \oddsidemargin -1.0in \begin{document} \[ {r_{nk}} = \frac{{{j_{nk}}}}{{{j_{nN}}}}R{\rm{ for }}1 \le k \le N - 1 \] \end{document}

For the finite frequency domain range [0, W_ρ] and a Hankel transform of order n, the sampling points are given by

(9)

ρ n m = j n m j n N W ρ for 1 ≤ m ≤ N − 1

M9 \documentclass[10pt]{article} \usepackage{wasysym} \usepackage[substack]{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage[mathscr]{eucal} \usepackage{mathrsfs} \usepackage{pmc} \usepackage[Euler]{upgreek} \pagestyle{empty} \oddsidemargin -1.0in \begin{document} \[ {\rho _{nm}} = \frac{{{j_{nm}}}}{{{j_{nN}}}}{W_\rho }{\rm{ for }}1 \le m \le N - 1 \] \end{document}

It is important to note that as in the case of the computation of the transformation matrix Y^nN, the first Bessel zero j_n₁ used in computing the discretization points is the first non-zero value.

The relationship $W ρ = j n N R$ M10 \documentclass[10pt]{article} \usepackage{wasysym} \usepackage[substack]{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage[mathscr]{eucal} \usepackage{mathrsfs} \usepackage{pmc} \usepackage[Euler]{upgreek} \pagestyle{empty} \oddsidemargin -1.0in \begin{document} \[ {W_\rho } = \frac{{{j_{nN}}}}{R} \] \end{document} , derived in [], 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.

Scaling Factor

The scaling factor α necessary for using the DHT to approximate the CHT depends on whether the function is space-limited or band-limited. Since it might be hard to determine if a function is space or band limited, the concept of effective limit is introduced. Therefore, a function defined as being “effectively limited in space by R” means that if r > R, then as r → ∞, f(r) → 0. 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 would be denoted by W_ρ. The conditions and corresponding scaling factors are listed in Table 1.

Table 1

Scaling factor under various conditions.

	Condition	Scaling Factor

1	Space-limited function	$α = R 2 j n N$ M11 \documentclass[10pt]{article} \usepackage{wasysym} \usepackage[substack]{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage[mathscr]{eucal} \usepackage{mathrsfs} \usepackage{pmc} \usepackage[Euler]{upgreek} \pagestyle{empty} \oddsidemargin -1.0in \begin{document} \[ \alpha = \frac{{{R^2}}}{{{j_{nN}}}} \] \end{document}
2	Frequency-limited function	$α = j n N W ρ 2$ M12 \documentclass[10pt]{article} \usepackage{wasysym} \usepackage[substack]{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage[mathscr]{eucal} \usepackage{mathrsfs} \usepackage{pmc} \usepackage[Euler]{upgreek} \pagestyle{empty} \oddsidemargin -1.0in \begin{document} \[ \alpha = \frac{{{j_{nN}}}}{{W_\rho ^2}} \] \end{document}

The detailed derivation of these scaling factors was shown in []. It can be observed that the scaling factors for the space-limited or frequency limited cases can be derived from each other via W_ρR = J_nN.

Implementation and architecture

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

Calculations of N Bessel zeros of the Bessel function of order n
Generation of N sample points (if using the DHT to approximate the continuous transform)
Discretization of the function (if needed)
Creation of the Y^nN transformation matrix
Performing the matrix-function multiplication

The steps are the same regardless of if the function is in the space or frequency domain and are summarized in Figure 1.

Figure 1

Steps for evaluation of DHT.

Furthermore, the Y^nN transformation matrix is used for both the forward and inverse transform. Steps 2–3 only need to be performed if the function (vector) to be transformed is not already given as a set of discrete points. In the case of a continuous function in either the space or frequency domain, it is important to use the sampling points as proposed in equations (8), (9) and then to discretize the continuous function by evaluating at these points. 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 listed in Table 2.

Table 2

Set of available functions.

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 []
freqSampler	freqSampler(R,zeros)	Creation of sample points in the frequency domain (eq. (9))
spaceSampler	spaceSampler(R,zeros)	Creation of sample points in the space domain (eq. (8))
YmatrixAssembly	YmatrixAssembly(n,N,zeros)	Creation of Y^nN matrix (eq. (4)) from the zeros

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

Quality control

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 []

(10)

e (v) = 20 log 10 [| f (v) − f * (ν) | max | f * (v) |]

M13 \documentclass[10pt]{article} \usepackage{wasysym} \usepackage[substack]{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage[mathscr]{eucal} \usepackage{mathrsfs} \usepackage{pmc} \usepackage[Euler]{upgreek} \pagestyle{empty} \oddsidemargin -1.0in \begin{document} \[ e\left( v \right) = 20{\log _{10}}\left[ {\frac{{\left| {f(v) - {f^*}(\nu )} \right|}}{{\max \left| {{f^*}(v)} \right|}}} \right] \] \end{document}

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*(v), scaled with the maximum value of the discretely estimated samples. Equation (10) can be used to evaluate the computation of either forward or inverse Hankel transform via the DHT/IDHT and compared with known continuous Hankel relationships. 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. This is performed by consecutive forward and then inverse transformation in order to verify that the transforms themselves do not add errors. For this evaluation, the average absolute error $1 N ∑ i = 1 N | f i − f i * |$ M14 \documentclass[10pt]{article} \usepackage{wasysym} \usepackage[substack]{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage[mathscr]{eucal} \usepackage{mathrsfs} \usepackage{pmc} \usepackage[Euler]{upgreek} \pagestyle{empty} \oddsidemargin -1.0in \begin{document} \[ \frac{1}{N}\sum\limits_{i = 1}^N {\left| {{f_i} - {f_i}^*} \right|} \] \end{document} was used.

The methodology of the testing is given in further detail in [] and also in the theory paper, [].

(2) Availability

Operating system

Windows XP and higher.

Programming language

Matlab

Additional system requirements

If using Matlab 7, minimum system requirements are 512MB of RAM (1024 MB recommended) and 460MB of hard disk space. System requirements for Matlab R2014 require 1024 MB (At least 2048 MB recommended) or RAM and 1 GB for MATLAB only, 3–4 GB for a typical installation. This code should run on older versions of Matlab although it was tested using Matlab R2014.

Dependencies

Matlab software by mathworks.

List of contributors

Ugo Chouinard, Natalie Baddour, Greg von Winckel

Software location

Archive (e.g. institutional repository, general repository) (required)

Name: Figshare

Persistent identifier: http://dx.doi.org/10.6084/m9.figshare.1453205

Licence: BSD

Publisher: Natalie Baddour

Date published: 18/06/15

Code repository (e.g. SourceForge, GitHub etc.) (required)

Name: GitHub

Persistent identifier: https://github.com/uchouinard/DiscreteHankelTransform

Licence: BSD

Publisher: Natalie Baddour

Date published: 07/12/16

Language

English

(3) Reuse potential

The Discrete Hankel Transform is applicable to many areas of scientific computation and potential reuse could be very high. Further details on applications can be found in [].

[B1] Leutenneger, M (). Hankel transform – File Exchange – MATLAB Central,

[B2] Wyatt, A (). Hankel Transform – File Exchange – MATLAB Central,

[B3] Baddour, N and Chouinard, U (2015). Theory and operational rules for the discrete Hankel transform J. Opt. Soc. Am. A Apr. 201532(4): 611–622, DOI: https://doi.org/10.1364/JOSAA.32.000611

[B4] Debnath, L and Bhatta, D (2007). Integral transforms and their applications. 2nd ed. Boca Raton: Chapman & Hall/CRC.

[B5] von Winckel, G (). Bessel Function Zeros – File Exchange – MATLAB Central [Online]. Available at: http://www.mathworks.com/matlabcentral/fileexchange/6794-bessel-function-zeros. [Accessed: 06-Jun-2015].

[B6] Guizar-Sicairos, M and Gutiérrez-Vega, J C (2004). Computation of quasi-discrete Hankel transforms of integer order for propagating optical wave fields J. Opt. Soc. Am. A Jan. 200421(1): 53–58, DOI: https://doi.org/10.1364/JOSAA.21.000053

[B7] Chouinard, U (2015). Numerical Simulations for the Discrete Hankel Transform, B.A.Sc. thesis. University of Ottawa, Ottawa, Canada, DOI: https://doi.org/10.13140/RG.2.1.1208.3364

Journal of Open Research Software

Software Metapapers