Angular spectrum calculations for arbitrary focal length with a scaled convolution

Satoru Odate; Chiaki Koike; Hidemitsu Toba; Tetsuya Koike; Ayako Sugaya; Katsumi Sugisaki; Katsura Otaki; Kiyoshi Uchikawa

doi:10.1364/OE.19.014268

1. Introduction

The angular spectrum (AS) method is known as one of the most precise calculations in the scalar propagation of a wave field [1,2]. This method can be simply computed very fast by using the fast Fourier transform (FFT) twice. However, there has been a shortcoming in the FFT-based AS computation due to the Nyquist sampling theorem [3]. The FFT-based AS calculation in case of long distance propagation causes aliasing errors. So far, many authors have studied the application of its calculation only in the short distance propagation. Recently, the calculation for long distance was successfully computed by expressing the AS formula as a convolution representation [4]. This convolution representation is indirectly to be performed the calculation by the FFT (e.g., the convolution theorem). Thus, the AS method overcame the restriction of the propagation distance.

Similarly to the propagation restriction, when one tries to compute the AS of an image formed by a lens, the phase distribution generated by the lens happens to cause aliasing errors too. In such computation, the focal length of the lens is bounded because of the sampling condition. M. Mansuripur nicely introduced the parabolic function so that the phase oscillation of the lens is suppressed [5]. To minimize the oscillation, the phase distribution was extracted from the parabolic function multiplied by a parameter. After that, he chose a value of the parameter to minimize the maximum absolute derivative value of the phase. As a result, the applicable range of the focal length was extended fairly but still depending on the sampling interval. By contrast, if one chooses to use the Fresnel formula instead of the AS method for the propagation calculation, that limitation can be resolved by applying a suitable convolution representation [6]. In addition to this, a focal pivot method based on the Fresnel formula can perform the image calculation formed by the lens in any defocus plane [7]. Thus, the Fresnel formula has basically no limitations regarding the sampling problems. Unfortunately, it is hardly applied to the image computation for a high-NA lens. One must utilize the AS formula in the case of the high-NA lens.

In this work, we report an AS method of the image computation formed by the lens with arbitrary focal length. To defeat the sampling problems with the AS method, suitable convolution formulas are derived and proposed. To the best of our knowledge, this work is the first study to avoid the aliasing error in the AS method in lens computations. Furthermore, a scaled convolution for each formula can magnify the lateral image in an observation plane. Accordingly, our AS method may compute the image generated from the lens of the arbitrary focal length in any distant observation plane with any small lateral resolution.

2. Aberration-free lens model

Figure 1 shows an aberration-free lens model of a focal length f. A monochromatic plane wave field parallel to the optical axis z transmits an aperture in ξ-η plane. A numerical aperture (NA) of the lens depends on the diameter d for a lens aperture as follows

Fig. 1 A schematic model of an aberration-free lens. ξ denotes vector (ξ,η). Details are shown in the text.

Download Full Size | PDF

NA = \frac{d / 2}{f} .

Note that the shorter f, the larger NA if d is constant. The wave field $A (x : z)$ propagating from the lens at a distance z is calculated by the AS method as

a (u) = \int \frac{f^{2}}{ξ \cdot ξ + f^{2}} T (ξ) A (ξ) \exp (- 2 π i u \cdot ξ) d ξ^{2},

A (x : z) = \int a (u) \exp (2 π i z \sqrt{λ^{- 2} - u \cdot u}) \exp (2 π i u \cdot x) d u^{2},

where $A (ξ)$ is the wave field on the aperture plane just before the lens (See Appendix A). Throughout this paper, boldfaced letters represent two-dimensional vectors. ξ and x denote vectors in ξ-η (image plane) and x-y plane (observation plane), respectively. u is a spatial-frequency vector against ξ vector. The function T over the aperture is defined by

T (ξ) = \exp [- \frac{2 π i}{λ} (\sqrt{ξ \cdot ξ + f^{2}} - f)] .

One can easily perform these calculations by means of the FFT.

Let us apply the sampling theorem to these equations in order to investigate the required conditions [4]. The resultant conditions are given by

z \leq \frac{M {(Δ ξ)}^{2}}{λ} \begin{matrix} , \end{matrix} f \geq \frac{M {(Δ ξ)}^{2}}{λ},

