Title Shifted angular spectrum method for off-axis num erical propagation

A novel method is proposed for simulating free-space propagation from an input source field to a destination sampling window laterally shifted from that in the source field. This off-axis type numerical propagation is realized using the shifted-Fresnel method (Shift-FR) and is very useful for calculating non-paraxial and large-scale fields. However, the Shift-FR is prone to a serious problem, in that it causes strong aliasing errors in short distance propagation. The proposed method, based on the angular spectrum method, resolves this problem. Numerical examples as well as the formulation are presented. ©2010 Optical Society of America OCIS codes: (050.1940) Diffraction; (090.1760) Computer holography; (070.0070) Fourier optics and signal processing. References and links 1. T. M. Kreis, M. Adams, and W. P. O. Jüptner, “Methods of digital holography: A comparison,” Proc. SPIE 3098, 224–233 (1997). 2. F. Zhang, I. Yamaguchi, and L. P. Yaroslavsky, “Algorithm for reconstruction of digital holograms with adjustable magnification,” Opt. Lett. 29(14), 1668–1670 (2004). 3. L. Yu, and M. K. Kim, “Pixel resolution control in numerical reconstruction of digital holography,” Opt. Lett. 31(7), 897–899 (2006). 4. 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). 5. R. P. Muffoletto, J. M. Tyler, and J. E. Tohline, “Shifted Fresnel diffraction for computational holography,” Opt. Express 15(9), 5631–5640 (2007), http://www.opticsinfobase.org/oe/abstract.cfm?URI=oe-15-9-5631. 6. M. Sypek, “Light propagation in the Fresnel region. New numerical approach,” Opt. Commun. 116(1-3), 43–48 (1995). 7. J. W. Goodman, Introduction to Fourier Optics, 2nd ed. (McGraw-Hill, 1996), chap. 3.10. 8. K. Matsushima, and T. Shimobaba, “Band-limited angular spectrum method for numerical simulation of freespace propagation in far and near fields,” Opt. Express 17(22), 19662–19673 (2009), http://www.opticsinfobase.org/oe/abstract.cfm?URI=oe-17-22-19662. 9. 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), http://www.opticsinfobase.org/josaa/abstract.cfm?URI=josaa-20-9-1755. 10. K. Matsushima, “Formulation of the rotational transformation of wave fields and their application to digital holography,” Appl. Opt. 47(19), D110–D116 (2008), http://www.opticsinfobase.org/ao/abstract.cfm?URI=ao-4719-D110. 11. 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). 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. K. Matsushima, “Computer-generated holograms for three-dimensional surface objects with shade and texture,” Appl. Opt. 44(22), 4607–4614 (2005), http://www.opticsinfobase.org/ao/abstract.cfm?URI=ao-44-22-4607. 15. K. Matsushima, and S. Nakahara, “Extremely high-definition full-parallax computer-generated hologram created by the polygon-based method,” Appl. Opt. 48(34), H54–H63 (2009), http://www.opticsinfobase.org/oe/abstract.cfm?URI=oe-17-22-19662. 16. J. W. Goodman, Introduction to Fourier Optics, 2nd ed. (McGraw-Hill, 1996), chap. 2.2. #131557 $15.00 USD Received 12 Jul 2010; accepted 8 Aug 2010; published 13 Aug 2010 (C) 2010 OSA 16 August 2010 / Vol. 18, No. 17 / OPTICS EXPRESS 18453


