Open access peer-reviewed chapter

# The Discrete Hankel Transform

Written By

Submitted: November 7th, 2018 Reviewed: January 14th, 2019 Published: February 18th, 2019

DOI: 10.5772/intechopen.84399

From the Edited Volume

## Fourier Transforms

Edited by Goran S. Nikoli? and Dragana Z. Markovi?-Nikoli?

Chapter metrics overview

View Full Metrics

## 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.

### Keywords

• Fourier-Bessel
• Hankel transform
• transform rules
• discrete transform
• polar coordinates

## 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ρof the function frof a real variable, r0, is defined by the integral [4]

Fρ=Hnfr=0frJnρrrdr,E1

where Jnzis the nth-order Bessel function of the first kind. If nis real and n>1/2, the transform is self-reciprocating and the inversion formula is given by

fr=0FρJnρrρdρ.E2

Thus, Hankel transforms take functions in the spatial rdomain and transform them to functions in the spatial frequency ρdomain frFρ. The notationis 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 0R, can be expanded in terms of a Fourier Bessel series [5] given by

fr=k=1fkJnjnkrR,E3

where the order, n, of the Bessel function is arbitrary and jnkdenotes the kth root of the nth Bessel function Jnz. The Fourier Bessel coefficients fkof the function frare given by

fk=2R2Jn+12jnk0RfrJnjnkrRrdr.E4

Eqs. (3) and (4) can be considered to be a transform pair where the continuous function fris forward-transformed to the discrete vector fkgiven in (4). The inverse transform is then the operation which returns frif given fk, 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 frcontains no frequencies higher than Wcycles per meter, then it is completely determined by a series of sampling points given by evaluating frat r=jnkWρwhere Wρ=2πW.

Proof:suppose that a function fris band-limited in the frequency Hankel domain so that its spectrum Fρis zero outside of an interval 02πW. The interval is written in this form since Wwould 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 s1or m1. Let us define Wρ=2πW. Since the Hankel transform Fρis defined on a finite portion of the real line 0Wρ, it can be expanded in terms of a Fourier Bessel series as

Fρ=k=1FkJnjnkρWρ.E5

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

Fk=2Wρ2Jn+12jnk0WρFρJnjnkρWρρ=2Wρ2Jn+12jnkfjnkWρE6

In (6), we have used the fact that frcan 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 0Wρcan be written as

Fρ=k=12Wρ2Jn+12jnkfjnkWρJnjnkρWρρ<Wρ0ρWρE7

Eq. (7) states that the samples fjnkWρdetermine the function frcompletely since (i) Fρis determined by Eq. (7), and (ii) fris known if Fρis known. Another way of looking at this is that band-limiting a function to 0Wρresults in information about the original function at samples rnk=jnkWρ. 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 Jnzare almostevenly spaced at intervals of πand that the spacing becomes exactly πin the limit as z. To determine the (bandlimited) function frcompletely, we need to sample the function at fjnkWρ=fjnk2πWand these samples are (eventually) multiples of π(2πW)=1(2W)apart, which is consistent with the standard Shannon sampling theorem which requires samples at multiples of 1(2W)[6].

### 3.2 Interpolation theorem for a band-limited function

It follows from Eq. (7) that frcan be found by inverse Hankel transformation to give

fr=0Wρk=12Wρ2Jn+12jnkfjnkWρJnjnkρWρJnρrρdρ=k=12Wρ2Jn+12jnkfjnkWρ0WρJnjnkρWρJnρrρdρ.E8

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

JnαzJ¯nβzzdz=zαJn+1αzJ¯nβzβJnαzJ¯n+1βzα2β2E9

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

fr=k=1fjnkWρ2jnkJn+1jnkJnrWρjnk2r2Wρ2E10

Eq. (10) gives the formula for interpolating the samples fjnkWρto reconstruct the continuous band-limited function fr. Each term used to reconstruct the original function frconsists of the samples of the function frat r=jnkWρmultiplied by a reconstructing function of the form

2jnkJn+1jnkJnrWρjnk2r2Wρ2.E11