where Δξ is a sampling interval on the aperture plane, the area of the sampling window is defined by $L \times L$ square, the number of samples is $M^{2}$ (Δξ = L/M) . Here we assumed that the phase of $A (ξ)$ is slowly varying function. When either of the inequalities is violated, the phase of the integrand of either Eq. (2) or (3) rapidly changes. As a result, an aliasing error occurs because the sampling interval is not sufficiently small [4,5]. Hence, the computations in Eqs. (2) and (3) can be performed only in short distance $z \leq M {(Δ ξ)}^{2} / λ$ with long focal length $f \geq M {(Δ ξ)}^{2} / λ$ .

Conversely, a convolution representation of Eq. (3) (Rayleigh-Sommerfeld representation) enables us to compute for long distance $z \geq M {(Δ ξ)}^{2} / λ$ [4], namely

A (x : z) = \int A (ξ) H (x - ξ) d ξ^{2},

H (ξ) = \frac{z \exp (2 π i f \sqrt{ξ \cdot ξ + z^{2}})}{\sqrt{ξ \cdot ξ + z^{2}}} [\frac{1}{2 π \sqrt{ξ \cdot ξ + z^{2}}} + i \frac{1}{λ}] .

When computing these formulas by the FFT (e.g., the convolution theorem), one must expand the area of the sampling window to four times its original size. The extended sampling area is padded with zeros so that the aliasing error is suppressed effectively [4].

Furthermore, by using to the derivation of Eqs. (6) and (7), let us perform the calculation for a short focal length $f \leq M {(Δ ξ)}^{2} / λ$ . This inequality indicates that the sampling interval Δξ on the aperture plane is not sufficiently small. In this case, it is necessary to rewrite Eq. (2) as a convolution representation. Specifically, we derived such formulas:

a (u) = \int g (p) t (u - p : f) d p^{2},

t (p : f) = \frac{\exp [- 2 π i f (\sqrt{λ^{- 2} - p \cdot p} - \frac{1}{λ})]}{λ^{- 2} - p \cdot p} [\frac{1}{2 π λ \sqrt{λ^{- 2} - p \cdot p}} + i \frac{f}{λ}],

g (p) = \int \frac{f^{2}}{ξ \cdot ξ + f^{2}} A (ξ) \exp (- 2 i π p \cdot ξ) d ξ^{2} .

p is a spatial-frequency vector against ξ vector. Here we utilized Weyl’s representation of a spherical wave (See appendix B for details). When computing these formulas by the FFT (e.g., the convolution theorem), one must expand the sampling area of $g (p)$ to four times its original size with zero-padding in order to suppress the aliasing error effectively. After all, we found the suitable formulas expressing by Eqs. (2), (3), (6), (8) depending on the sampling conditions (i.e., f ≥ M(Δξ) ²/λ, z ≤ M(Δξ)²/λ, z ≥ M(Δξ)²/λ, f ≤ M(Δξ)²/λ, respectively). Thus, one can make it possible to compute the wave field $A (x : z)$ generated by the lens of arbitrary focal length for an arbitrary propagation distance.

3. Scaled convolution

The sampling intervals in the aperture and the observation plane are basically the same spacing (Δξ = Δx) if one performs the AS computation by means of the FFT. Therefore, such method is not adequate for practical computation of the spot image near the focal plane. Let us modify the formulas of the propagations so that the lateral scale of the image can be magnified by reducing the sampling interval Δx in the observation plane. Equation (6) of the convolution representation for long distance propagation is directly expressed as

A (x : z) = \int a (u) [\int H (ξ) \exp (- 2 π i u \cdot ξ) d ξ^{2}] \exp (2 π i u \cdot x) d u^{2},

which can be implemented for computation by the FFT adoptions of three times. On the other hand, Eq. (3) for short distance propagation computation is carried out by the FFT also. Consequently, let us consider the Fourier transform for both propagation computations:

f (x - x_{0}) = \int F (u) \exp [2 π i u \cdot (x - x_{0})] d u^{2} .

The discrete integral which the center of the sampling area is set to be origin of the axes is given by

f [(m - M) Δ x - x_{0}] = \sum_{p = 0}^{2 M - 1} F_{p - M} \exp [2 π i (p - M) \cdot ((m - M) Δ x - x_{0}) Δ u],

where mand q are two-dimensional vectors of integers, $M = (M, M)$ . Here we extended the sampling area to be four times for the purpose of transforming this equation to the convolution representation. To this end, we converted Eq. (13) into

