## Abstract

In this paper, we address the problem of calculating Fresnel diffraction integrals using a finite number of uniformly spaced samples. General and simple sampling rules of thumb are derived that allow the user to calculate the distribution for any propagation distance. It is shown how these rules can be extended to fast-Fourier-transform-based algorithms to increase calculation efficiency. A comparison with other theoretical approaches is made.

© 2014 Optical Society of America

## 1. INTRODUCTION

The propagation and diffraction of light is a classical problem in optics that remains important today. Maxwell’s equations provide a rigorous framework for examining the behavior of light, describing it as an electromagnetic wave. The behavior of both the electric and magnetic field components needs to be specifically considered: for example, when analyzing high numerical aperture lenses, such as those found in microscopes, or when modeling the propagation of surface plasmons. There are many occasions, however, when a full vectorial description of the electromagnetic field is not required to model a given physical problem, and a simpler scalar model can be used in its place. Using a scalar model to describe an optical field is a useful approximation. In the optical regime, it provides highly accurate predictions for the diffraction of light from apertures (where the aperture size is much larger than the wavelength), see for example the experimental results in [1]. Other relevant theoretical and numerical work is presented in [2,3]. In this paper, we assume that a scalar model is sufficiently accurate for our purposes.

Within the scalar description of light propagation, it is useful to distinguish two subgroups: nonparaxial and paraxial diffraction integrals. Both types of diffraction integrals can be derived from Maxwell’s equations, see [4,5]. However, the paraxial treatment requires that more approximations are made in the derivation process. Hence, the paraxial model seems further removed from the real underlying physical processes of diffraction, with the result that the nonparaxial diffraction integrals are considered to be more accurate. We refer to the treatment of diffraction presented in Chapters 3 and 4 of [4] and accordingly catagorize the Kirchhoff or Raleigh–Sommerfeld I and II as nonparaxial diffraction integrals. The Fresnel transform is a paraxial diffraction integral.

In this paper, we examine the numerical calculation of the Fresnel diffraction integral using a finite number of samples. This is because numerical techniques are of considerable importance for simulating the propagation of light and for calculating solutions for general diffraction problems that do not have an analytical solution. The Fresnel transform plays an indispensable role in digital holography [6–12], iterative phase retrieval techniques [13–15], and transport of intensity methods for estimating the phase distribution, which has an important role in biological imaging applications [16,17]. The Fresnel transform is also widely used in the design of diffractive optical elements for beam shaping [18,19], in metrology applications, where decorrelation of diffracted speckle fields are used to determine deformations such as stress, strain, and cracks in an optically rough surface under test, [20–22], or in the detection of displacement and tilting [23,24]. In all of these applications, numerical calculation of the Fresnel transform is necessary and hence it is of significant practical and theoretical importance to fully understand how the calculation should be performed.

Perhaps the first issue one encounters when implementing a numerical calculation is to choose an appropriate number of uniformly spaced samples to represent the complex scalar wavefield in the input plane. There are a number of manuscripts in the literature that consider this problem by examining how many samples are required so that the Nyquist sampling rate (or a related criterion) is adhered to when sampling the kernel or the Fresnel chirp function, see for example [25–31]. While this approach can yield accurate numerical solutions, it can often be unnecessarily numerically intensive [28]. The size of the input and output extents of the calculated signal is found to also affect the optimal sampling rate, however, there appears to be some divergence of views on the ideal sampling rate in the literature, see for example [29–31]. Often it is not clear how inaccurate the numerical calculation becomes as these sampling rules are violated. For these reasons, it is necessary to reexamine this issue of an appropriate sampling rate for the calculation of diffraction integrals. Although our analysis is performed for the Fresnel transform, the results generalize to the nonparaxial integrals discussed in [29–31].

Here, we approach this sampling question in a completely different manner, developing an intuitive physical description of the numerical diffraction process. This description is useful both for analyzing general diffraction problems and for determining exactly the expected errors introduced by the sampling process. In contrast to the other manuscripts in the literature, we do not consider the role of the Fresnel kernel or chirp function in the calculation at all when determining an appropriate sampling rate. Rather, we define a sampling rate based on the spatial and spatial frequency extent of the input field and on the location of the replicas in the diffraction plane. These replicas arise due to the sampling operation (discretization of the complex signal at Plane 1, see Fig. 1 and Section 3). The separation between these replicas and their properties depends on how the Fresnel transform is calculated. Provided that the output extent of the signal we are calculating is less than the separation between neighboring replicas, we consider our signal to be properly sampled.

The argument made in this paper rests on the following considerations:

- (1) A signal that has an effective finite support in the spatial and spatial frequency domains is defined using concepts from the space–bandwidth product described in [27,28,32–35]. Such a signal contains a defined (and user chosen) maximum spatial frequency, ${f}_{x}$, and a finite input extent. The real physical process of diffraction causes this signal to have a definite output spatial extent in the diffraction plane, which we later define as
*SE*. - (2) Calculating the diffraction distribution from a finite number of samples instead of the continuous case considered in point 1 fundamentally changes the nature of the problem. The numerically calculated signal will reproduce the analytical results from point 1 leading to an output field of extent
*SE*. In addition, however, an infinite set of replicas are generated in the diffraction plane. The properties of these replicas depend on how the calculation is performed. - (3) Two techniques for numerically calculating the diffraction integral are considered. Direct calculation of the diffraction integral produces replicas that are separated from each other by $\lambda z/\delta X$, where $\lambda $ is the wavelength, $z$ is the propagation distance, and $\delta X$ is the sampling interval in the input spatial domain. No spatial frequencies present in the signal in the input domain are removed by this calculation and spatial frequencies present that are not sampled according to the Nyquist criterion will nevertheless propagate and contribute to the field in the diffraction domain. The spectral-based calculation produces an infinite set of replicas in the diffraction plane that are separated from each other by $1/\delta v$, where $\delta v$ is sampling interval in the spatial frequency domain. The spatial frequency content of the diffracted signal is usually limited in the spatial frequency domain during the spectral calculation.
- (4) The sampling guidelines we propose follow from points 1–3 above:
*SE*must be less than either $1/\delta v$ or $\lambda z/\delta X$. When this is assured, the diffracted signal is considered well sampled. Again, we note that the approach taken above is different from many treatments in the literature since the behavior of the Fresnel chirp function does not feature in our considerations.