### 3.3 Interpretation in terms of a jinc

Eq. (8) states

fr=k=12Wρ2Jn+12jnkfjnkWρ0WρJnjnkρWρJnρrρdρE12

In other research work [8], the generalized shift operator Rr0indicating a shift of r0acting on the function frhas been defined by the formula

Rr0fr=0FρJnρr0Jnρrρdρ.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 jnkWρ. More explicitly,

0WρJnjnkρWρJnρrρdρ=RjnkWρgeneralizedshift ofjnkWρ   0ΠWρρJnρrρdρinverse Hankel transform ofΠWρρ=jnkWρ2jnk2rWρ2Jn+1jnkJnrWρE14

where

ΠWρρ=10ρWρ0otherwiseE15

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/21/2and 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

0ΠWρρJ0ρrρdρ=0WρJ0ρrρdρ=WρrJ1Wρr=Wρ2J1WρrWρr.E16

The function 2J1rris 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

Rj0kWρ2J1WρrWρr=2j0kJ1j0kj0k2rWρ2J0rWρ.E17

Hence, for a zeroth-order Fourier Bessel transform, Eq. (12), the expansion for frreads

fr=k=1fj0kWρ1J12j0k2j0kJ1j0kJ0rWρj0k2rWρ2=Rj0kWρ2J1WρrWρrE18

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 frso that fris zero outside of an interval 0R. It then follows that it can be expanded on 0Rin terms of a Fourier Bessel series so that

fr=k=1fkJnjnkrR,E19

where the Fourier Bessel coefficients can be found from

fk=2R2Jn+12jnk0RfrJnjnkrRrdr=2R2Jn+12jnkFjnkR.E20

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

fr=k=12R2Jn+12jnkFjnkRJnjnkrRr<R0rRE21

From Eq. (21), it is evident that the samples FjnkRdetermine the function frand hence its transform Fρcompletely. Another way of looking at this is that space limiting a function to 0Rimplies discretization in spatial frequency space, at frequencies ρnk=jnkR.

### 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

Fρ=0frJnρrrdr=k=12R2Jn+12jnkFjnkR0RJnjnkrRJnρrrdrE22

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

Fρ=k=1FjnkR2jnkJn+1jnkJnρRjnk2ρR2E23

Eq. (23) gives the formula for interpolating the samples FjnkRto give the continuous function Fρ.

## 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 ispossible 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 ρ<Wρ=2πWcan be written as

Fρ=k=12Wρ2Jn+12jnkfjnkWρJnjnkρWρ.E24

Evaluating the previous Eq. (24) at the sampling points ρnm=jnmWρjnN(for any integer N) gives for m<N

FjnmWρjnN=k=12Wρ2Jn+12jnkfjnkWρJnjnkWρjnmWρjnN   m<N.E25

For m<N, then ρnm=jnmWρjnN<Wρ, and Eq. (25), summing over infinite k,is exact. For mN, then ρnm=jnmWρjnNWρ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 effectivelyspace limited, this means that there exists an integer Nfor which fjnkWρ0for 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 limkfjnkWρ=0, which means that for any arbitrarily small ε, there must exist an integer Nfor which fjnkWρ<ε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 Nthat ensures the effective space limit. Terminating the series at k=Nis the same as assuming that fr0for r>R=jnNWρ. Noting that at k=N, the last term in (25) is JnjnNjnmjnN=Jnjnm=0, then after terminating at N, Eq. (25) becomes for m=1..N1

FjnmWρjnN=k=1N12Wρ2Jn+12jnkfjnkWρJnjnkjnmjnN.E26

In this case, the truncated sum in Eq. (26) does not represent Fρnmexactlydue to the truncation at Nterms, 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 fjnkWρ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

fr=m=12R2Jn+12jnmFjnmRJnjnmrR.E27

More specifically, suppose that we follow the logic from the previous section that the function frthat was bandlimited but also “effectively space-limited” due the truncation of the series in Eq. (25) at N. In that case then R=jnNWρand the band-limit implies that Fρ=0for ρ>Wρ. Following this logic and using R=jnNWρ, then Eq. (27) becomes

