## Abstract

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

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

## 2. The angular spectrum method and its inherent problem

#### 2.1 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 $g(x,y,0)$ is given by:

*λ*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

*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*:

Note that the impulse response and transfer function in the case of the Fresnel approximation are given, respectively, as:

#### 2.2 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 $g(x,y,0)$ 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 $g(x,y,0)$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.

#### 2.3 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 is computed by three methods. Here, the sampling interval and the number of samplings are $\Delta x=2\lambda $ and ${N}_{x}=1024$, respectively, in the source and destination, and therefore the size of the sampling window is ${S}_{x}=2048\lambda $. The width of the rectangular aperture and the propagation distance are ${W}_{x}={S}_{x}/2$ and $z=50{S}_{x}$, respectively.

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 about$10{S}_{x}$.

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.

## 3. The aliasing error of the sampled transfer function

#### 3.1 One-dimensional wave fields

For simplicity, one-dimensional wave fields that are a function of *x* are discussed in this section. The transfer function of Eq. (6) for the AS is a kind of chirp function with respect to *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 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:

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(a) have disappeared and the SNR does not decrease in long distance propagation as shown in Fig. 6(b).

#### 3.2 Two-dimensional wave fields

To avoid aliasing errors, the region of 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 $H(u,v;z)=\mathrm{exp}[i\varphi (u,v)]$ is given by:

*v*and

*S*is the size of the sampling window in the

_{y}*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:

*λ*

^{−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 ${u}_{\text{limit}}$ and ${v}_{\text{limit}}$, and therefore, the eccentricities are given by $2\Delta uz$ and $2\Delta vz$, 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 $\Delta {x}^{-1}$ and $\Delta {y}^{-1}$, 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:

where*S*and

_{x}*S*are again the sizes of the sampling window in the

_{y}*x*and

*y*directions, respectively. Note that the sampling intervals in the linearized convolution are once again given by $\Delta u={(2{S}_{x})}^{-1}$ and $\Delta v={(2{S}_{y})}^{-1}$. As a result, if the condition ${S}_{x}\text{and}{S}_{y}\ll 2z$ is satisfied, the common region is a simple rectangle. In this case, the band-limited transfer function is written as:

Amplitude images of the output fields calculated using the original and the band-limited AS are shown in Fig. 9 . These fields are diffracted by square and circular apertures. Numerical errors such as high-frequency noise found in the original AS are not present in the band-limited AS.

## 4. Discussion on the minimum bandwidth required for exact numerical propagation

We can improve the angular spectrum method by limiting bandwidth of the sampled transfer function. However, because limiting bandwidth of the transfer function is equivalent to limiting bandwidth of the propagation field, the question arises of what bandwidth of the propagation field is required for accurate calculation.

A model for estimating the minimum bandwidth necessary for exact numerical propagation is shown in Fig. 10
. An aperture with size ${W}_{x}$ is placed at the center of a sampling window with size ${S}_{x}$. The highest spatial frequency, observed at the upper end of the destination sampling window, may be given by the field emitted from the point at the lower end of the source window. Therefore, the maximum frequency is given by ${u}_{\mathrm{max}}=\mathrm{sin}\theta /\lambda $ using the angle *θ* shown in Fig. 10 [16].

We introduced a band-limit of $2{u}_{\text{limit}}$into the source field and the transfer function in order to avoid numerical errors. However, if the cutting frequency is less than the maximum frequency, i.e. ${u}_{\text{limit}}<{u}_{\mathrm{max}}$, we might introduce another physical error into the diffraction field, because the source field loses a part of the frequency band necessary for exact diffraction.

Supposing that the destination plane is placed at a distance *z* from the aperture, the required minimum bandwidth is obtained from the geometry of the model as follows:

In Fig. 11 , the SNR is depicted as a function of the bandwidth of the propagation field in the AS. For the normal sampling of $\Delta u={\left(2{S}_{x}\right)}^{-1}$, the SNR has the maximum value for a bandwidth around $2{u}_{\text{limit}}$, and it decreases rapidly for bandwidths below ${u}_{\text{need}}$. In the case of the over sampling of $\Delta u={\left(4{S}_{x}\right)}^{-1}$, since the numerical error owing to the sampling problem is not caused, the SNR does not decrease for bandwidths above $2{u}_{\text{limit}}$, and there is no clear peak in the curve. However, the SNR decreases for bandwidths below ${u}_{\text{need}}$ as in the normal sampling. This means that the error in bandwidths below ${u}_{\text{need}}$ 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.

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

## Acknowledgments

This work was partially supported by the JSPS.KAKENHI (21500114).

## 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). [CrossRef]

**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). [CrossRef] [PubMed]

**3. **I. Yamaguchi and T. Zhang, “Phase-shifting digital holography,” Opt. Lett. **22**(16), 1268–1270 (1997). [CrossRef] [PubMed]

**4. **T. Nakatsuji and K. Matsushima, “Free-viewpoint images captured using phase-shifting synthetic aperture digital holography,” Appl. Opt. **47**(19), D136–D143 (2008). [CrossRef] [PubMed]

**5. **K. Matsushima, “Formulation of the rotational transformation of wave fields and their application to digital holography,” Appl. Opt. **47**(19), D110–D116 (2008). [CrossRef] [PubMed]

**6. **K. Matsushima, “Computer-generated holograms for three-dimensional surface objects with shade and texture,” Appl. Opt. **44**(22), 4607–4614 (2005). [CrossRef] [PubMed]

**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). [CrossRef] [PubMed]

**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). [CrossRef] [PubMed]

**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). [CrossRef]

**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). [CrossRef]

**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). [CrossRef] [PubMed]

**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). [CrossRef] [PubMed]

**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). [CrossRef] [PubMed]

**15. **L. Yu and M. K. Kim, “Pixel resolution control in numerical reconstruction of digital holography,” Opt. Lett. **31**(7), 897–899 (2006). [CrossRef] [PubMed]

**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). [CrossRef] [PubMed]

**17. **R. P. Muffoletto, J. M. Tyler, and J. E. Tohline, “Shifted Fresnel diffraction for computational holography,” Opt. Express **15**(9), 5631–5640 (2007). [CrossRef] [PubMed]

**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). [CrossRef]

**19. **M. Sypek, “Light propagation in the Fresnel region. New numerical approach,” Opt. Commun. **116**(1-3), 43–48 (1995). [CrossRef]

**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). [CrossRef] [PubMed]

**21. **F. Shen and A. Wang, “Fast-Fourier-transform based numerical integration method for the Rayleigh-Sommerfeld diffraction formula,” Appl. Opt. **45**(6), 1102–1110 (2006). [CrossRef] [PubMed]

**22. **J. W. Goodman, Introduction to Fourier Optics, 2nd ed. (McGraw-Hill, 1996), chap. 2.2.