Modified shifted angular spectrum method for numerical propagation at reduced spatial sampling rates

The shifted angular spectrum method allows a reduction of the number of samples required for numerical off-axis propagation of scalar wave fields. In this work, a modification of the shifted angular spectrum method is presented. It allows a further reduction of the spatial sampling rate for certain wave fields. We calculate the benefit of this method for spherical waves. Additionally, a working implementation is presented showing the example of a spherical wave propagating through a circular aperture. © 2014 Optical Society of America OCIS codes: (050.1940) Diffraction; (070.0070) Fourier optics and signal processing; (070.7345) Wave propagation. References and links 1. D. M. Paganin, Coherent X-Ray Optics (Oxford University, 2006). 2. 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). 3. K. Matsushima, “Shifted angular spectrum method for off-axis numerical propagation,” Opt. Express 18(17), 18453–18463 (2010). 4. C. Guo, Y. Xie, and B. Sha, “Diffraction algorithm suitable for both near and far field with shifted destination window and oblique illumination,” Opt. Lett. 39(8), 2338–2341 (2014). 5. X. Yu, T. xiong, P. Hao, and W. Wei, “Wide-window angular spectrum method for diffraction propagation in far and near field,” Opt. Lett. 37(23), 4943–4945 (2012). 6. T. Kozacki, K. Falaggis, and M. Kujawinska, “Computation of diffracted fields for the case of high numerical aperture using the angular spectrum method,” Appl. Opt. 51(29), 7080–7088 (2012). 7. L. Onural, “Sampling of the diffraction field,” Appl. Opt. 39(32), 5929–5935 (2000). 8. L. Onural, “Exact analysis of the effects of sampling of the scalar diffraction field,” J. Opt. Soc. Am. A 24(4), 359–367 (2007). 9. M. Hillenbrand, D. P. Kelly, and S. Sinzinger, “Numerical solution of nonparaxial scalar diffraction integrals for focused fields,” J. Opt. Soc. Am. A 31(8), 1832–1841 (2014). 10. R. P. Muffoletto, J. M. Tyler, and J. E. Tohline, “Shifted Fresnel diffraction for computational holography,” Opt. Express 15(9), 5631–5640 (2007). 11. P. Lobaz, “Memory-efficient reference calculation of light propagation using the convolution method,” Opt. Express 21(3), 2795–2806 (2013). 12. K. Falaggis, T. Kozacki, and M. Kujawinska, “Computation of highly off-axis diffracted fields using the bandlimited angular spectrum method with suppressed Gibbs related artifacts,” Appl. Opt. 52(14), 3288–3297 (2013). 13. S. Odate, C. Koike, H. Toba, T. Koike, A. Sugaya, K. Sugisaki, K. Otaki, and K. Uchikawa, “Angular spectrum calculations for arbitrary focal length with a scaled convolution,” Opt. Express 19(15), 14268–14276 (2011). 14. P. Lobaz, “Reference calculation of light propagation between parallel planes of different sizes and sampling rates,” Opt. Express 19(1), 32–39 (2011). 15. K. Yamamoto, Y. Ichihashi, T. Senoh, R. Oi, and T. Kurita, “Calculating the Fresnel diffraction of light from a shifted and tilted plane,” Opt. Express 20(12), 12949–12958 (2012). 16. T. Shimobaba, T. Kakue, M. Oikawa, N. Okada, Y. Endo, R. Hirayama, and T. Ito, “Nonuniform sampled scalar diffraction calculation using nonuniform fast Fourier transform,” Opt. Lett. 38(23), 5130–5133 (2013). #217174 $15.00 USD Received 16 Jul 2014; revised 8 Sep 2014; accepted 17 Sep 2014; published 17 Oct 2014 (C) 2014 OSA 20 October 2014 | Vol. 22, No. 21 | DOI:10.1364/OE.22.026265 | OPTICS EXPRESS 26265 17. 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). 18. A. Ritter, P. Bartl, F. Bayer, K. C. Gödel, W. Haas, T. Michel, G. Pelzer, J. Rieger, T. Weber, A. Zang, and G. Anton, “A simulation framework for coherent and incoherent X-ray imaging and its application in Talbot-Lau dark-field imaging,” Opt. Express 22(19), 23276–23289 (2014).