fr=m=1N12Wρ2jnN2Jn+12jnmFjnmWρjnNJnjnmWρjnNrE28

where we truncated the series in Eq. (28) at Nby using the fact that Fρ=0for ρWρto deduce that FjnmWρjnN=0for mN. Evaluating (28) at rnk=jnkRjnN=jnkWρgives for k=1..N1

fjnkWρ=m=1N1 2Wρ2jnN2Jn+12jnmFjnmWρjnNJnjnmjnkjnN.E29

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

FjnmWρjnN=k=1N12Wρ2Jn+12jnkfjnkWρJnjnkjnmjnN.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 0Rand finite frequency space 0Wρis given by

rnk=jnkWρ=jnkRjnNρnk=jnkR=jnkWρjnNk=1N1.E31

The relationship WρR=jnNcan 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ρR=jnN.

It is pointed out in [10] that the zeros of Jnzare almostevenly 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

Jnz2πzcoszn+12π2E32

Eq. (32) becomes a better approximation to Jnzas z. The zeros of the cosine function are at odd multiples of π/2. Therefore, an approximation to the Bessel zero, jnkis given by

jnk2k+n12π2.E33

Using this approximation, then WρR=jnNbecomes

2πWR=jnN2N+n12π2E34

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

2WRN+n2E35

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 Tand the samples are spaced 1/(2 W) seconds apart (where Wis the bandwidth), there will be a total of 2WTsamples 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 Tif, and only if, all the samples outside this interval are exactly zero. Then we can say that any function limited to the bandwidth Wand the time interval Tcan be specified by giving N=2WTnumbers”. Following this line of thinking, Eq. (35) states that for an nth-order Hankel transform, any function limited to the bandwidth Wand the space interval Rcan be specified by giving N=2WRn/2numbers and it will certainly be true that specifying N=2WRnumbers 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

2Jn+12jnkJnjnkjnmjnN.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 1m,iN1

k=1N14JnjnmjnkjnNJnjnkjnijnNJn+12jnk=jnN2Jn+12jnmδmiE37

where jnmrepresents the mth zero of the nth-order Bessel function Jnx, and δmiis the Kronecker delta function, defined as

δmn=1ifm=n0otherwise.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 Nand is true for N>30within the limits of computational error ±107. For smaller values of N, Eq. (37) holds with the worst case for the smallest value of Ngiving ±103.

## 5. Transformation matrices

### 5.1 Transformation matrix

With inspiration from the notation in [11], and an additional scaling factor of 1/jnN, we define an N1×N1transformation matrix with the (m,k)th entry given by

Ym,knN=2jnNJn+12jnkJnjnmjnkjnN1m,kN1.E39

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

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

k=1N1Yi,knNYk,mnN=k=1N14JnjnijnkjnNJnjnkjnmjnNjnN2Jn+12jnkJn+12jnm=δim.E40

Eq. (40) states that the rows and columns of the matrix Ym,knNare orthonormal and can be written in matrix form as

YnNYnN=Ι,E41

where Iis the N1dimensional identity matrix and we have written the N1square matrix Ym,knNas YnN. The forward and inverse truncated and discretized transforms given in Eqs. (26) and (29) can be expressed in terms of YnN. The forward transform, Eq. (26), can be written as

Fρnm=jnNWρ2k=1N1Ym,knNfrnk.E42

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

frnk=Wρ2jnNm=1N1Yk,mnNFρnm.E43

### 5.2 Another choice of transformation matrix

Following the notation in [12], we can also define a different N1×N1transformation matrix with the (m,k)th entry given by

Tm,knN=2JnjnmjnkjnNJn+1jnmJn+1jnkjnN1m,kN1.E44

In Eq. (44), the superscripts nand Nrefer to the order of the Bessel function and the dimension of the space that are being considered, respectively. The subscripts mand krefer to the (m,k)th entry of the matrix. From (39), it can be seen that Tm,knN=Tk,mnNso that TnNis a real, symmetric matrix. The relationship between the Tm,knNand Ym,knNmatrices is given by