In Sections 2–5, we develop upon points 1–4 above. In Section 6, we turn our attention to numerically efficient algorithms for calculating the Fresnel transform which use the fast Fourier transform (FFT). Two algorithms are specifically addressed: the direct method (DM) and the spectral method (SM) which are fast (FFT-based) implementations of the direct and spectral calculations discussed in Sections 4 and 5. Our theoretical results are compared to results of other analysis presented in the literature. Finally, in Section 7, we discuss an important case where the sampling guidelines, provided in Section 4 and 5, are not appropriate: namely, when the diffraction field is sampled and we wish to calculate an inverse Fresnel transform on the result. This latter case has important implications for digital holography. We then summarize the main findings of this paper.

## 2. AN ANALYTICAL SOLUTION

In what follows, for notational simplicity, we will perform a 1D analysis. In Fig. 1, we sketch the nature of our problem. The field in Plane 1, denoted with uppercase letters $U(X)$, is assumed to be known and our task is to find an expression for the diffracted field, ${u}_{z}(x)$ in Plane 2. From the analysis presented in [4], we now write the formal definition of the Fresnel transform:

In some diffraction problems it is possible to solve the integral in Eq. (1) analytically. These solutions are important for several reasons: they provide significant insight into the diffraction process, and they also serve as a means of assessing the accuracy of a numerical approach. This can be achieved by comparing the results of a numerical calculation directly against those produced by the correct analytical solution. Some relevant properties of the Fresnel transform are discussed in [27,36–41]. We first examine the analytical solution of Eq. (1) for the following function:

which can also be expressed as*SE*, in Plane 2 of This represents our first important result. It indicates that a signal that has effectively finite support and a given spatial frequency component, ${f}_{x}$, will increase in extent as it propagates, due to two separate factors: ${f}_{x}$, and ${\alpha}_{z}$. We have chosen this particular analytical form because of the relative simplicity of the result which allows one to visualize how the power associated with specific spatial frequencies “walks-off” the optical axis under propagation, while remaining localized within the Gaussian envelope, see Section 18.6 of [42] and [43]. The higher the spatial frequency is, the more it “walks-off” the optical axis as it propagates. Although this is a solution for a relatively simple diffraction problem, its results can be extended to other more complex cases.

Nearly all diffraction problems that we wish to numerically evaluate have a finite extent in the input plane, similar to Eq. (2). These signals will in general have some spatial structure associated with them in this input plane. A problem arises, however, when we try to extend the results from Eq. (10) to these more general signals. This is because the result in Eq. (10) arises from a function, Eq. (2), that has an explicit spatial frequency component. How do we choose a suitable spatial frequency value, ${f}_{x}$, for a general signal? This question is examined in more detail in Section 5, where we show how a suitable value for ${f}_{x}$ can be chosen by examining the distribution of the signal’s power in its Fourier domain, see also [27,34,42,44]. We assume that if more than a specified amount (or percentage) of the power of the signal’s Fourier transform lies below some spatial frequency, this upper frequency value becomes ${f}_{x}$. The total power of the signal can be calculated in the spatial domain by integrating the intensity of the input wave field over the entire input plane.

## 3. INPUT TO THE NUMERICAL CALCULATION

In the following sections, we describe two different numerical approaches that can be used to calculate the Fresnel transform. We refer to these approaches as the direct calculation and the spectral calculation. These approaches are different from the fast numerical calculation techniques (DM and SM) discussed later in Section 6, since they do not make use of the FFT algorithm and the output spatial coordinate, $x$, can be chosen freely. The purpose of this section is to clearly state the input to our numerical calculation technique, which will take the form of two vectors: a spatial vector and a vector of complex values representing our signal. Once these values are chosen, we can then compare the relative performance of each algorithm.

To begin our numerical analysis, we assume that we have a vector of complex values, $\mathbf{U}$, that correspond to the actual values of the function, $U(X)$ at uniformly spaced sampling points (separated from each other by $\delta X$), defined with the spatial vector

where and where $N$ is the total number of samples. In the rest of the paper, we use bold-faced letters to indicate a vector. These two vectors, with a given $\lambda $ and $z$, we define the input parameters for the numerical calculation.## 4. DIRECT CALCULATION

Since we are now using only a finite number of samples to define the input field in Plane 1, the diffraction integral in Eq. (1) reduces to a finite summation

We now choose to examine ${u}_{z}^{S}(x)$ for a specific example

whereWe sample $U(X)$, taking $N=100$, over the range $-L\le X\le L$, where $L=0.1\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{mm}$, giving a $\delta X$ of 2 μm. The other parameters were chosen to be $z=1\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{cm}$, $\lambda =505.7\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{nm}$, and $\mathrm{\Gamma}=30\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{lines}/\mathrm{mm}$. In Figs. 2 and 3, we present the result of this calculation over the ranges $-SE/2\le X\le \lambda z/\delta X+0.5\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{mm}$, and $0\le X\le \delta X\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{mm}$, respectively. In Fig. 2, we present the magnitude distribution while the phase distribution is plotted in Fig. 3.

