1. Introduction

Diffraction calculation methods [1–7] are essential to model the free-space wave propagation, which is vital to reconstruct the object image in digital holography [8,9] and to synthesize the holograms in computer holography [10,11]. There are several different categories of these calculation methods, such as diffraction between parallel planes [4,12–17] and between non-parallel planes [18–20]. In the first category, there are still two sub-categories: on-axis diffraction [12–14] and off-axis diffraction [4,15–17]. The so-called off-axis diffraction refers to the situation where the origin of the observation window at the output plane is shifted with that of the source window at the input plane. This is very useful for simulating non-paraxial fields, which travel in an off-axis direction, and large-scale fields, such as object fields of an extremely high definition computer-generated hologram [21].

Among all these diffraction calculation methods, angular spectrum method (ASM), as a rigorous diffraction theory in free space, plays an important role. Recently, for the on-axis diffraction, band-extended angular spectrum method (BEASM) was proposed to achieve accurate calculation of the diffraction field in a wide propagation range [13]. On the other hand, for the off-axis diffraction, there is a lack of such a method. Shifted-Fresnel method (Shift-FR) has the capability of off-axis diffraction calculation and scaling the sampling interval thanks to the scaled fast Fourier transform (FFT) [4]. However, aliasing error would occur with Shift-FR in near-field propagation [15]. Although aliasing-reduced Shift-FR (AR-Shift-FR) was proposed to alleviate this problem [15], the calculation accuracy is still limited due to the paraxial approximation of Fresnel diffraction and the band limiting used to reduce the noise, as discussed in Section 4.

The shifted band-limited angular spectrum method (Shift-BLASM) has been proposed to avoid the problem of the off-axis near-field propagation [16]. It successfully achieved accurate calculation of near-field diffraction. However, the Shift-BLASM suffers from low accuracy in far field because the effective bandwidth and sampling number in the spatial frequency domain shrink rapidly with the increase of the propagation distance. To overcome the problem, in this work, a shifted-BEASM is proposed based on the band extended technique [13]. In cases of off-axis propagation, asymmetry caused by the offset of the output sampling window produces severe sampling problems. We analyze the asymmetry in this study, which is most likely important in various situations of diffraction calculation. An off-axis parameter is added to the transfer function of ASM, and the band-extended method is also applied therein. Through rearranging the sampling points in the spatial frequency domain, extension of the effective bandwidth is achieved; all the sampling points are used for the calculation. Therefore, the shifted-BEASM is suitable for both near- and far-field off-axis diffraction calculations with a high accuracy. Theoretical analysis and numerical examples are presented to verify the validity.

2. Shortage of sampling points in shift-BLASM

As shown in Fig. 1, there is an off-axis offset between the origins of the source window and the observation window. The y-axis is omitted for simplification of the following discussion. The diffraction field $u{^{\prime}_o}(x )$ can be calculated from the input field ${u_i}(x )$ as

 

Fig. 1. Schematic diagram of off-axis diffraction

Download Full Size | PDF

In Shift-BLSAM, the problem is solved by shifting the origin of the output field to the center of the observation window, indicated by the yellow dot in Fig. 1, in order to offset the field shift: $x = \hat{x} + {x_0}$. In this way, using Eq. (1), the output field can be rewritten as [16]

It is evident that the offset ${x_0}$ is independent of the propagation distance z. However, if we directly use Eqs. (3) and (4) for the diffraction calculation, aliasing errors would occur due to the violation of the sampling theorem. Supposing the input field ${u_i}(x )$ is uniformly sampled inside the source window and the sampling number and interval are N and ${\Delta _x}$, respectively, the size of the sampling window is properly $N{\Delta _x}$. According to the sampling theorem, to sample ${H_n}({{f_x}} )$ correctly, the following must be satisfied,

where

It should be noted that, if $f_{\textrm{limit}}^{(+ )} > 1/2{\Delta _x}$ or $f_{\textrm{limit}}^{(- )} < - 1/2{\Delta _x}$, they should be replaced by $f_{\textrm{limit}}^{({\pm} )} ={\pm} 1/2{\Delta _x}$ because ${f_x}\epsilon [{ - 1/2{\Delta_x},1/2{\Delta_x}} )$ according to the sampling parameters.

It is clear that effective bandwidth ${B_{\textrm{limit}}}$ would shrink with increasing z and decreasing (${x_0} \pm N{\Delta _x}$) values. Because the sampling interval $\Delta f$ is fixed, the effective sampling number of ${H_n}({{f_x}} )$, which is ${N_{\textrm{limit}}} = {B_{\textrm{limit}}}/\Delta f$, also decrease in this case. Here, the effective bandwidth is represented as ${B_{\textrm{limit}}} = {f_{\textrm{ul}}} - {f_{\textrm{ll}}}$, where ${f_{\textrm{ul}}}$ and ${f_{\textrm{ll}}}$ are upper and lower limits of the spatial frequency, respectively, given by relations (6). As a result, increase of the propagation distance always reduces the calculation accuracy in the Shift-BLASM. To show the reduction of ${B_{\textrm{limit}}}$ and ${N_{\textrm{limit}}}$ intuitively, actual values for specific parameters are shown in Fig. 2. The parameters used to draw the curves are listed in Table 1.

 

Fig. 2. (a) Schematic diagram of shrink of the effective bandwidth with increasing the propagation distance in Shift-BLASM. Actual change of (b) ${B_{\textrm{limit}}}$ and (c) ${N_{\textrm{limit}}}$ is depicted as functions of the propagation distance with parameters in Table 1.

Download Full Size | PDF

Table 1. Parameters used for examples in Fig. 1.

View Table

As shown in Fig. 2, the effective bandwidth and sampling number shrink rapidly with increasing the propagation distance. The components outside the effective band are forced to be zero, as in (a); the corresponding sampling points at those positions are wasted in practice.

3. Shifted band-extended angular spectrum method (Shift-BEASM)

To make all of the input sampling points contribute to the output field, we uniformly rearrange the sampling points within an extended band, as shown in Fig. 3(a). In this way, the sampling interval of ${H_n}({{f_x}} )$ is always reduced because the whole bandwidth is no more than $1/{\Delta _x}$ and all of the 2N sampling points are used for the calculation.

 

Fig. 3. (a) Schematic diagram of the effective bandwidth shrinking with the propagation distance in the Shift-BEASM. (b) The effective bandwidth ${B_{\textrm{extend}}}$ and (c) ${N_{\textrm{extend}}}$, calculated from the parameters in Table 1.

Download Full Size | PDF

To formulate the Shift-BEASM, we introduce a reduction parameter R, where R $> $ 1, to the sampling interval of ${H_n}({{f_x}} )$ as follows:

By replacing $\Delta f$ with $\Delta {f_{\textrm{new}}}$ in (5), relations (6) are rewritten as:

where

It is important to choose an appropriate value of the reduction parameter R in order to use as much effective bandwidth of the input field as possible in the proposed Shift-BEASM. However, it is not easy to determine the reduction parameter because the effective bandwidth of the transfer function also depends on the reduction parameter. When the sampling interval is reduced to $\Delta {f_{\textrm{new}}}$, the transfer function ${H_n}({{f_x}} )$ is simply calculated within a bandwidth:

Substituting Eq. (10) in Eq. (12), we have

Because it is difficult to solve Eq. (13) for R directly, let us assume that the propagation distance is much larger than the maximum of $|{{x_0} + RN{\Delta_x}} |$ and $|{{x_0} - RN{\Delta_x}} |$. In this case, Eq. (13) is simplified down to

Because R>1, the propagation distance z should be larger than $2N\Delta _x^2/\lambda $. Regarding other two cases in Eq. (9), we also get the same R value as an optimum reduction parameter.

Here note that the assumption of the long propagation distance is reasonable in the Shift-BEASM. The Shift-BLASM suffers from narrow bandwidth and shortage of the sampling points of the transfer function in far fields. This is the problem that we want to solve in the Shift-BEASM. When R=1, Eq. (15) gives $z = 2N\Delta _x^2/\lambda $. Thus, when $z < 2N\Delta _x^2/\lambda $, Shift-BEASM degenerates into the Shift-BLASM in practice. It should be noted that, if $f_{\textrm{extend}}^{(+ )} > 1/2{\Delta _x}$ or $f_{\textrm{extend}}^{(- )} < - 1/2{\Delta _x}$, they should be also replaced by $f_{\textrm{extend}}^{({\pm} )} ={\pm} 1/2{\Delta _x}$. In this case, $\Delta {f_{\textrm{new}}}$ should be renewed again because the whole bandwidth ${B_{\textrm{extend}}}$ changes. Note that 2N points are used for the analysis in this paper because we intend to make full use of the input sampling resources. More sampling points may give even better results properly.

The effective bandwidth ${B_{\textrm{extend}}}$ in this technique is shown in Fig. 3(b). The parameters are the same as those in Fig. 2. When comparing with Fig. 2(b), the effective bandwidth increases over an order of magnitude. In addition, let us emphasize that the effective sampling number ${N_{\textrm{extend}}} = {B_{\textrm{extend}}}/\Delta {f_{\textrm{new}}}$ is a constant of 2N, i.e., 2048 for the whole propagation distance in this example, as shown in Fig. 3(c). This is almost 100 times larger than that of the Shift-BLAMS in the far-field region indicated by a rectangle in Fig. 2(c). Conclusively, the proposed Shift-BEASM has a much wider bandwidth and more sampling points than the Shift-BLASM; these properties most likely benefit the diffraction calculation.

The procedures of the Shift-BEASM is summarized as follows:

The above steps 5 and 6 cannot be realized by conventional FFT-based methods. This is because the input field ${u_i}(x )$ is sampled with the interval of ${\Delta _x}$, while the spectrum $U({{f_x}} )$ must be sampled with the interval of $\Delta {f_{\textrm{new}}}$. In addition, the sampling positions of $U({{f_x}} )$ are no longer symmetrical with respect to the coordinate origin. Because the ordinary FFT cannot achieve such a transform, a shifted and scaled FFT is needed for the Fourier transform in this technique. There are two fast methods that can be used for the required Fourier transform: chirp-z transform [22] and nonuniform FFT (NUFFT) [23].

Chirp-z transform uses convolution to realize the shifted and scaled FFT and thus requires three ordinary FFTs. Alternatively, NUFFT is realized by the combination of interpolation and FFT, and has flexibility to give arbitrary sampling positions. Suppose the number of sampling points is M, the computational complexity of NUFFT is O(MlogM) and therefore is comparable with ordinary FFT. When taking the computational efficiency into account, NUFFT is suited for performing the Fourier transforms in steps 5 and 6. Please note that with NUFFT zero padding to the input field in step 5 is not necessary, alternatively, we can directly transform the N-point sampled input field to the 2N-point sampled Fourier spectrum by adaptive-sampling ASM [24].

4. Numerical verification of the shift-BEASM

The setup for the numerical experiment is shown in Fig. 4(a). Off-axis diffractions by a one-dimensional (1-D) rectangular aperture and two-dimensional (2-D) circular aperture are demonstrated for numerical verification. The aperture is illuminated by an oblique plane wave ($\lambda = 500\; \textrm{nm}$) at the incident angle of $\theta = 5^\circ $. The off-axis offset ${x_0}$ is set as $z \cdot \textrm{tan}\theta $ along the incident direction.

 

Fig. 4. (a) Setup for the numerical verification of the proposed Shift-BEASM: Amplitude and phase of the input fields in (b) 1-D and (c) 2-D diffraction.

Download Full Size | PDF

4.1 One-dimensional diffraction calculation

The complex amplitude of the input field, sampled with N = 1024 and ${\Delta _x} = 1\; \mathrm{\mu}\textrm{m}$, is shown in Fig. 4(b). The width of the aperture is 0.768 mm (=$3N{\mathrm{\Delta }_x}/4)$.

In this demonstration, we use the type 3 of Greengard and Lee’s NUFFT to realize the proposed Shift-BEASM [23,25]. Many methods of NUFFT have been proposed to perform flexible sampling in both space and spatial frequency domains; some of them have been used for diffraction calculation to realize different functions [13,26]. The NUFFT used here is defined as

Besides the proposed Shift-BEASM, the diffraction field was also calculated with two other methods: Shift-BLASM [16] and AR-Shift-FR [15] for comparison. Figure 5 shows amplitude of diffraction fields calculated by these three methods with three propagation distances: z$= 2$ mm, $100$ mm and $1000$ mm. To evaluate accuracy of these methods, diffraction fields calculated by the Rayleigh-Sommerfeld (R-S) integral are also shown as the ground truth [27].

 

Fig. 5. Amplitude of diffraction fields calculated by the AR-Shift-FR, Shift-BLASM, Shift-BEASM (this work) and R-S integral.

Download Full Size | PDF

We find out the following from Fig. 5: (1) the AR-Shift-FR does not work well in the case of large offsets because the paraxial condition is violated in this case, (2) the Shift-BLASM fails in far-field diffraction calculation, and (3) the proposed Shift-BEASM gives very similar results to those by the R-S integral in the far-field diffraction as well as the near-field diffraction.

To quantitatively evaluate the calculation accuracy, the following SNR metric is used [18]:

Here, superscript * denotes complex conjugation and ${s_{\textrm{obs}}}$ is an integral region that agrees with the observation window. The SNR evaluates similarity of complex amplitude distribution of wave field $u(x )$ to reference field $\; {u_{\textrm{ref}}}(x )$ calculated by the R-S integral. Thus, the higher the SNR, the better the accuracy.

Figure 6(a) shows the curves of SNR as a function of the propagation distance in the range of 2 mm to 1000 mm. The parameters are the same as those in Fig. 4. The proposed Shift-BEASM obviously gives better results than the Shift-BLASM at a distance approximately more than 100 mm. Shift-BLASM works well only at a limited small distance. For the AR-Shift-FR, the SNR is always lower than that of the other two techniques and decreases with increasing the propagation distance. This is most likely because the off-axis offset increases as the diffraction distance increases; the paraxial approximation may be unsuitable more and more with increasing the distance.

 

Fig. 6. SNRs varying with the propagation distance at (a) $\theta = 5^\circ $, (b) $\theta = 1^\circ $, and (c) $\theta = 10^\circ $.

Download Full Size | PDF

To verify the reason of low accuracy of the AR-Shift-FR, SNR curves at a small incident angle of $\theta = 1^\circ $ are calculated as shown in Fig. 6(b). The curves of the Shift-BEASM and Shift-BLASM are almost unchanged, while accuracy of the AR-Shift-FR is improved considerably. On the contrary, when increasing the incident angle, as shown in Fig. 6(c) where $\theta = 10^\circ $, SNR of the AR-Shift-FR reduces drastically. In addition, the SNR monotonically decreases with increasing the propagation distance in this case. Because the offset is given by ${x_0} = z \cdot \textrm{tan}\theta $, the fact suggests that a large offset violates the assumption of the Fresnel approximation, i.e., the paraxiality in the AR-Shift-FR.

The results of Fig. 6 show that the accuracy of the proposed Shift-BEASM is, unlike the AR-Shift-FR, independent of the offset and illumination angle. The technique always gives the highest accuracy of the three techniques. As a result, we conclude that the proposed Shift-BEASM is only one technique to make it possible to calculate shifted diffraction fields in a wide range from the near to far field. The offset does not affect the calculation efficiency because it is just a parameter of the transfer function. Since the numbers of sampling points of AR-Shift-FR, Shift-BLASM and Shift-BEASM are the same, the computational complexity is almost the same. Therefore, the calculation time of all examples are comparable and are all at the order of milliseconds (6 ms, 9 ms, 13 ms) with a laptop (CPU 2.9 GHz, RAM 8G).

4.2 Two-dimensional diffraction calculation

The 2-D Shift-BEASM is a little more complicated, but can be easily formulized by a similar procedure in Section 3. When the propagation distance is relatively large, x and y coordinates can be treated independently [16]. In the 2-D case, the Shift-BEASM is formulated as follows.

For the 2-D case, the transfer function in Eq. (4) is written as

Here, note that $f_{x\_\textrm{extend}}^{({\pm} )}$ is a function of ${f_y}$, and $f_{y\_\textrm{extend}}^{({\pm} )}$ is a function of ${f_x}$. Therefore, it is very difficult to determine the sampling range and interval in the $({{f_x},{f_y}} )$ domain.

As mentioned in section 3, to determine the sampling parameters, it is important to determine the values of ${R_x}$ and ${R_y}$. However, as shown in Eq. (24), there are three cases for each dimension. Thus, there are totally nine cases for a combination of (${x_0}$, ${y_0}$) as shown in Fig. 7. Although some cases, e.g. case (i) and (ix), have a symmetry, we have to treat at least five cases to provide the general solution. Instead of this complicated work, let us introduce one more simple assumption to solve the problem without dividing the case.

 

Fig. 7. Cases of the combination of (${x_0}$, ${y_0}$); and the symmetry: (i) and (ix), (ii) and (viii), (iii) and (vii), (iv) and (vi) have a symmetry, respectively.

Download Full Size | PDF

Here, we analyze the problem in case (ix) for example. In this case, the bandwidths are represented as follows:

To solve Eq. (26) for ${R_x}$ and ${R_y}$, we again use the long propagation-distance condition used in the derivation of the 1-D solution in section 3. In addition, let us introduce another assumption that is practical in many applications; we assume that ${\Delta _x} = {\Delta _y} = \Delta \ge 2\lambda $. Because spatial frequencies ${f_x}$ and ${f_y}$ are given within intervals [$- \frac{1}{{2{\Delta _x}}},\; \frac{1}{{2{\Delta _x}}}$] and [$- \frac{1}{{2{\Delta _y}}},\; \frac{1}{{2{\Delta _y}}}$], respectively, in this case, we can approximate the numerators in Eq. (26) as $\sqrt {1 - {{({\lambda {f_x}} )}^2}} \approx 1$ and $\sqrt {1 - {{({\lambda {f_y}} )}^2}} \approx 1$. Therefore, Eq. (26) is simplified to

As a result, we can get a simple solution:

This agrees with the reduction parameter in the 1-D case. The same result is also derived in all other cases under the second assumption. Using the values of ${R_x}$ and ${R_y}$, the sampling intervals are calculated from Eq. (23) and the ranges are calculated from Eqs. (24) and (25) by using the same approximation; $\sqrt {1 - {{({\lambda {f_i}} )}^2}} \approx 1\; ({i\; \textrm{denotes}\; x\; \textrm{or}\; y} )$. The following calculation process is similar to that in the 1-D case demonstrated in section 3.

We demonstrate 2-D diffraction by the Shift-BEASM using a circular aperture. Geometry and setup for the 2-D case are shown in Fig. 4(a) and (c). Here note that ${x_0}$= z tan $\theta $ and ${y_0} = 0$. The sampling number inside the source window is ${N_x} \times {N_y}$=1024 ${\times} 1024$, and the sampling intervals are ${\Delta _x} = {\Delta _y}$=$\; 1\; \mathrm{\mu}\textrm{m}$. Therefore, the size of the source window is ${S_x} \times {S_y}$=1.024 mm ${\times} 1.024\; \textrm{mm}$. The diffraction fields calculated at three propagation distances using the proposed Shift-BEASM are shown in Fig. 8.

 

Fig. 8. Diffraction fields calculated for the circular aperture using the proposed Shift-BEASM.

Download Full Size | PDF

5. Discussions and conclusions

We successfully confirmed that the proposed Shift-BEASM makes it possible to calculate off-axis diffraction with a high accuracy in both near and far fields. However, there is one important question worthwhile to be discussed. The bandwidth required to calculate exact shifted diffraction is analyzed on the basis of geometry [16]. In the Shift-BLASM, the limited bandwidth exactly agrees with the required bandwidth. Why does the Shift-BEASM improve the accuracy? This can be explained as follow. When the propagation distance is relatively large, the required bandwidth is much narrow, and thus, sampling points within the bandwidth become very few under the condition of a fixed sampling interval. In such cases, sharp truncation of the spatial spectrum and the small sampling number would be harmful for the calculation accuracy [2]. Therefore, we avoid these two problems by band extension and improve the accuracy, as verified by the numerical demonstrations.

In addition to high accuracy over a wide range of the propagation distance, the proposed Shift-BEASM has a considerable advantage that the observation window can be flexibly scaled. This property is derived from the used NUFFT, as demonstrated in [24]. Scaling of the observation window is provided by the AR-Shift-Fresnel but not achieved by the Shift-BLASM.

Furthermore, the proposed Shift-BEASM focuses on the off-axis diffraction calculation in a wide propagation range. Previously, several methods were developed for different diffraction calculation tasks, such as Rayleigh-Sommerfeld convolution (RSC) [7] and BEASM [13], those methods focus on the on-axis case. Since off-axis diffraction is more generalized and complex, the proposed Shift-BEASM goes one step further in the numerical diffraction calculation. Readers may wonder that why off-axis case is complex than on-axis case. At first glance, only a linear phase factor [Eq. (4)] is added in the off-axis case which just shifts the diffraction pattern and the shift value can be determined by the linear phase. It seems like that we can calculate the corresponding on-axis pattern and then apply a shift for the off-axis pattern. It is true that there is a shift. However, besides the shift, there is a “stretch effect” of the diffraction pattern along the shifted direction [17]. The stretch effect is obvious with large linear phase and it exists as long as the linear phase is not zero. Unfortunately, we cannot quantitatively know how much the stretch is. Therefore, we cannot obtain the accurate off-axis diffraction pattern from the corresponding on-axis calculation result. To show this effect, a numerical demonstration is carried out. A 1-D rectangular aperture is used as the input object and the propagation distance z is set to be 200 mm and the illumination wavelength is 500 nm. First, we calculate the on-axis diffraction pattern by the proposed Shift-BEASM with offset ${x_0} = 0$ and $\theta = 0^\circ $; and the result is shown in Fig. 9(a). Then, we calculate the off-axis diffraction pattern by RSI with offset ${x_0} = z \cdot \textrm{tan}\theta $ and $\theta = 10^\circ $; and the result is shown in Fig. 9(b). There is obvious difference between the on-axis pattern and the off-axis pattern. The SNR (28.1 dB) is very low which means the accuracy is low. Please note we only evaluate the amplitude field to calculate the SNR. On the other hand, we calculate the off-axis diffraction pattern by the proposed Shift-BEASM and the result is shown in Fig. 9(c). The SNR (47.7 dB) is much higher which means the accuracy is high. Because we do not know how to stretch the on-axis pattern for the off-axis pattern, even for the simple apertures, the off-axis diffraction calculation is necessary which can directly calculate the off-axis pattern. For general input objects, it is impossible to establish the stretching mapping between the on-axis pattern and the off-axis pattern; and the direct off-axis diffraction calculation would be the only and optimal option. Furthermore, the proposed Shift-BEASM for off-axis diffraction calculation is more likely a framework. Although in the demonstration we use an oblique plane wave as the illumination and the off-axis offset ${x_0}$ is set as $z \cdot \textrm{tan}\theta $, it is just because this demonstration case is intuitive and easy to understand and follow. Actually, the illumination can be in any form as long as it can be sampled correctly, such as spherical wave. In such cases, the off-axis diffraction pattern can only be obtained by the off-axis diffraction calculation.

 

Fig. 9. Amplitude field given by (a) on-axis Shift-BEASM calculation, (b) off-axis RSI calculation and (c) off-axis Shift-BEASM calculation.

Download Full Size | PDF

In conclusion, the shifted band-extended angular spectrum method was proposed for accurate off-axis diffraction calculation in a large propagation distance range. Extension of the effective bandwidth of the transfer function and full utilization of the sampling points are contributed to the high accuracy over the wide range. This method would be useful for applications that require accurate modeling of off-axis wave propagation, such as large-scaled computer holography, digital holography, and design of diffractive optical elements.