\begin{array}{l} f [(m - M) Δ x / α - x_{0}] = \exp [\frac{π i}{γ M} (m \cdot m - 2 M \cdot m + 4 M^{2})] \sum_{p = 0}^{2 M - 1} F_{p - M} \\ \times \exp (\frac{π i}{γ M} (p \cdot p - 2 M \cdot p) - \frac{2 π i}{L} (p - M) \cdot x_{0}) \exp [- \frac{π i}{γ M} {(m - p)}^{2}], \end{array}

where the sampling intervals are defined as $Δ u = 1 / L$ and $Δ u Δ x = 1 / γ M$ [8]. Notice that the parameter γ indicates the ability to vary the total number of data for the Discrete Fourier Transform in one period. The parameter $γ \geq 1$ satisfies the following relation:

\frac{Δ ξ}{Δ x} = γ .

Here we assume that the identity $Δ u Δ x = 1 / M$ is fulfilled because the computation of Eq. (2) is performed by FFT. Let us note that the larger γ, the higher the magnification of the lateral image in the observation plane. Thus, one can simply achieve the computation of Eq. (14) by the FFT adoptions of three times:

\begin{array}{l} f [(m - M) Δ x / γ - x_{0}] = \exp [\frac{π i}{γ M} (m \cdot m - 2 M \cdot m + 4 M^{2})] \\ \times FFT (IFFT [F_{p - M} \exp [\frac{π i}{γ M} (p \cdot p - 2 M \cdot p) - \frac{2 π i}{L} (p - M) \cdot x_{0}]] \\ \times IFFT [\exp (- \frac{π i}{γ M} {(p - M)}^{2})]) . \end{array}

To avoid the aliasing error, the zero-padding in the extended area must be performed. After the computation, one must cut and reduce the sampling window to the original size. Finally, it is possible to magnify the lateral image by applying those methods to both Eqs. (3) and (6).

So far, we limited the range of the parameter $γ \geq 1$ . If one chooses to use the parameter range $γ < 1$ , the original sampling data are reduced and lost in the computation. Consequently, aliasing error can occur. Thus, the scaled convolution is permitted to vary the lateral scale of the image not for the scale-down but for the scale-up. In order to enable us to do the scale-down of the image, let us note that Eqs. (2) and (3) are invariants under transformations of $u \to γ u$ , $ξ \to ξ / γ$ , $x \to x / γ$ , $z \to z / γ$ , $λ \to λ / γ$ . If one applies the angular spectrum $a (γ u)$ to the computation of Eq. (3), the calculated scalar wave field can be regarded as $A (x / γ : z / γ)$ . The computation of the $a (γ u)$ can be simply achieved when one considers the coordinates of the angular spectrum as u after computing the scaled convolution of Eq. (3). Thus, one can obtain the scale-down image of the scalar wave field in the observation plane.

Let us reconsider the computation method of Eq. (8). $g (p)$ has the non-zero values only in a small sampled area because the number of samples is small ( $M \leq L^{2} / (f λ)$ ). Therefore, the aliasing error may occur in the computation of $a (u)$ . To avoid this, one should expand the distribution of $g (p)$ . The scaled convolution applied to Eq. (10) can accomplish such expansion as we mentioned above. Note that the calculated distribution of $A (x : z)$ is magnified, as the scaled convolution expands the $g (p)$ .

Table 1 shows the summary of the suitable formulas for the image calculations depending on the conditions of the propagation distance z and the focal length f. It also indicates original data requiring for the zero-padding in each scaled convolution.

Table 1. Summary of the suitable formulas and the conditions

View Table

4. Numerical verifications

To verify the scaled convolution, let us consider the magnification of the Airy disc image. When computing the image near the focal plane, the propagation distance z and focal length f should take almost the same values. According to Table 1, one should use the Type-II or-III for this computation. Figure 2 shows the magnifications of the Airy disc images using Type-II and -III. We computed the Airy disc images of the aberration-free lenses having NA 0.2 and 0.9 at the focal planes by using Type-II and -III, respectively. The images where the scale convolutions are not applied are shown for comparison (1x magnification). As expected, we can clearly confirm the magnifications of these Airy disc images. Furthermore, we can hardly see the fine features of Airy disc images, when we expand the images of 1x magnifications by digitally scale-up. Thus, we verified the capability of the scale convolution for the magnification.

Fig. 2 Airy disc images using the scaled convolutions (λ = 0.6 um, F = 2⁸ um, M = 2¹²-1). (a) Type-II: NA = 0.2, magnification 1x. (b) Type-II: NA = 0.2, magnification 20x. (c) Type-III: NA = 0.9, magnification 1x. (d) Type-III: NA = 0.9, magnification 625x.