It is instructive to calculate *SE* for this numerical example. For this figure, we initially choose a value of ${f}_{x}=\mathrm{\Gamma}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{lines}/\mathrm{mm}$, however as we shall see, Eq. (14) contains power in spatial frequencies higher than $\mathrm{\Gamma}$ due to the rectangular function that limits the finite spatial extent of the signal. To calculate the value of ${\alpha}_{z}$, we need to relate ${\alpha}_{i}$ to the finite extent of the rectangular function defined in Eq. (15), i.e., we set $4{\alpha}_{i}=2L$. The result is plotted in Fig. 2. From Fig. 2, we can also see that there is a replica present, centered at the location $x=\lambda z/\delta X=2.5\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{mm}$. Comparing the phase of the central and first-order replicas, we find that the first-order phase distributions have an additional linear phase and constant phase. We note there are an infinite number of replicas in Plane 2 as a result of this calculation, which can be verified by simply plotting ${u}_{z}^{S}(X)$ over a wider range of $x$. The following relationship holds and can be derived following the approach outlined in Appendix A:

When the DC is used to calculate the Fresnel integral, all spatial frequencies propagate and contribute to the complex distribution in Plane 2, even those spatial frequencies that are sampled below the so-called Nyquist sampling rate. The presence of spatial frequencies at spatial locations greater than $SE/2$ is evident from Fig. 2, where an oscillating signal extends outward from $x\ge SE/2$ (distribution plotted in blue). This oscillation is due to spatial frequencies in the input signal that are higher than $\mathrm{\Gamma}$. A poor choice of ${f}_{x}=\mathrm{\Gamma}$ leads to a value for $SE/2$ (as presented in Fig. 2) that does not reflect the proper output extent of the diffracted signal. A proper choice for ${f}_{x}$ should be made on the basis of how the input signal’s power is distributed in its Fourier domain. In the next section, we examine in more detail how to choose an appropriate value for ${f}_{x}$ based on concepts from the space–bandwidth product of a signal.

It is possible, however, to limit the spatial frequency extent of ${u}_{z}^{S}(x)$ by subjecting it to a filtering operation in the Fourier domain and hence controlling *SE*. Here we remove all frequencies higher than $35\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{lines}/\mathrm{mm}$ [the continuous field ${u}_{z}^{S}(x)$ is sampled at a very high rate; these samples are Fourier transformed using an FFT operation, and the resulting spatial frequency distribution is multiplied by a rectangular aperture; this filtered spatial frequency distribution is then subjected to an inverse FFT operation]. In Fig. 2, the results of such a filtering operation can be seen by examining the black plot. This distribution has a much clearer finite spatial extent and is more in agreement with the value for *SE*, Eq. (10), calculated for a value of ${f}_{x}=\mathrm{\Gamma}=30\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{lines}/\mathrm{mm}$. The distribution (plotted in black) now has a spatial extent that agrees with the results presented in Section 2, however, the distribution is not as accurate when compared to the correct analytical solution. The filtering operation has removed some of the structural details evident in the replica distribution (blue) located at $x=2.5\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{mm}$.

Nevertheless, if the spatially filtered result (the distribution plotted in black) presented in Fig. 2 is considered satisfactory for a given application, we note that an almost identical result can be achieved by reducing the number of samples by a factor of 3.5, $\approx 28$ samples. Reducing the number of samples for a given input extent reduces the separation between the replicas; the spatial frequency filtering operation controls *SE*. In Section 5, we find that a sensible choice for ${f}_{x}$ is $60\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{lines}/\mathrm{mm}$ (and removing spatial frequencies higher than this) provides a reasonable balance of maintaining high accuracy while ensuring there is little cross talk between replicas.

We suggest the relation in Eq. (17) as a sampling guideline to be used when the spatial and spatial frequency extent of the input signal are known. When we use this sampling guideline for the specific example discussed, Eq. (14), we find that for a range of different values for $z$ and $\mathrm{\Gamma}$, an absolute minimum of about 98% of the signal’s power lies within $-SE/2<x<SE/2$, see Eq. (18), when ${f}_{x}=60\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{lines}/\mathrm{mm}$. Although the *SE* criterion has been specifically examined for Eq. (14), we have found that it holds for a range of different signal types, including random signals of finite spatial frequency extent.

An important point relating to the first sampling rule is that as $z$ becomes smaller, the replicas move closer together. This effect can be compensated for by increasing $N$ (while maintaining the same input extent of $2L$) and thereby reducing $\delta X$. However, after a certain point this calculation becomes extremely numerically intensive and hence unfeasible for practical calculations.

In closing this section we calculate the total power associated with the input signal, ${P}_{i}$, with the following expression

which for the numerical example presented here amounts to a value ${P}_{i}=0.0001$ arbitrary units at Plane 1.## 5. SPECTRAL CALCULATION

This approach makes use of the Fourier transform in order to perform the Fresnel transform. It is appropriate here to define the forward and inverse Fourier transform for some complex signal $f(x)$:

From [4], we find that Eq. (1) can also be written in the following manner:

Equations (21) and (1) are mathematically equivalent statements. In the following subsections, we examine the different steps that occur when we numerically implement Eq. (21) using a finite number of samples. To that end it, will be helpful to introduce a discrete form of the Fourier transform where again the superscript “$S$” indicates that this continuous function is generated from a finite set of complex data points, in our case from the vector $\mathbf{U}$, which we explicitly defined in Section 3:We note that the distribution described by Eq. (22) also contains an infinite number of replicas whose centers are separated from each other in the Fourier domain by a distance $\mathrm{\Delta}B=1/\delta X$. We shall refer to these replicas as Fourier replicas so as to distinguish them from the replicas discussed in Section 4.

