A Wavelet Algorithm for Fourier-Bessel Transform Arising in Optics

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


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] ( 2 2 + 1 ) ( ) ←→ − 2 0 ( ) , 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][6][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

International Journal of Engineering Mathematics
In the case of the finite Hankel transform only a direct transform has an integral form. Without loss of generality its expression is (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 ( ) ] ( ) (or ] ( ) ] ( )) is a quick oscillating function if (or ) 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][15][16][17][18][19][20][21][22][23][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][28][29][30][31].
The organization of the paper is as follows: Section 2 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.

Wavelets and CAS Wavelets.
Wavelets constitute a family of functions constructed from dilation and translation of a single function ( ) 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]: where form a wavelet orthonormal basis for 2 ( ). CAS wavelets ( ) = ( , , , ) involve four arguments , , , and , where = 0, 1, . . . , 2 − 1, is assumed to be any nonnegative integer, is any integer, and is normalized time. CAS wavelets are defined as [33] ( ) where CAS ( ) = cos (2 ) + sin (2 ) .
An efficient algorithm has been presented for the Fourier-Bessel transform.

Outline of Algorithm
The function ( ) representing physical fields either are zero or have an infinitely long decaying tail outside a disk of finite radius . Hence, in most practical applications either the signal ( ) has a compact support or for a given > 0 there exists > 0 such that | ∫ ∞ ( ) ] ( ) | < . Therefore, in either case, known as the finite Hankel transform (FHT), is a good approximation of the HT as given by (2). Writing ( ) = ( ) in (8), we get̂] We may expand ( ) as follows: where = ⟨ ( ), ( )⟩. By truncating infinite series (10) at levels = 2 − 1 and = , we obtain an approximate representation for ( ) as where the matrices and are given by Substituting (11) in (9), we get Taking = 1 and = 1, (13) reduces tô Now, we relabel and write (14) aŝ where ] 's are the th-place integral in (14). The integrals arising in (14) are evaluated by using the formulae Re ] > −1 (16) (see [34]) and are calculated with the help of Simpson's onethird 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.

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 = 1, = 1, and = 60 in (15). We observed that the accuracy of the method is very high even at such a low level of truncation.   Simpson's One-Third Rule. See Figures 1 and 2.
Example 2. The following example was solved numerically by [35].
Simpson's One-Third Rule. See Figures 9 and 10.
Composite Simpson's One-Third Rule. See Figures 13 and 14.

Example 3 (sombrero function). A very important and often
used function is the Circ function that can be defined as [22] Circ This function is quite common in optical problems where it is used, for instance, to represent a circular pupil of radius .
International Journal of Engineering Mathematics       The Fourier-Bessel transform of (20) is the well-known "sombrero function." The zeroth-order Hankel transform of Circ( / ) is the sombrero function [29], given by The exact and numerical transforms differ very slightly but the differences are hardly visible.
Simpson's One-Third Rule. See Figures 17 and 18.
a well-known result. The pair ( ( ), 0 ( )) arises in optical diffraction theory [36]. The function ( ) is the optical transfer function of an aberration-free optical system with a circular aperture, and 0 ( ) is the corresponding spread function.   Barakat and Sandler [26] evaluated 0 ( ) numerically using Filon quadrature philosophy but the associated error is appreciable for < 1, whereas our method gives almost zero error in that range.

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.

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.