# A Wavelet Algorithm for Fourier-Bessel Transform Arising in Optics

Authors: Nagma Irfan; A. H. Siddiqi
Journal: International Journal of Engineering Mathematics
 (2015)

---

## Abstract

The aim of the paper is to propose an efficient and stable algorithm that is quite accurate and fast for numerical evaluation of the Fourier-Bessel transform of orderν,ν>-1, using wavelets. The philosophy behind the proposed algorithm is to replace the part tf(t) of the integral by its wavelet decomposition obtained by using CAS wavelets thus representing Fν(p) as a Fourier-Bessel series with coefficients depending strongly on the input function tf(t). The wavelet method indicates that the approach is easy to implement and thus computationally very attractive.


---

## Body

## 1. Introduction

The Fourier-Bessel transform (also designated as Hankel transform) is a very useful tool of mathematical physics [1]. It is a very useful instrument in a wide range of physical problems which have an axial symmetry. It is particularly important in optics and two-dimensional image processing, it naturally occurs in image reconstruction from projections or from reflected pulses, and it is a useful tool in the analysis and synthesis of three-dimensional wave fields. The present development is essentially motivated by optics application. The influence of the Laplacian on a function in cylindrical coordinates is equal to the product of the squared parameter of the transformation and the transform of the function [2](1)d2dr2+1rddrfr⟷-p2F0p,d2dr2+1rddr-1r2fr⟷-p2F1p.There are two types of the Hankel transform. The first one is defined on the semi-infinite interval. In this case the direct and inverse transforms of theνth kind are represented as a symmetric pair. When we are dealing with problems that show circular symmetry, Hankel transforms may be very useful [3, 4]. Laplace’s partial differential equation in cylindrical coordinates can be transformed into an ordinary differential equation by using the Hankel transform. Because the Hankel transform is the two-dimensional Fourier transform of a circularly symmetric function, it plays an important role in optical data processing [5–7]. In optics, the Hankel transform appears in many contexts, not the least of which is the propagation of cylindrically symmetric laser beams. Most classical optical systems like mirrors or lenses are axially symmetric devices. Hankel transform also proved to be extremely useful in problems associated with seismology, geophysics [8, 9], electroscattering, acoustics, hydrodynamics, image processing [10], time dependent Schrodinger equation, and so forth.Mathematical Background. The Fourier-Bessel transform may be defined by the following expression: (2)Fνp≡χνft≡∫0∞tftJνptdt,(3)ft=χν-1Fνp≡∫0∞pFνpJνptdp,where Jν is the νth-order Bessel function of the first kind.In the case of the finite Hankel transform only a direct transform has an integral form. Without loss of generality its expression is(4)Fνp=∫01tftJνptdt(see [11]). Practical calculation of direct and inverse Hankel transform is connected with two problems. The first problem is based on the fact that not every transform in the real physical situation has analytical expression for result of inverse Hankel transform. The second one is the determination of functions as a set of their values for numerical calculations. The classical trapezoidal rule, Cotes rule, and other rules connected with the replacement of the integrand by sequence of polynomials have high accuracy if integrand is a smooth function. But tf(t)Jν(pt) (or pFν(p)Jν(pt)) is a quick oscillating function if t (or p) is large. There are two general methods of the effective calculation in this area. The first is the fast Hankel transform [12]. The specification of that method is transforming the function to the logarithmical space and fast Fourier transform in that space. This method needs a smoothing of the function in log space. The second method is based on the separation of the integrand into product of slowly varying component and a rapidly oscillating Bessel function [13]. But it needs the smoothness of the slow component for its approximation by lower-order polynomials.To overcome these difficulties, various different techniques are available in the literature. Several papers have been written on the numerical evaluation of the HT in general and the zeroth order in particular [14–24]. There are two general methods of the effective calculation in this area. The first is the fast Hankel transform [25]. The specification of that method is transforming the function to the logarithmical space and fast Fourier transform in that space. This method needs a smoothing of the function in log space. The second method is based on the separation of the integrand into product of slowly varying component and a rapidly oscillating Bessel function [26]. But it needs the smoothness of the slow component for its approximation by lower-order polynomials. From variety of algorithm, a potential user would probably find it difficult to select any one algorithm that might be best for a particular application. For an overview of these algorithms and their numerical complexity, the reader is referred to [27–31].The organization of the paper is as follows: Section2 gives a brief description of the CAS wavelets followed by the derivation of the algorithm in Section 3. The efficiency and stability of the algorithm are shown by applying it to four test functions with known analytical transform in Section 4. At the end, a brief conclusion and future work are given in Section 5.


