Title Band-Limited Angular Spectrum Method for Numeric al Simulation of Free-Space Propagation in Far a nd Near Fields

A novel method is proposed for simulating free-space propagation. This method is an improvement of the angular spectrum method (AS). The AS does not include any approximation of the propagation distance, because the formula thereof is derived directly from the Rayleigh-Sommerfeld equation. However, the AS is not an all-round method, because it produces severe numerical errors due to a sampling problem of the transfer function even in Fresnel regions. The proposed method resolves this problem by limiting the bandwidth of the propagation field and also expands the region in which exact fields can be calculated by the AS. A discussion on the validity of limiting the bandwidth is also presented. ©2009 Optical Society of America OCIS codes: (050.1940) Diffraction; (090.1760) Computer holography; (070.0070) Fourier optics and signal processing References and links 1. J. W. Goodman, and R. W. Lawrence, “Digital image formation from electronically detected holograms,” Appl. Phys. Lett. 11(3), 77–79 (1967). 2. U. Schnars, and W. Jüptner, “Direct recording of holograms by a CCD target and numerical reconstruction,” Appl. Opt. 33(2), 179–181 (1994). 3. I. Yamaguchi, and T. Zhang, “Phase-shifting digital holography,” Opt. Lett. 22(16), 1268–1270 (1997). 4. T. Nakatsuji, and K. Matsushima, “Free-viewpoint images captured using phase-shifting synthetic aperture digital holography,” Appl. Opt. 47(19), D136–D143 (2008). 5. K. Matsushima, “Formulation of the rotational transformation of wave fields and their application to digital holography,” Appl. Opt. 47(19), D110–D116 (2008). 6. K. Matsushima, “Computer-generated holograms for three-dimensional surface objects with shade and texture,” Appl. Opt. 44(22), 4607–4614 (2005). 7. L. Ahrenberg, P. Benzie, M. Magnor, and J. Watson, “Computer generated holograms from three dimensional meshes using an analytic light transport model,” Appl. Opt. 47(10), 1567–1574 (2008). 8. H. Kim, J. Hahn, and B. Lee, “Mathematical modeling of triangle-mesh-modeled three-dimensional surface objects for digital holography,” Appl. Opt. 47(19), D117–D127 (2008). 9. T. M. Kreis, M. Adams, and W. P. O. Jüptner, “Methods of digital holography: A comparison,” SPIE Proc. 3098, 224–233 (1997). 10. N. Delen, and B. Hooker, “Free-space beam propagation between arbitrarily oriented planes based on full diffraction theory: a fast Fourier transform approach,” J. Opt. Soc. Am. A 15(4), 857–867 (1998). 11. K. Matsushima, H. Schimmel, and F. Wyrowski, “Fast calculation method for optical diffraction on tilted planes by use of the angular spectrum of plane waves,” J. Opt. Soc. Am. A 20(9), 1755–1762 (2003). 12. S. J. Jeong, and C. K. Hong, “Pixel-size-maintained image reconstruction of digital holograms on arbitrarily tilted planes by the angular spectrum method,” Appl. Opt. 47(16), 3064–3071 (2008). 13. S. De Nicola, A. Finizio, G. Pierattini, P. Ferraro, and D. Alfieri, “Angular spectrum method with correction of anamorphism for numerical reconstruction of digital holograms on tilted planes,” Opt. Express 13(24), 9935– 9940 (2005). 14. F. Zhang, I. Yamaguchi, and L. P. Yaroslavsky, “Algorithm for reconstruction of digital holograms with adjustable magnification,” Opt. Lett. 29(14), 1668–1670 (2004). 15. L. Yu, and M. K. Kim, “Pixel resolution control in numerical reconstruction of digital holography,” Opt. Lett. 31(7), 897–899 (2006). #115525 $15.00 USD Received 10 Aug 2009; revised 27 Sep 2009; accepted 4 Oct 2009; published 15 Oct 2009 (C) 2009 OSA 26 October 2009 / Vol. 17, No. 22 / OPTICS EXPRESS 19662 16. D. Wang, J. Zhao, F. Zhang, G. Pedrini, and W. Osten, “High-fidelity numerical realization of multiple-step Fresnel propagation for the reconstruction of digital holograms,” Appl. Opt. 47(19), D12–D20 (2008). 17. R. P. Muffoletto, J. M. Tyler, and J. E. Tohline, “Shifted Fresnel diffraction for computational holography,” Opt. Express 15(9), 5631–5640 (2007). 18. Y. M. Engelberg, and S. Ruschin, “Fast method for physical optics propagation of high-numerical-aperture beams,” J. Opt. Soc. Am. A 21(11), 2135–2145 (2004). 19. M. Sypek, “Light propagation in the Fresnel region. New numerical approach,” Opt. Commun. 116(1-3), 43–48 (1995). 20. E. A. Sziklas, and A. E. Siegman, “Mode calculations in unstable resonators with flowing saturable gain. 2: Fast Fourier transform method,” Appl. Opt. 14(8), 1874–1889 (1975). 21. F. Shen, and A. Wang, “Fast-Fourier-transform based numerical integration method for the RayleighSommerfeld diffraction formula,” Appl. Opt. 45(6), 1102–1110 (2006). 22. J. W. Goodman, Introduction to Fourier Optics, 2nd ed. (McGraw-Hill, 1996), chap. 2.2.


