Nyquist Sampling Conditions of Some Diffraction Algorithms with Adjustable Magnification

Diffraction algorithms with adjustable magnification are dominant in holographic projection and imaging. However, the algorithms are limited by the Nyquist sampling conditions, and simulation results with inappropriate parameters sometimes appear with aliasing. At present, many diffraction algorithms have been proposed and improved, but there is a need for an overall analysis of their sampling conditions. In this paper, some classical diffraction algorithms with adjustable magnification are summarized, and their sampling conditions in the case of plane wave or spherical wave illumination are analyzed and compared, which helps to select the appropriate diffraction algorithm according to the specific parameter conditions of the simulation to avoid aliasing.


Introduction
In the applications of large field of view diffraction projection, such as advanced vehicle front-lighting [1], holographic projection [2,3], and 3D sensing [4], the generation of holograms is a key problem [5]. This process inevitably requires the simulation of the forward and backward diffraction processes of light between the initial plane and the diffraction plane; therefore, a diffraction algorithm with adjustable magnification is necessary. Figure 1 shows a typical diffraction process. The complex amplitude distribution of the light field on the initial plane is denoted as U 0 (x 0 , y 0 , 0). The complex amplitude distribution at the diffraction distance ∆z on the diffraction plane is denoted as U(x, y, ∆z). The diffraction process can be described by two equations [6].
Since the Fourier transform can be calculated via fast Fourier transform (FFT) [6], Equation (3) is widely used in the study of diffraction problems. The single fast Fourier transform (S-FFT) [7] is a classical method for calculation of the diffraction process using FFT. This method uses FFT once to calculate the diffraction process. If we want to obtain diffraction results satisfying the sampling theorem, the number of sampling points, N, the initial plane size, 0, the diffraction plane size, , the sampling interval, ∆x0, and the scaling factor, m, of relative to 0 should satisfy the following condition: When the diffraction distance is small, the target plane size will be close to 0. In order to ensure that is not too small under the premise that 0 is fixed, a large number of sampling points, N, is necessary, which will increase the amount of calculations. However, in a large field of view diffraction projection system, the initial plane size is often much smaller than the diffraction plane size, and the use of S-FFT has great limitations. Therefore, some more flexible numerical algorithms have been developed to break the Since the Fourier transform can be calculated via fast Fourier transform (FFT) [6], Equation (3) is widely used in the study of diffraction problems. The single fast Fourier transform (S-FFT) [7] is a classical method for calculation of the diffraction process using FFT. This method uses FFT once to calculate the diffraction process. If we want to obtain diffraction results satisfying the sampling theorem, the number of sampling points, N, the initial plane size, L 0 , the diffraction plane size, L, the sampling interval, ∆x 0 , and the scaling factor, m, of L relative to L 0 should satisfy the following condition: When the diffraction distance is small, the target plane size L will be close to 0. In order to ensure that L is not too small under the premise that L 0 is fixed, a large number of sampling points, N, is necessary, which will increase the amount of calculations. However, in a large field of view diffraction projection system, the initial plane size is often much smaller than the diffraction plane size, and the use of S-FFT has great limitations. Therefore, some more flexible numerical algorithms have been developed to break the limit and make flexible calculations with adjustable magnification between the sizes of the initial plane and the diffraction plane.
One of them is the double Fresnel transform (DBFT) algorithm, which was proposed by Fucai Zhang et al. for the reconstruction of recorded digital holograms [8]. Later, an improved algorithm with similar principles was also proposed [9]. The algorithm involves two reconstruction steps implemented by conventional S-FFT, such that the magnification between the sizes of the initial plane and the diffraction plane can be adjusted flexibly by adjusting the diffraction distances of the two steps.
In the book Numerical Simulation of Optical Wave Propagation, the author Jason D. Schmidt mentioned an algorithm containing a scaling factor [10] (scaled angular spectrum method, SASM). Furthermore, Richard P. Muffletto et al. proposed the shifted Fresnel diffraction algorithm [11] (SFD), which can compute the Fresnel propagation of light between parallel planes with different sampling intervals. The common feature of the two algorithms is that the magnification between the sampling intervals of the initial plane and the diffraction plane can be set independently before simulation. Although the first three diffraction algorithms are improved compared with the ordinary S-FFT in terms of sampling conditions, they are still limited. If all algorithms are not applicable to the set simulation parameters, the algorithms [12][13][14] with relatively smaller sampling limit can also be used at the expense of computing efficiency. Their mathematical principle is similar to nonuniform fast Fourier transform (NUFFT) in signal processing theory [15]. One of them is the matrix product angular spectrum method (MPASM). This algorithm was proposed by Wanli Zhao et al. in 2020 [14]. In this method, a DFT calculation based on the matrix product is realized.
The ideas of these four algorithms to break the sampling limit of S-FFT are classical and representative. Therefore, their specific sampling conditions are further analyzed in this paper. Our conclusion can also be referred to when analyzing other algorithms with similar ideas.
In simulations of the diffraction process of the light field, the Nyquist sampling condition needs to be considered. That is, the condition that the sampling frequency is not less than twice of the maximum spatial frequency in the light field distribution needs to be satisfied, otherwise aliasing will occur.
Generally, there are two kinds of aliasing results. One is that the image is divided into several parts. The other is the appearance of horizontal and vertical high brightness stripes. They darken the whole image due to brightness normalization, which is also one of the bases for judging whether the image is aliased. These two different aliasing phenomena are caused by the failure of different steps in the calculation to meet the sampling conditions.
In most diffraction algorithms, the maximum spatial frequency of the light field generally depends on the phase factors in the diffraction equation, especially the quadratic phase factors. The spatial frequency of the amplitude is generally much smaller than that of the phase, which is often not considered [6]. In the general quadratic phase factors, the spatial frequency becomes larger upon approaching the field edge. Therefore, when analyzing the sampling frequency, the spatial frequency of the phase distribution at the four corners of the rectangular field is generally taken as the maximum spatial frequency for analysis, so as to ensure that every 2π change in the phase will pass through at least two sampling points. Different diffraction algorithms involve different phase factors; therefore, their sampling conditions are also different. There is an urgent need for overall analysis of their sampling conditions. Compared with the S-FFT algorithm with a simple calculation process, the calculation processes of these algorithms are more complex. Therefore, each step of the calculation process needs to be analyzed when analyzing the sampling conditions, and the mutual influence between the sampling limits of different parts should also be considered. In addition, compared with simple plane wave illumination, spherical wave illumination with an additional quadratic phase factor also makes the analysis of sampling conditions more complex.
Four classical diffraction algorithms are summarized in this paper: the double Fresnel transform [8] (DBFT), the scaled angular spectrum method [10] (SASM), the shifted Fresnel diffraction [11] (SFD), and the matrix product angular spectrum method [14] (MPASM). Their sampling conditions are analyzed and compared, which is helpful when selecting the appropriate diffraction algorithm according to specific simulation conditions to avoid aliasing. As a result, when using them, people often encounter aliasing due to inappropriate parameters, but do not know where and how to adjust them. Next, this paper briefly introduces the principles and analyzes the sampling conditions for each algorithm, and finally summarizes and compares them to provide important references.