Introduction
A variety of numerical methods for free space propagation are required for computational holography, such as digital holography (DH) or computer-generated holograms (CGH).DH can be regarded as analog to digital (A/D) conversion of optical wave fields, whereas CGH is digital to analog (D/A) conversion.Evolution of these conversions may contribute to the feasibility of digital signal processing (DSP) of light, in which numerical simulation of freespace propagation plays a role in the signal processing.
The first category of free-space propagation is propagation between parallel planes [1].In this category, various methods, such as the single Fourier-transform-based Fresnel method (SFT-FR) [1][2][3][4][5] and the convolution-based Fresnel method (CV-FR) [1,6], are continually being proposed.However, the angular spectrum-based method (AS) [7] is potentially the most powerful, because it is rigorously derived from the Rayleigh-Sommerfeld integral.However, the traditional AS cannot serve as an all-round method in a numerical implementation due to sampling problems.The author recently proposed the band-limited AS (BL-AS) to avoid the sampling problems [8].This is a simple, yet effective, improvement of the AS that magnifies the range of the effective propagation distance of the AS.
Another category of free space propagation is propagation between non-parallel planes [9][10][11][12].The AS also plays an important role in this category.The rotational transformation of wave fields [9,10], formulated as an expansion of the AS, makes it possible to calculate wave fields in arbitrarily tilted planes from a given source field.This method is used for clear imaging of deeply tilted surfaces [10,13] and forming surface sources of light in the polygonbased method for creating CGHs of surface-modeled objects [14].
Recently, a new category of free space propagation, off-axis numerical propagation, was added to the field.This involves propagation between parallel planes, but with the sampling window of the output destination field shifted from that of the input source field.This is very useful in cases where a field is not paraxial and travels in an off-axis direction.In such cases, the sampling window needs to be expanded in conventional methods, as described in Section 2. This generally leads to a much higher computational effort.The techniques of off-axis numerical propagation can remove this difficulty.Off-axis numerical propagation is also useful for the propagation of large-scale wave fields that cannot be stored simultaneously in main memory due to the size of the data.An extremely high-definition CGH, reconstructing true fine 3D images, can also be calculated using off-axis numerical propagation [15].
The most notable method for off-axis numerical propagation is the shifted Fresnel method (Shift-FR) [5].This excellent technique is derived from the SFT-FR using a scaled FFT and includes the capabilities of off-axis propagation and a variable sampling interval.However, the Shift-FR has a serious problem in that it causes strong aliasing when applied to short distance propagation for given sampling intervals and field sizes, as discussed in Section 4.
The novel method for off-axis numerical propagation proposed in this paper is a generalization of the BL-AS introduced in [8].Therefore, the proposed method, called the shifted angular spectrum method (Shift-AS), does not cause any aliasing errors in short distance propagation, unlike the Shift-FR.The formulation thereof and numerical examples are presented in the following sections.

Off-axis numerical propagation
In numerical propagation of wave fields, the destination field is sometimes required to be in a region set apart from the optical axis in the destination plane.This situation is depicted in Fig. 1.The source sampling area is placed around the origin of the coordinate system, whereas the region of interest in the destination plane is located apart from the origin.In this case, the destination field is generally calculated by the following procedure: Step 1. Extend the source sampling window so that the sampling window includes the region of interest after numerical propagation.
Step 2. Pad the extended source sampling window with zeros.
Step 3. Numerically propagate the source wave field onto the destination plane.
Step 4. Cut out the region of interest from the destination sampling window.
Here, in Step 3, the source wave field is propagated using conventional numerical methods that do not change the sampling interval, such as the CV-FR, AS, or BL-AS.The diffracted field can be calculated using this procedure.However, the procedure usually requires a huge computational effort, especially in cases where the region of interest is far from the optical axis.
On the other hand, as shown in Fig. 1(b), a numerical method can compute directly the diffracted field in a sampling window set apart from the optical axis without extending the sampling window.This type of propagation is referred to as off-axis numerical propagation in this paper.The most notable method for off-axis numerical propagation is the Shift-FR [5], which makes it possible to change the sampling distance and the center of the sampling area in the destination plane by applying a scaled FFT to the SFT-FR method.Although this is an excellent method, it causes a serious problem in short-distance propagation as discussed in Section 4.

The angular spectrum method for shifted coordinates
The coordinate system used for the formulation is shown in Fig. 2. The source field ( , , 0) g x y is given in the source plane ( , , 0) x y , while the destination field 0 ˆ( , , ) g x y z is given in the destination plane 0 ˆ( , , ) x y z , in that the origin of the lateral coordinates x and ŷ is shifted from the source coordinates as follows: Without any shifts ( 0 0 0 x y = = ), the spectrum of the destination field is given by a convolution form of the Rayleigh-Sommerfeld integral (see Section 2.1 in [8]).
g x y z g x y z g x y z g x y h x x y y z dx dy where 0 ( , , ) h x y z is the propagation kernel (the impulse response) of the Rayleigh-Sommerfeld formula.This convolution is rewritten using the convolution theorem as: where u and v are Fourier frequencies with respect to x and y , respectively.The spectrum of the source field is given as: Here, F represents the Fourier transform.The transfer function 0 ( , ; ) H u v z is given by: [ ] ( ) where λ denotes the wavelength.

