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.

 

Fig. 1. Illustration of the diffraction process. An infinite plane wave is incident on the aperture in Plane 1, denoted by the spatial variable X. Consider how to calculate the complex amplitude at $P(x,z)$. Only the light that lies within the aperture opening can contribute to the complex amplitude at $P(x,z)$. We envisage the light within the aperture opening as consisting of an infinite number of secondary point sources (PS), each with an intensity and phase value (only a few PS are shown here). Since the aperture is illuminated with an ideal plane wave, each secondary PS will have an identical intensity and phase value in Plane 1. Nevertheless, each PS will accumulate a different phase value $P(x,z)$ due to the different path lengths, compared to the distance $l_1$ and $l_4$.

Download Full Size | PDF

The argument made in this paper rests on the following 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:

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

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

We sample $U(X)$, taking $N=100$, over the range $-L\le X\le L$, where $L=0.1\mathrm{mm}$, giving a $\delta X$ of 2 μm. The other parameters were chosen to be $z=1\mathrm{cm}$, $\lambda =505.7\mathrm{nm}$, and $\mathrm{\Gamma }=30\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\mathrm{mm}$, and $0\le X\le \delta X\mathrm{mm}$, respectively. In Fig. 2, we present the magnitude distribution while the phase distribution is plotted in Fig. 3.

 

Fig. 2. Result of the diffraction calculation. The presence of a replica is clearly visible at $x=2.5\mathrm{mm}$. The red dotted line indicates the location of $SE/2$, illustrating the extent of the diffracted field. The distribution plotted in black has been spatially filtered and contains no spatial frequencies higher than $35\mathrm{lines}/\mathrm{mm}$. Note the finite extent of the signal in comparison to the blue plot. (For presentation purposes, the central order replica in blue is not plotted over the range $-SE/2\le x\le SE/2$.)

Download Full Size | PDF

 

Fig. 3. Result of the diffraction calculation for the phase distributions. Both plots are over a range of $\delta X$. The plot on the left is the phase of the zero-order replica (offset to zero). The plot on the right is for the first replica, where there is a constant phase offset of $\pi /4$ and a linear phase with slope $1/\delta X$.

Download Full Size | PDF

It is instructive to calculate SE for this numerical example. For this figure, we initially choose a value of $f_x=\mathrm{\Gamma }\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\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\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\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\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\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\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

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:

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

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\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\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\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

Comparing 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:

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

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

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).

 

Fig. 4. Reproducing the numerical results from Voelz and Roggemann [30].

Download Full Size | PDF

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].