OSA's Digital Library

Optics Express

Optics Express

  • Editor: C. Martijn de Sterke
  • Vol. 19, Iss. 15 — Jul. 18, 2011
  • pp: 14268–14276
« Show journal navigation

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, and Kiyoshi Uchikawa  »View Author Affiliations


Optics Express, Vol. 19, Issue 15, pp. 14268-14276 (2011)
http://dx.doi.org/10.1364/OE.19.014268


View Full Text Article

Acrobat PDF (1074 KB)





Browse Journals / Lookup Meetings

Browse by Journal and Year


   


Lookup Conference Papers

Close Browse Journals / Lookup Meetings

Article Tools

Share
Citations

Abstract

Nyquist sampling theorem in an image calculation with angular spectrum method restricts a propagation distance and a focal length of a lens. In order to avoid these restrictions, we studied suitable expressions for the image computations depending on their conditions. Additionally, a lateral scale in an observation plane can be magnified freely by using a scaled convolution in each expression.

© 2011 OSA

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

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

,2

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

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

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]

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

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]

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

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

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

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]

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

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

]. 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
Fig. 1 A schematic model of an aberration-free lens. ξ denotes vector (ξ,η). Details are shown in the text.
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

NA=d/2f.
(1)

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)=f2ξξ+f2T(ξ)A(ξ)exp(2πiuξ)dξ2,
(2)
A(x:z)=a(u)exp(2πizλ2uu)exp(2πiux)du2,
(3)

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[2πiλ(ξξ+f2f)].
(4)

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

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]

]. The resultant conditions are given by

zM(Δξ)2λ,fM(Δξ)2λ,
(5)

where Δξ is a sampling interval on the aperture plane, the area of the sampling window is defined by L×L square, the number of samples is M2 (Δξ = 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

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

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

]. Hence, the computations in Eqs. (2) and (3) can be performed only in short distance zM(Δξ)2/λ with long focal length fM(Δξ)2/λ.

Conversely, a convolution representation of Eq. (3) (Rayleigh-Sommerfeld representation) enables us to compute for long distance zM(Δξ)2/λ [4

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]

], namely

A(x:z)=A(ξ)H(xξ)dξ2,
(6)
H(ξ)=zexp(2πifξξ+z2)ξξ+z2[12πξξ+z2+i1λ].
(7)

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

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]

].

Furthermore, by using to the derivation of Eqs. (6) and (7), let us perform the calculation for a short focal lengthfM(Δξ)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)=g(p)t(up:f)dp2,
(8)
t(p:f)=exp[2πif(λ2pp1λ)]λ2pp[12πλλ2pp+ifλ],
(9)
g(p)=f2ξξ+f2A(ξ)exp(2iπpξ)dξ2.
(10)

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ξ) 2/λ, z Mξ)2/λ, z Mξ)2/λ, f Mξ)2/λ, 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)=a(u)[H(ξ)exp(2πiuξ)dξ2]exp(2πiux)du2,
(11)

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(xx0)=F(u)exp[2πiu(xx0)]du2.
(12)

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

f[(mM)Δxx0]=p=02M1FpMexp[2πi(pM)((mM)Δxx0)Δu],
(13)

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

f[(mM)Δx/αx0]=exp[πiγM(mm2Mm+4M2)]p=02M1FpM×exp(πiγM(pp2Mp)2πiL(pM)x0)exp[πiγM(mp)2],
(14)

where the sampling intervals are defined as Δu=1/L and ΔuΔx=1/γM [8

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]

]. Notice that the parameter γ indicates the ability to vary the total number of data for the Discrete Fourier Transform in one period. The parameter γ1 satisfies the following relation:

ΔξΔx=γ.
(15)

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:

f[(mM)Δx/γx0]=exp[πiγM(mm2Mm+4M2)]×FFT(IFFT[FpMexp[πiγM(pp2Mp)2πiL(pM)x0]]×IFFT[exp(πiγM(pM)2)]).
(16)

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γ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γu,ξξ/γ, xx/γ, zz/γ, λλ/γ. 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 (ML2/(fλ)). Therefore, the aliasing error may occur in the computation of a(u). To avoid this, one should expand the distribution ofg(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

Table 1. Summary of the suitable formulas and the conditions

table-icon
View This Table
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.

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
Fig. 2 Airy disc images using the scaled convolutions (λ = 0.6 um, F = 28 um, M = 212-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.
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.

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 = 214, 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
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 = 214.
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.

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 = 210, 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
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 = 214.
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.

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

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

].

Appendix A

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

dΩ=sinθdθdφ,

where the solid angle is expressed by the wedge product on the spherical coordinates (θ, ϕ) [9

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

]. This solid angle can be represented in tangent (ξ, η) and sine (α, β) coordinates as

dΩ=1cosθdαdβ=cos3θdξdη,

where we defined tangent and sine coordinates as follows:

ξ=ftanθcosφη=ftanθsinφ
α=fsinθcosφβ=fsinθsinφ.

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

dαdβ=cos4θdξdη.

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

dαdβ=cos2θdξdη.

Then, the relation for the scalar field is given by

dαdβ=f2ξξ+f2dξ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 ofT(ξ):

t(p:f)=exp[ik(ξξ+f2f)]exp(2iπpξ)dξ2=exp(ikf)kiξξ+f2exp[ikξξ+f2]exp(2iπpξ)dξ2.

Then, this formula can be converted as follows:

t(p:f)=exp(ikf)kds2dξ21λ2ss×exp[2πifλ2ss]exp[2iπ(p+s)ξ],

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

exp(ikξξ+f2)ξξ+f2=iλ2ssexp[2πi(sξ+fλ2ss)]ds2.

By calculating a derivative at k, we obtain

t(p:f)=exp[2πif(λ2pp1λ)]λ2pp[12πλλ2pp+ifλ].

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

a(u)=exp[ik(ξξ+f2f)]A(ξ)exp(2πiuξ)dξ2=g(p)t(up:f)dp2,

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

g(p)=A(ξ)exp(2iπpξ)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).

OCIS Codes
(220.2560) Optical design and fabrication : Propagating methods
(260.1960) Physical optics : Diffraction theory
(070.7345) Fourier optics and signal processing : Wave propagation

ToC Category:
Physical Optics

History
Original Manuscript: May 27, 2011
Revised Manuscript: June 24, 2011
Manuscript Accepted: June 25, 2011
Published: July 11, 2011

Citation
Satoru Odate, Chiaki Koike, Hidemitsu Toba, Tetsuya Koike, Ayako Sugaya, Katsumi Sugisaki, Katsura Otaki, and Kiyoshi Uchikawa, "Angular spectrum calculations for arbitrary focal length with a scaled convolution," Opt. Express 19, 14268-14276 (2011)
http://www.opticsinfobase.org/oe/abstract.cfm?URI=oe-19-15-14268


Sort:  Author  |  Year  |  Journal  |  Reset  

References

Cited By

Alert me when this paper is cited

OSA is able to provide readers links to articles that cite this paper by participating in CrossRef's Cited-By Linking service. CrossRef includes content from more than 3000 publishers and societies. In addition to listing OSA journal articles that cite this paper, citing articles from other participating publishers will also be listed.

Figures

Fig. 1 Fig. 2 Fig. 3
 
Fig. 4
 

« Previous Article  |  Next Article »

OSA is a member of CrossRef.

CrossCheck Deposited