Introduction
The propagation of scalar wave fields from one plane to a second parallel plane through free space is an important problem in the field of wave optics [1].The angular spectrum method or the Rayleigh-Sommerfeld diffraction integrals of first and second kind provide analytical exact solutions to this problem.The angular spectrum method and the Rayleigh-Sommerfeld diffraction integral of the first kind, which are equivalent, solve this problem for given wave field amplitudes on the input plane.The Rayleigh-Sommerfeld diffraction integral of the second kind solves this problem for the derivative of the wave field amplitude in propagation direction.The Rayleigh-Sommerfeld diffraction integral can be evaluated directly in the spatial domain while the angular spectrum method requires a Fourier transform to obtain the angular spectrum of the wave field.There are several methods to obtain numerical solutions for all these methods, which are suited for different conditions and problems.In the following, we will focus on the angular spectrum method, especially with regard to the problem described in the next paragraphs.
Concerning X-ray wave optics, the wavelength can be several orders of magnitude smaller than the wavelength in visbile light optics.Thus, the spatial sampling rates needed for a correct spatial sampling of wave fields for numerical propagation in the X-ray optics regime tend to be very high.This is especially true for the off-axis propagation of such wave fields at large oblique angles, where the needed sampling rates can be even higher by several orders of magnitude compared to the case of propagation parallel to the optical axis.A drawback of the angular spectrum method is that it is only suited for numerical near field propagation up to a certain limit distance.In the case of X-ray optics this limitation might be a minor problem compared to the problem of high sampling rates due to small wavelenghts.The angular spectrum method and explicitly the band limited angular spectrum method [2] are well suited for numerical propagation in X-ray optics or generally in wave optics with similar conditions.
Another drawback of a plain numerical implementation of the angular spectrum method is that the input sampling window and the output sampling window have the same size and center.For off axis propagation the size of this sampling window has to be increased to be able to contain the field both in the input as well as in the output plain.In [3], the shifted angular spectrum method is presented as a solution to this problem.It allows for off-axis numerical propagation of scalar wave-fields by shifting the output sampling window relative to the input sampling window by using the Fourier shift theorem.With this method, the number of samples needed in the numerical propagation can be reduced as the input and output sampling windows do not neccesarily have the same center.Still, the sampling condtions concerning the spatial sampling rate are the same as in the original band limited angular spectrum method.
To further reduce the spatial sampling rate that is needed, especially for off axis propagation, a modification of the shifted angular spectrum method ist presented.As a prerequisite, the angular frequencies contained in the wave field have to be within limits and the limits are not centered around zero in the domain of angular frequencies.The method is then able to shift the spectrum of angular frequencies prior to a discrete Fourier transform.The sampling can then be done with a lesser sampling rate compared to the original signal without producing aliasing errors.The angular frequency shift can be considered in the subsequent operations of the method, especially in the transfer function.
In section 2, the modification of the shifted angular spectrum is introduced in detail.Section 3 discusses the angular sampling conditions of the modified transfer function of the method and compares these to other methods found in literature, also with regard to the maximum allowed propagation distance.Subsequently, the spatial sampling of the wave field and conditions, which allow to reduce the spatial sampling rate in general and more specific for spherical waves, are discussed in section 4. The shown reduction in spatial sampling rate can be of several orders of magnitude.Finally, the numerical propagation of a spherical wave through a circular aperture is presented in section 5. There, it can be shown that with the original shifted angular spectrum method it is not possible to reduce the sampling rate below a certain rate while it is possible to further reduce the sampling rate with the modification presented in section 2.
There are several methods which also try to reduce the sampling demands for the angular spectrum method.In [4], a similar approach applicable to oblique illumination with plane waves is used.Compared to that, the method presented in section 2 is directly applicable for more general wave fields under the conditions given in section 4, and thus for example illumination with spherical waves, which is also shown in section 5.In [5,6], methods are presented which reduce the amount of zero padded samples while spatial sampling rate plays no role in these approaches.A combination of these methods with the method presented in section 2 could be possible.In [7][8][9], different conditions and methods are analyzed that might assume approximations or could result in the presence of aliasing errors that can be accepted which also might reduce the sampling demands.Further methods for numerical propagation, with shifted or tilted destination window or which try to reduce the computational demands, can be found in [10][11][12][13][14][15][16][17].