Tm,knNJn+1jnmJn+1jnk=Ym,knN.E45

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

k=1N14JnjnmjnkjnNJnjnkjnijnNJn+12jnmJn+12jnkjnN2=k=1N1Tm,knNTk,inN=δmi.E46

Eq. (40) states that the rows and columns of the matrix TnNare orthonormal so that TnNis an orthogonal matrix. Furthermore, TnNis also symmetric. Eq. (46) can be written in matrix form as

TnNTnN=TnNTnNT=I.E47

Therefore, the TnNmatrix is unitary and furthermore orthogonal since the entries are real.

Using the symmetric, orthogonal transformation matrix TnN, the forward transform from Eq. (26) can be written in as

Fρnm=R2jnNk=1N1Tm,knNJn+1jnmJn+1jnkfrnk=jnNWρ2k=1N1Tm,knNJn+1jnmJn+1jnkfrnkE48

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

frnk=jnNR2m=1N1Tk,mnNJn+1jnkJn+1jnmFρnm=Wρ2jnNm=1N1Tk,mnNJn+1jnkJn+1jnmFρnm.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 Ym,knNor Tm,knN. Ym,knNis closer to the discretized version of the continuous Hankel transform that we hope the discrete version emulates. However, Tm,knNis 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:

Fm=k=1N1Ym,knNfkorFm=k=1N1Tm,knNfk.E50

Here, the transform is of anyN1dimensional vector fkto anyN1dimensional vector Fmfor the integers m,kwhere 1m,k<N1. This can be written in matrix form as

F=YnNforF=TnNfE51

where Fis anyN1dimensional column vector and fis also anycolumn vector, defined in the same manner.

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

fk=m=1N1Yk,mnNFmorfk=m=1N1Tk,mnNFm.E52

This can also be written in matrix form as

f=YnNForf=TnNF.E53

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

Proof

We show the proof for the YnNformulation, but it proceeds similarly if YnNis replaced with TnN. Substituting Eq. (52) into the right hand side of (50) gives

k=1N1Ym,knNfk=k=1N1Ym,knNp=1N1Yk,pnNFp.E54

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

p=1N1k=1N1Ym,knNYk,pnNδmpFp=p=1N1δmpFp=FmE55

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 TnN), 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 TnNmatrix formulation. To see this, consider any two vectors given by the transform g=TnNG, h=TnNHthen

gTh=TnNGTTnNH=GTTnNTTnN=IH=GTH.E56

The YnNmatrix formulation does not directly preserve inner products:

gTh=YnNGTYnNH=GTYnNTYnNH.E57

However, under the YnNformulation, the inner product between gkJn+1jnkand hkJn+1jnkis preserved. To see this, we calculate the inner product between them as

k=1N1gkJn+1jnkhkJn+1jnk=k=1N11Jn+12jnkp=1N1Yk,pnNGpq=1N1Yk,qnNHq=p=1N1q=1N11Jn+12jnpk=1N14jnkjnpjnNJnjnkjnqjnNjnN2Jn+12jnkJn+12jnqδpqHqGpE58

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

k=1N1gkJn+1jnkhkJn+1jnk=p=1N1HpJn+1jnpGpJn+1jnp.E59

In other words, under a DHT using the YnNmatrix, 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 TnN, the TnNbased DHT is energy preserving, meaning that

F¯TF=f¯Tf.E60

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

For the formulation with YnNas the transformation kernel, the equivalent expression is

F¯TF=YnNf¯TYnNf=f¯TYnNTYnNf.E61

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

fkScaled=fkJn+1jnkandFpScaled=FpJn+1jnp,E62

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

FScaled¯TFScaled=fScaled¯TfScaled.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, YnNis used but all expressions apply equally if YnNis replaced with TnN.

### 8.1 Transform of Kronecker-Delta function

