Band-Limited Angular Spectrum Method for

Numerical Simulation of Free-Space

Propagation in Far and Near Fields
Kyoji Matsushima1,* and Tomoyoshi Shimobaba2
1
Department of Electrical and Electronic Engineering, Kansai University, 3-3-35 Yamate-cho, Suita,
Osaka 564-8680, Japan
2
Graduate School of Engineering, Chiba University, 1-33 Yayoi-cho, Inage-ku,
Chiba 263-8522, Japan
matsu@kansai-u.ac.jp

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
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 Rayleigh-
Sommerfeld 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.

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

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

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

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:

exp ( i 2π r ' λ −1 ) z  1 1
g ( x, y, z ) = ∫∫ g ( x ', y ', 0)  +  dx ' dy ', (1)
r' r '  2π r ' iλ 

#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 19664
 1/ 2
where r ' = ( x − x ') + ( y − y ') + z 2  and λ is the wavelength. This diffraction integral is
2 2
 
rewritten by a convolution form using the propagation kernel (the impulse response)
h( x, y, z ) as follows:

g ( x, y, z ) = g ( x, y, 0) ∗ h( x, y, z ), (2)
where the symbol ∗ represents the two-dimensional convolution with respect to x and y. The
propagation kernel is given by

exp ( i 2π r λ −1 ) z  1 1
h ( x, y , z ) =  + , (3)
r r  2π r iλ 

where r = ( x 2 + y 2 + z 2 )
1/ 2
. The convolution of Eq. (2) is rewritten using the convolution
theorem as:
G (u, v; z ) = G (u, v; 0) H (u, v; z ). (4)

Here, the spectrum, G (u, v; 0) , and the transfer function, H (u, v; z ) , are given by:

G (u, v;0) = ∫∫ g ( x, y, 0) exp [ −i 2π (ux + vy )] dxdy

(5)
= F { g ( x, y, 0)} ,

and
H (u, v; z ) = F {h( x, y, z )}
(6)
= exp [i 2π wz ] ,

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:
( λ −2 − u 2 − v 2 )1/ 2 ⋯ u 2 + v 2 ≤ λ −2
w = w(u, v) =  (7)
0 ⋯ otherwise
Consequently, the AS is formulated as follows:

g ( x, y, z ) = F −1 {G (u, v; 0) exp [ i 2π w(u, v) z ]} . (8)

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

1  x 2 + y 2 
hFR ( x, y; z ) = exp i 2πλ −1  z +  ,
iλ z   2z  
(9)
 2 
H FR (u, v; z ) = exp iπ z  − λ ( u 2 + v 2 )   .
  λ 

In the CV-FR, the transfer function H FR (u , v; z ) is used instead of H (u, v; z ) .

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.

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

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.

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 ∆x = 2λ and N x = 1024 , respectively, in the source
and destination, and therefore the size of the sampling window is S x = 2048λ . The width of
the rectangular aperture and the propagation distance are Wx = S x / 2 and z = 50 S x ,
respectively.

#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 19666
 Fig. 3. Numerical simulation for verifying accuracy of the AS.

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.

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.

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
about10S 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

#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 19667
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:
H (u; z ) = exp [iφ (u ) ] ,
(10)
φ (u ) = 2π ( λ −2 − u 2 ) z,
1/ 2

the local signal frequency of the function H (u; z ) is given by [22]:

1 ∂φ
fu =
2π ∂u
uz (11)
= .
2 1/ 2
 λ −2 − u 

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 ∆u = (2S x ) −1 and not by ∆u = S x−1 , because the source sampling area is doubled to
linearize the discrete convolution, as mentioned in Section 2.2.

Fig. 5. Example of a sampled transfer function of the AS. Only the real part of the transfer
−1 −1
function H (u; z ) is depicted in the sampling interval ∆u = (2 S x ) = (2 N x ∆x ) , where
N x = 1024 , ∆x = 2λ , and z = 50 S x .

The Nyquist theorem requires the following relation to avoid the aliasing error in the
sampled transfer function.
∆u −1 ≥ 2 fu . (12)

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:
1
u ≤ 1/ 2
 (2∆u z ) 2 + 1 λ (13)
≡ ulimit .

As a result, the transfer function must be clipped within a bandwidth of 2ulimit .

#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 19668
  u 
H ′(u; z ) = H (u; z )rect  , (14)
 2ulimit 
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(a) have disappeared and the SNR does not
decrease in long distance propagation as shown in Fig. 6(b).

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.

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 ) = exp[iφ (u , v)] is given by:

1 ∂φ uz
fu = = , (15)
2π ∂u λ − u 2 − v 2 1/ 2
−2
 

1 ∂φ vz
fv = = , (16)
2π ∂v  λ −2 − u 2 − v 2 1/ 2
 

where φ (u, v) = 2π ( λ −2 − u 2 − v 2 )
1/ 2
z . Relations to avoid the aliasing error are given by:

∆u −1 ≥ 2 f u and ∆v −1 ≥ 2 f v , (17)

#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 19669
 Fig. 7. Regions in the (u , v ) space to avoid aliasing errors of the sampled transfer function.
−1 −1
The sampling intervals in the Fourier space are ∆u = ∆v = (2 S x ) = (2 S y ) for the sampling
windows S x = N x ∆x and S y = N y ∆y in the real space, where N x = N y = 1024 and
∆x = ∆y = 2λ .

where ∆v = (2S y ) −1 is the sampling interval for the spatial frequency v and Sy 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:

u2 v2
+ −2 ≤ 1, (18)
ulimit λ
2

u2 v2
−2
+ ≤ 1, (19)
λ 2
vlimit

where ulimit is defined in Eq. (13) and vlimit is given by:

−1/ 2
vlimit = (2∆v z ) 2 + 1 λ −1 . (20)

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 ulimit and vlimit ,
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 ∆x −1 and ∆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:
u ≤ ulimit if S x ≪ 2z , (21)

v ≤ vlimit if S y ≪ 2 z , (22)

where Sx and Sy 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

#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 19670
given by ∆u = (2S x ) −1 and ∆v = (2S y ) −1 . As a result, if the condition S x and S y ≪ 2z is
satisfied, the common region is a simple rectangle. In this case, the band-limited transfer
function is written as:

 u   v 
H ′(u, v; z ) = H (u , v; z )rect   rect  . (23)
 2 limit 
u  2vlimit 
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.

Fig. 8. A schematic illustration of the approximated rectangular region.

Fig. 9. Amplitude images calculated by the original and band-limited AS. (a) Diffraction by a
square aperture with dimensions S 2 × S 2 and z = 100 Sx. (b) Diffraction by a circular
x y

aperture with diameter Sx/2 and z = 200 Sx.

#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 19671
 Fig. 10. Model for the minimum bandwidth required for exact calculation of field propagation.

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 Wx 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
umax = sin θ / λ using the angle θ shown in Fig. 10 [16].
We introduced a band-limit of 2ulimit 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. ulimit < umax , 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:
−1/ 2
 2 z  2 
uneed ≡ 2umax = 2   + 1 λ −1 . (24)
 Wx + S x  

When the size of the aperture Wx is the same as the size of the sampling window S x , uneed
has the maximum value, which is the same as 2ulimit defined in Eq. (13), otherwise uneed is
always less than 2ulimit . As a result, the bandwidth of the proposed method always satisfies
the condition of the minimum bandwidth.
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 ∆u = ( 2S x ) , the SNR has the maximum value for a
−1

bandwidth around 2ulimit , and it decreases rapidly for bandwidths below uneed . In the case of
the over sampling of ∆u = ( 4S x ) , since the numerical error owing to the sampling problem
−1

is not caused, the SNR does not decrease for bandwidths above 2ulimit , and there is no clear
peak in the curve. However, the SNR decreases for bandwidths below uneed as in the normal

#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 19672
sampling. This means that the error in bandwidths below uneed 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.

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 Wx = S x / 2 and z = 50 S x . The other
parameters are the same as in Fig. 4.

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

#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 19673