## 2. Properties of CAS Wavelets

### 2.1. Wavelets and CAS Wavelets

Wavelets constitute a family of functions constructed from dilation and translation of a single functionϕ(t) called the mother wavelets. When the dilation parameter is 2 and the translation parameter is 1 we have the following family of discrete wavelets [32]:(5)ψknt=2k/2ϕ2kt-n,where ψkn form a wavelet orthonormal basis for L2(R).CAS waveletsψnm(t)=ψ(k,n,m,t) involve four arguments n, k, m, and t, where n=0,1,…,2k-1, k is assumed to be any nonnegative integer, m is any integer, and t is normalized time. CAS wavelets are defined as [33](6)
ψ
n
m
t
=
2
k
/
2
CAS
m
2
k
t
-
n
,
for
n
2
k
≤
t
<
n
+
1
2
k
,
0
,
otherwise
,where(7)CASmt=cos⁡2mπt+sin⁡2mπt.An efficient algorithm has been presented for the Fourier-Bessel transform.


## 2.1. Wavelets and CAS Wavelets

Wavelets constitute a family of functions constructed from dilation and translation of a single functionϕ(t) called the mother wavelets. When the dilation parameter is 2 and the translation parameter is 1 we have the following family of discrete wavelets [32]:(5)ψknt=2k/2ϕ2kt-n,where ψkn form a wavelet orthonormal basis for L2(R).CAS waveletsψnm(t)=ψ(k,n,m,t) involve four arguments n, k, m, and t, where n=0,1,…,2k-1, k is assumed to be any nonnegative integer, m is any integer, and t is normalized time. CAS wavelets are defined as [33](6)
ψ
n
m
t
=
2
k
/
2
CAS
m
2
k
t
-
n
,
for
n
2
k
≤
t
<
n
+
1
2
k
,
0
,
otherwise
,where(7)CASmt=cos⁡2mπt+sin⁡2mπt.An efficient algorithm has been presented for the Fourier-Bessel transform.


## 3. Outline of Algorithm

The functionft representing physical fields either are zero or have an infinitely long decaying tail outside a disk of finite radius R. Hence, in most practical applications either the signal ft has a compact support or for a given ε>0 there exists R>0 such that |∫R∞tf(t)Jν(pt)dt|<ε.Therefore, in either case,(8)F^νp=∫0RtftJνptdt=∫01tftJνptdt,by scalingknown as the finite Hankel transform (FHT), is a good approximation of the HT as given by (2). Writing tf(t)=g(t) in (8), we get(9)F^νp=∫01gtJνptdt.We may expand gt as follows:(10)gt=∑n=0∞∑m∈Ζcnmψnmt,where cnm=g(t),ψnmt.By truncating infinite series (10) at levels n=2k-1 and m=M, we obtain an approximate representation for gt as(11)gt≈∑n=02k-1∑m=-MMcnmψnmt=CTψt,where the matrices C and ψ are given by (12)C=c0-M,c0-M+1,…,c0,M,c1-M,…,c1M,…,c2k-1-M,…,c2k-1MT,ψt=ψ0-M,ψ0-M+1,…,ψ0,M,ψ1-M,…,ψ1M,…,ψ2k-1-M,…,ψ2k-1MT.Substituting (11) in (9), we get(13)F^νp≈CT∫01ψtJνptdt.Taking M=1 and k=1, (13) reduces to (14)F^νp≈CT∫01ψ0-1tJνptdt,∫01ψ0,0tJνptdt,∫01ψ0,1tJνptdt,…,∫01ψ1-1tJνptdt,…,∫01ψ1,1tJνptdtT.Now, we relabel and write (14) as(15)F^νp≈c10,c11,…,c42Iν0,In1,…,In11T,where Iνl’s are the lth-place integral in (14).The integrals arising in (14) are evaluated by using the formulae(16)∫0aJνtdt=2limL→∞⁡∑z=0LJν+2z+1a,Re⁡ν>-1(see [34]) and are calculated with the help of Simpson’s one-third rule, Simpson’s three-eighth rule, composite Simpson’s one-third rule, and composite Simpson’s three-eighth rule, respectively. In numerical analysis, Simpson’s rule and composite Simpson’s rule are method for numerical integration, the numerical approximation of definite integrals.