Download Full Size | PDF

Also, let us verify that one must use to select Eq. (2) or (8) depending on the focal length f. When we consider an aberration-free lens having f = 11.456 mm, NA = 0.4, L = 10 mm and D = 9.999 mm in λ = 0.6 um with M = 2¹⁴, threshold value $M {(Δ ξ)}^{2} / λ$ of the inequality becomes 10.173 mm. The suitable formulas in these conditions require choosing Type-I or -II. We calculated the intensity distributions of 1 D profile in an observation plane far from the focal plane in distance −10 mm. Figure 3 shows the resultant intensity distributions, when we apply this to all formulas from Type-I to -IV. As expected, Type-I shows the correct intensity distribution. The other distributions excluding Type-I happen to cause aliasing error.

Fig. 3 The intensity distributions in the observation planes far from the focal planes in distance −10 mm using (a) Type-I. (b) Type-II. (c) Type-III. (d) Type-IV, when aberration-free lens having f = 11.456 mm, NA = 0.4, L = 10 mm and D = 9.999 mm in λ = 0.6 um with M = 2¹⁴.

Download Full Size | PDF

Let us continue the verification for high-NA lens. When we consider an aberration-free lens having f = 0.128 mm, NA = 0.9, L = 0.778 mm and D = 0.529 mm in λ = 0.6 um with M = 2¹⁰, threshold value $M {(Δ ξ)}^{2} / λ$ of the inequality becomes 0.988 mm. The suitable formulas in these conditions require choosing Type-III or -IV. We calculated the intensity distributions of 1 D profile in an observation plane far from the focal plane in distance 0.1 mm. Figure 4 shows the resultant intensity distributions, when we apply to use all formulas from Type-I to -IV. As expected, Type-III shows correct intensity distribution. The other distributions excluding Type-III happen to cause aliasing error. So, we confirmed that Eq. (8) which is derived in this paper is the proper expression for the image computation generated by the aberration-free lens having short focal length. Thus, our AS method accomplish the image computation formed by the lens of the arbitrary focal length in any distant observation plane with any small lateral resolution.

Fig. 4 The intensity distributions in the observation planes far from the focal planes in distance 0.1 mm using (a) Type-I. (b) Type-II. (c) Type-III. (d) Type-IV, when aberration-free lens having f = 0.128 mm, NA = 0.9, L = 10 mm and D = 0.779 mm in λ = 0.6 um with M = 2¹⁴.

Download Full Size | PDF

5. Summary and Discussion

We proposed and derived the suitable AS formulas for the image computation of the lens. We verify that our AS method can perform the computation of the image generated by the aberration-free lens of arbitrary focal length at any distant observation plane with any small lateral resolution. Let us note that one must consider a vector wave field calculation instead of the scalar wave field for the case when the high-NA lens calculation. To extend the AS method for the vector field, one must modify their formulas. Such modification can be easily accomplished [5].

Appendix A

Let us consider the energy flow over the solid angle dΩ:

d Ω = \sin θ d θ \land d φ,

where the solid angle is expressed by the wedge product on the spherical coordinates (θ, ϕ) [9]. This solid angle can be represented in tangent (ξ, η) and sine (α, β) coordinates as

d Ω = \frac{1}{\cos θ} d α \land d β = \cos^{3} θ d ξ \land d η,

where we defined tangent and sine coordinates as follows:

\begin{array}{l} ξ = f \tan θ \cos φ \\ η = f \tan θ \sin φ \end{array}

\begin{array}{l} α = f \sin θ \cos φ \\ β = f \sin θ \sin φ . \end{array}

By comparing the identity, we obtain the relation of the surface elements between tangent and sine coordinates,

d α \land d β = \cos^{4} θ d ξ \land d η .

In case of the scalar field, we regard its relation as

d α \land d β = \cos^{2} θ d ξ \land d η .

Then, the relation for the scalar field is given by

d α \land d β = \frac{f^{2}}{ξ \cdot ξ + f^{2}} d ξ \land d η .

Thus, we found the weight function in tangent coordinates for the AS calculation.

Appendix B

Let us calculate $t (p : f)$ , that is, the Fourier transform of $T (ξ)$ :

\begin{array}{l} t (p : f) = \int \exp [- i k (\sqrt{ξ \cdot ξ + f^{2}} - f)] \exp (- 2 i π p \cdot ξ) d ξ^{2} \\ \begin{matrix}  \end{matrix} = \exp (i k f) \frac{\partial}{\partial k} \int \frac{i}{\sqrt{ξ \cdot ξ + f^{2}}} \exp [- i k \sqrt{ξ \cdot ξ + f^{2}}] \exp (- 2 i π p \cdot ξ) d ξ^{2} . \end{array}