The discrete counterpart of the Dirac-delta function is the Kronecker-delta function, δkk0. We recall that the DHT as defined above is a matrix transform from a N1dimensional vector to another. The vector δkkois interpreted as the vector as having zero entries everywhere except at position k=k0(k0fixed so δkk0is a vector), or in other words, the k0th column of the N1sized identity matrix. The DHT of the Kronecker-delta can be found from the definition of the forward transform via

Hδkko=k=1N1Ym,kδkko=Ym,k0nNE64

The symbol His used to denote the operation of taking the discrete Hankel transform. This gives us our first DHT transform pair of order ndimension N1, and we denote this relationship as

δkkoYm,k0n,NE65

Here, fkFmdenotes a transform pair and Ym,k0nNis k0th column of the matrix YnN.

### 8.2 Inverse transform of the Kronecker Delta function

From Eq. (65), we can deduce the vector fkthat transforms tothe Kronecker-delta vector δmmofunction. Namely, we take the forward transform of

fk=Yk,m0n,N.E66

As before, Yk,m0nNrepresents the m0th column of the transformation matrix YnN. From the forward definition of the transform, Eq. (50), the transform of Yk,m0n,Nis given by

Fm=k=1N1Ym,knNfk=k=1N1Ym,knNYk,m0nN=δmm0,E67

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

Yk,m0n,Nδmmo.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

fxa=F1eiaωf̂ω=12πeiaωf̂ωeiωt.E69

In Eq. (69), f̂ωis the Fourier transform of fx, F1denotes an inverse Fourier transform and eiaω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

fk,koshift=p=1N1Yk,pnNYp,konNFpshift in Hankeldomain,E70

where 1k,koN1. For a single, fixed value of ko, then fk,koshiftis another N1vector, with the notation fk,koshiftimplying a k0-shifted version of fk. This generalizes the notion of the shift, usually denoted fkko, which inevitably encounters difficulty when the subscript kkofalls outside of the range 1N1. We note that if allpossible shifts koare considered, then fk,koshiftis a N1square matrix (in other words, a twodimensional structure), whereas the original un-shifted fkis an N1vector. For the discrete Fourier transform, when the shifted subscript kkofalls 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 Jnxyand Jnxis the reason that the continuous Hankel transform does not have a convolution-multiplication rule [13]. Thus, the notation fkkowould 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 kand ko.

The shifted function fk,koshiftcan also be expressed in terms of the original un-shifted function fk. Using the definition of Fmfrom Eq. (50) and a dummy change of variable, then Eq. (70) becomes

fk,koshift=p=1N1Yk,pnNYp,konNFp=p=1N1Yk,pnNYp,konNm=1N1Yp,mnNfm.E71

Changing the order of summation gives

fk,koshift=p=1N1Yk,pnNYp,konNFp=m=1N1p=1N1Yk,pnNYp,konNYp,mnNshift operatorfm.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

Sk,ko,mnN=p=1N1Yk,pnNYp,konNYp,mnN.E73

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

fk,koshift=m=1N1Sk,ko,mnNfm.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 fk,koshift. From the definition, the DHT of the shifted function can be found from

k=1N1Ym,knNfk,koshift=k=1N1Ym,knNp=1N1Yk,pnNYp,konNFp.E75

Changing the order of summation gives

p=1N1k=1N1Ym,knNYk,pnN=δmpYp,konNFp=p=1N1δmpYp,konNFp=Ym,konNFm.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:

fk,koshiftYm,konNFm.E77

Note that Eq. (77) does notimply a summation over the mindex. For a fixed value of koon the left hand side, the corresponding transformed value of Fmis multiplied by the mkoth entry of the YnNmatrix.

### 8.5 Modulation

We consider the forward DHT of a function “modulated” in the space domain fk=Yk,konNgk. Here, the interpretation of fk=Yk,konNgkis that the kth entry of the vector gis multiplied by the kkoth entry of YnNfor a fixed value of ko. No summation is implied so this is not a dot product; both fkand Yk,konNgkare N1vectors. Again, we implement the definition of the forward transform

k=1N1Ym,knNfk=k=1N1Ym,knNYk,konNgk,E78

and write gkin terms of its inverse transform

gk=p=1N1Yk,pnNGp.E79