g x y z g x x y y z g x y h x x x y y y z dx dy
This can be rewritten using the convolution theorem and shift theorem of the Fourier transform as: Consequently, we obtain the destination field as follows: { } Nevertheless, simple numerical implementations of the above procedure are most likely to cause considerable sampling problems as described in the next section.

One-dimensional wave fields
For simplicity, one-dimensional wave fields that are a function of x are discussed in this section.In this case, the transfer function of Eq. ( 7) is redefined as follows: [ ] ( ) The local signal frequency of the transfer function is given by [16]: where the local signal frequency denotes neither physical frequencies of time nor space, but the frequency of peaks and valleys in the function 0 ( ; ) H u z in a certain period u.Therefore, supposing that the transfer function is sampled at intervals of u ∆ , the Nyquist theorem requires the following relation to avoid sampling problems: By substituting Eq. ( 10) into Eq.( 11), the Nyquist condition is rewritten as: Finally, the conditions are summarized in the following three cases: ( ) ( ) limit limit 0 ( ) ( ) limit limit 0 ( ) is the size of the source sampling window.Here, it should be noted that the size of the sampling window is given by 1 (2 ) , and not by source sampling window is doubled to linearize the discrete convolution, as described in Section 2.2 in [8].The constant ( ) limit u ± is given by Where 0 0 x = , this agrees with  defined in Eq. ( 13) in [8].(2 ) (2 ) , where 1024 , and 0 20 x z S = .

Table 1. Constants used for the band-limit to avoid aliasing errors
An example of the sampled function of 0 ˆ( ; ) H u z is depicted in Fig. 3. Here, the center of the destination sampling window is shifted by an amount half the size of the sampling window in this example, i.e., 0 / 2 x x S = + is assumed.Since aliasing errors occur below ( ) limit u − − or above ( ) limit u + , in this case the sampled transfer function and the source field should be bandlimited within ( ) ( ) limit limit u u u − + − < < to avoid aliasing errors.Consequently, the transfer function for shifted coordinates is given by: where rect( ) ξ is a rectangular function with unity width.The constants 0 u and width u are given in Table 1.

Two-dimensional wave fields
In two-dimensional cases, the Nyquist condition is given as two-dimensional regions in the ( , ) u v plane.The local signal frequencies of the function where . The Nyquist conditions are given by: (2 ) is the sampling interval and S y is the size of the source sampling window in the y direction.Note that the sampling window is also doubled in the y direction to linearize the discrete convolution.
Both conditions Eqs. ( 17) and ( 18) should be satisfied simultaneously to avoid sampling problems in two-dimensional wave fields.Since these conditions are very similar, let us formulate the region that satisfies only condition Eq. ( 17) in the ( , ) u v space as follows: Case (i): Case (ii): 1 and 0 or 1 and 0 , Case (iii): 0 2 2 2 2 ( ) ( ) limit limit 0 and 1 and 1.
Here, the region satisfying the other condition [Eq.( 18)] is simply obtained by switching the symbols x and u to y and v in the above relations, respectively.Note that the region under the conditions 0 0 x = or 0 0 y = , that is, the special case of (ii) is given by: These agree with relations Eqs. ( 18) and ( 19) in [8].
Case (i) : Relations Eqs. ( 19)-(21) for shifting 0 x give a region specified by the combination of vertical ellipsoidal regions with a major diameter 2λ −1 in the (u, v) plane, while the relations for shifting 0 y give horizontal ellipsoidal regions.Figure 4 shows these ellipsoidal regions in cases where 0 0 x ≠ and 0 0 y = .In this figure, the overlap regions that satisfy both Nyquist conditions Eqs. (17) and (18) are schematically depicted as red-hatched areas.The sampled transfer function should be limited in the overlap region to avoid sampling problems.This is equivalent to limiting the bandwidth of the source field within the overlap region.
Here, note that the inside ellipse in Case (i) is switched to the outside ellipse in Case (iii).This switching occurs when 0 x changes its sign in Case (ii).Although Eq. (20) seems, in a sense, to be written as a separate form, the mathematical expression of Case (ii) is the same, irrespective of the sign of 0 x .