Modification of the shifted angular spectrum method
The method proposed in this work is based on the shifted angular spectrum method [3].The shifted angular spectrum method is given by with ϕ in (x, y) being the wave field on the input plane and ϕ out (x, y) being the wave-field on the output plane.The origin of the output plane is translated by the vector (x 0 , y 0 , z) against the origin of the input plane.x and y denote the cartesian coordinates on each plane whereas k x and k y are the corresponding angular frequencies.The transfer function for the shifted angular spectrum method is given by with k = 2π λ being the angular wavenumber and λ being the wavelength of the wave field.The Fourier transform φ(k x , k y ) of the wave field ϕ(x, y) is given by and its inverse We assume that the wave field can be written as a product of a wave field ϕ * (x, y) and a plane wave ϕ(x, y) = ϕ * (x, y) exp(iq x x + iq y y). ( This means, that the angular spectrum of the wave field ϕ(x, y) is given by the angular spectrum of the wave field φ * (k x , k y ) shifted in the domain of angular frequencies by the vector (q x , q y ).
If we take the Fourier transform of Eq. ( 1) and take Eqs. ( 3) and ( 6) into account, we arrive at φ * out (k We then first introduce k x with and accordingly k y .We substitute this into Eq.( 7) and rename k x by k x and k y by k y , to arrive at We apply the inverse Fourier transform of Eq. ( 4) to this equation and insert the inverse of Eq. ( 5) Finally, we gain an equation for the modified shifted angular spectrum method: For q x = 0 and q y = 0 we obtain the shifted angular spectrum method as a special case of the modified shifted angular method.The resulting Eq. 12 of this calculation can be interpreted in the following way.First, the angular spectrum is obtained by a Fourier transform of the product of input wave field and a plane wave.Due to the Fourier shift theorem, the angular spectrum of this product is shifted by the angular frequencies q x and q y compared to the angular spectrum of the input wave field alone.The angular spectrum of this product is then multiplied with a modified transfer function of the shifted angular spectrum method.The modified transfer function takes into account, that the angular frequencies are shifted by q x and q y .The resulting angular spectrum, which is still shifted, will then be transformed by the inverse Fourier transform.The result of this transform is the product of the output wave field and the plane wave used in the product with the input wave field.Therefore, this result has to be multiplied by the inverse of the plane wave that was used in the first place to undo the shift in angular frequencies and to obtain the output wave field.

Sampling limits and band limiting
In the following we will first discuss the x component, the y component can be discussed in analogy.We assume that the width of the sampling window in the spatial domain is given by at least twice the value of L x .The angular sampling interval ∆k x for the domain of angular frequencies k x is then given by and thus the angular sampling rate r k x in the domain of angular frequencies is given by Following the discussion in [2,3], the local angular signal frequency ω x (k x ) of the transfer function is given by and therefore and analogous To avoid aliasing the Nyquist conditions for both dimensions have to be satisfied.In a numerical implementation this can be used to band-limit the transfer function by taking the value of p k,x 0 ,y 0 (k x + q x , k y + q y ) if both conditions are true otherwise a zero value is used.

Discussion on the sampling limit
By using the conditions in Eq. ( 18) the transfer function of the modified shifted angular spectrum method can be band limited.To make sure that the numerical propagation is done in a physical way, the angular frequencies contained in the wave field have to be within the limits of the band.In our formalism, regarding Eqs. ( 15) and ( 16), the angular frequencies of the wave field are given by k x + q x and k y + q y and are equivalent to k x and k y in the unmodified angular spectrum method.Therefore, the discussion on the sampling limit given in [3] is also true for the case of the modified angular spectrum method.This has consequences on the maximum propagation distance that can effectively be done with this method.For a given bandwidth of the input wave field a given size of the sampling window and a given offset between the sampling windows, the propagation can only be done up to a certain distance limit without loss of accuracy.Detailed discussions on this limit for the angular spectrum and shifted angular spectrum method, which is also true for this modified method, can be found for example in [4][5][6]17].To enable propagation for distances above the distance limit one has to choose another method for example [4][5][6].Or, one has to increase the size of the sampling window, which according to Eq. 18, will increase the acceptable bandwidth of the input wave field, but at the cost of an increased number of samples.This has also been discussed in [2,5,6,17].Especially, for the case of X-ray optics, where wavelengths are very small compared to the sampling window size, a sufficient sampling rate can be harder to obtain than a sufficient sampling window size.

4.
Reducing the spatial sampling rate with the modified shifted angular spectrum method

General discussion
At first, the modified method seems to be more complex than the original shifted angular spectrum method and the result is the same.But, for certain wave fields ϕ in , that have an angular spectrum, which is limited between boundaries k min x , k min y and k max x , k max y , this might open the opportunity to reduce the spatial sampling rate r x and r y in discrete numerical computations of the angular spectrum method.To gain a correct sampling of the wave field the sampling rate, in terms of an angular frequency, has to be twice as high as the highest angular frequency of the field.Therefore, the sampling rates in terms of a spatial frequency are given by Concerning the modified shifted angular spectrum method, the spatial sampling of the wave field prior to its Fourier Transform has to be discussed.For these wave fields angular shifts may be found, which effectively set the boundaries for the shifted angular spectrum symmetric to the origin in the angular frequency domain.This will lower the highest angular frequency and therefore the spatial sampling rate needed for the modified shifted angular spectrum method.The reduced rates are then given by

Spherical waves
Assume a spherical wave within a sampling window (L x , L y ) centered around zero starting at (∆x, ∆y, −∆z).The values of the angular shifted wave field in the sampling plane are given by We obtain the local angular frequencies of the spherical wave by Again for q x = 0 and q y = 0 we obtain the special case of the shifted angular spectrum method.
To find values for q x and q y we use Eq. ( 21).For this we need the minimum and maximum angular frequencies within the sampling window for q x = 0 and q y = 0. We will first discuss the function k x (x, y), which is strictly monotonic with regard to x for all y.Regarding the y-axis the function has a minimum at y = ∆y for all x.This effectively creates two half spaces y < ∆y and y > ∆y where the function again is strictly monotonic with regard to y for all x.Therefore, we can find the minimum and maximum values k min x and k max x within the sampling window by evaluating all four corners of the sampling window.Additionally, we have to include the points (x = ± L x 2 , y = ∆y) within the search for minimum and maximum, if |∆y| < L y 2 .In analogy, we find k min y and k max y .Figure 1 shows the quotient between reduced rate and rate for a spherical wave in a quadratic sampling window L x = L y = L for different values of L z and different offsets ∆x z and ∆y z .It is clearly observable that a reduction of the sampling rate with many orders of magnitude can be achieved.This is especially true for smaller L z and bigger ∆x z and ∆y z .Although these reductions might be somehow lessened for more realistic wave fields containing additional contributions, e.g. from a sample.Nonetheless, this could remove high sampling demands for largely off axis propagation in spherical setups.

Numerical propagation of spherical waves through circular apertures
The subfigures in Fig. 2 show the resulting intensities of a numerical off-axis propagation of a spherical wave through a circular aperture.The propagation has been done using our framework cxi for coherent X-ray imaging [18].The photon energy of the wave-field is given at 25 keV, which corresponds to a wavenumber of about 126•10 9 1  m .The origin of the spherical wave is displaced against the center of the aperture by 30 µm in x-direction and -10 mm in z-direction.The field is propagated from the aperture plane to a plane at a distance of 10 mm.Each subfigure shows the intensity of the propagated wave-field.In all subfigures, the complete sampling window used in the propagation is shown and has a width of 120 µm x 120 µm.This also in-  The high spatial sampling is chosen to be just above the sampling rate specified by the Nyquist theorem for the angular spectrum method and the shifted angular spectrum method.The results in Figs.2(a), 2(c) and 2(e) show that all three methods give the same correct results for the high spatial sampling rate.The low spatial sampling rate is chosen to be just above the sampling rate specified by the Nyquist theorem for the modified shifted angular spectrum method.For the low spatial sampling rate the angular spectrum method (see Fig. 2(b)) and the shifted angular spectrum method (see Fig. 2(d)) are not able to reproduce the correct result due to an insufficient spatial sampling rate and thus aliasing errors.The centers of the propagated waves in Figs.2(b) and 2(d) are not at (-30 µm, 0) as expected.Additionally, artifacts can be observed especially at the left hand borders of the circular beams in Figs.2(b) and 2(d).In contrast, the modified shifted angular spectrum method is still able to produce the same correct result as before (see Fig. 2(f)) at a sampling rate which is lower by a factor of 0.4.It should be noted here that the reduction of the sampling rate is not as high as one might expect from the theoretical consideration done in the section before.This is most likely due to higher angular frequencies, introduced by the edge of the circular aperture.
In Fig. 3 we see a similar propagation test, where the origin of the spherical wave is displaced by 100 µm on the x axis perpendicular to the optical axis (z axis).In Fig. 3(a), the shifted angular spectrum method shows good results for a sampling window containing 968×5711 samples.The spatial sampling rate was chosen to be above the sampling rate specified by the Nyquist theorem for the shifted angular spectrum method.In Fig. 3(b) the spatial sampling rate was decreased slightly below the Nyquist sampling rate by a factor of 0.92, yielding 968×5227 samples.This is enough to be able to see artifacts appearing in Fig. 3(b) when using the original shifted angular spectrum method.These artifacts are due to insufficient sampling and thus aliasing errors.In Fig. 3(c) the propagation is done with the modified shifted angular spectrum method.The spatial sampling rate is chosen to be above the sampling rate specified by the Nyquist theorem for the shifted angular spectrum method.The resulting spatial sampling rate is lower by a factor of 0.17 compared to Fig. 3(a) and a factor 0.19 compared to Fig. 3(b) yielding 968×968 samples.The modified shifted angular spectrum method provides a correct result for this case at a spatial sampling rate that is at least five times lower, than the sampling rate that is needed for correct propagation with the original shifted angular spectrum method.
In Fig. 4 the same propagation as in Fig. 3 is done.But, the origin of the spherical wave is displaced even further by 1000 µm perpendicular to the optical axis.Again we see that a reduction of the spatial sampling rate for the original shifted angular spectrum method (see Fig. 4(a)) by a factor of 0.99 produces artifacts in Fig. 4(b).Additionally, the spatial sampling rate for the original shifted angular spectrum method has to be increased by a factor of 8.6 compared to the former propagation with a displacement of 100 µm compared to 1000 µm.While the modified shifted angular spectrum method in Fig. 4(c) works at the same spatial sampling rate as in Fig. 3(c) before, and therefore with a spatial sampling rate that is at least 50-times lower than in the working case of the original angular spectrum method in Fig. 4(a) for a displacement of 1000 µm.

Conclusion
In this work, a modification to the shifted angular spectrum method for numerical off-axis propagation is introduced.The sampling properties in the domain of angular frequencies are the same as in the original shifted angular spectrum method and band limiting can be done in the same way.The method is applicable to near field propagation up to a certain distance limit, above this limit accuracy decreases but due to band limiting aliasing errors can be avoided.Within these limits, the modification allows to reduce the spatial sampling rate up to many orders of magnitude for of axis propagation of spherical waves to the original shifted angular spectrum method.Especially, for the case of X-ray optics or comparable conditions with small wavelengths compared to the sampling window size, the limitation of propagation distance might not be as severe as the requirements on spatial sampling rate.Under this conditions, numerical off-axis propagation of a spherical wave through a circular aperture has been shown.With the original shifted angular spectrum method the spatial sampling rate has to be increased with increasing oblique angle to avoid aliasing errors.In comparison, with the modified angular spectrum method the spatial sampling rate could be reduced up to 50 times for the tested angles.Theoretical, considerations show that reductions up to several orders of magnitude are possible.
In conclusion, the limitations of the modified shifted angular spectrum method with regard to sampling in the domain of angular frequencies and thus the maximum propagation distance are the same as for the original shifted angular spectrum method.Within these limits, the modified shifted angular spectrum method is able to work at reduced spatial sampling rates compared to the original shifted angular spectrum method.For certain cases with wavelengths being very small compared to the sampling window, for example X-ray optics, the reduction of the spatial sampling rate might lessen the sampling demands more than sampling demands have to be increased due to an increased size of the sampling window to allow greater distances.

Fig. 1 .
Fig.1.Ratio of reduced sampling rate r * of the modified angular spectrum method to the sampling rate r for the shifted angular spectrum method.