Then Eq. (78) becomes

k=1N1Ym,knNfk=k=1N1Ym,knNYk,konNgk=k=1N1Ym,knNYk,konNp=1N1Yk,pnNGp.E80

Interchanging the order of summation gives

p=1N1k=1N1Ym,knNYk,konNYk,pnNshift operatorGp=Gm,koshift.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:

Yk,konNgkGm,koshift.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

fk=ghk=k0=1N1gkohk,koshift.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

fk=k0=1N1gkohk,koshift=k0=1N1q=1N1Yko,qnNGqgkop=1N1Yk,pnNYp,konNHphk,koshift=q=1N1p=1N1k0=1N1Yp,konNYko,qnN=δpqYk,pnNHpGq.E84

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

ghk=k0=1N1gkohk,koshift=q=1N1p=1N1δpqYk,pnNHpGq=p=1N1Yk,pnNHpGpE85

The right hand side of Eq. (85) is clearly the inverse transform of the product of the transforms HpFp. This gives us another transform pair

ghk=k0=1N1gkohk,koshiftHmGm.E86

It follows from Eq. (85) that interchanging the roles of gand hwill yield the same result, meaning

k0=1N1gk,koshifthko=p=1N1Yk,pnNGpHp.E87

Therefore, it follows that

hgk=k0=1N1gk,koshifthko=k0=1N1gkohk,koshift=ghk.E88

### 8.7 Multiplication

We now consider the forward transform of a product in the space domain fk=gkhkso that

k=1N1Ym,knNgkhk=k=1N1Ym,knNq=1N1Yk,qnNGqgk p=1N1Yk,pnNHphk.E89

Rearranging gives

k=1N1Ym,knNgkhk=q=1N1Gqp=1N1k=1N1Ym,knNYk,qnNYk,pnNshift operatorHp=q=1N1GqHm,qshift=GHm.E90

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

gkhkq=1N1GqHm,qshift=GHm.E91

Interchanging the roles of Gand Hin Eq. (91) demonstrates that convolution in the transform domain also commutes:

GHm=q=1N1GqHm,qshift=q=1N1Gm,qshiftHq=HGm.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 frevaluated at the discrete points rnk(given by Eq. (31)) in the space domain (1kN1), its nth-order Hankel-transform function Fρevaluated at the discrete points ρnm(given in Eq. (31)) in the frequency domain (1mN1), can be approximately given by

Fm=αk=1N1Ym,knNfkF=αYnNfE93

where αis a scaling factor to be discussed below, and Fm=Fρnm, fk=frnk.

Similarly, given a continuous function Fρevaluated at the discrete points ρnmin the frequency domain (1mN1), its nth-order inverse Hankel transform frevaluated at the discrete points rnk(1kN1), can be approximately given by

fk=1αm=1N1Ym,knNFmf=1αYnNFE94

For both the forward and inverse transforms, αis a scaling factor and α=R2jnNor equivalently α=jnNWρ2, where Ris 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 fr0for all r>R. In other words, the function can be made as close to zero as desired by selecting an Rthat 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=jnN, derived in the previous sections, holds between the ranges in space and frequency. Choosing Ndetermines the dimension (size) of the DHT and determines jnN. The determination of jnN(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ρ,jnNcan be chosen but the third must follow from WρR=jnN. 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=WρRπN+n2.

### 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 0Rand a Hankel transform of order n, these sampling points are given in the space domain as rnkand frequency domain by ρnm, given in Eq. (31) and repeated here for convenience

rnk=jnkWρ=jnkRjnNρnm=jnmR=jnmWρjnNk,m=1N1E95

It is important to note that as in the case of the computation of the transformation matrix YnN, the first Bessel zero jn1used 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 ρnmget closer together). Similarly, making the frequency domain larger (larger Wρ) implies a tighter sampling density (smaller step size) in the spatial domain. Although jnmare not equispaced, they are nearly so for higher values of mand for purposes of developing quick intuitions on ideas such as sampling density, if is convenient to approximately think of jnkk+n2π.

### 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 NBessel zeros of the Bessel function of order n