Introduction
The study of the propagation of wave fields in homogeneous and isotropic mediums has a long history.However, many researchers are still reporting new methods for numerical simulation of free-space propagation.In recent years, developments in computational holography such as digital holography or computer-generated holograms (CGH), have obviously driven this investigation.In digital holography, object fields are recorded using an image sensor and a numerical reconstruction is obtained as amplitude or phase images of the wave field that is numerically propagated to a plane around the object surface from the sensor [1][2][3].This propagation process can be regarded as digital signal processing of light.Therefore, various reconstructions that are impossible in ordinary digital imaging are possible in digital holography.For example, free viewpoint images [4] and a clear image of a deeply tilted surface [5] have been presented in digital holography using specific numerical propagation methods for the purpose.In the field of CGHs, numerical propagation also plays an important role in synthesizing object waves numerically [6][7][8].
There are three main methods for propagating wave fields between parallel planes [9].These methods, commonly based on the fast Fourier transform (FFT), are referred to in this paper as the single Fourier-transform-based Fresnel method (SFT-FR), the convolution-based Fresnel method (CV-FR), and the angular spectrum method (AS).There is also another category of numerical propagation, that is, propagation between non-parallel planes [5,[10][11][12][13], Methods for this type of propagation are usually given as an extension or modification of the parallel propagation methods.
The SFT-FR is simple and the fastest of the three methods as it uses only a single FFT.However, it has a serious disadvantage, in that the sampling window and intervals are proportionate to the propagation distance.To overcome this problem and add the ability to control the sampling intervals, the multi-step Fresnel method [14][15][16] and shifted Fresnel method (Shift-FR) [17] have been proposed as improvements to the original SFT-FR.
The CV-FR and AS are convolution-based methods, and as such, require at least two FFTs in the computation.The sampling window and intervals are not dependent on the propagation distance, i.e., they are always the same as the source field.The CV-FR is suitable for paraxial wave fields, while the AS is applicable to both non-paraxial fields and paraxial fields, since the formulas of the AS are derived directly from the Kirchhoff or Rayleigh-Sommerfeld diffraction theory without approximation.The CV-FR can be regarded as a kind of Fresnel approximation and a subset of the AS.Therefore, the behavior of the CV-FR is similar to that of the AS except in non-paraxial fields.Hence, we focus on the AS and barely consider the CV-FR in subsequent sections.
The current AS is, however, not an all-round method.It is suitable only for near field regions, whereas the SFT-FR and its family of algorithms are suitable for numerical propagation in far fields [18].The reason that the AS is not applicable for far field propagation is the sampling problem in its transfer function.Since the CV-FR is also not suitable for far-field calculation for the same reason, a multi-step method has been proposed for the CV-FR [19].This method, however, causes a different error especially in long distance propagation, because the cascaded sampling windows used are equivalent to diffraction by cascaded apertures.In addition, the computation time of the method is equal to the product of the number of the steps and the computation time of the original CV-FR.Total computation time is, therefore, much greater than that of the original method, especially in multiple step cases.This technique is applicable to the AS without any modification, but causes the same problem.Consequently, the current AS is not suitable for long distance propagation.The field size usually increases with increasing the propagation distance, whereas the size of the output sampling window of the AS does not change during propagation.Therefore, to achieve long distance propagation without sampling problem, the input sampling window must be extended so as to cover the whole filed in the output diffraction plane.This, however, requires a high computational effort.If the whole field in the diffraction plane is necessary for a specific purpose, there is no method other than the extension of the sampling window.However, in some cases, only a small region of the diffraction field within the aperture of an optical element is necessary for numerical simulation.In this case, the current AS has the big disadvantage.These sampling and window size problems of the AS already pointed out in early works using the AS [20].
We note that direct integration methods of Rayleigh-Sommerfeld diffraction formula are suitable for this kind of simulation [21].However, the methods require three FFTs for the computation, whereas the AS can be executed by two FFTs.
In this paper, we propose an improved AS that features suitability for long distance propagation as well as short distance propagation.This new method resolves the sampling problem in the AS and avoids the aliasing error of the transfer function by limiting the bandwidth and truncating unnecessary high-frequency signals in the input source field.Computation time of the proposed method is the same as the original AS.The suitability of limiting the bandwidth is also discussed in relation to the minimum bandwidth required for exact numerical propagation.

Formulation of the angular spectrum method
The coordinate system used in formulation is shown in Fig. 1.Source fields given in the source plane (x, y, 0) are propagated to the destination plane parallel to the source plane.The AS is equivalent to the Rayleigh-Sommerfeld formula [21].Here, the Rayleigh-Sommerfeld solution for an input monochromatic source field ( , , 0) g x y is given by: ( ) where and λ is the wavelength.This diffraction integral is rewritten by a convolution form using the propagation kernel (the impulse response) ( , , ) h x y z as follows: where the symbol * represents the two-dimensional convolution with respect to x and y.The propagation kernel is given by ( ) where ( ) The convolution of Eq. ( 2) is rewritten using the convolution theorem as: Here, the spectrum, ( , ; 0) G u v , and the transfer function, ( , ; ) H u v z , are given by: and where F represents the Fourier transform.The symbols u, v, and w are Fourier frequencies in the x, y, and z directions, respectively.These frequencies are not independent, i.e., frequency w is a function of u and v: ( ) Consequently, the AS is formulated as follows: Note that the impulse response and transfer function in the case of the Fresnel approximation are given, respectively, as: ( ) In the CV-FR, the transfer function FR ( , ; ) H u v z is used instead of ( , ; ) H u v z .

Discrete linear convolution in the convolution-based methods
All wave fields, spectrums, and the transfer function are sampled using an equidistant grid in a numerical simulation.Fourier transforms of Eqs. ( 5) and ( 8) are also replaced by the FFT.
However, a discrete Fourier transform of the sampled input field of ( , , 0) g x y involves periodicity of the fields in both the Fourier and real space, and therefore, the convolution with the transfer function using the FFT is a circular convolution.The terminology, circular convolution, is of the field of signal processing.Since circular convolution is for periodic functions, an error is caused by aperiodic functions of ( , , 0) g x y and ( , , ) h x y z in the edge of the computation window of ( , , ) g x y z .If the spatial extent of the output field ( , , ) g x y z is sufficiently small compared with the computation window and invasion of the fields in neighbor periods can be ignored, the influence of the circular convolution may be within a permissible level.However, in cases where the field diffuses strongly and runs over the computation window in the output plane, numerical propagation by convolution-based methods such as the AS or CV-FR no longer gives accurate results.To convert a circular convolution to a linear convolution, the area of the sampling window of the input field needs to be doubled along both the x-and y-axes as shown in Fig. 2, and the additional sampling points must be padded with zeros.Furthermore, the output window should be clipped and reduced to the original size once again.The manner of this clipping depends on the position of the input window in the given coordinates, as shown in Fig. 2. If the origin of the coordinates system is placed at the center of sampling array as shown in (a), the clipped output window is the same as that of the input sampling window, whereas if the origin is at the left lower corner, the output window is right upper quadrant of the sampling array and the origin is again left lower corner of the clipped output window.

Numerical errors of the angular spectrum method in long distance propagation
Computational results of long distance field propagation using the AS are not accurate, even if the convolution is linearized.To verify accuracy of the AS, one-dimensional diffraction by a rectangular aperture shown in Fig. 3   Amplitude distributions computed by the AS, Shifted-FR, and numerical integration of the diffraction integral, are shown in Fig. 4(a).Here, the numerical integration of Eq. ( 1) is calculated by using the trapezium rule.To ensure that errors of numerical integration are sufficiently smaller than that in other methods, the step size of numerical integration, which strongly affects the accuracy, is reduced to a value that the resultant wave field is no longer changed below.The Shift-FR gives almost the same results as the numerical integration, whereas the results obtained by the AS are very noisy.A comparison of accuracy between the AS and Shift-FR is shown as a function of the propagation distance in Fig. 4(b).Accuracy is defined as the SNR of the wave fields calculated by the methods as compared with the reference fields [11].Here the reference field is provided by numerical integration.Although the AS obtains good results in the near field regions, the accuracy declines for distances greater than about10 x S .We note that this numerical error of the AS can be avoided by extending the sampling window and computing the whole field.The exact result can be obtained by clipping the interest region of the whole field.However, this procedure requires a huge computational effort especially in long distance propagation of two dimensional wave fields.u, i.e., the signal frequency increases as u increases.Note that "signal frequency" as used here refers neither to physical frequencies of time nor space, but means the frequency of peaks and valleys of the function ( ; ) H u z in a certain period of u.Supposing that the one-dimensional transfer function of the AS is redefined as follows: [ ] ( ) the local signal frequency of the function ( ; ) H u z is given by [22]: When the transfer function is sampled at intervals of u ∆ , an aliasing error may be introduced in the sampled transfer function, as shown in Fig. 5.Note that the sampling interval is given by 1 (2 ) and not by , because the source sampling area is doubled to linearize the discrete convolution, as mentioned in Section 2.2.(2 ) (2 ) , where , and 50 x S = z .
The Nyquist theorem requires the following relation to avoid the aliasing error in the sampled transfer function.
This relation provides the means to determine the sampling intervals of the transfer function when band limited wave fields are numerically propagated.However, in practical applications such as digital holography, the sampling intervals are generally fixed and cannot be freely chosen.Therefore, the frequency range for which the transfer function causes no aliasing error must be determined by Eq. ( 12) as follows: .
As a result, the transfer function must be clipped within a bandwidth of limit 2u .
where rect( ) ξ is a rectangular function with width unity.
The wave field calculated by the AS using the band-limited transfer function of Eq. ( 14) is shown in Fig. 6(a).The errors seen in Fig. 4

Two-dimensional wave fields
To avoid aliasing errors, the region the two-dimensional transfer function of Eq. ( 6) must be limited, by applying the same procedure as in the one-dimensional case.The local signal frequency of the function ( , ; ) exp[ ( , )] is given by: where ( ) Relations to avoid the aliasing error are given by:  where 1 (2 ) is the sampling interval for the spatial frequency v and S y is the size of the sampling window in the y direction.Both relations in ( 17) must be satisfied in order to avoid the aliasing error in two-dimensional wave fields.
Consequently, the sampled transfer function is limited within the region expressed by: where limit u is defined in Eq. ( 13) and limit v is given by: Although both relations give ellipsoidal regions with major diameter 2λ −1 in the (u, v) plane as shown in Fig. 7, Eq. ( 18) depicts a vertical ellipse, whereas Eq. ( 19) produces a horizontal one.The transfer function and the spectrum of the wave field must be limited within the common region of these ellipsoidal regions.Minor radii of the ellipsoidal regions for Eqs. ( 18) and ( 19) are given by limit u and limit v , and therefore, the eccentricities are given by 2 u ∆ z and 2 v ∆ z , respectively.If either eccentricity of the ellipses is sufficiently greater than zero, the ellipse is oblate enough to regard the region as a rectangle within the sampling window of 1 x − ∆ and 1 y − ∆ , as shown in Fig. 8. Therefore, the ellipsoidal region can be approximated by a simple form in this case.When we adopt a value of 1/2 as the criterion of the eccentricity, the approximated regions and the criteria for applying the approximation are given by: limit if 2 , where S x and S y are again the sizes of the sampling window in the x and y directions, respectively.Note that the sampling intervals in the linearized convolution are once again sampling.This means that the error in bandwidths below need u is not attributed to numerical problems but physical property of field propagation.
These results verify the validity of the idea of a minimum bandwidth of field propagation and the suitability of limiting bandwidth proposed in this paper.

Conclusion
We proposed a new method for the exact calculation of field propagation in the free-space.This method is based on the angular spectrum method, but resolves the problem of numerical errors produced in far field propagation by the original AS.Since the errors can be attributed to the sampling problem of the sampled transfer function of the original AS, the bandwidth of the propagation field is limited to the range in which the sampling problem does not occur.As a result, the band-limited AS proposed in this paper is applicable to far field propagation as well as near field propagation.In addition, we verify the validity of limiting the bandwidth of the propagation field by considering the minimum bandwidth required for exact calculation of field propagation.

Fig. 1 .
Fig. 1.Definition of the coordinate system and the geometry of the model.

Fig. 2 .
Fig. 2. Conversion of a circular convolution into a linear convolution in cases where the origin of the coordinates system is placed at the center (a) and the left lower corner (b) of the sampling array.
is computed by three methods.Here, the sampling interval and the number of samplings are 2 in the source and destination, and therefore the size of the sampling window is 2048 $15.00 USD Received 10 Aug 2009; revised 27 Sep 2009; accepted 4 Oct 2009; published 15 Oct 2009 (C) 2009 OSA

Fig. 3 .
Fig. 3. Numerical simulation for verifying accuracy of the AS.

Fig. 4 .
Fig. 4. Comparison of the accuracy of the AS and Shift-FR.(a) One-dimensional amplitude distribution in the destination plane.(b) SNR of the AS and Shift-FR.

Fig. 5 .
Fig. 5. Example of a sampled transfer function of the AS.Only the real part of the transfer function ( ; ) H u z is depicted in the sampling interval limit (a) have disappeared and the SNR does not decrease in long distance propagation as shown in Fig. 6(b).

Fig. 6 .
Fig. 6.Accuracy of the band-limited AS proposed in this work.(a) One-dimensional amplitude distribution in the destination plane.(b) SNR of the AS and Shift-FR.Parameters used in the calculation are the same as in Fig. 4.

Fig. 7 .
Fig.7.Regions in the( , )  u v space to avoid aliasing errors of the sampled transfer function.The sampling intervals in the Fourier space are

Fig. 10 .
Fig. 10.Model for the minimum bandwidth required for exact calculation of field propagation.

Fig. 11 .
Fig. 11.SNR as a function of bandwidth in the AS in cases of the normal sampling (black solid line) and over sampling (red broken line).Here / 2 x x W S = and 50 x S = z .The other parameters are the same as in Fig. 4.