Methods
Generally, in the diffraction equation, if the phase factors in the Fourier transform term do not satisfy the sampling conditions, the amplitude distribution of the simulated diffraction result will be aliased, and if the phase factors outside the Fourier transform term do not satisfy the sampling conditions, the phase distribution of the simulated diffraction result will be aliased. In display applications, while only the amplitude distribution of the diffracted result is considered, the phase factors outside the Fourier transform term are often ignored. Yet, for a specific diffraction algorithm, it is also necessary to specifically analyze which phase factors are related to the amplitude and which are related to the phase in the target plane before making a choice.
When the initial plane is illuminated by a spherical wave, a quadratic phase factor is introduced into it. The phase factor is related to the distance from the point light source to the initial plane. Generally, the phase factor is similar to the phase factors in the diffraction equation; thus, it can be combined with them to form new phase factors. Then, it is necessary to analyze the sampling conditions of the combined new phase factor. The results obtained by analyzing these phase factors separately and then taking their intersection are inaccurate.
On the basis of the discussions above, we briefly introduction the four algorithms and discuss their sampling conditions separately.

SFD
The SFD algorithm considers the condition that the center of the initial plane and the center of the diffraction plane are not coaxial and the sampling intervals of the two sampling planes are different. Taking the center point as a relative reference, each sampling point can be written as: where p 0 , q 0 , p, and q are within [−N/2, N/2 − 1] and N is the number of horizontal (vertical) sampling points. Only the case where the two sampling planes are square is considered here, and the sampling points of the two planes are the same. Then, according to Equation (3), we obtain: where Equation (7) contains a convolution operation and can be converted into a Fourier transform operation: Then, FFT can be used, and the calculation speed will be greatly accelerated. A total of three FFT operations are required. This algorithm has been used in holographic projections without lenses [16] and 3D holography [17].
In addition, Tomoyoshi Shimobaba's group further proposed an improved algorithm ARSS (aliasing reduced Fresnel diffraction with scale and shift operations) to address the drawback of SFD that aliasing will occur with a short diffraction distance [18]. By introducing a rectangular window, the algorithm eliminates the aliasing phenomenon with a short diffraction distance. They used the ARSS in a virtual converged spherical wave computergenerated holographic system [19], in which the method of cyclic iteration was used to SFD determines the scaling factor before deriving the diffraction process, and then extracts the convolution operation, so that FFT can be used for fast calculation.
The quadratic phase factors of a(p 0 , q 0 ) and b(p 0 , q 0 ) in Equation (9) are mainly analyzed. Here, in order to keep the simulation conditions consistent with other algorithms, only the case where the centers of the initial plane and diffraction plane are coaxial is considered. The scaling factor is still assumed to be m. Next, as shown in Figure 1, we consider the phase factor introduced by the spherical wave for illumination, which will also affect the sampling conditions. The phase factor introduced by the illumination spherical wave P div is expressed as: It should be noted that the spherical wave here is not necessarily a divergent spherical wave illuminated by a point light source. This is only true when r is > 0. When r is < 0, the spherical wave is convergent. When r tends to +∞, the incident light tends to the plane wave. Here, we only analyze the case where r is > 0 or r tends to +∞.
In order to satisfy the sampling theorem, the sampling frequency should be higher than twice the maximum spatial frequency. Generally, the phase frequency is much higher than the amplitude frequency. Therefore, only the phase factor in the diffraction equation is considered. Given that the maximum spatial frequency is located at the edge, where the sampling frequency shall satisfy the sampling theorem, for a(p 0 , q 0 ), we have: The analysis of the quadratic phase factor in b(p 0 , q 0 ) shows that: Next, further analysis is conducted to observe the results of F {a(p 0 , q 0 )}, F {b(p 0 , q 0 )}, and F {a(p 0 , q 0 )}F {b(p 0 , q 0 )}. Here, taking ∆z = 500 mm as an example, the analysis process is shown in Figure 2.
In Figure 2, the aliasing phenomenon occurs at the edge of F {b(p 0 , q 0 )}, but the final diffraction result is normal, because there is an invalid blank area around F {a(p 0 , q 0 )}. In addition, in the process of multiplying F {b(p 0 , q 0 )}, the aliased part in F {b(p 0 , q 0 )} does not affect the effective area at the center of F {a(p 0 , q 0 )}; hence, it does not affect the final result. According to this discussion, the analysis of the quadratic phase factor in b(p 0 , q 0 ) can be adjusted. It is not necessary to satisfy the sampling conditions at the farthest edge of the field, as long as the sampling conditions are satisfied at the edge of the effective area of F {a(p 0 , q 0 )}. Therefore, we need to calculate the range of the effective area in F {a(p 0 , q 0 )}. From the form of a(p 0 , q 0 ), we can get that F {a(p 0 , q 0 )} is actually the propagation result of the initial plane illuminated by a spherical wave of radius r through the distance ∆z/(1 − m), which will scale the size of image by [r(m − 1) − ∆z]/[r(m − 1)] times. In addition, after a single FFT calculation, the sampling range is scaled by [N∆x 0 2 (m − 1)]/(λ∆z) times. Given these two scales, the scale of the range of the effective area relative to the original image should be N∆x 0 2 [r(m − −1) − ∆z]/(λr∆z), according to which we can scale N in the sampling conditions obtained from the analysis of the quadratic phase factor in b(p 0 , q 0 ) to obtain a new sampling condition. It should be noted here that, in order to ensure that the amplitude of F {b(p 0 , q 0 )} does not appear aliasing at the boundary of the effective area, b(p 0 , q 0 ) is required to satisfy the sampling conditions in the diagonal direction at the corners of the Sensors 2023, 23, 1662 6 of 18 effective area, and its sampling interval is √ 2 times of the horizontal or vertical interval. Then, we obtain: The analysis of the quadratic phase factor in b(p0, q0) shows that: Next, further analysis is conducted to observe the results of {a(p0, q0)}, {b(p0, q0)}, and {a(p0, q0)}{b(p0, q0)}. Here, taking ∆z = 500 mm as an example, the analysis process is shown in Figure 2. In Figure 2, the aliasing phenomenon occurs at the edge of {b(p0, q0)}, but the final diffraction result is normal, because there is an invalid blank area around {a(p0, q0)}. In addition, in the process of multiplying {b(p0, q0)}, the aliased part in {b(p0, q0)} does not affect the effective area at the center of {a(p0, q0)}; hence, it does not affect the final result. According to this discussion, the analysis of the quadratic phase factor in b(p0, q0) can be adjusted. It is not necessary to satisfy the sampling conditions at the farthest edge of the field, as long as the sampling conditions are satisfied at the edge of the effective area of {a(p0, q0)}.
Therefore, we need to calculate the range of the effective area in {a(p0, q0)}. From the form of a(p0, q0), we can get that {a(p0, q0)} is actually the propagation result of the initial plane illuminated by a spherical wave of radius r through the distance ∆z/(1 − m), which will scale the size of image by [r(m − 1) − ∆z]/[r(m − 1)] times. In addition, after a single FFT calculation, the sampling range is scaled by [N∆x0 2 (m − 1)]/(λ∆z) times. Given these two scales, the scale of the range of the effective area relative to the original image should be N∆x0 2 [r(m − −1) − ∆z]/(λr∆z), according to which we can scale in the sampling conditions obtained from the analysis of the quadratic phase factor in b(p0, q0) to obtain a new sampling condition. It should be noted here that, in order to ensure that the amplitude of {b(p0, q0)} does not appear aliasing at the boundary of the effective area, b(p0, q0) is required to satisfy the sampling conditions in the diagonal direction at the corners of the Next, we analyze the phase factor in F {a(p 0 , q 0 )}F {b(p 0 , q 0 )}. U(p 0 , q 0 ) and the linear phase factor in a(p 0 , q 0 ) have little influence on F {a(p 0 , q 0 )}; thus, they can be ignored. Then, the phase factor of a(p 0 , q 0 ) can be regarded as exp{iπ[(1 − m)r + ∆z](x 0 2 + y 0 2 )/(λr∆z)}, and the phase factor of F {a(p 0 , q 0 )} can be calculated as exp{iπ(λr∆z)(f Meanwhile, the phase factor in F {b(p 0 , q 0 )} can be calculated as exp{−iπλ∆z(f x 2 + f y 2 )/m}. Multiplying the two phase factors, we obtain exp{iπλ∆z(r/[(m − 1)r − ∆z] − 1/m)(f x 2 + f y 2 )}, which meets the sampling conditions at the boundary of the effective area. Then, we obtain: Since Equation (13) is more restrictive than Equation (11), the final adjusted sampling conditions are: Next, we performed simulations to verify this sampling condition for the SFD algorithm. Two verifications were performed. One is to adjust the propagation distance ∆z with other parameters fixed, and the other is to adjust ∆x 0 and N with other parameters and initial plane size (L 0 = N∆x 0 ) fixed. For the first verification, the simulation parameters are set as shown in Table 1. According to Equation (15) and the simulation parameters, the sampling condition is satisfied when 450 mm ≤ ∆z ≤ 750 mm.

Parameters Values
Scaling factor (m) 6 Number of sampling points (N) 1080 Distance from point source to initial plane (r) 150 mm Figure 3 shows the amplitude distributions on the diffraction plane with different ∆z. It can be seen that the calculated sampling range is consistent with the simulation results. There is a hologram on the initial plane. In theory, with the condition of r = 150 mm given in Table 1, a clear image in focus will be obtained on the diffraction plane 600 mm behind it. All verifications in this paper use the hologram as the initial plane. When ∆z decreases to 450 mm, horizontal and vertical bright stripes begin to appear in the middle of the image, and the image brightness starts to decrease. With a further decrease in ∆z, the phenomenon of aliasing becomes more serious. Multiple images overlap and interlace, and the image brightness lowers. When ∆z increases to 750 mm, bright stripes parallel to the edges begin to appear around the image, and the image brightness starts to decrease. With the further increase in ∆z, the aliasing becomes more serious, and the image brightness lowers.  For the second verification, the simulation parameters are set as shown in Table 2. According to Equation (15) and the simulation parameters, the sampling condition is satisfied when ∆x0 ≤ 15.086 μm.  Figure 4 shows the amplitude distributions on the diffraction plane with different ∆x0 and N. The calculated sampling range is also consistent with the simulation results. When ∆x0 increases to 15.086 μm, multiple images overlap and interlace, horizontal and vertical bright stripes begin to appear, and the image brightness starts to decrease. With the fur- For the second verification, the simulation parameters are set as shown in Table 2. According to Equation (15) and the simulation parameters, the sampling condition is satisfied when ∆x 0 ≤ 15.086 µm.  Figure 4 shows the amplitude distributions on the diffraction plane with different ∆x 0 and N. The calculated sampling range is also consistent with the simulation results. When ∆x 0 increases to 15.086 µm, multiple images overlap and interlace, horizontal and vertical bright stripes begin to appear, and the image brightness starts to decrease. With the further increase in ∆x 0 , the phenomenon of aliasing becomes more serious, and the image becomes illegible. The results prove that the sampling conditions are correct from another perspective. In addition, this sampling condition is also applicable to the Fresnel-Bluestein algorithm, the mathematical essence of which is similar to that of SFD [20].

SASM
In the book Numerical Simulation of Optical Wave Propagation [10], the author Jason D. Schmidt mentioned an algorithm, which assumes that both sides are coaxial, and also introduces the scaling factor m, but uses FFT once less than SFD. The author uses angular spectrum propagation to define this method because its calculation form is similar to the angular spectrum form of the Fresnel diffraction integral. This method still uses a paraxial approximation, so it is not strictly equivalent to the true angular spectrum form.
In order to simplify the derivation process of SASM, the following operators are introduced: where ∆z is the diffraction distance and r0 and r denote the coordinate vectors at the initial plane and the diffraction plane, respectively. Then, a size scaling factor m of the diffraction plane relative to the initial plane is introduced, and the identity transformation can be performed:

SASM
In the book Numerical Simulation of Optical Wave Propagation [10], the author Jason D. Schmidt mentioned an algorithm, which assumes that both sides are coaxial, and also introduces the scaling factor m, but uses FFT once less than SFD. The author uses angular spectrum propagation to define this method because its calculation form is similar to the angular spectrum form of the Fresnel diffraction integral. This method still uses a paraxial approximation, so it is not strictly equivalent to the true angular spectrum form.
In order to simplify the derivation process of SASM, the following operators are introduced: (16) where the operators Q[c, r], F [r, f], and F −1 [r, f] in the equation indicate multiplication by the phase factor, Fourier transform, and inverse Fourier transform, respectively. The scalar forms of the vector symbol are r = (x, y) and f = (f x , f y ), representing the spatial coordinate and the frequency domain coordinate, respectively. Next is the derivation of SASM. First, the Fresnel diffraction equation can be written as: where ∆z is the diffraction distance and r 0 and r denote the coordinate vectors at the initial plane and the diffraction plane, respectively. Then, a size scaling factor m of the diffraction plane relative to the initial plane is introduced, and the identity transformation can be performed: With Equation (17), we can further obtain: Finally, the convolution theorem can be used to obtain U(r 2 ) as: In this algorithm, the scaling factor is determined before the diffraction process is derived, and the convolution operation is extracted such that FFT can be used for fast calculations.
This is similar to SFD in essence. The difference is that, for the Fourier transform of b in SFD and h in SASM, SFD calculates the Fourier transform, while SASM directly uses the phase factor of the calculation result. The use of FFT in the calculation process is an important factor influencing the sampling conditions; therefore, there are still differences between them. That is, SASM does not have to satisfy the second sampling condition in Equation (13) of SFD, which is obtained from the analysis of b(p 0 , q 0 ).
The final sampling conditions are as follows: Next, we performed simulations to verify these sampling conditions for the SASM algorithm. The two verifications were similar to those of SFD. For the first, the simulation parameters were also set as shown in Table 1. According to Equation (21) and the simulation parameters, the sampling condition is satisfied when 316 mm ≤ ∆z ≤ 750 mm. Figure 5 shows the amplitude distributions of the diffraction plane with different ∆z. It can be seen that the calculated sampling range is consistent with the simulation results. Within the calculated range, the image is only scaled. Beyond this range, aliasing occurs. When ∆z is large, the result of SASM is the same as that of SFD. When ∆z is reduced to 316 mm, the edge parts of the image start to separate from the center part, and the image is divided into nine parts, which is different from the common effect of aliasing containing replications, and is caused by the violation of the Nyquist sampling theorem by the second Fourier transform in Equation (20). With a further decrease in ∆z, the segmented part becomes larger and farther away from the center.
For the second verification, the simulation parameters were set as shown in Table 2. According to Equation (21) and the simulation parameters, the sampling condition is satisfied when ∆x 0 ≤ 43.944 µm. Figure 6 shows the amplitude distributions on the diffraction plane with different ∆x 0 and N. The calculated sampling range is also consistent with the simulation results. When ∆x 0 increases to 43.944 µm, the edge of the image begins to blur, and the image brightness starts to decrease. With a further increase in ∆x 0 , the phenomenon of aliasing becomes more serious; the edge parts of the image separate from the center part evidently such that the image is divided into nine parts and the whole image becomes blurred. It can be seen that the sampling conditions of SASM are not as strict as those of SFD due to the lack of the sampling limit for b(p 0 , q 0 ). For the second verification, the simulation parameters were set as shown in Table 2. According to Equation (21) and the simulation parameters, the sampling condition is satisfied when ∆x0 ≤ 43.944 μm. Figure 6 shows the amplitude distributions on the diffraction plane with different ∆x0 and N. The calculated sampling range is also consistent with the simulation results. When ∆x0 increases to 43.944 μm, the edge of the image begins to blur, and the image brightness starts to decrease. With a further increase in ∆x0, the phenomenon of aliasing becomes more serious; the edge parts of the image separate from the center part evidently such that the image is divided into nine parts and the whole image becomes blurred. It can be seen that the sampling conditions of SASM are not as strict as those of SFD due to the lack of the sampling limit for b(p0, q0).  Figure 7 is the schematic of DBFT. 0( 0, 0) and ( , ) represent the complex amplitude distribution of light field in the initial plane and the diffraction plane, respectively, while 1( 1, 1) is the complex amplitude distribution of the light field in a virtual intermediate plane that does not physically exist, which is indicated by the dotted lines. The diffraction process of the light field from 0( 0, 0) to 1( 1, 1), and from 1( 1, 1) to ( , ) is taken into account, such that the size of the diffraction plane can be adjusted by changing the position of the virtual intermediate plane.  For the second verification, the simulation parameters were set as shown in Table 2. According to Equation (21) and the simulation parameters, the sampling condition is satisfied when ∆x0 ≤ 43.944 μm. Figure 6 shows the amplitude distributions on the diffraction plane with different ∆x0 and N. The calculated sampling range is also consistent with the simulation results. When ∆x0 increases to 43.944 μm, the edge of the image begins to blur, and the image brightness starts to decrease. With a further increase in ∆x0, the phenomenon of aliasing becomes more serious; the edge parts of the image separate from the center part evidently such that the image is divided into nine parts and the whole image becomes blurred. It can be seen that the sampling conditions of SASM are not as strict as those of SFD due to the lack of the sampling limit for b(p0, q0).  Figure 7 is the schematic of DBFT. 0( 0, 0) and ( , ) represent the complex amplitude distribution of light field in the initial plane and the diffraction plane, respectively, while 1( 1, 1) is the complex amplitude distribution of the light field in a virtual intermediate plane that does not physically exist, which is indicated by the dotted lines. The diffraction process of the light field from 0( 0, 0) to 1( 1, 1), and from 1( 1, 1) to ( , ) is taken into account, such that the size of the diffraction plane can be adjusted by changing the position of the virtual intermediate plane.  Figure 7 is the schematic of DBFT. U 0 (x 0 , y 0 ) and U(x, y) represent the complex amplitude distribution of light field in the initial plane and the diffraction plane, respectively, while U 1 (x 1 , y 1 ) is the complex amplitude distribution of the light field in a virtual intermediate plane that does not physically exist, which is indicated by the dotted lines. The diffraction process of the light field from U 0 (x 0 , y 0 ) to U 1 (x 1 , y 1 ), and from U 1 (x 1 , y 1 ) to U(x, y) is taken into account, such that the size of the diffraction plane can be adjusted by changing the position of the virtual intermediate plane.

DBFT
The sampling intervals on the three planes are ∆x 0 = ∆y 0 = L 0 /N, ∆x 1 = ∆y 1 = L 1 /N, and ∆x = ∆y = L/N, where L represents the sampling range and N represents the number of sampling points. The sampling range of the frequency domain (L u , L v ) has the following mathematical relation with the spatial domain sampling interval in the discrete Fourier transform (DFT):  The sampling intervals on the three planes are ∆ 0 = ∆ 0 = 0/ , ∆ 1 = ∆ 1 = 1/ , and ∆ = ∆ = / , where represents the sampling range and represents the number of sampling points. The sampling range of the frequency domain ( , ) has the following mathematical relation with the spatial domain sampling interval in the discrete Fourier transform (DFT): Therefore, for the three planes, the following relations can be obtained: Then, we have the following relation: The relation between 0 and is: So far, the derivation of DBFT has been completed. The core idea to solve the sampling problem is as follows: in single FFT, the sampling interval in the spatial domain and that in the frequency domain are reciprocal; hence, the size is not large enough. However, with double FFTs, a positive proportional relation will be obtained after two reciprocal relations, as shown in Equation (25), which is applicable to the diffraction process.
It is worth noting that the absolute value is used in Equation (25), which indicates that the distance parameter can take a negative value, i.e., the intermediate virtual plane does not have to be between the initial plane and the diffraction plane but can also be placed on the left side of the initial plane or the right side of the diffraction plane. For example, the double sampling Fresnel (DSF) algorithm with similar mathematical essence places the intermediate plane at the point light source of the spherical wave for illumination [21].
From 0( 0, 0) to 1( 1, 1), The phase factor outside the Fourier transform term only affects the phase distribution of 1( 1, 1); thus, it is not considered. Only the phase factor exp[(ik/2z1)(x0 2 + y0 2 )] Therefore, for the three planes, the following relations can be obtained: Then, we have the following relation: The relation between L 0 and L is: So far, the derivation of DBFT has been completed. The core idea to solve the sampling problem is as follows: in single FFT, the sampling interval in the spatial domain and that in the frequency domain are reciprocal; hence, the size L is not large enough. However, with double FFTs, a positive proportional relation will be obtained after two reciprocal relations, as shown in Equation (25), which is applicable to the diffraction process.
It is worth noting that the absolute value is used in Equation (25), which indicates that the distance parameter can take a negative value, i.e., the intermediate virtual plane does not have to be between the initial plane and the diffraction plane but can also be placed on the left side of the initial plane or the right side of the diffraction plane. For example, the double sampling Fresnel (DSF) algorithm with similar mathematical essence places the intermediate plane at the point light source of the spherical wave for illumination [21].
From U 0 (x 0 , y 0 ) to U 1 (x 1 , y 1 ), The phase factor outside the Fourier transform term only affects the phase distribution of U 1 (x 1 , y 1 ); thus, it is not considered. Only the phase factor exp[(ik/2z 1 ) (x 0 2 + y 0 2 )] inside the Fourier transform term is considered. Thus, we obtain the sampling condition |z 1 | ≥ (∆x 0 2 N)/λ. From U 1 (x 1 , y 1 ) to U(x, y), Similarly, we obtain the sampling condition |z 2 | ≥ (∆x 1 2 N)/λ. Next, further analysis is conducted. The spherical wave phase factor can be combined with the phase factor in the diffraction equation. In addition, the phase factor of U 1 (x 1 , y 1 ) can also be combined with the phase factor in the Fourier transform term in the second step. The phase factor of U 1 (x 1 , y 1 ) is composed of two parts, one is the phase factor outside the Fourier transform term, and the other is the phase factor of the new term obtained by the Fourier operation of the part inside the Fourier transform term. The influence of the amplitude term U 0 (x 0 , y 0 ) in the Fourier transform term can be ignored; hence, the phase factor of the new term obtained after calculation is exp{−iπr(x 1 2 + y 1 2 )/[λz 1 (r + z 1 )]}. In addition, similar to the SFD analysis process, the phase factor inside the Fourier transform term in U(x, y) only has to satisfy the sampling conditions at the boundary of effective area of U 1 (x 1 , y 1 ). The calculation shows that the scale of the range of the effective area relative to the original image is N∆x 0 2 (r + z 1 )/(λrz 1 ). Therefore, the sampling conditions of the combined phase factor can be analyzed.
The final sampling conditions can be obtained as Next, we performed simulations to verify this sampling condition for the DBFT algorithm. The simulation parameters were set as shown in Table 1. According to Equation (28) and the simulation parameters, the sampling condition is satisfied when −150 mm ≤ z 1 ≤ −63 mm for any ∆z = (z 1 + z 2 ) > 0 mm. Figure 8 shows the amplitude distributions of the diffraction plane with different z 1 . It can be seen that the calculated sampling range is consistent with the simulation results. When z 1 is reduced to −150 mm, bright stripes parallel to the edge begin to appear around the image, and the image brightness starts to decrease. With a further reduction in z 1 , the aliasing becomes more serious and the image brightness becomes lower. When z 1 is increased to −63 mm, the edge part of the image starts to separate from the center part, and the image is divided into nine parts. With a further increase in z 1 , the segmented part becomes larger and further away from the center. When z 1 exceeds 97 mm, horizontal and vertical bright stripes appear in the middle of the image and the image brightness decreases.
Similarly, we obtain the sampling condition |z2|≥(∆x1 2 N)/λ. Next, further analysis is conducted. The spherical wave phase factor can be combined with the phase factor in the diffraction equation. In addition, the phase factor of 1 ( 1,1) can also be combined with the phase factor in the Fourier transform term in the second step. The phase factor of 1( 1, 1) is composed of two parts, one is the phase factor outside the Fourier transform term, and the other is the phase factor of the new term obtained by the Fourier operation of the part inside the Fourier transform term. The influence of the amplitude term 0( 0, 0) in the Fourier transform term can be ignored; hence, the phase factor of the new term obtained after calculation is exp{−iπr(x1 2 + y1 2 )/[λz1(r + z1)]}.
In addition, similar to the SFD analysis process, the phase factor inside the Fourier transform term in ( , ) only has to satisfy the sampling conditions at the boundary of effective area of 1 ( 1, 1). The calculation shows that the scale of the range of the effective area relative to the original image is N∆x0 2 (r + z1)/(λrz1). Therefore, the sampling conditions of the combined phase factor can be analyzed.
The final sampling conditions can be obtained as (28) Next, we performed simulations to verify this sampling condition for the DBFT algorithm. The simulation parameters were set as shown in Table 1. According to Equation (28) and the simulation parameters, the sampling condition is satisfied when −150 mm ≤ z1 ≤ −63 mm for any ∆z = (z1 + z2) > 0 mm. Figure 8 shows the amplitude distributions of the diffraction plane with different z1. It can be seen that the calculated sampling range is consistent with the simulation results. When z1 is reduced to −150 mm, bright stripes parallel to the edge begin to appear around the image, and the image brightness starts to decrease. With a further reduction in z1, the aliasing becomes more serious and the image brightness becomes lower. When z1 is increased to −63 mm, the edge part of the image starts to separate from the center part, and the image is divided into nine parts. With a further increase in z1, the segmented part becomes larger and further away from the center. When z1 exceeds 97 mm, horizontal and vertical bright stripes appear in the middle of the image and the image brightness decreases. Next, in order to keep the same form with the sampling conditions of other algorithms for easy comparison, z 1 and z 2 in the sampling conditions are expressed in terms of m and ∆z according to Equation (25). Here, two cases are considered. One is z 2 /z 1 = m and the other is z 2 /z 1 = −m. For the former, the final sampling condition after conversion is: For the latter, the final sampling condition after conversion is: It can be seen that the sampling condition of the latter is better, which is consistent with that of SASM. Therefore, in order to avoid aliasing as much as possible, the virtual intermediate plane should be placed on the left side of the initial plane rather than between the initial plane and the diffraction plane when the algorithm is used for large field of view diffraction calculations. Next, two verifications similar to those of other algorithms were performed. Here, m is fixed, which means that the ratio of z 2 /z 1 is fixed. In addition, z 1 is always negative while z 2 is always positive.
For the first verification, the simulation parameters were also set as shown in Table 1. According to Equation (30) and the simulation parameters, the sampling condition is satisfied when 316 mm ≤ ∆z ≤ 750 mm.
For the second verification, the simulation parameters were set as shown in Table 2. According to Equation (30) and the simulation parameters, the sampling condition is satisfied when ∆x 0 ≤ 43.944 µm.
It can be seen from Figures 9 and 10 that not only are the final sampling conditions of DBFT and SASM consistent, but their verification results with the same parameters are also consistent, which shows that although the overall ideas of the two algorithms are completely different, their mathematical essence is similar.
For the latter, the final sampling condition after conversion is: It can be seen that the sampling condition of the latter is better, which is consistent with that of SASM. Therefore, in order to avoid aliasing as much as possible, the virtual intermediate plane should be placed on the left side of the initial plane rather than between the initial plane and the diffraction plane when the algorithm is used for large field of view diffraction calculations. Next, two verifications similar to those of other algorithms were performed. Here, m is fixed, which means that the ratio of z2/z1 is fixed. In addition, z1 is always negative while z2 is always positive.
For the first verification, the simulation parameters were also set as shown in Table  1. According to Equation (30) and the simulation parameters, the sampling condition is satisfied when 316 mm ≤ ∆z ≤ 750 mm.
For the second verification, the simulation parameters were set as shown in Table 2. According to Equation (30) and the simulation parameters, the sampling condition is satisfied when ∆x0 ≤ 43.944 μm.
It can be seen from Figures 9 and 10 that not only are the final sampling conditions of DBFT and SASM consistent, but their verification results with the same parameters are also consistent, which shows that although the overall ideas of the two algorithms are completely different, their mathematical essence is similar.

MPASM
It can be seen that there are still sampling restrictions for the above three algorithms. If all algorithms are not applicable, the matrix product angular spectrum method (MPASM) can be considered. This algorithm was proposed by Wanli Zhao et al. in 2020 [14]. In this method, the calculation of the true angular-spectrum form without paraxial approximation based on matrix product is realized.
To get rid of the limitation ∆fx = 1/(N∆x), this algorithm uses the matrix product instead of FFT to calculate the diffraction process. Moreover, the frequency domain range can be adjusted by adjusting the number of frequency domain sampling points to satisfy

MPASM
It can be seen that there are still sampling restrictions for the above three algorithms. If all algorithms are not applicable, the matrix product angular spectrum method (MPASM) can be considered. This algorithm was proposed by Wanli Zhao et al. in 2020 [14]. In this simulation results are shown in Figures 11 and 12. It can be seen that when s is small, the sampling conditions of the algorithm are relatively strict, and aliasing is easy to occur. With an increase in s, the sampling condition of the algorithm becomes looser, and the upper limits of ∆z and ∆x 0 become higher.

Results
The characteristics and sampling conditions of the algorithms are analyzed and summarized in Table 3, where ∆x0 is the sampling interval of the initial plane, is the number of horizontal (vertical) sampling points, L0 is the size of initial plane, λ is the wavelength of light, m is the scaling factor of the diffraction plane relative to the initial plane, ∆z is the

Results
The characteristics and sampling conditions of the algorithms are analyzed and summarized in Table 3, where ∆x0 is the sampling interval of the initial plane, is the number of horizontal (vertical) sampling points, L0 is the size of initial plane, λ is the wavelength of light, m is the scaling factor of the diffraction plane relative to the initial plane, ∆z is the The parameters were consistent with the simulation of SFD (Tables 1 and 2).

Results
The characteristics and sampling conditions of the algorithms are analyzed and summarized in Table 3, where ∆x 0 is the sampling interval of the initial plane, N is the number of horizontal (vertical) sampling points, L 0 is the size of initial plane, λ is the wavelength of light, m is the scaling factor of the diffraction plane relative to the initial plane, ∆z is the diffraction distance, and r is the distance between the initial plane and the center of spherical wave for illumination (for a plane wave, r tends to +∞). It can be seen from the results that for the simulation of large field of view diffraction, three FFT-based algorithms (SFD, SASM, and DBFT) have upper limits on the diffraction distance, ∆z, and sampling pitch, ∆x 0 . Both ∆z and ∆x 0 need to satisfy the sampling conditions to avoid aliasing. If one of them does not satisfy the sampling conditions, a smaller value of the other one will not help. In general, these three algorithms have the same restriction on the diffraction distance, ∆z. For the limits on sampling pitch ∆x 0 , SASM and DBFT have the same restriction, while SFD has a stricter restriction. Lastly, for the MPASM algorithm based on matrix calculation, we can break through any limits by increasing the parameter s to avoid aliasing, but the computational efficiency is generally much lower than other algorithms.

Discussion
In the simulation of large field of view diffraction, in the case that the centers of the initial plane and diffraction plane are coaxial, SASM and DBFT are more recommended from the aspect of aliasing because of looser sampling restrictions. Before the simulation or after the aliasing occurs during the simulation, the specific parameters can be substituted into the formula in Table 3 to check whether they satisfy the sampling conditions. This paper recommends that the check and adjustment of ∆z take precedence over ∆x 0 , because the limit on ∆z is independent of ∆x 0 , while the limit on ∆x 0 is related to ∆z. In addition, the restrictions on ∆z of various FFT-based algorithms are basically consistent. Therefore, if ∆z does not satisfy the sampling conditions and the parameters ∆z, m, and r involved are not adjustable, a smaller value of ∆x 0 will not help; therefore, the FFT-based algorithms should be abandoned, and the MPASM algorithm should be used for calculation to avoid aliasing. However, this generally reduces computational efficiency. If ∆z satisfies the condition but there is still aliasing, it means that ∆x 0 does not satisfy the condition and needs to be reduced to below the limit. If it is difficult to reduce ∆x 0 , other parameters can be adjusted to increase the upper limit of ∆x 0 according to the specific situation, while ensuring that ∆z still satisfies the sampling condition. If both ∆z and ∆x 0 satisfy the sampling conditions, there will be no aliasing for FFT-based algorithms.
In addition, since the use of FFT is the source of sampling restrictions, fewer uses of FFT in the algorithm lead to looser sampling restrictions. Therefore, when using FFT-based algorithms, the number of times FFT is used in the algorithm should be paid attention to and minimized. For example, the Fourier transform of some expressions without information of the initial plane can be calculated manually. The algorithm should directly deduce the analytical expression instead of using FFT to calculate the Fourier transform process. Generally, sampling conditions consistent with those of SASM and DBFT can be obtained by reducing the use of FFT to two times. Further use of FFT will make the sampling restrictions stricter.
Lastly, for FFT-based algorithms, this paper summarizes the following problems that are easy to ignore or may be confusing, which can provide reference for the process of analyzing the sampling conditions of other similar algorithms: 1. When analyzing whether the phase distribution in FFT satisfies the Nyquist sampling conditions, sometimes there are two phase factors multiplied. If the two phase factors have similar forms, they should be combined and analyzed to obtain the correct results. If the two phase factors are analyzed separately, the restrictions obtained will be stricter. For example, the phase factor in U 1 (x 1 , y 1 ) can be combined with exp[(ik/2z 2 )(x 1 2 + y 1 2 )] in Equation (27); 2. When analyzing whether the phase distribution in FFT satisfies the Nyquist sampling conditions, it is necessary to pay attention to its corresponding amplitude distribution synchronously. If there is a blank area in the amplitude distribution, then, even if the phase distribution in this area breaks the Nyquist sampling conditions, it will not affect the results. Accordingly, the sampling restrictions can be further loosened. The process shown in Figure 2 is an example; 3. Sometimes it is necessary to analyze complex phase factors, such as F {a(p 0 , q 0 )} in Equation (9). The theoretical expression of its phase distribution cannot be obtained because a(p 0 , q 0 ) in Equation (8) is very complex. However, the factor with less influence in a(p 0 , q 0 ) can be ignored to simplify the problem. For example, the influence of the amplitude term and the primary phase term in a(p 0 , q 0 ) on the phase distribution of F {a(p 0 , q 0 )} can be explored. If the influence is small, we can ignore them to simplify a(p 0 , q 0 ) to a quadratic phase factor, and then calculate the simple approximate expression of F {a(p 0 , q 0 )}.
Lastly, we qualitatively explored the definition of diffraction images at various propagation distances near the focal length at different sampling rates to study the influence of sampling choice on 3D factors. In general, for all propagation distances, including the focal length, when the sampling rate is higher than the original resolution of the initial hologram, the results of different sampling rates do not change. When the sampling rate is lower than the original resolution of the initial hologram, a lower sampling rate will lead to lower diffraction image definition. We plan to go deeper into this issue in another article.