#### A. Inverse Fourier Transform Operation

Replacing the function $FT\{U(X)\}(v)$ in Eq. (21) with Eq. (22) we get

In practice, however, the integration limits of the inverse Fourier transform in Eq. (23) are limited to a finite extent of the Fourier domain and to a finite number of samples. Usually when one implements the fast-Fourier-transform-based SM calculation (see the next section) the finite extent taken when calculating this integral is $\mathrm{\Delta}B$, which contains $N$ samples. It is precisely this limiting of the Fourier extent and the secondary sampling procedure that takes place here that causes the fast-Fourier-transform-based DM and SM calculations to have fundamentally different properties with regard to the output extent (DM-scaled output extent, while the SM remains constant as a function of propagation distance). We refer the reader to the following reference where the contributions of Fourier replicas to the output of the Fresnel transform are specifically examined [45].

#### B. Choosing the Fourier Extent: Power Considerations

Since ${\tilde{U}}^{S}(v)$ is actually a continuous function, we are free to choose any sampling interval, $\delta v$, and any Fourier domain extent. According to Parseval’s theorem, a Fourier transform conserves power, and we should expect the total power in the Fourier domain to be equal to ${P}_{i}$. A feature of the sampling operation is that now this power is contained within each replica window, i.e., within the range $-\mathrm{\Delta}B/2\le v\le \mathrm{\Delta}B/2$. We now observe that this latter statement is another way of saying that errors are introduced by the sampling operation.

A signal can have a finite extent in only one domain, see [36]. Since $U(X)$ is bounded, the power of its Fourier transform pair $\tilde{U}(v)$ must be distributed over the entire Fourier domain. This means that when a signal of finite extent is sampled (thereby creating “Fourier replicas” of infinite extent), it is inevitable that power from neighboring replicas (in the Fourier domain) leak into each other. And do so in such a manner that the total power in each replica window is ${P}_{i}$. This power leakage is in fact the aliasing effect and can be made arbitrarily small by increasing $\mathrm{\Delta}B$ (by increasing the number of samples $N$ for a constant input extent). Defining the spatial extent and a spatial frequency extent in this manner leads to the idea of a space–bandwith product. This concept has been discussed in some detail by other authors, see for example [27,28,34,35,44], and a region in the Fourier domain can be defined such that ${P}_{v}/{P}_{i}>R$, where

and where $R$ is some ratio, perhaps 0.9. As $FR$ in Eq. (24) increases and $R$ approaches 1, the distributions ${\tilde{U}}^{S}(v)$ and $\tilde{U}(v)$ become increasingly alike over the range $-\mathrm{\Delta}B/2\le v\le \mathrm{\Delta}B/2$, indicating a reduction in errors introduced by aliasing. From these considerations, we can define a region in the Fourier domain that contains $100R=90\%$ of the signal’s power.We can examine this question in more detail by returning to the numerical example discussed in Section 4. From these simulation values, we find that the distance between the centers of neighboring replicas in the Fourier domain is $\mathrm{\Delta}B=1/\delta X=500,000\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{lines}/\mathrm{m}$. When we examine how the signal’s power is distributed in the Fourier plane using Eq. (24), we find that $100R\approx 96\%$ of the signal’s power (in the Fourier domain) lies within the range $-60,000\le v\le 60,000\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{lines}/\mathrm{m}$. This indicates that we could achieve reasonable calculation accuracy using only that part of the signal that lies within the range $-60,000\le v\le 60,000\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{lines}/\mathrm{m}$. It is instructive to perform this calculation and compare the results with the DC depicted in Figs. 2 and 3. Reducing the contributing spatial frequencies can be used to reduce computational load and as we shall see in Section 6.3, can reduce the aliasing effect or cross talk between neighboring replicas by deliberately removing higher spatial frequencies and hence reducing *SE*.

#### C. Spectral Calculation: Sampling Guideline

In the previous subsection, we considered how to define an extent in the Fourier domain based on power considerations. We now would like to define a vector of complex values in the Fourier plane that describe our signal. Keeping $N=100$, we can define a spatial frequency step size $\delta v=FR/N$ and then a spatial frequency vector as

and using this vector we sample the continuous fieldComparing the terms within the curly brackets of Eqs. (13) and (27), we observe that these are very similar mathematical operations: (i) both take a finite vector of samples, (ii) multiply these by a chirp function, and (iii) then perform a Fourier transform operation on the result. This particular insight then allows us to develop a similar sampling rule for the SC approach that we have developed for the direct calculation, which can be summarized with the following equation:

where*SE*is given by Eq. (10).

This result indicates that there are fundamental differences between the replicas that arise in the direct and spectral calculations. The most important is that there is a different spacing between them given by ($\lambda z/\delta X$) and ($1/\delta v$) for the DC and SC cases, respectively.

## 6. FAST FOURIER TRANSFORM ALGORITHMS

In this section, we describe how the DC and the SC calculations discussed in the previous section can be implemented using FFT-based algorithms. Implementing the DC and SC calculations means that we retain control over the output spatial variable $x$ and are able to calculate the value of this continuous function for any given value of $x$. The downside to this flexibility is the length of time required for the calculation and hence the desire to somehow use the FFT algorithm. We note in advance that the intention here is to sketch one particular solution to the problem. We expect that it is possible to implement more numerically efficient techniques for achieving the same goal.

#### A. Zero-Padding: Increasing the Number of Samples in the Output Window

The FFT algorithm takes a vector of $N$ complex values and returns a vector of $N$ complex values. In Section 3 of this paper, we discussed the given input to the numerical problem we have at hand, defining both the spatial vector $\mathbf{X}$ and the associated set of complex values $\mathbf{U}$. The physical spacing between each sample is given by $\delta X$, and accordingly, the spatial frequency extent, i.e., the spacing between successive replicas in the spatial frequency domain, is given by $\mathrm{\Delta}B=1/(\delta X)$. We recall that in the spatial frequency domain, each of the samples are spatially separated from each other by $\delta v=\mathrm{\Delta}B/N$. We now zero-pad the vector $\mathbf{U}$ so that have a new vector

where the subscript “zp” refers to zero-padding.The vector ${\mathbf{U}}_{\mathbf{zp}}$ now has $M$ samples where $M>N$. This zero-padding operation does not change the value of $\delta X$, and hence the spatial frequency extent $\mathrm{\Delta}B$ also remains constant. We recall that the FFT algorithm maps $M$ samples to $M$ samples, with the effect that we have increasing the number of samples in the spatial frequency domain over the same $\mathrm{\Delta}B$. The spacing between samples in the spatial frequency domain is now given as $\delta {v}_{\mathrm{zp}}$, where $\delta {v}_{\mathrm{zp}}<\delta v$.

#### B. Inserting Zeros: Increasing the Extent of the Output Window

Here, we outline a different technique where we can increase the output extent of the FFT algorithm by inserting zeros in our input complex vector. This technique has been discussed in [9,46], and we refer the reader to Appendix B in [9] for more detail on the technique. Suffice to say that it is possible to define a vector

where the subscript “iz” refers to inserting zeros and the complex values ${U}_{1}$, ${U}_{2}$ are the individual elements of the vector $\mathbf{U}$. As we insert zeros into $\mathbf{U}$, we increase the output range of the FFT algorithm, however, we do not change the Fourier distribution in any manner. Now it is possible, however, to examine the distribution over a larger spatial frequency extent and to examine the properties of the resulting replicas as we shall see in the following numerical example.#### C. Another Numerical Example

In this section, we present the results of a second numerical calculation. We repeat the calculation presented by Voelz and Roggemann in Section 2.D of [30] and refer the reader to Figs. 3(a)–3(h) in their publication. These authors examine the diffraction pattern formed when a plane wave is incident on a 1D slit. In this section of the paper [30], two different algorithms are compared, which the authors term the impulse response (IR) method and the transfer function (TF) method. The TF approach there is the same as the SM technique discussed here. Their IR approach is, however, conceptually different from our DM technique in that the IR output has a similar extent to the SM.

We now consider Fig. 3(h) from [30], where the SM is shown to produce incorrect results reproducing the plot in Fig. 4(a). There are distinctive high-frequency spikes present in the distribution. The plot in Fig. 4(a) has a larger output window than that presented in Fig 2(h) of [30]. To increase the output window size in our numerical results, we have inserted zeros into the Fourier plane distribution (i.e., we insert zeros into the vector ${\tilde{\mathbf{U}}}^{\mathbf{S}}$ as outlined in Section 6.2 to give ${\tilde{\mathbf{U}}}_{\mathbf{iz}}^{\mathbf{S}}$), and we are now able to view the neighboring replicas. The simulation values indicate that $1/\delta v=0.5$, and hence that the first order replica for the SM technique should be located at this spatial location. This is indeed the case as can be observed from Fig. 4(a).

Now, using the zero-padding technique discussed in Section 6.1, we zero-pad the Fourier distribution ${\tilde{\mathbf{U}}}_{\mathbf{iz}}^{\mathbf{S}}$ (while keeping the inserted zeros present). Now, with a finer step size between neighboring samples, we can recognize that in fact higher frequencies are also present in the signal, as can be seen by comparing Figs. 4(a) and 4(b).

We wish to interpret the numerical results here in terms of the sampling guidelines that were suggested in Sections 4 and 5 of this paper. We first note that the DM produces the correct result for the SM calculation shown in Figs. 4(a) and 4(b) in this paper. If we set ${f}_{x}$ in Eq. (10) to $1/(2\delta X)$, i.e., the Nyquist sampling frequency, we find from Eq. (10) that $SE=5.2$, and so the physical extent of the signal lies mostly within the output sampling window $\lambda z/\delta X=5$. This is in contrast to the results of the SM approach, where $1/\delta v<SE$ and hence the neighboring replicas will overlap with each other causing the erroneous numerical artifacts that we see in Figs. 4(a) and 4(b).

We now wish to examine how we can control the SM calculation so that the correct diffraction distribution can be calculated. We begin by considering how we might make $\delta v$ smaller and hence increase the separation between the SM replicas. This can be achieved by zero-padding the input complex vector $\mathbf{U}$ to give ${\mathbf{U}}_{\mathbf{zp}}$ (we point out to the reader that this zero-padding operation occurs in the spatial domain). In Fig. 4(c), we have halved the value of $\delta v$ and the neighboring replica shifts and is no longer visible in the plot. There is still some overlap between the replicas and hence the residual distortions in the distribution. By further reducing the size of $\delta v$, the replicas can be further separated and these distortions can be made arbitrarily small.

Alternatively, we can change the physical spreading of each replica by removing higher spatial frequencies. This is what we present in Fig. 4(d), where we allow the Fourier distribution of the signal to be filtered by a rectangular aperture so that the maximum value for ${f}_{x}=1/(8T)$. We return the value of $\delta v$ to the value chosen in Figs. 4(a) and 4(b) and the replica again appears located at $x=0.5$. We note, however, that by removing the higher spatial frequencies, the cross talk or aliasing between neighboring replicas has been significantly reduced as we expect from Eq. (10).

## 7. CONCLUSION

In this paper, we have discussed how the Fresnel transform may be calculated using a finite number of samples. This topic has been examined by several other authors [25–31], where sampling guidelines were derived by considering the role of the quadratic phase factor or chirp function in the Fresnel kernel of the diffraction integral. Here, we have adopted a different approach to this problem.

The first observation we make is that an analytical signal of finite extent in the input plane and that contains specific spatial frequencies will increase in extent as the signal propagates. In Section 2, we examined a specific function, Eq. (2), which when inserted into the Fresnel diffraction integral can be integrated to give a known analytical solution. Using this analytical solution, we derived a relationship between the spatial extent (SE) of a diffracted signal, its finite extent, and spatial frequency content in the input plane. We then attempted to generalize these properties to other signals. In Section 3, we outlined the nature of our numerical calculation.

In Section 4, we found that using a finite number of samples to calculate the Fresnel transform produces an infinite set of replicas in the output plane that are separated from each other by a distance $\lambda z/\delta X$. In this direct calculation (DC), higher order replicas also have additional linear and constant phase terms associated with them. Provided that the SE of the signal that we wish to calculate is less than $\lambda z/\delta X$, we consider our signal to be well sampled.

In Section 5, we used a Fourier-transform-based approach to calculate the Fresnel transform, the SC. The nature of this calculation is different from the DC in that a filtering and resampling operation will, in practice, occur in the Fourier plane. These steps produce a different set of replicas in the diffraction plane that are separated from each other by $1/\delta v$. Higher order replicas do not have any additional phase terms associated with them. Again, our sampling guidelines indicate that SE should be less than the separation between the replicas.

In Section 6, we examined how to the DC and SC calculations could be implemented using FFT-based algorithms. Our numerical results were compared to other literature [30].

There is at least one occasion when the guideline outlined in Section 2 will not correctly predict the output extent of the propagated analytical solution. Since this situation occurs in practice in Fresnel-based digital holography systems, it is worthy of mention. In such systems, an object of interest is illuminated with a finite spot of coherent laser light. Let us assume that the finite illuminating spot is small, smaller than the extent of the CCD recording device. The scattered light propagates to the CCD plane where it is mixed with a known flat plane reference wave and the intensity is recorded. This is, in fact, a digital hologram. If phase shifting interferometric techniques are used (and at least three phase-shifted holograms recorded) it is possible to separate the complex diffracted field ${u}_{z}(x)$ from the DC terms and twin image. When a filtered and sampled version of ${u}_{z}(x)$ is back propagated numerically, one finds that the extent of the signal does not increase but rather reduces in extent approaching the size of the finite illuminating spot, see for example [9,36,37,43,47–49]. Although this may seem obvious, it does not arise from the considerations presented in Section 2. It also has some important consequences, namely that it is possible to sample the field in the diffraction plane at a rate far lower than the Nyquist limit and yet still be able to recover spatial frequencies higher than the sampling rate of the CCD digital array. The reconstruction of these sub-Nyquist sampled signals by inverse Fresnel propagation has been discussed in [45], where it was found that the DC is the optimal technique for reconstruction. This is related to the secondary filtering and resampling operations that occur in the Fourier plane when the SC is being implemented.

Although we have concentrated on representing our signal using a finite set of data points, other techniques have been shown to be useful for modeling propagation [50–53]. It is often of interest to extend the results here to the cases of pulse diffraction, see [54–56].

## APPENDIX A

Let us assume that we have a continuous function $U(X)$. Its Fresnel transform is denoted ${u}_{z}(x)$, as defined in Eq. (1), i.e.,

We now wish to perform a Fresnel transform on a sampled version of $U(X)$. To aid us in this task, we use a Dirac delta train where $m$ is an integer and $\delta $ is the Dirac delta function [57]. This can also be expressed as [58]## ACKNOWLEDGMENTS

DPK is now Junior-Stiftungs professor of Optics Design and is supported by funding from the Carl-Zeiss-Stiftung (FKZ: 21-0563-2.8/121/1). DPK is very grateful to Thomas Meinecke for helpful discussions regarding numerical propagation of speckle fields.

## REFERENCES

**1. **M. Hrynevych, “Diffraction effects in Michelson stellar interferometry,” Ph.D. thesis (University of Sydney, 1992).

**2. **W. H. Southwell, “Validity of the Fresnel approximation in the near field,” J. Opt. Soc. Am. **71**, 7–14 (1981). [CrossRef]

**3. **G. W. Forbes, “Validity of the Fresnel approximation in the diffraction of collimated beams,” J. Opt. Soc. Am. A **13**, 1816–1826 (1996). [CrossRef]

**4. **J. Goodman, *Introduction to Fourier Optics*, 2nd ed. (McGraw-Hill, 1966).

**5. **M. Born and E. Wolf, *Principles of Optics: Electromagnetic Theory of Propagation, Interference and Diffraction of Light*, 7th ed. (Cambridge University, 1999).

**6. **G. Pedrini, P. Fröning, H. J. Tiziani, and M. E. Gusev, “Pulsed digital holography for high-speed contouring that uses a two-wavelength method,” Appl. Opt. **38**, 3460–3467 (1999). [CrossRef]

**7. **P. Picart, J. Leval, D. Mounier, and S. Gougeon, “Some opportunities for vibration analysis with time averaging in digital Fresnel holography,” Appl. Opt. **44**, 337–343 (2005). [CrossRef]

**8. **P. Picart and J. Leval, “General theoretical formulation of image formation in digital Fresnel holography,” J. Opt. Soc. Am. A **25**, 1744–1761 (2008). [CrossRef]

**9. **D. P. Kelly, B. M. Hennelly, N. Pandey, T. J. Naughton, and W. T. Rhodes, “Resolution limits in practical digital holographic systems,” Opt. Eng. **48**, 095801 (2009). [CrossRef]

**10. **K. M. Molony, B. M. Hennelly, D. P. Kelly, and T. J. Naughton, “Reconstruction algorithms applied to in-line Gabor digital holographic microscopy,” Opt. Commun. **283**, 903–909 (2010). [CrossRef]

**11. **D. S. Monaghan, D. P. Kelly, N. Pandey, and B. M. Hennelly, “Twin removal in digital holography using diffuse illumination,” Opt. Lett. **34**, 3610–3612 (2009). [CrossRef]

**12. **T. Meinecke, N. Sabitov, and S. Sinzinger, “Information extraction from digital holograms for particle flow analysis,” Appl. Opt. **49**, 2446–2455 (2010). [CrossRef]

**13. **P. F. Almoro, P. N. Gundu, and S. G. Hanson, “Numerical correction of aberrations via phase retrieval with speckle illumination,” Opt. Lett. **34**, 521–523 (2009). [CrossRef]

**14. **J. R. Fienup, “Phase retrieval algorithms: a comparison,” Appl. Opt. **21**, 2758–2769 (1982). [CrossRef]

**15. **M. Guizar-Sicairos and J. R. Fienup, “Understanding the twin-image problem in phase retrieval,” J. Opt. Soc. Am. A **29**, 2367–2375 (2012). [CrossRef]

**16. **M. R. Teague, “Deterministic phase retrieval: a Green’s function solution,” J. Opt. Soc. Am. **73**, 1434–1441 (1983). [CrossRef]

**17. **C. J. R. Sheppard, “Three-dimensional phase imaging with the intensity transport equation,” Appl. Opt. **41**, 5951–5955 (2002). [CrossRef]

**18. **V. Arrizón, M. Testorf, S. Sinzinger, and J. Jahns, “Iterative optimization of phase-only diffractive optical elements based on a lenslet array,” J. Opt. Soc. Am. A **17**, 2157–2164 (2000). [CrossRef]

**19. **R. Kleindienst, L. Moeller, and S. Sinzinger, “Highly efficient refractive Gaussian-to-tophat beam shaper for compact terahertz imager,” Appl. Opt. **49**, 1757–1763 (2010). [CrossRef]

**20. **I. Eriksson, P. Haglund, J. Powell, M. Sjodahl, and A. F. H. Kaplan, “Holographic measurement of thermal distortion during laser spot welding,” Opt. Eng. **51**, 030501 (2012). [CrossRef]

**21. **O. V. Angelsky, A. P. Maksimyak, P. P. Maksimyak, and S. G. Hanson, “Optical correlation diagnostics of rough surfaces with large surface inhomogeneities,” Opt. Express **14**, 7299–7311 (2006). [CrossRef]

**22. **J. W. Goodman, *Speckle Phenomena in Optics* (Roberts and Company, 2007).

**23. **D. P. Kelly, B. M. Hennelly, and J. T. Sheridan, “Magnitude and direction of motion with speckle correlation and the optical fractional fourier transform,” Appl. Opt. **44**, 2720–2727 (2005). [CrossRef]

**24. **D. P. Kelly, J. E. Ward, U. Gopinathan, B. M. Hennelly, F. T. O’Neill, and J. T. Sheridan, “Generalized Yamaguchi correlation factor for coherent quadratic phase speckle metrology systems with an aperture,” Opt. Lett. **31**, 3444–3446 (2006). [CrossRef]

**25. **D. Mas, J. Garcia, C. Ferreira, L. M. Bernardo, and F. Marinho, “Fast algorithms for free-space diffraction patterns calculation,” Opt. Commun. **164**, 233–245 (1999). [CrossRef]

**26. **D. Mas, J. Garcia, C. Ferreira, L. M. Bernardo, and F. Marinho, “Fast numerical calculation of Fresnel patterns in convergent systems,” Opt. Commun. **227**, 245–258 (2003). [CrossRef]

**27. **B. M. Hennelly and J. T. Sheridan, “Generalizing, optimizing, and inventing numerical algorithms for the fractional Fourier, Fresnel, and linear canonical transforms,” J. Opt. Soc. Am. A **22**, 917–927 (2005). [CrossRef]

**28. **D. P. Kelly, B. M. Hennelly, W. T. Rhodes, and J. T. Sheridan, “Analytical and numerical analysis of linear optical systems,” Opt. Eng. **45**, 088201 (2006). [CrossRef]

**29. **F. Shen and A. Wang, “Fast-Fourier-transform based numerical integration method for the Rayleigh–Sommerfeld diffraction formula,” Appl. Opt. **45**, 1102–1110 (2006). [CrossRef]

**30. **D. G. Voelz and M. C. Roggemann, “Digital simulation of scalar optical diffraction: revisiting chirp function sampling criteria and consequences,” Appl. Opt. **48**, 6132–6142 (2009). [CrossRef]

**31. **T. Kozacki, K. Falaggis, and M. Kujawinska, “Computation of diffracted fields for the case of high numerical aperture using the angular spectrum method,” Appl. Opt. **51**, 7080–7088 (2012). [CrossRef]

**32. **A. W. Lohmann, “Image rotation, Wigner rotation, and the fractional Fourier transform,” J. Opt. Soc. Am. A **10**, 2181–2186 (1993). [CrossRef]

**33. **A. W. Lohmann, R. G. Dorsch, D. Mendlovic, Z. Zalevsky, and C. Ferreira, “Space–bandwidth product of optical signals and systems,” J. Opt. Soc. Am. A **13**, 470–473 (1996). [CrossRef]

**34. **D. Mendlovic and A. W. Lohmann, “Space–bandwidth product adaptation and its application to superresolution: fundamentals,” J. Opt. Soc. Am. A **14**, 558–562 (1997). [CrossRef]

**35. **D. Mendlovic, A. W. Lohmann, and Z. Zalevsky, “Space–bandwidth product adaptation and its application to superresolution: examples,” J. Opt. Soc. Am. A **14**, 563–567 (1997). [CrossRef]

**36. **F. Gori, “Fresnel transform and sampling theorem,” Opt. Commun. **39**, 293–297 (1981). [CrossRef]

**37. **L. Onural, “Sampling of the diffraction field,” Appl. Opt. **39**, 5929–5935 (2000). [CrossRef]

**38. **A. Stern and B. Javidi, “Analysis of practical sampling and reconstruction from Fresnel fields,” Opt. Eng. **43**, 239–250 (2004). [CrossRef]

**39. **A. Sokołowski, “Consequences of sampling an image and Fourier planes in a numerical light propagation model based on the Helmholtz–Kirchhoff approximation,” J. Opt. Soc. Am. A **23**, 2764–2767 (2006). [CrossRef]

**40. **L. Onural, “Exact analysis of the effects of sampling of the scalar diffraction field,” J. Opt. Soc. Am. A **24**, 359–367 (2007). [CrossRef]

**41. **A. Sokołowski and T. Więcek, “Consequences of sampling an image and Fourier planes in a numerical light propagation model based on the Helmholtz–Kirchhoff approximation. ii. comparison between two numerical algorithms,” J. Opt. Soc. Am. A **27**, 1688–1693 (2010). [CrossRef]

**42. **A. W. Lohmann, *Optical Information Processing* (Universitätsverlag Ilmenau, 2006).

**43. **D. P. Kelly, J. J. Healy, and J. T. S. B. M. Hennelly, “Quantifying the 2.5d imaging performance of digital holographic systems,” J. Eur. Opt. Soc. Rapid Pub. **6**, 11034 (2011). [CrossRef]

**44. **A. W. Lohmann, R. G. Dorsch, D. Mendolovic, Z. Zalevsky, and C. Ferreira, “Space–bandwidth product of optical signals and systems,” J. Opt. Soc. Am. A **13**, 470–473 (1996). [CrossRef]

**45. **D. P. Kelly, N. Sabitov, T. Meinecke, and S. Sinzinger, “Some considerations when numerically calculating diffraction patterns,” in *Digital Holography and Three-Dimensional Imaging* (Optical Society of America, 2011), p. DTuC5.

**46. **A. Stern and B. Javidi, “Improved-resolution digital holography using the generalized sampling theorem for locally band-limited fields,” J. Opt. Soc. Am. A **23**, 1227–1235 (2006). [CrossRef]

**47. **A. Stern and B. Javidi, “Sampling in the light of wigner distribution,” J. Opt. Soc. Am. A **21**, 360–366 (2004). [CrossRef]

**48. **A. Stern and B. Javidi, “Sampling in the light of Wigner distribution: errata,” J. Opt. Soc. Am. A **21**, 2038 (2004). [CrossRef]

**49. **D. P. Kelly and D. Claus, “Filtering role of the sensor pixel in Fourier and Fresnel digital holography,” Appl. Opt. **52**, A336–A345 (2013). [CrossRef]

**50. **L. Onural, “Diffraction from a wavelet point of view,” Opt. Lett. **18**, 846–848 (1993). [CrossRef]

**51. **M. Liebling, T. Blu, and M. Unser, “Fresnelets: new multiresolution wavelet bases for digital holography,” IEEE Trans. Image Process. **12**, 29–43 (2003). [CrossRef]

**52. **M. Liebling, “On Fresnelets, interference fringes, and digital holography,” Ph.D. thesis (Swiss Federal Institute of Technology Lausanne, 2004).

**53. **N. Chacko, M. Liebling, and T. Blu, “Discretization of continuous convolution operators for accurate modeling of wave propagation in digital holography,” J. Opt. Soc. Am. A **30**, 2012–2020 (2013). [CrossRef]

**54. **M. Gu and X. S. Gan, “Fresnel diffraction by circular and serrated apertures illuminated with an ultrashort pulsed-laser beam,” J. Opt. Soc. Am. A **13**, 771–778 (1996). [CrossRef]

**55. **D. P. Kelly, B. M. Hennelly, A. Grün, and K. Unterrainer, “Numerical sampling rules for paraxial regime pulse diffraction calculations,” J. Opt. Soc. Am. A **25**, 2299–2308 (2008). [CrossRef]

**56. **R. J. Mahon and J. A. Murphy, “Diffraction of an optical pulse as an expansion in ultrashort orthogonal Gaussian beam modes,” J. Opt. Soc. Am. A **30**, 215–226 (2013). [CrossRef]

**57. **R. Bracewell, *The Fourier Transform and its Applications* (McGraw-Hill, 1965).

**58. **A. Stern, “Sampling of linear canonical transformed signals,” Signal Process. **86**, 1421–1425 (2006). [CrossRef]