Then, this formula can be converted as follows:

\begin{array}{l} t (p : f) = \exp (i k f) \frac{\partial}{\partial k} \int d s^{2} d ξ^{2} \frac{- 1}{\sqrt{λ^{- 2} - s \cdot s}} \\ \begin{matrix}  \end{matrix} \times \exp [- 2 π i f \sqrt{λ^{- 2} - s \cdot s}] \exp [- 2 i π (p + s) \cdot ξ], \end{array}

where we utilized the Weyl’s representation of a spherical wave:

\frac{\exp (- i k \sqrt{ξ \cdot ξ + f^{2}})}{\sqrt{ξ \cdot ξ + f^{2}}} = \int \frac{i}{\sqrt{λ^{- 2} - s \cdot s}} \exp [- 2 π i (s \cdot ξ + f \sqrt{λ^{- 2} - s \cdot s})] d s^{2} .

By calculating a derivative at k, we obtain

t (p : f) = \frac{\exp [- 2 π i f (\sqrt{λ^{- 2} - p \cdot p} - \frac{1}{λ})]}{λ^{- 2} - p \cdot p} [\frac{1}{2 π λ \sqrt{λ^{- 2} - p \cdot p}} + i \frac{f}{λ}] .

Finally, Eq. (4) can be described as convolution representation:

a (u) = \int \exp [- i k (\sqrt{ξ \cdot ξ + f^{2}} - f)] A (ξ) \exp (- 2 π i u \cdot ξ) d ξ^{2} = \int g (p) t (u - p : f) d p^{2},

where we defined the Fourier transform of $A (ξ)$ as $t (p)$ :

g (p) = \int A (ξ) \exp (- 2 i π p \cdot ξ) d ξ^{2} .

Acknowledgments

We wish to thank Takashi Gemma, Yucong Zhu, Zhiqiang Liu, Takashi Mori and Hiroshi Ohki for fruitful discussions.

References and links

1. L. Mandel and E. Wolf, Optical Coherence and Quantum Optics (Cambridge University Press, 1995).

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

3. K. Matsushima and T. Shimobaba, “Band-limited angular spectrum method for numerical simulation of free-space propagation in far and near fields,” Opt. Express 17(22), 19662–19673 (2009). [CrossRef] [PubMed]

4. 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]

5. M. Mansuripur, “Certain computational aspects of vector diffraction problems,” J. Opt. Soc. Am. A 6(6), 786–805 (1989). [CrossRef]

6. D. Mas, J. Garcia, C. Ferreira, L. M. Bernardo, and F. Marinho, “Fast algorithms for free-space diffraction patterns calculation,” Opt. Commun. 164(4-6), 233–245 (1999). [CrossRef]

7. S. Nishiwaki, “Calculations of optical field by fast Fourier transform analysis,” Appl. Opt. 27(16), 3518–3521 (1988). [CrossRef] [PubMed]

8. 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]

9. M. Nakahara, Geometry, Topology and Physics, 2nd ed. (Taylor & Francis, 2003).

Type	I	II	III	IV
Propagation Distance z	$z \leq \frac{M {(Δ ξ)}^{2}}{λ}$	$z \geq \frac{M {(Δ ξ)}^{2}}{λ}$	$z \leq \frac{M {(Δ ξ)}^{2}}{λ}$	$z \geq \frac{M {(Δ ξ)}^{2}}{λ}$
Focal length f	$f \geq M {(Δ ξ)}^{2} / λ$	$f \geq \frac{M {(Δ ξ)}^{2}}{λ}$	$f \leq M {(Δ ξ)}^{2} / λ$	$f \leq \frac{M {(Δ ξ)}^{2}}{λ}$
Suitable formulas	(2), (3)	(2), (6)	(8), (3)	(8), (6)
Zero padding	$a (u)$	$A (ξ)$	$g (p)$ , $a (u)$	$g (p)$ , $A (ξ)$

Angular spectrum calculations for arbitrary focal length with a scaled convolution

Abstract

1. Introduction

2. Aberration-free lens model

3. Scaled convolution

4. Numerical verifications

5. Summary and Discussion

Appendix A

Appendix B

Acknowledgments

References and links

Cited By

Figures (4)

Tables (1)

Equations (29)

Optics Express