Approximated rectangular region to avoid aliasing errors
Regions that avoid sampling problems are considerably more complicated in two-dimensional wave fields.However, in most cases where 0 and , the two-dimensional region can be approximated as a simple rectangular region by combining the one-dimensional relations of Eq. ( 13).An example of the approximated region is illustrated schematically in Fig. 5. Finally, the sampled transfer function for two-dimensional wave fields is approximately given by: where the additional constants are listed in Table 1 with:

Numerical verification of the shifted angular spectrum method
The optical setup used for numerical verification is shown in Fig. 6(a).We assume that a plane wave travels at an incident angle of θ .Supposing that the plane wave is diffracted by a circular aperture, we calculate the diffracted field in the plane at position 0 z .Here, the source field including the plane wave and the aperture is sampled at intervals of   The diffracted field is calculated using three different methods: the band-limited angular spectrum (BL-AS), shifted Fresnel (Shift-FR), and shifted angular spectrum method (Shift-AS) proposed in this paper.Since the BL-AS is not a method for off-axis numerical propagation, the number of sampling points in the source field is doubled both on the x-and yaxes so that the diffracted field is included within the destination sampling window.In addition, the approximated transfer function given by Eq. ( 23) is used for calculating the Shift-AS.
The calculated destination field is shown in Fig. 7(a)-7(c) for three positions of z 0 .The shift of the destination sampling window is x 0 = 0.4 cm both for z 0 = 5 cm in (a) and 10 cm in (b), whereas x 0 = 1 cm for z 0 = 40 cm in (c).In the case of z 0 = 5 cm in (a), the position of the destination field is shifted slightly along the x-axis because of the slanted plane wave.In this case, the field calculated by the Shift-AS agrees with that by the BL-AS, whereas the field by the Shift-FR disagrees with these fields because of the aliasing error.The field calculated by the Shift-FR is also not normal in z 0 = 10 cm in (b), while the BL-AS and Shift-AS again give similar fields.
In the case of z 0 = 40 cm in (c), the destination field lies outside the extended sampling area of the BL-AS.However, the Shift-AS and Shift-FR can calculate the field within the sampling area owing to their off-axis property.The Shift-FR, in this case, gives the same result as that using the Shift-AS, due to the long distance of the propagation.

Discussion on the limit frequency
We limit the spatial frequencies of the source fields to a band specified by ( ) limit u ± .This limit frequency can be physically interpreted as the local spatial frequency of the field emitted from a point within the source sampling window, as shown in Fig. 8.The highest spatial frequency is given by the field emitted from the point at the lower end of the source sampling window to the point at the upper end of the destination window.The lowest frequency is given in the same manner.Here, angles max x S u − = ∆ , the highest and lowest spatial frequencies are given by: As a result, the limits imposed on the frequencies by the procedure to avoid sampling problems are in agreement with the highest and lowest frequencies.Therefore, we can interpret the condition specified by ( ) limit u ± as the frequency band that is required to physically propagate the source wave field onto the shifted destination area.

Conclusion
A numerical technique called the shifted angular spectrum method is proposed for off-axis numerical propagation.The Shift-AS is a generalization of the band-limited angular spectrum method, and therefore, the method does not have any defined restrictions on the propagation distance.This new method makes it possible to calculate exact wave fields without any aliasing errors, such as those caused by the shifted Fresnel method.

Fig. 1 .
Fig. 1.Calculation of the diffracted field using (a) conventional methods and (b) methods for off-axis propagation, where the region of interest is apart from the optical axis in the destination pane.

Fig. 4 .
Fig. 4. Schematic illustrations of frequency regions to avoid aliasing errors of the sampled transfer function.Here, 0 0 x ≠ and 0 0 y = .

#(Fig. 5 .
Fig. 5. Schematic illustrating the approximated rectangular region for the band-limit of a twodimensional wave field.Here, the region is depicted with 0 x x S > + and 0 The amplitude and phase images of the sampled source field are shown in Fig. 6(b).

Fig. 6 .Fig. 7 .
Fig. 6.Setup for (a) numerical simulation of off-axis propagation, and (b) the sampled source field used in the simulation.

Fig. 8 .
Fig. 8. Model for discussing the highest and lowest spatial frequencies required for off-axis numerical propagation.
θ and min θ in Fig.8are given by the geometry of the source and destination window as follows: