1. Introduction

Diffraction matters in optics, either as a physical phenomenon or as a tool used in engineering [1]. Particularly, diffraction in free space is an important branch. To study it, the mathematical description was used and many theories were proposed, such as Rayleigh-Sommerfeld theory and angular spectrum theory [1,2]. The goal of such theories is to predict an accurate field distribution after the light propagation. It is well known that the diffraction field is an integral transform of the input field. Thus, theoretical works were carried out for a high-accuracy transform, such as a modified Fresnel propagation [3]. However, we can only get the analytical diffraction field of a few simple input fields. On the other hand, in many cases, the input fields have already been discretely sampled, such as the holograms [4,5] or x-ray diffraction patterns [6,7] captured by digital sensors. Therefore, to get the diffraction field, such as the reconstruction of holograms or x-ray diffraction patterns, discrete computation is necessary in practice. In such cases, the discretization process should be implemented correctly based on the sampling theorem.

From the analytical solution to the discrete computation, three points need attention. First, high accuracy should be maintained. Second, computational efficiency should be high. Third, computational flexibility is desired. The first two points are basic requirements of diffraction calculations. The last one requires the flexibility to observe the diffraction field with adjustable propagation distance, adjustable position and size of the observation window, and adjustable sampling interval. To achieve such goals, many methods have been proposed. Band-limited angular spectrum method (ASM) [8], band-extended ASM [9] were proposed for accurate calculation in a large propagation distance range. Shift Fresnel diffraction [10], two-step Fresnel diffraction [11], and scaled ASM [12] were proposed for an adjustable sampling interval of the diffraction field. Shift band-limited ASM [13] and oblique ASM [14] were proposed for an off-axis diffraction calculation. A generalized Debye integral was proposed for a larger focus range [15]. Semi-analytical propagation techniques [16] and homeomorphic Fourier transform [17] were proposed to improve the computational efficiency. A brief review can be found in the introduction of our recent paper [18].

Here, we would like to emphasize that the numerical diffraction calculation via discrete computation is a research field where both physics and mathematics should be considered. Previous studies only focused on the discrete computation but ignored the physical picture which should be carefully considered. The conventional pure discrete sampling analysis would cause error in the non-convolution-based method. Another representative example is the zero padding operation in the numerical diffraction calculation. Most studies gave the reason of zero padding from the discrete computation perspective: to avoid the circular convolution errors [8,9,19]. Of course, zero padding has this function. However, there are other more important reasons why zero padding is necessary, which are given from the perspective of physics in this paper. Furthermore, with the physical picture, an adaptive zero padding method can be achieved based on needs for a higher efficiency. In short, a unified frequency sampling strategy from a joint physics-mathematics analysis could benefit numerical diffraction calculations a lot.

In this work, from the perspective of frequency response (physics) and discrete computation (mathematics), we propose different sampling strategies for five common-used diffraction calculation methods: single-Fourier transform based Fresnel transform (SFR), ASM/Fresnel transfer function approach (Fres-TF), and Rayleigh-Sommerfeld convolution (RSC)/Fresnel impulse response approach (Fres-IR). Note that Fres-TF is the paraxial version of ASM; and Fres-IR is the paraxial version of RSC. For each method, the characteristic of its frequency response is analyzed associated with the sampling theorem. After understanding the physics-consistent sampling request, the corresponding frequency sampling strategies are then proposed. Such strategies either guarantee the calculation correctness, which was affected in the conventional methods, or improve the calculation efficiency, which was limited by the wasted sampling resources. Besides, a further discussion related to the band-limited ASM (BLASM) [8] and band-extended ASM (BEASM) [9] is given. Here, frequency response refers to the spatial distribution of the diffraction field determined by the spatial frequency of the input field, which will be further explained with schematic diagrams in the following.

This paper is organized as following: part 2 briefly introduces the physical diffraction methods; part 3 discusses the properties of discrete computation; part 4 analyzes the frequency response of above mentioned methods and proposes corresponding strategies; parts 5 analyzes the reason why BLASM fails in far field and introduces the solution of BEASM; part 6 presents a discussion and concludes this work.

2. Physical diffraction methods

The diffraction model is shown in Fig. 1. The input field $u(x )$ is known and the task of diffraction calculation is to obtain the diffraction field $\tilde{u}(X )$.

 

Fig. 1. Schematic diagram of free-space diffraction.

Download Full Size | PDF

Here, we first summarize the mathematical expressions of five common-used methods in Table 1. Note that all the constants are neglected in the expressions. The first three methods are all based on the Fresnel approximation. Thus, analytically, they are equivalent. Since they are based on Fresnel approximation, they can be only used in the paraxial case. The last two methods, ASM and RSC, are also analytically equivalent with each other because the transfer function of ASM and the impulse response function of RSC are a pair of Fourier transformations. Because there is no paraxial approximation used in these two methods, they can be used in both paraxial and non-paraxial cases.

Table 1. Summary of several common-used physical diffraction methods

View Table | View all tables in this article

Physically, we can observe the diffraction field with any propagation distance z; inside an observation window with any size and located at any lateral position. However, after the discretization, all these flexibilities are limited in the discrete computation. More seriously, the accuracy of the diffraction field may be affected if the discrete computation is carried out in an incorrect way. Therefore, it is important to know the properties of the discrete computation.

3. Properties of discrete computation

As seen from Table 1, all these five methods need Fourier transform to get the results. In practice, Fourier transform is implemented by the fast Fourier transform (FFT) algorithm for a high computational efficiency.

Note that the latter four methods are all convolution-based. Take RSC as example, the diffraction field $\tilde{u}(X )$ is the linear convolution between the input field $u(x )$ and the impulse response function $h(x )$. To get a fast calculation, the linear convolution is carried out in the spatial frequency domain by multiplying their Fourier transformations because the computational complexity of multiplication is much cheaper than that of convolution. Then, inverse Fourier transform is implemented to the product for calculating the diffraction field. The theoretical basis for this operation is the convolution theorem. However, in discrete computation, the convolution theorem is for circular convolution not the linear convolution [20]. Therefore, to use the convolution theorem for a high computational efficiency, zero padding is required because circular convolution between zero-padded signals is equivalent to linear convolution between the original signals. That is to say padding zeros to $u(x )$ and $h(x )$ is necessary.

Now let us examine the properties of the discrete computation carefully. Suppose the input field $u(x )$ is sampled by N points with the sampling interval of ${\Delta_x}$; and the number of padded zeros is ${N_p}$; i.e. $N + {N_p}$ points for the zero-padded input field ${u_p}(x )$. According to the transform property of FFT, the sampling interval of its Fourier spectrum ${U_p}({{f_x}} )$ is

Equation (1) means that the sampling interval in the spatial frequency domain is determined by the sampling parameters in the space domain. This is the most important property of the FFT-based discrete computation. In the following, we would like to analyze each method in detail and provide corresponding frequency sampling strategy based on this property.

4. Frequency sampling strategies

4.1 SFR

Among the above five methods, SFR is the only one that is not convolution-based. Therefore, in the conventional method, zero padding is not required. However, it has been proven that the conventional method would cause aliasing error of high-frequency components [21]. In addition to the solution given in [21], here, we would like to provide another approach for an aliasing-free SFR calculation. As indicated in Eq. (2), there are two chirp functions used in SFR calculation. We call $\textrm{exp}({\pi {x^2}/\lambda z} )$ inner chirp and $\textrm{exp}({\pi {X^2}/\lambda z} )$ outer chirp.

To sample them correctly, the sampling theorem should be satisfied. The input field $u(x )$ is sampled by N points with the sampling interval of ${\Delta _x}$, so is the inner chirp. Therefore, we have

Equation (8) means only the transform result within this range is correct. According to the transform property of FFT, the whole transform result would spread on the range of $[{ - \lambda z/2{\Delta _x},\; \; \lambda z/2{\Delta _x}} ]$, but only the part inside the range given by Eq. (8) is correct because outside parts have aliasing error due to the incorrect sampling. It is evident that the correct range increases with the propagation distance.

Next, we need to calculate the product of the Fourier transformation and the outer chirp. Similarly, for the individual outer chirp,

Readers may ask that why N points before the Fourier transform and $\hat{N}$ points after the Fourier transform. The answer is that we can do the Fourier transform of any number of points by FFT. Note that when calculate Eq. (9), the maximum value of the left-hand side should be taken when $X = \lambda z/2{\Delta _x} - N{\Delta _x}/2$, because this is the largest value of correct sampling range.

Combining Eq. (9) and Eq. (10), we have

For the final complex amplitude of the diffraction field, there should be

Substituting Eq. (10) into Eq. (13), we can get

If $\hat{N} < N$, we can set $\hat{N} = N$ since there are N points as the input of the calculation.

Summarizing the above analysis, we can get the following calculation strategy:

 

Fig. 2. Flow chart of the proposed aliasing-free calculation strategy for SFR.

Download Full Size | PDF

To verify the proposed strategy, we carry out a numerical demonstration. A linear chirp grating is used as the input field. The reason why using such an input field is because it has rich frequency components so that it can clearly show the response of different frequencies. Furthermore, an input field with rich frequency components is more in line with the actual situation. The analytical expression of the input field and simulation parameters are listed in Table 2; and the propagation distance is z = 10 mm. Note that the expression of the input field is in a concise form of $|{\textrm{cos}({400\pi {x^2}} )} |\textrm{exp}[{i\varphi (x )} ]$. Phase $\varphi (x )= 0$ for those x values where $\textrm{cos}({400\pi {x^2}} )\ge 0$ and $\varphi (x )= \pi $ for those x values where $\cos ({400\pi {x^2}} )< 0$. For simplicity, $|{\textrm{cos}({400\pi {x^2}} )} |\textrm{exp}[{i\varphi (x )} ]$ is written as $\textrm{cos}({400\pi {x^2}} )$. This work is studied from the perspective of holography and diffractive optics. Commonly used sampling intervals are a few microns. Therefore, we choose a $2\; \mathrm{\mu }\textrm{m}$ sampling interval which is appropriate.

Table 2. Input field and simulation parameters used in numerical calculations

View Table | View all tables in this article

As shown in Fig. 3, where the results by the conventional method [22], proposed strategy, and analytical Fresnel integral (as reference) are given, the proposed strategy gives an accurate result by avoiding the aliasing error caused by the incorrect sampling. In contrast, the conventional method gives the aliased result because only considering the sampling of inner chirp cannot satisfy the sampling theorem. For the amplitude, there is aliasing error of the result by the conventional method because only the range given in Eq. (8) is correct. By cutting off the aliasing part, no error occurs in the result by the proposed strategy. For the phase, worse situation happens in the conventional method because the outer chirp cannot be correctly sampled. Most of the output field is wrong. On the contrary, all the output field given by the proposed strategy is accurate thanks to the $\hat{N}$-point FFT. For a clear comparison, the results by the conventional method within the corresponding range, which is $[{ - L/2,L/2} ]$ with $L = \lambda z/{\Delta _x} - N{\Delta _x} = 1.5\;\textrm{mm}$, are stretched to align with the results by the proposed strategy. By a detailed frequency analysis, the proposed strategy can give an accurate result, either the amplitude or the phase, whose correlation to that by analytical Fresnel integral is 1. Note that this strategy guarantees the transform correctness of $u(x )$ within the bandwidth of $B = [{N{\Delta _x}/2\lambda z - 1/2{\Delta _x},\; \; 1/2{\Delta _x} - N{\Delta _x}/2\lambda z} ]$, which is lower than the full bandwidth of $[{ - 1/2{\Delta _x},\; \; 1/2{\Delta _x}} ]$. The strategy for correct full-bandwidth transfer can be found in [21].

 

Fig. 3. Calculation results by conventional method, proposed strategy and analytical Fresnel integral. The x-axis unit is mm.

Download Full Size | PDF

All the above analyses are carried out mainly from the perspective of discrete computation. They should be consistent with the physics. The physical model is shown in Fig. 4, where the frequency response is intuitively illustrated. Before going into the detail, we would like to mention that the free-space propagation has the linear space-invariant property. Thus, conventional method only considers the central point of the input field. The input field is sampled by N points with the interval ${\Delta _x}$; and the Fourier transform is implemented by FFT. Therefore, the diffraction range of the output field is determined as $\lambda z/{\Delta _x}$. Within such a range, all the diffracted light from the central point can be correctly calculated. Problems come with the diffraction light from the edge points. If the spatial frequency of the input field located at the edge point is $1/2{\Delta _x}$, its diffracted light would exceed the received range, as indicated by the symbol ① in Fig. 4. Due to the circular property of FFT, this part would go into the received range from the other side, as indicated by the symbol ②, which would cause errors. Therefore, the range without errors is

 

Fig. 4. Physical model used to explain the SFR method.

Download Full Size | PDF

4.2 ASM and Fres-TF

ASM and Fres-TF model diffraction in the spatial frequency domain, where the transfer function is used, as indicated in Eq. (16). Because the transfer function of Fres-TF is the paraxial approximation version of that of ASM [18], we focus on the analysis of ASM in this section. Fres-TF case can be easily obtained under a paraxial approximation in a similar way.

As a convolution-based method, ASM usually requires zero padding for accurate calculation [8]. Traditionally, the number of padded zeros is equal to the sampling number of the input field, which is ${N_p} = N$. However, we find that it is not always necessary to pad such many zeros. The number of padded zeros increases from 0 to N with the propagation distance, which would be given with detailed analysis in the following. Due to the symmetrical operation of the ASM calculation indicated in Eq. (18), one FFT and one IFFT, the sampling interval of the output field is the same as the input field. And in most cases, the size of the output field is set to be the same as the input field. For ASM, the key point is to correctly discretize the transfer function.

According to the transform property of FFT as given in Eq. (1), the Fourier spectrum of the input field is sampled by ${N_p} + N$ points with the interval

Thus, the transfer function $H({{f_x}} )$ is also sampled with these parameters. By applying the sampling theorem to the phase kernel ${\varphi }({{f_x}} )= \frac{{2\pi }}{\lambda }z\sqrt {1 - {{({\lambda {f_x}} )}^2}} $ of $H({{f_x}} )$,

Substituting the highest achievable spatial frequency ${f_{x\_\textrm{max}}} = 1/2{\Delta _x}$ into Eq. (19) and after some simplifications, we can get

Equation (20) means that the transfer function can only be correctly sampled in such a distance range. Generally speaking, the number of padded zeros is no more than the sampling number of the input field, i.e. ${N_p} \le N$; because when the propagation distance is larger than the critical distance ${z_c}$, which is

We can see that the upper limit of propagation distance given by Eq. (20) varies with the number of padded zeros. Does this mean the number of padded zeros required at different distances is also different? The answer is yes. To quantitatively study this problem, we need to analyze it from the perspective of frequency response of the discrete ASM diffraction model. As shown in Fig. 5(a), physically, the diffraction range increases with the propagation distance. However, in ASM calculation, the output field has the same size as the input field. Thus, we need to deal with the exceeded part of the diffraction field. Numerically, the Fourier transform in Eq. (16) is implemented by the FFT algorithm, which has the circular property. The exceeded part would get into the output field from the other side, as shown in Fig. 5(a), which would cause calculation error. It is evident that the size of wrong part would increase with the propagation distance. This effect is usually called circular convolution errors. To get rid of it, zero padding is used to enlarge the calculation window, as shown in Fig. 5(b). Since the exceeded part would cause error, we can control the number of padded zeros to prevent the exceeded part from entering the output field. From Fig. 5(b), the number of padded zeros can be calculated as

 

Fig. 5. Schematic diagram to show the effect of circular convolution errors (a) and adaptive zero padding method to avoid it (b).

Download Full Size | PDF

 

Fig. 6. Number of padded zeros varies with the propagation distance.

Download Full Size | PDF

To prove that Eq. (22) can guarantee the correctness of ASM, we demonstrate a numerical calculation. The input information are listed in Table 2 and the propagation distance is set as z = 3 mm. Figure 7(a) shows the results with ${N_p} = 190$ padded zeros; and Fig. 7(b) shows the results with N = 500 padded zeros. The cross correlation between the amplitude fields by these two methods is 1, so is the cross correlation between the phase fields, which proves the effectiveness of the proposed method.

 

Fig. 7. Calculation results by the proposed adaptive zero padding ASM: (a) amplitude and (c) phase; by the conventional ASM: (b) amplitude and (d) phase. The x-axis unit is mm.

Download Full Size | PDF

The procedure of the proposed adaptive zero padding ASM, as flow chart shown in Fig. 8, can be concluded as following:

 

Fig. 8. Flow chart of the proposed adaptive zero padding ASM.

Download Full Size | PDF

4.3 RSC and Fres-IR

Different from ASM (FRes-TF), RSC (Fres-IR) models the diffraction in the space domain, where the impulse response function is used, as indicated in Eq. (24). Because the impulse response function of Fres-IR is the paraxial approximation version of that of RSC [18], here, we focus on the analysis of RSC. Fres-IR case can be easily obtained under a paraxial approximation in a similar way.

Also, as a convolution-based method, RSC needs zero padding for accurate calculation. However, zero padding is not only being used to get rid of circular convolution errors, but also to increase the transfer bandwidth. Let us go into detail of RSC. Due to the symmetrical operation of the RSC calculation indicated in Eq. (24), two FFTs and one IFFT, the sampling interval of the output field is the same with that of the input field. And in most cases, the size of the output field is set to be the same as the input field. For RSC, the key point is to correctly discretize the impulse response function, including the sampling number and sampling interval. There are two options: (1) modeling $h(x )$ by $N$ points and padding ${N_p}$ zeros; or (2) modeling $h(x )$ by $({N + {N_p}} )$ points. From the perspective of discrete computation, i.e. to get rid of the circular convolution errors, option (1) should be used with ${N_p} = N$ [20]. However, the appropriate option should be determined through the frequency response of the RSC model.

RSC is a spherical-wave model [18]. The impulse response function $h(x )$ describes a weighted spherical wave from one point of the input field; and the diffraction field is the coherent superposition of these spherical waves from all the points of the input field. Due to the linear space-invariant property of free-space diffraction, we just analyze the response from the central point of the input field and then extend to the edge points. As shown in Fig. 9(a), if option (1) is used, the response range of the central point (green line) cannot cover the full calculation window. Neither can the edge points. Therefore, the response from point A of the input field cannot reach point B of the output field, which limits the transfer bandwidth. In other words, at any point (except the central one), the output field is not the superposition of the response from all the points of the input plane. To achieve such a full-bandwidth transfer, option (2) is used, as shown in Fig. 9(b). Because the impulse response function is modeled by $N + {N_p} = 2N$ points, the response range of the central point of the input field can cover the full calculation window. However, the diffraction light from the edge point would exceed the calculation window and the exceeded part would get into the calculation window from the other side, which would cause errors in those areas. Fortunately, the output field in concern has a smaller size than the calculation window; and the calculation accuracy of the output field can be guaranteed, as shown in Fig. 9(b).

 

Fig. 9. Schematic diagram to show zero padding. (a) Option (1) is used with a limited transfer bandwidth; (b) option (2) is used with a full transfer bandwidth.

Download Full Size | PDF

There is one point needing attention. The diffraction angle $\beta = \textrm{arctan}({N{\Delta _x}/z} )$ shown in Fig. 9(b) should be no more than the highest potential diffraction angle $\theta = \textrm{arcsin}({\lambda /2{\Delta _x}} )$, therefore, we have

After some simplification, we get

To prove the above analysis, we carry out one numerical demonstration with the parameters shown in Table 2 and the propagation distance is set as z = 50 mm. Figure 10(a) shows the results by option (1); and Fig. 10(b) shows the results by option (2). It is clear that there are no circular convolution errors of the result by option (1) inside the whole calculation window, but the high-frequency components are lost, which cannot give an accurate physics-consistent output field as shown in Fig. 10(c) given by Rayleigh-Sommerfeld integral (RSI) [2]. In contrast, the result by option (2), a full-bandwidth calculation, has more high-frequency details which is consistent with the real physics. As a price, only the range of the same size as the input field is correct. Readers may notice that there is still a little difference between the result by option (2) and that by RSI. This is because only N discrete points of the input field contribute to the result by option (2); but the result by RSI is analytically calculated which would give a more accurate result. However, such a difference would not have much impact on calculation accuracy. Through the above analysis, we reveal that in RSC, zero padding not only avoids circular convolution errors, but also expands the transfer bandwidth, as comparison between Fig. 10(a) and Fig. 10(b).

 

Fig. 10. Calculated amplitude field by option (1): (a); option (2): (b); and RSI: (c). Figure 10(b1) is the zoomed part of Fig. 10(b) with a renormalization. The x-axis unit is mm.

Download Full Size | PDF

The procedure of option (2), as flow chart shown in Fig. 11, can be concluded as following:

 

Fig. 11. Flow chart of option (2) of the common-used RSC.

Download Full Size | PDF

4.4 Summary

In this section, we have analyzed several popular diffraction calculation methods. For each method, we provide a detailed frequency response analysis and give the corresponding sampling strategy. Readers may find that ASM is only applicable for near-field calculation and RSC is only for far-field calculation. To make ASM suitable for far field, some methods have been proposed [8,9,24,25]. One popular method of them is BLASM [8]. Compared with ASM, BLASM can work in a larger distance range. However, the performance of BLASM at long distance is also limited [9,25]. In the next section, we would give the reason why this happens although all the effective bandwidth is maintained in BLASM; and provide a solution of BEASM [9].

5. BLASM and BEASM

