Open Access
Add to BibSonomyAbstract
Unlike the FFT, the Quasi Discrete Hankel Transform (QDHT) is not sampled on a uniform grid; in particular the field may no longer be sampled on axis. We demonstrate how the generalised sampling theorem may be applied to optical problems, evaluated with the QDHT, to efficiently and accurately reconstruct the optical field at any point. Without sacrificing numerical accuracy this is demonstrated to be typically 50× faster than using an equivalent 2D FFT.
©2010 Optical Society of America
1. Motivation
Two dimensional beam propagation, utilising spectral techniques such as the split-step method [1], has been well studied. For radially symmetric beams the treatment of their propagation may be simplified by replacing the 2D Fourier transform with a Hankel transform. Although the FFT is well known for being extremely efficient; reducing the number of sample points, by assuming radial symmetry, means that we may use a less efficient, but 1D, Hankel transform technique and ultimately achieve a saving. Either transform may be further optimised by sampling close to the Nyquist limit therefore limiting the number of sample points to the minimum required to adequately describe the beam.
The Hankel transform is commonly evaluated with the Quasi Fast Hankel Transform (QFHT) [2], whereby the Hankel transform is expressed as a convolution and evaluated using the FFT technique. Although the QFHT is an very fast technique it has been shown that the Quasi Discrete Hankel Transform (QDHT) [3] is more efficient for a given accuracy [4].
Fig. 1. Sample points for the FFT (o) and zeroth order QDHT (x). Separation of QDHT points has been exaggerated for illustrative purposes.
Figure 1 shows the sample points of the QDHT relative to a slice, through the origin, of the uniform 2D FFT. The QDHT sampling begins non-uniformly, asymptotically approaching a uniform grid separation. Although there are no sample points for r < 0 the grey points marked on Fig. 1 represent the negative positions due to the radial symmetry. Unlike the 2D FFT, the zero order QDHT does not explicitly represent the field at the axis. Although this does not present a problem for the numerical propagation of the beam, it is often the axial intensity we are most interested in [5,6]. Not knowing this point could cause us to draw incorrect conclusions about the nature of some localised process of interest.
As a concrete example, Fig. 2 shows the intensity for two linear Bessel-Gauss beams propagating in the direction z. The initial field is given by
where α 1 and α 2 are the principle transverse frequencies of the separate Bessel-Gauss beams. These beams have been sampled just above the numerical space-bandwidth product and are therefore entirely defined within numerical limits. The analytical axial intensity [7] is compared to the first sample point; clearly the first sample provides an extremely poor representation for the axial value!
By increasing the sampling density it is clearly possible to move the first point closer to the axis. This may result in an undesirable increase in the number of sample points required to propagate the beam, clumsily sabotaging the efficiency of the technique. Ultimately this may cause the radial scheme to be less efficient than a 2D FFT due to the oversampling in real space.
To be able to properly observe the axial behaviour, while maintaining the efficient nature of the Nyquist limited QDHT, we desire a technique that is able to determine the optical field at points which do not coincide with sample points, in particular the axial value. We may achieve this by applying the generalised sampling theorem [8].
2. Theory
The n th order Hankel transform, for a function F(kr), may be expressed as the finite integral [4],
provided F(kx) = 0 for kx > K max. For optical beams, the transform order, n, depends on the azimuthal phase variation term, exp(-inθ). For radially symmetric beams with no azimuthal phase, the n = 0 Hankel transform is used.
The size of the numerical spatial domain, R, and spatial frequency domain, K, for the QDHT are related by the product;
Here N is the number of sample points and j n,N+1 is the (N + 1)th root of the Bessel function Jn. If a function posses no value above r = R max and no frequency component above kr = K max then it is entirely defined, within the computational limits of the QDHT, by N samples when
no additional precision is gained by increasing N above this value [3, 4]. This radial space-bandwidth product limit is analogous to the familiar Shannon-Nyquist sampling theorem. Of course, no function may strictly be simultaneously band-limited in both space and frequency space. In practice, however, we may say that the function is space-bandwidth limited by numerical noise when the function possesses no value above numerical precision outside these bounds.
It can be shown [9,10] that functions sampled on the grid rn,p, where rn,p = jn,p/K, may be reconstructed at an arbitrary point r as;
with the sampling kernel
From this we may reconstruct the function f(rn,p) as f 0(r) at any arbitrary point in r, including r = 0 provided the function is sampled on the non-uniform grid rn,p. This reconstruction method is generalised for any Hankel transform of order n but is only considered here in relation to the QDHT for which n ∈ ℤ.
Additionally, we may be interested in points which do not lie at r = 0 but do not coincide with a point on rn,p either. For a single point it may be reasonable to simply adjust the number of sample points such that the desired point now lies on rn,p. In general however this may not be appropriate. Consider for example the comparison of a radially propagated beam with data sampled on grid defined for use with the FFT. In this situation the uneven separation of rn,p is not compatible with the even spacing of the FFT.
The application of Eq. (5) is not limited to the study of beams for which the axis is the main point of interest. Because this equation is generalised to Hankel transforms of order n it is possible to study the propagation of beams possessing angular momentum. For such vortex beams the field on axis will always be zero, however it is often the near-axial points which are of interest. Using Eq. (5) it would be possible to model these beams efficiently using a minimal number of sample points and then reconstruct the near-axial region – or any region of interest – to a high degree of accuracy. It is also possible to compare the propagation of beams with differing azimuthal phase where the sample points would lie on the grids r n1,p and r n2,p where n 1 ≠ n 2. This generalised interpolation should allow mixing of various techniques such as resampling linear combinations of vortex beams on to grids suitable for finite differences or 2D FFT grids without the addition of significant numerical artifacts.
3. Analysis
To analyse the reconstruction error we test five beam-like functions with different rates of spectral convergence [11]. The reconstruction errors of these functions, for increasing sample densities, are shown in Fig. 3. The forms, and rates of convergence in frequency space, of the test functions are given in Table 1. These are compared against the same functions reconstructed from a 2D Cartesian grid of size N 2 FFT using the cardinal series [9] in Fig. 3(b). Such reconstruction is compatible with the grid spacing required for the FFT.
The error was determined by comparing the function f 0(r), reconstructed from the sampled function f(rn,p), against the analytic form of the test function f(r). From this we determine the maximum error,
This error was estimated by sampling r over a large number of points linearly spaced over the region 0 ≤ r ≤ R max, for Fig. 3(a) 1000 points were used to determine ε ∞ and for Fig. 3(b) the function was reconstructed onto 2D grid of 1000 × 1000 points.
The reconstruction error is introduced when the function is not effectively space-bandwidth limited. We therefore expect the error due to reconstruction to fall, with increasing sample density, at the same rate as the function. The larger noise floor in Fig. 3(a) when compared to the 2D FFT is caused by the limited precision to which the Bessel function and its roots may be calculated. To reach a nominal accuracy of 10-10 in the reconstruction of a Gaussian function requires N = 95 points using Eq. (5) and N FFT = 176 along one dimension (30976 points in total) using the cardinal series. For the sech function the same precision is reached when N = 314 and N FFT = 566.
The rate of convergence for a function in the real or spectral domain may be generalised into two forms, geometric and algebraic. Functions exhibiting algebraic convergence, converge as f(r) ~ r -q for any positive q; geometric convergence is typified by the rate f(r) ~ exp(-rp). This can be subdivided into Super- (p > 1), Simple- (p = 1) and Sub-Geometric (p < 1).
Table 1. Typical beam-like functions.
The Gaussian, Bessel-Gauss and Laguerre-Gauss functions all exhibit a super-geometric rate of convergence in both space and frequency space and are consequently excellent candidates for representation by spectral techniques of this form. This results in a rapid convergence in the reconstruction error on a noise floor. The off-axis maximum in the Bessel-Gauss reconstruction error may be simply attributed to the off-axis nature of the Bessel-Gauss spectrum. Although the sech function converges at a geometric rate in both domains, the reconstruction error still quickly converges on a noise floor.
Finally we consider a step function which is effectively a pupil, due to the radial symmetry. The spectrum of a step function converges at an algebraic rate. This poor rate of convergence means that the function is never well sampled. As a result we see that the error also decreases algebraically ~ O(N -2). To reduce the error to 10-10, significantly greater precision than is often of experimental interest, would require 105 sample points for the QDHT. Therefore the reconstruction technique, and indeed spectral techniques in general, are less suited to functions of this form. The 2D reconstruction of this function using the cardinal series results in a constantly large maximum error [Fig. 3(b)] as the pupil function is poorly represented by a Cartesian grid.
To assess the accuracy and speed of this reconstruction technique, when combined with a beam propagation model, we compare the axial value obtained in four ways. First we use a 2D FFT technique for which the axial value will already lie on the sample grid; second we use the QDHT in combination with the reconstruction technique described above; third we use the QDHT and assume r = 0 ≃ r 0,1; finally we use the oversampled QDHT so that r 0,1 is much closer to r = 0.
To assess the accuracy and efficiency of this technique when combined with a simple beam propagation model; a Gaussian beam profile is propagated over one Rayleigh length from its waist using 4000 steps for each technique. The reconstructed value on axis, f 0(r = 0,z) is compared to the analytic value f(r = 0,z) and the maximum error for the axial point is determined from
Due to the radial symmetry of the Gaussian beam, the n = 0 QDHT transform is used. To make the most efficient use of each transform the total number of sample points N, given in Table 2, was chosen so that the beam is sampled just above the numerical limit of the space-bandwidth product through the entire propagation. For the 2D FFT N was padded from the observed Nyquist minimum of N = 134 points to 135, a product of small prime numbers (3,3,3,5) to make optimal use of the FFTW routine [12], the 2D FFT will therefore have N 2 points. This is approximately twice as fast as padding to the nearest power of two (256). As such the 2D FFT method is presented in the best possible light.
Table 2. Accuracy of field on axis and relative speed for each method.
The speed for each technique, shown in Table 2, has been normalised relative to the 2D FFT which is 100 time units. The QDHT, without reconstruction, is the fastest but there is a considerable error in the value of the r = 0 point, consistent with Fig. 2. The QDHT technique, with reconstruction at r = 0, exhibits negligible additional time penalty but vastly improves the accuracy of the axial value. Increasing the number of sample points in the QDHT causes the time to increase as N 2 and only reduces the error by N -2, consequently the oversampled field only sees a small improvment in accuracy. Clearly to approach the accuracy of the QDHT technique with reconstruction, or indeed the naive 2D FFT, oversampling the QDHT is an unwise proposition.
4. Conclusions
The QDHT represents an accurate and efficient technique to supplement the numerical modelling of beam propagation. One shortcoming of the technique is the location of the sampling points; in particular the lack of a sample point at r = 0 for the zero order QDHT. We may accurately reconstruct this, or any other, point by applying a specific form of the generalised sampling theorem. Providing the beam is well sampled in the space-bandwidth sense the error in this technique is only limited by the precision to which the Bessel roots may be calculated. Reconstructing the beam at r = 0 is considerably faster than using a 2D FFT model and is both faster and more accurate than oversampling the beam.
Acknowledgements
A. Norfolk is generously supported by an EPSRC doctoral training award. E. Grace is generously supported by a fellowship from the RAE/EPSRC.
References and links
1. M. D. Feit and J. A. Fleck Jr., “Light propagation in graded-index optical fibers,” Appl. Opt. 17, 3990–3998 (1978). [CrossRef] [PubMed]
2. A. E. Siegman, “Quasi fast Hankel transform,” Opt. Lett. 1, 13–15 (1977). [CrossRef] [PubMed]
3. L. Yu, M. Huang, M. Chen, W. Chen, W. Huang, and Z. Zhu, “Quasi-discrete Hankel transform”, Opt. Lett. 23(6), 409–411 (1998). [CrossRef]
4. M. Guizar-Sicairos and J. C. Gutièrrez-Vega, “Computation of quasi-discrete Hankel transforms of integer order for propagating optical wave fields,” J. Opt. Soc. Am. A 21, 53–8 (2004). [CrossRef]
5. P. Srisungsitthisunti, O. K. Ersoy, and X. Xu, “Beam propagation modeling of modified volume Fresnel zone plates fabricated by femtosecond laser direct writing,” J. Opt. Soc. Am. A 26(1), 188–194 (2009). [CrossRef]
6. A. W. Norfolk and E. J. Grace, “New beat length for writing periodic structures using Bessel beams,” Opt. Commun. 283(3), 447–450 (2010). [CrossRef]
7. D. Ding and X. Liu, “Approximate description for Bessel, Bessel-Gauss, and Gaussian beams with finite aperture.” J. Opt. Soc. Am. A 16, 1286–1293 (1999). [CrossRef]
8. H. P. Kramer, “A Generalized Sampling Theorem,” J. Math. & Phys. 38, 68–72 (1959).
9. A. J. Jerri and D. W. Kreisler, “Sampling Expansions with Derivatives for Finite Hankel and Other Transforms,” SIAM J. Math. Anal. 6(2), 262–267 (1975). [CrossRef]
10. L. L. Campbell, “A comparison of the sampling theorems of Kramer and Whittaker,” J. Soc. Indust. Appl. Math. 12, 117–130 (1964). [CrossRef]
11. J. P. Boyd, Chebyshev and Fourier Spectral Methods (Dover Publications, 2001).
12. M. Frigo and S. G. Johnson, “The design and implementation of FFTW3,” Proc. IEEE 93(2), 216–231 (2005). [CrossRef]