• Generate of Nsample points (if using the DHT to approximate the continuous transform)

• Sample the function (if needed)

• Create the YnNtransformation matrix

• Perform the matrix-function multiplication

The steps are the same regardless if the function is in the space or frequency domain. Furthermore, the YnNtransformation 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 YnN. In order to perform the steps listed above, several Matlab functions have been developed. These functions are shown in Table 1.

Function nameCalling sequenceDescription
besselzerobesselzero(n,k,kind)Calculation of kBessel zeros of the nth-order Bessel function of kind—developed in [17]
freqSamplerfreqSampler(R,zeros)Creation of sample points in the frequency domain (Eq. (95))
spaceSamplerspaceSampler(R,zeros)Creation of sample points in the space domain (Eq. (95))
YmatrixAssemblyYmatrixAssembly(n,N,zeros)Creation of YnNmatrix from the zeros

### Table 1.

Set of available functions.

Additionally, the matlab script GuidetoDHT.mis 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 nwere evaluated.

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

ev=20log10fvfνmaxfvE96

This error function compares the difference between the exact function values fv(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 1Ni=1Nfifiwas 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.

## References

1. 1. Baddour N. Operational and convolution properties of two-dimensional Fourier transforms in polar coordinates. Journal of the Optical Society of America. A. 2009;26(8):1767-1777
2. 2. Baddour N. Operational and convolution properties of three-dimensional Fourier transforms in spherical polar coordinates. Journal of the Optical Society of America. A. 2010;27(10):2144-2155
3. 3. Baddour N, Chouinard U. Theory and operational rules for the discrete Hankel transform. JOSA A. 2015;32(4):611-622
4. 4. Piessens R. The Hankel transform. In: The Transforms and Applications Handbook. Vol. 2. Boca Raton: CRC Press; 2000. pp. 9.1-9.30
5. 5. Schroeder J. Signal processing via Fourier-Bessel series expansion. Digital Signal Processing. 1993;3(2):112-124
6. 6. Shannon CE. Communication in the presence of noise. Proceedings of the IEEE. 1998;86(2):447-457
7. 7. Watson GN. A Treatise on the Theory of Bessel Functions. Cambridge, UK: Cambridge University Press; 1995
8. 8. Levitan BM. Generalized displacement operators. In: Encyclopedia of Mathematics. Heidelberg: Springer. p. 2002
9. 9. Shannon CE. Communication in the presence of noise. Proceedings of the IRE. 1949;37(1):10-21
10. 10. Arfken GB. Mathematical Methods for Physicists. 6th ed. Boston: Elsevier; 2005
11. 11. Johnson HF. An improved method for computing a discrete Hankel transform. Computer Physics Communications. 1987;43(2):181-202
12. 12. Guizar-Sicairos M, Gutiérrez-Vega JC. Computation of quasi-discrete Hankel transforms of integer order for propagating optical wave fields. Journal of the Optical Society of America. A. 2004;21(1):53-58
13. 13. Baddour N. Application of the generalized shift operator to the Hankel transform. Springerplus. 2014;3(1):246
14. 14. Belhadj M, Betancor JJ. Hankel convolution operators on entire functions and distributions. Journal of Mathematical Analysis and Applications. 2002;276(1):40-63
15. 15. de Sousa Pinto J. A generalised Hankel convolution. SIAM Journal on Mathematical Analysis. 1985;16(6):1335-1346
16. 16. Malgonde SP, Gaikawad GS. On a generalized Hankel type convolution of generalized functions. Proceedings of the Indian Academy of Sciences–Mathematical Sciences. 2001;111(4):471-487
17. 17. G. von Winckel, “Bessel Function Zeros—File Exchange—MATLAB Central.” [Online]. Available from:http://www.mathworks.com/matlabcentral/fileexchange/6794-bessel-function-zeros. [Accessed: 06-Jun-2015]
18. 18. Chouinard U. Numerical simulations for the discrete Hankel transform [B.A.Sc. thesis]. Ottawa, Canada: University of Ottawa; 2015

Written By