BLASM works by limiting the effective bandwidth of the transfer function $H({{f_x}} )$. When the propagation distance is larger than the critical distance given in Eq. (21), $H({{f_x}} )$ is aliased which would cause high-frequency noise [8,9]. By cutting off the aliased part, BLASM works in a larger distance range. It has been proven that the reserved bandwidth without aliasing is exactly what is needed for the output field. Here, we give a brief introduction of BLASM and detailed analysis can be found in [8].

Let us revisit Eq. (19) and set ${N_p} = N$. By regarding ${f_x}$ as the variable, we can get

 

Fig. 12. Schematic diagram of the bandwidth required in BLASM.

Download Full Size | PDF

A numerical demonstration is carried out to illustrate this effect. A rectangle aperture with the width of 0.8 mm is used as the input field. The sampling number and interval of the input field are N = 500 and ${\Delta _x}$ = 2 µm, respectively. The illumination wavelength is $\lambda = 500\;\textrm{nm}$ and the propagation distance is z = 300 mm. Figure 13(a) shows the result by BLASM. As analyzed in Section 4.3, RSC is applicable in far-field calculation. Therefore, we give the result by RSC in Fig. 13(b). From the comparison, we can see that the performance of BLASM is poor. The only difference between BLASM and RSC is the transfer function. In BLASM, the transfer function is directly modeled with $H({{f_x}} )$ and followed with a band limiting operation. In RSC, the transfer function is obtained by the Fourier transformation of the impulse response function $h(x )$. To show the difference, we draw these two transfer functions in Figs. 13(d) and (e). It is evident that BLASM has a hard truncation and RSC has a soft truncation. Therefore, it is the hard truncation that makes BLASM fail in relative far-field calculation.

 

Fig. 13. Calculation results by BLASM (a); by RSC (b) and by BEASM (c). Imaginary part of the transfer function of BLASM (d); RSC (e); and BEASM (f). The x-axis unit for (a-c) is mm and for (d-f) is $\textrm{m}{\textrm{m}^{ - 1}}$.

Download Full Size | PDF

To overcome this problem, we have previously proposed a BEASM for accurate calculation in a wide propagation distance range [9]. In BEASM, the bandwidth in the calculation is larger than that is physically needed. Therefore, for the needed bandwidth, there is no hard truncation any more, which gives an accurate calculation. For the same calculation task, the result by BEASM is shown in Fig. 13(c); and the transfer function of BEASM is shown Fig. 13(f). It is clear that the BEASM works better than BLASM, which verifies our analysis.

6. Discussion and conclusion

It is obvious that the Nyquist sampling criteria is used in this work. It is worth mentioning that there may be two kinds of numerical diffraction calculations. Diffraction calculation of an analytical signal and diffraction calculation of a discrete signal. For the analytical case, the first step is to discretize the signal with proper sampling parameters. After discretization, the numerical diffraction can be followed. The space bandwidth product (SBP) of the signal is a prior knowledge for the discretization and calculation [26]. For the discrete case, the sampling parameters, i.e., the sampling number and sampling interval, are first given, not determined by analytical signals, such as the diffraction calculation used in digital holography. Thus, we need to ensure the calculation correctness by applying the Nyquist criteria for a robust result. Of course, if the SBP of the discrete signal can be evaluated, more efficient calculation could be possible by considering the sampling of the evaluated SBP not the upper bound given by the sampling parameters. In this work, we study the discrete case.

To summarize, we have analyzed several popular diffraction calculation methods from a joint physics-mathematics perspective. For each method, the frequency response has been carefully studied, by which we found (1) the conventional single-Fourier transform based Fresnel diffraction (SFR) cannot deal with the full bandwidth of the input field; (2) there is no need to always pad N zeros in angular spectrum method (ASM); (3) zero padding in Rayleigh-Sommerfeld convolution (RSC) not only avoids circular convolution errors, but also expands the transfer bandwidth. Based on these findings, we have proposed a new aliasing-free SFR method and an adaptive zero-padding ASM. In addition, through the comparison between band-limited ASM (BLASM) and RSC, we found it is the hard truncation that makes BLASM fail in long propagation distance and introduced band-extended ASM as an effective solution. All the Matlab codes of this work can be found in [27].

The key contribution of this work is the idea of frequency response analysis based on the joint perspective of physics and mathematics. With this idea, for any diffraction calculation task, readers can make a do-it-yourself analysis and choose or develop a method for accurate and efficient computation. More flexible calculation methods, such as shifting and scaling, can also be analyzed with this idea.