## 4. Numerical Results

In this section, we test the proposed algorithm (15) by evaluating the approximate Hankel transforms of 4 well-known test functions with known analytical Hankel transforms. Note that in all the examples the truncation is done at level N=1, k=1, and L=60 in (15). We observed that the accuracy of the method is very high even at such a low level of truncation.Example 1.
Letf(r)=rνsin⁡(πr2/4), 0≤r<1; then(17)Fνp=12π2-ν-1pνUν+1π2,p-Uν+2π2,p(obtained  from  [34]  by  putting  a=π/4, b=1), where Uν(w,p) is a Lommel function of two variables,(18)=12p∑η=0L-1ηπ2p2ηJν+2η+1p-π2pJν+2η+2pasL⟶∞.The comparison of the approximation Hν(p) (dotted line) with the exact Hankel transform Fν(p) (solid line) is shown in Figures 1, 3, and 5 and the error E(p)=Hν(p)-Fν(p) in Figures 2, 4, and 6.
Simpson’s One-Third Rule. See Figures 1 and 2.
Simpson’s Three-Eighth Rule. See Figures 3 and 4.
Composite Simpson’s One-Third Rule. See Figures 5 and 6.
Composite Simpson’s Three-Eighth Rule. See Figures 7 and 8.Figure 1
The exact transform,Fν(p) (solid line), and the approximate transform, Hν(p) (dotted line), where ν=0.Figure 2
Comparison of the errors.Figure 3
The exact transform,Fν(p) (solid line), and the approximate transform, Hν(p) (dotted line), where ν=0.Figure 4
Comparison of the errors.Figure 5
The exact transform,Fν(p) (solid line), and the approximate transform, Hν(p) (dotted line), where ν=0.Figure 6
Comparison of the errors.Figure 7
The exact transform,Fν(p) (solid line), and the approximate transform, F(p) (dotted line), where ν=0.Figure 8
Comparison of the errors.Example 2.
The following example was solved numerically by [35].
For(19)fr=e-r,its HT isF0p=11+p23/2.We solve the above problem by the proposed algorithm and observe that our method gives a result comparable to [35]. Note that F0p and F^0p are indicated by F0p (solid line) and Fp (dotted line) in Figures 9, 11, 13, and 15 and the error Ep=F(p)-F0(p) is shown in Figures 10, 12, 14, and 16.
Simpson’s One-Third Rule. See Figures 9 and 10.
Simpson’s Three-Eighth Rule. See Figures 11 and 12.
Composite Simpson’s One-Third Rule. See Figures 13 and 14.
Composite Simpson’s Three-Eighth Rule. See Figures 15 and 16.Figure 9
The exact transform,F0(p) (solid line), and the approximate transform, F(p) (dotted line).Figure 10
Comparison of the errors.Figure 11
The exact transform,F0(p) (solid line), and the approximate transform, F(p) (dotted line).Figure 12
Comparison of the errors.Figure 13
The exact transform,Fν(p) (solid line), and the approximate transform, Hν(p) (dotted line).Figure 14
Comparison of the errors.Figure 15
The exact transform,F0(p) (solid line), and the approximate transform, F(p) (dotted line).Figure 16
Comparison of the errors.Example 3 (sombrero function).
A very important and often used function is the Circ function that can be defined as [22](20)Circ⁡ta=1,t≤a,0,t>a.This function is quite common in optical problems where it is used, for instance, to represent a circular pupil of radius a. The Fourier-Bessel transform of (20) is the well-known “sombrero function.”
The zeroth-order Hankel transform ofCirc⁡(t/a) is the sombrero function [29], given by(21)S0p=a2J1apap.The exact and numerical transforms differ very slightly but the differences are hardly visible.
Simpson’s One-Third Rule. See Figures 17 and 18.
Simpson’s Three-Eighth Rule. See Figures 19 and 20.
Composite Simpson’s One-Third Rule. See Figures 21 and 22.
Composite Simpson’s Three-Eighth Rule. See Figures 23 and 24.Figure 17
The exact transform,F0(p) (solid line), and the approximate transform, F(p) (dotted line).Figure 18
Comparison of the errors.Figure 19
The exact transform,F0(p) (solid line), and the approximate transform, F(p) (dotted line).Figure 20
Comparison of the errors.Figure 21
The exact transform,F0(p) (solid line), and the approximate transform, F(p) (dotted line).Figure 22
Comparison of the errors.Figure 23
The exact transform,F0(p) (solid line), and the approximate transform, F(p) (dotted line).Figure 24
Comparison of the errorsExample 4.
Letft=2/π[arccos⁡t-t(1-t2)1/2], 0≤t≤1; then,(22)F0p=2J12p/2p2,0≤p≤∞,a well-known result. The pair (f(t),F0(p)) arises in optical diffraction theory [36]. The function f(t) is the optical transfer function of an aberration-free optical system with a circular aperture, and F0(p) is the corresponding spread function.
Barakat and Sandler [26] evaluated F0p numerically using Filon quadrature philosophy but the associated error is appreciable for p<1, whereas our method gives almost zero error in that range.
Simpson’s One-Third Rule. See Figures 25 and 26.
Simpson’s Three-Eighth Rule. See Figures 27 and 28.
Composite Simpson’s One-Third Rule. See Figures 29 and 30.
Composite Simpson’s Three-Eighth Rule. See Figures 31 and 32.Figure 25
The exact transform,F0(p) (solid line), and the approximate transform, F(p) (dotted line).Figure 26
Comparison of the errors.Figure 27
The exact transform,F0(p) (solid line), and the approximate transform, F(p) (dotted line).Figure 28
Comparison of the errors.Figure 29
The exact transform,F0(p) (solid line), and the approximate transform, F(p) (dotted line).Figure 30
Comparison of the errors.Figure 31
The exact transform,F0(p) (solid line), and the approximate transform, F(p) (dotted line).Figure 32
Comparison of the errors.


## 5. Conclusion

Since the basis functions used to construct the wavelets are orthogonal and have compact support, it makes them more useful and simple in actual computations. Also, since the numbers of mother wavelet’s components are restricted to one, they do not lead to the growth of complexity of calculations. Our choice of wavelets makes them more attractive in their applications in the applied physical problems as they eliminate the problems connected with the Gibbs phenomenon taking place in [30]. A good agreement between the obtained solution and some well-known results has been obtained. Four test examples are provided to show the advantage of using wavelets. This method is capable of greatly reducing the size of calculations while still maintaining high accuracy of the numerical solution.Proposed wavelet method is very simple and attractive. The implementation of the current approach in analogy to existing methods is more convenient and the accuracy is high. The numerical example and the compared results support our claim. The difference between the exact and approximate solutions for each example was plotted graphically to determine the accuracy of numerical solutions.


### 5.1. Future Work

Since computational work is fully supportive of compatibility of the proposed algorithm, the same may be extended to other physical problems also. A very high level of accuracy explicitly reflects the reliability of this scheme for such problems. We would like to stress that the approximate solution includes not only time information but also frequency information due to the localization property of wavelet basis; with some change we can apply this method with the help of other wavelet bases.


## 5.1. Future Work

Since computational work is fully supportive of compatibility of the proposed algorithm, the same may be extended to other physical problems also. A very high level of accuracy explicitly reflects the reliability of this scheme for such problems. We would like to stress that the approximate solution includes not only time information but also frequency information due to the localization property of wavelet basis; with some change we can apply this method with the help of other wavelet bases.



---

*Source: 789675-2015-08-31.xml*