The numerical implementation of Fresnel diffraction is widely used in optics [1], especially in holography [2,3]. It can efficiently calculate how the wavefield propagates over space, forming the backbone of many applications in computer-generated holography [4] and holographic microscopy [5].

Currently, numerical Fresnel diffraction can be said to come in two main forms. The spatial domain form computes Fresnel diffraction directly as a spatial convolution. The first method is highly flexible and efficient for small distances, used for techniques such as lookup table methods [6,7], layer-based methods [8], and methods related to the wavefront recording plane method [9,10]. However, it does not scale well for larger propagation distances and may exhibit some artifacts as a result of the finite-sized zone plate.

The other main form uses the frequency domain. This computes diffraction as a convolution via Fourier space, which is more efficient than the spatial domain approach for larger distances. However, it is inherently less flexible, often only working for certain propagation distance ranges, needs zero padding, can result in aliasing for some distances, may need multiple fast Fourier transforms (FFTs), and requires extensions whenever there are differences in resolution, pixel pitch, or lateral position between the source and destination planes. That is why they come in many variants, such as shifted Fresnel diffraction [11] or aliasing-reduced Fresnel diffraction [12]. Recent approaches also support more flexible angular spectrum methods, including efficient off-axis diffraction algorithms [13,14], scalable angular spectrum models [15], and the homeomorphic Fourier transform [16], but do not resolve all of the aforementioned limitations.

We propose a novel approach for computing Fresnel diffraction using Gabor frames. Because the point-spread function (PSF) is sparse in Gabor space, we can achieve linear computational complexity independently of propagation distance for any signal [17]. The transfer function exhibits both properties from spatial and frequency domain methods. Because of the analytical derivations, the method has inherently no limitations w.r.t. the parameters of the source and destination planes; any combination of pixel pitches, resolutions, and propagation distances is supported. The method requires no zero-padding. Because of its phase-space-based modeling, it is more flexible for integration with computer-generated holography (CGH) algorithms relying on, e.g., ray tracing [18] or ray-wavefront conversion [19].

This summarizes our specific contributions: •we derive an analytical model for Fresnel diffraction with chirplets and Gabor frames;•we utilize this to devise an algorithm for calculating accurate discrete Fresnel diffraction with freely chosen parameters;•we propose an efficient approximation for high-speed algorithm evaluation, leading to a linear complexity algorithm;•the method is implemented on GPU with C++/CUDA and compared to conventional algorithms;•multiple experiments are performed, analyzing the PSF shape and measuring the computation speed and numerical accuracy.

The remainder of this paper is organized as follows. In Section 2, we introduce the theory for numerical diffraction, phase space, and frames and derive an analytical model using chirplets for Fresnel diffraction with Gabor frames. Section 3 details the proposed algorithm and details an efficient approximation for fast implementations. The algorithm is evaluated in Section 4, covering the PSF analysis, comparing the accuracy with different Fresnel diffraction algorithms, and for different parameters. Finally, we conclude and discuss future work in Section 5.

The Fresnel diffraction operator expresses the propagation of a wavefield from a source plane z=z1 to a parallel destination plane z=z2. The continuous Fresnel integral is given by [1] $$\begin{gathered} U(x,y;z2)=eikdikd∬−∞+∞U(ξ,η;z1) \\ ·exp(ik2d((x−ξ)2+(y−η)2))dξdη, \end{gathered}$$(1)where i is the imaginary unit, d=z2−z1, k=2πλ is the wavenumber with wavelength λ, and U is the complex-valued wavefield. Note that the Fresnel diffraction operator is separable along x and y. This means we can only consider the operation along x, giving us $$\begin{gathered} U(x;z2)=∫−∞+∞U(ξ;z1)exp(ik2d(x−ξ)2)dξ \\ =∫−∞+∞U(ξ;z1)K(x−ξ)dξ, \end{gathered}$$(2)where K(·) is the Fresnel convolution kernel. We ignore the leading constant for simplicity, which does not affect our calculations.

This integral can be computed numerically in the spatial domain by means of a spatial convolution. However, one should not discretize the integral expression in a straightforward way, since this would lead to aliasing. K(·) has a well-defined unique instantaneous frequency ν in every point (provided the assumption of stationary instantaneous phase [20]), given by ν(x)=12π∂∂x∠(K(x))=x2λd,(3)where ∠(·) is the phase function (or angle, “complex” argument). This instantaneous frequency may exceed the Nyquist rate of the sampled signal, especially for short propagation distances d. This can be addressed by only evaluating K(·) in the region where |ν(x)|≤12p, i.e., by staying below the Nyquist limit, and setting it to zero otherwise. Note that p is the pixel pitch of the hologram. See Ref. [21] for a more general discussion on the shape of the masking region. Formally, we can use $$\begin{gathered} U(x;z2)=∫−∞+∞U(ξ;z1)Π(2pλd(x−ξ))K(x−ξ)dξ \\ =U(x;z1)*Π(2pxλd)K(x), \end{gathered}$$(4)where Π(·) is the rectangular function and * is the convolution operator. The spatial convolution algorithm has an O(n2) complexity and is thus generally unsuitable for efficiently modeling diffraction. Although the complexity is generally high, the PSF support will be limited for short distances d, leading to fast run times that can surpass frequency domain algorithms. This has been leveraged, e.g., in Ref. [9]. Conversely, a very small d will result in a PSF with a tiny number of samples, leading to larger errors. In addition, the sharp transition to zero will cause ringing artifacts.

The second main class of techniques computes the Fresnel transform in the frequency domain with the FFT. For this, the Fresnel kernel is expressed analytically in Fourier space as K^(ξ)=F{K(x)}=iπdkexp(idξ22k)(5)using the Fourier transform F{·} so that we can compute the convolution as U(x;z2)=F−1{F{U(x;z1)}K^(ξ)}.(6)The main advantage of this approach is the O(nlogn) computational complexity, independently of the parameters of K(·). This form is mainly used for long-distance propagation with pixel pitch preservation, i.e., the pixel pitches of the source and destination planes match. It is one of the most widely used approaches for computing Fresnel diffraction. Other variants exist using one or three FFTs, where the pixel pitch may change as a function of d [22].

However, several considerations must be made when computing the discrete version of the Fresnel diffraction operator in the frequency domain. Aliasing can occur for large distances d and zero padding is required. The frequency domain approach is less flexible than the spatial domain convolution whenever there are differences in the sampling pitch p and the number of samples between the source and destination planes or whenever there is a lateral offset. To this end, several extensions have been proposed to (partially) address these limitations, as mentioned in the introduction [11,12,23]. Multiple shifted propagations can be used whenever the destination plane resolution is a whole multiple of the source plane, but this cannot efficiently calculate arbitrary resolution changes.

We propose a third way for encoding Fresnel diffraction, relying on a joint space-frequency (or time-frequency) representation, also known as the Wigner transform or phase space.

One important way of finding a suitable signal representation is the notion of sparsity; this is a concept in signal processing, where a small number of coefficients can accurately represent some class of signals in a well-chosen transform space. This property is useful for many applications, such as data compression, filtering, compressed sensing, and accelerating calculations. Classical signal analysis using the Fourier transform for analyzing and representing many kinds of generally non-stationary holograms tends to be suboptimal [24]. This fact is apparent when looking at the PSF in Fresnel diffraction, cf. Fig. 1(a): although the PSF can be dense simultaneously in both spatial x and frequency ν domains, it is sparse in phase space.

We can generally model many paraxial diffraction operations as linear canonical transforms [25]. Many of these operations can be expressed as linear phase space transforms, including the Fourier transform (as a 90° rotation), fractional Fourier transform (arbitrary rotation), magnification (area-preserving coordinate scaling), and, in particular, Fresnel diffraction (shearing). For a mathematically rigorous description of, among others, reversible phase space transforms, we refer to Refs. [26–29].

The shearing in phase space [cf. Fig. 1(b)] shows how Fresnel diffraction is a very structured convolution operation amenable to efficient modeling and processing beyond a pure Fourier representation. One common reversible phase space representation is the Wigner transform. As it is quadratic with respect to the signal, it introduces cross-terms for any signal that does not form a single straight line in phase space [27,28,30]—such as holograms, which typically consist of a superposition of millions of PSFs, encoded in phase space as lines or sigmoids [20]. Furthermore, it is often highly overcomplete (typically requiring N2 coefficients for an N-sample input signal). It is thus not a suitable candidate to leverage the phase space sparsity of holograms.

Therefore, we aim to select a better representation suitable for efficiently and accurately modeling Fresnel diffraction for any combination of propagation distances, lateral offsets, sampling pitches, and wavelengths. This requires a representation satisfying several requirements: •locality or phase space compactness, to have most of the signal energy concentrated over a small area in phase space, promoting sparsity for arbitrary diffraction conditions;•efficiently “tiling” phase space, representing any complicated signal effectively with as few elements as possible;•efficient computation of the forward and inverse transforms;•being analytically and numerically compatible with Fresnel diffraction, i.e., preserving the aforementioned properties after Fresnel diffraction over any distance;•linearity with respect to the signal, to avoid cross-terms.

We argue that chirped Gabor frames are good candidates for modeling Fresnel diffraction, which we will introduce in the next section. This is because, locally, any hologram in the Fresnel diffraction regime can be well approximated using chirplets.

We define a chirplet as a function composing the exponential function with a quadratic function with complex-valued coefficients. Formally, we define the set of chirplets G as G={exp(αt2+βt+γ)|α,β,γ∈C∧ℜ(α)<0},(7)where we consider elements where the real part ℜ(α)<0 so that they are square-integrable and tending to zero for t→±∞. They can be viewed as an extension of Gaussians, where the imaginary part of the coefficients ℑ(α), ℑ(β), and ℑ(γ) are additional degrees of freedom, cf. Table 1.Table 1.

Relation between Chirplet Parameters and Attributes

Chirplets can be viewed as idealized chirped pulses with a Gaussian envelope, cf. Fig. 2.

Because chirplets are an extension of Gaussians, they inherit the optimality of Gaussian beam modes in that they are shape-preserving, have a well-defined Fourier representation, and spread minimally as they propagate. They can be viewed as a wave-like extension of a ray bundle, similar to field tracing [31], which are efficient generalizations for ray tracing in geometric optics for numerical diffraction.

The set of chirplets G is closed under multiplications and convolutions. Consider two chirplets u=exp(at2+bt+c), u^=exp(a^t2+b^t+c^); we have that ∀u,u^∈G, u·u^=exp((a+a^)t2+(b+b^)t+(c+c^))∈G,(8)and $$\begin{gathered} u*u^=∫−∞+∞u(τ)u^(t−τ)dτ \\ =−πa+a^exp(aa^a+a^t2+ab^+a^ba+a^t−(b−b^)24(a+a^)+c+c^) \end{gathered}$$(9)are also part of G. In addition, chirplets can be integrated over R analytically as follows: ∀u∈G:∫−∞+∞u(t)dt=π−aexp(c−b24a).(10)This means we can freely multiply, convolve, or compute inner products with chirplets, making them suitable basis functions.

The Fresnel convolution kernel K can be viewed as a “degenerate” chirplet since it almost satisfies the condition of belonging to G, except for having a purely imaginary α, i.e., ℜ(α)=0 rather than ℜ(α)<0. This poses no problem in practice because multiplying or convolving any proper chirplet with K will still result in a chirplet. Formally, ∀u∈G:u·K∈G∧u*K∈G,(11)following immediately from Eqs. (8) and (9), seamlessly integrating Fresnel diffraction with chirplets. For notational simplicity, we will define K=exp(iKx2), where K=πp2λd is a dimensionless constant closely related to the Fresnel number.

Suppose we would like to express the relationship between the signals in the parallel input and output planes, spaced by a distance d, respectively represented with chirplets. For any u=exp(at2+bt+c) in the input plane and u^=exp(a^t2+b^t+c^) the output plane; we can now directly express the relationship of their coefficients with the inner product $$\begin{gathered} ⟨u^,u*K⟩=∫−∞+∞u^(t)·(u*K)(t)¯dt \\ =πaa^+K(a+a^)exp(−14ab^2+a^b2+(b−b^)2aa^+K(a+a^)+c+c^), \end{gathered}$$(12)where the overbar is the complex conjugate operation. In addition to this mapping, one should be able to represent arbitrary signals with a collection of chirplets. This can be achieved with Gabor frames, which provide a mathematical framework for a reversible, linear, and local transform.

Frames generalize vector space bases to sets that may be linearly dependent, providing a redundant and stable representation of a signal [29,32]. Although an “overcomplete” basis means that the number of transform coefficients will be larger than the number of input samples to ensure reversibility, this allows for more flexibility in designing the transform properties. Formally, a family of vectors {φj}j=1M in a separable Hilbert space H is called a frame for H, if there exist constants 0<A≤B<∞ such that [32] ∀h∈H,A||h||2≤∑j=1M|⟨h,φj⟩|2≤B||h||2.(13)For every frame {φj}j=1M, there exists a canonical dual frame {φ^j}j=1M, satisfying ∀h∈H,h=∑j=1M⟨h,φj⟩φ^j=∑j=1M⟨h,φ^j⟩φj,(14)namely, serving as an inverse operator.

We will focus on a particular class of frames: stationary Gabor systems. Such a Gabor frame can be represented as a lattice (m,n)∈Z2 in phase space, expressing a collection of translations T and modulations M of a window g: {MmaTnbg}m,n∈Z,(15)where a,b>0 are fixed, and the operators Tp and Mq operating on some f∈H are defined as Tp:(Tpf)(x)=f(x−p)andMq:(Mqf)(x)=e2πiqxf(x).(16)Many windows g are admissible, but they usually are chosen to be symmetric, with a maximum in the middle and tapering away from the center, e.g., rectangular windows, B-spline windows, or Hann windows, among many others. We chose a Gaussian window g to be compatible with chirplets. Only values for a and b where ab≤1 are admissible. Otherwise, the phase space lattice becomes too “sparse,” making inversion impossible. As the lattice gets denser for smaller values of ab, the dual window g^ for the canonical dual frame gets increasingly compact and regular; see Fig. 3. This causes a trade-off between the number of basis functions and their efficient computability. In practice, oversampling factors of ≈43 (or ab=34=0.75) provide sufficiently regular dual windows, which is superior to the 2× oversampling factor required for the zero-padding in frequency-domain-based Fresnel diffraction. In other words, the memory requirements of the proposed Gabor diffraction algorithm are lower than that of conventional Fresnel diffraction, only needing an ≈33% increase of coefficients in the intermediate representation as opposed to the memory doubling required for zero-padding for frequency domain Fresnel diffraction.

Since Gaussian Gabor frame elements are special cases of chirplets [with ℑ(α)=0], we have all the ingredients to express Fresnel diffraction directly via the Gabor transform instead of the spatial or frequency domain. The practical methodology will be explained in the next section.

The methodology for Fresnel diffraction of Gabor frames is summarized in Fig. 4. The general approach is similar in spirit to the frequency domain Fresnel transform, but instead of going back and forth with FFTs and applying operations in Fourier space, we use the Gabor transform.

First, we can independently process the rows and columns of the input hologram since the Fresnel transform is separable along x and y. This can also be done in parallel to accelerate the calculations. Every row and column are transformed into an invertible linear combination of chirplets; this is achieved by using the Gabor transform with the dual frame. Then, these chirplets are transformed using the analytical mapping expressed in Eq. (12); how this can be achieved in practice is elaborated in Section 3.B. The resulting Gabor coefficients now describe the hologram in the destination plane. We obtain the output hologram after a final inverse Gabor transform using the dual frame again.

We use the dual frame twice in the forward and inverse Gabor transforms. This is because the analytical mapping is expressed using chirplets with a Gaussian window, requiring initial coefficients obtained with the dual frame rather than the standard frame. This means we can use the same (precomputed) dual window for both transform parts.

The Gabor transform can be computed with a freely selectable number of modulations, independent of the sample count n, making it a transform of linear complexity. What remains to be done is to find a linear-complexity implementation of the analytical chirplet mapping given by Eq. (12), resulting in O(n) complexity Fresnel diffraction.

Although the mapping would indicate at first glance a relationship exists between all input and output coefficients, this mapping transformation is sparse: the large majority of mappings will be near-zero and thus negligible. Note that this is independent of the input signal, which does not have to be sparse. Suppose that we have u=exp(−s(x−δ)2+imx) as the input chirplet and u^=exp(−s^(x−δ^)2+im^x) as the output chirplet. Namely, the support is encoded by (s,s^), the translation by (δ,δ^), and the modulation by (m,m^). From Eq. (12), we can see that the inner product is characterized by an exponential function of the chirplet parameters. The magnitude will thus only be affected by the real part of the expression within the exponential, that is, $$\begin{gathered} ℜ(log⟨u^,u*K⟩)∝E(δ^,m^)≔(4K2(s+s^)2+4s2s^2)−1 \\ ·(K2(s+s^)(4ss^(δ−δ^)2+(m−m^)2) \\ +ss^(4K(ms^+m^s)(δ−δ^)+m2s^+m^2s)), \end{gathered}$$(17)as we can ignore the leading constant since it will remain constant for varying δ^ or m^. The function E(δ^,m^) is concave everywhere since we have d2dδ^2E=−2K2ss^(s+s^)K2(s+s^)2+s2s^2<0(18)and d2dm^2E=−K2(s+s^)+ss^22K2(s+s^)2+2s2s^2<0,(19)meaning that the mapping has a unique maximum in phase space and will quickly decay to zero everywhere else. In other words, every input chirplet will only significantly affect a small number of output chirplets, thereby enabling efficient numerical implementations.

The location of this phase space maximum is highly structured across chirplets. This can be deduced from the linear canonical transform model previously discussed in Section 2.B, expressing the Fresnel diffraction as a phase space shear. This principle is illustrated in Fig. 5.

This means one can accurately compute the value of an output Gabor coefficient by only considering a small sheared neighborhood in the appropriate location in phase space corresponding to a few input Gabor coefficients. The size of this χ×ω region-of-interest (ROI) neighborhood can be chosen based on the user’s computational complexity versus speed requirements.

Consider the source and destination plane represented by the discrete Gabor transform with M modulations and T translations. Suppose that we want to compute a particular destination Gabor coefficient with coordinates (td,md), where td∈{0,...,T−1} and md∈{0,...,M−1}. Given the chosen size χ×ω of the ROI, which source Gabor coefficient array coordinates (ts,ms) will have the highest contributions?

We observe that the shear is purely horizontal since the Fresnel transform is an all-pass filter. This indicates that the most relevant source frequencies will be centered around md. We let ms vary in the range of ms∈{⌊md−ω2⌋,…,⌊md+ω2⌋},(20)where ⌊·⌋ is the floor operator (rounding down). For each row ms, we have a different spatial displacement matching the phase space shear. The shearing displacement for that row is given by ss=π(M2−ms)KMT,(21)where T is the number of translations so that we only consider source Gabor coefficient coordinates in that row satisfying ts∈{⌊td+ss−χ2⌉,…,⌊td+ss+χ2−1⌉},(22)where ⌊·⌉ is the rounding operator. In this case, we treat coefficients for coordinates outside the source Gabor coefficient array as zero. Experiments measuring the impact of this method are reported in Section 4.D.

We will cover multiple aspects of the newly proposed propagation algorithm for the experiments section. We will look at the properties of the PSF, test the propagation operator on a sample image, and measure the computation speed and accuracy.

As discussed in Section 2.A and Ref. [1], the continuous Fresnel kernel K(x)=exp(iKx2) cannot be straightforwardly discretized. Different approaches can be taken depending on assumptions made about sampling conditions, propagation distance, wavelength, and desired signal properties.

For this experiment, we examine the PSF resulting from a single point at a distance of d=25mm from a one-dimensional 1024-sample hologram with pixel pitch p=6μm and wavelength of λ=633nm. Three different algorithms were tested, cf. Fig. 6: the spatial domain method using Eq. (4), the proposed Gabor domain method, and the frequency domain method using Eq. (6) with zero padding.

The spatial domain method exactly matches K(x) in the valid (i.e., nonaliased) region but goes sharply to zero, which will cause visual artifacts, especially apparent at edges (see also Section 4.B). The frequency domain method will smoothly drop toward zero in the spatial domain, being defined everywhere over the hologram. However, it will not match K(x) precisely in the spatial domain, showing highly oscillatory small perturbations, so-called ringing, especially in the low-frequency regions of the signal. The frequency domain method will also become less accurate at long distances due to the required bandpass filtering [12]. The proposed method combines good properties from the spatial and frequency methods. In the middle of the PSF corresponding to low frequencies, the signal exactly matches K(x) like the spatial method. As we go outwards, there is a transient region similar to the frequency method, but its extent is smaller due to the locality of the Gabor frame; it will drop faster and eventually become identically zero at larger distances, unlike the frequency domain method. Nonetheless, the overall properties of the proposed PSF are closer to the FFT version.

The desired properties of the discrete Fresnel diffraction operator are partially controllable depending on the properties of the Gabor frame, offering simultaneous accurate approximations to the PSF in the center for low frequencies and smooth transients at the edges. This aspect is controllable depending on the chosen space and frequency resolution of the Gabor transform by altering the shape of the Gabor lattice expressed by the number of translations T and modulations M.

For this series of experiments, we compare the relative accuracy of a propagated image for the three Fresnel diffraction methods: the spatial domain, frequency domain, and proposed Gabor domain method. We use the 1024×1024pixels hologram depicting a bird, whose properties were set to p=4μm and λ=532nm. We chose three different propagation distances d=5,10,20mm and compared the PSNR of the reconstructions for the three methods in Table 2. Example reconstructions are shown in Fig. 7.Table 2.

Comparing the PSNR between Different Fresnel Diffraction Methods as a Function of Propagation Distance

This confirms that the methods are visually highly similar but not identical due to the different PSF approximations. The frequency domain method and the Gabor method are more similar to each other, especially for shorter propagation distances, where the spatial method has a smaller PSF with fewer samples, leading to larger errors.

Another experiment we performed was to use different settings for the source and destination holograms, demonstrating the flexibility of the proposed algorithm. It can compute holograms where the pixel pitch and hologram resolution are different using off-axis propagation without any modifications or zero-padding. We used a 1024×1024pixels source hologram with a pixel pitch of pin=4μm and λ=532nm, propagating over a distance of d=5mm to a destination hologram with a resolution of 1536×1536pixels and a pixel pitch of pout=5μm, with an addition lateral displacement of 200 μm along x. This yielded a PSNR of 38.7 dB. The reconstructions are shown in Fig. 8.

The frequency domain algorithms do not support arbitrary changes in resolution. It is achievable with a spatial domain algorithm; however, it will no longer be a conventional discrete convolution, since the input samples have a different sampling density than the output samples.

This flexibility means that it is possible to compute small parts of the wavefield with any sampling density, allowing for the efficient computation of particular regions of interest, such as specific simulated viewpoints for holographic displays. Since only a subset of the Gabor coefficients of the input hologram is needed in that case, it may be possible to further accelerate the algorithm for view-dependent reconstructions by only computing the relevant coefficients of the input hologram.

We implemented the proposed algorithm on GPU for the calculation time experiments, comparing it to the frequency domain method utilizing the FFT. We used a machine with an AMD Ryzen Threadripper PRO 5995WX CPU, 256 GB RAM, NVIDIA RTX A6000 GPU, and a Windows 11 OS. The code was implemented in C++20 with CUDA 12.1, enabling CUDA compute capability 8.6 and using 32-bit floating-point precision.

We calculated the Fresnel diffraction operator of an N×512 pixels hologram, where N was a variable parameter. The properties of the hologram were set to p=8μm, λ=633nm, and d=1cm. We set the ROI window size for the proposed Gabor propagation algorithm to 5×5 coefficients. The results are graphed in Fig. 9.

The parameter N was varied from 210=1024 up to 219=524,288, and the experiment was run 20 times for each value of N, reporting average calculation time. The proposed algorithm is consistently faster, reporting speed gain ratios rising from approximately 1.7× up to 2.0× for increasing N. Although the FFT algorithm has O(nlogn) complexity, it is highly optimized and streamlined for GPUs, so exceeding it is a challenge.

The computational advantage of our proposed O(n) complexity algorithm should be more pronounced in highly memory-bound environments (such as with FPGA or ASIC) thanks to the phase space locality of the algorithm. It can be more efficient when the output resolution is not an integer multiple of the input resolution. Efficient approximations of the complex exponential and more sophisticated GPU implementations could also further extend the gap with FFT-based algorithms.

Depending on the application’s precision and speed requirements, one should be careful what specific programming instructions are used for computing the different coefficients [34,35]; this is especially important for lower-precision fixed-point integer representations, which are beneficial for significant further acceleration.

In this section, we evaluate the accuracy of the version of the approximated method by varying the size of the ROI neighborhood χ×ω detailed in Section 3.B. We do this by computing the PSNR of a random complex-valued signal with N=512, with hologram properties of p=6μm, λ=633nm, and d=1cm. The results were compared with the reference, non-approximated version of the algorithm. The results are reported in Table 3.Table 3.

Quantitative Assessment Using the PSNR^a

These show that it is not useful to increase only χ or ω, as the accuracy will eventually saturate. Rather, it is better to increase both parameters since the magnitude appears to decay radially like a Gaussian in phase space, confirming the appropriateness of the model derived in Eqs. (18) and (19). Thus, the PSNR will increase roughly proportionally to the ROI area for a given ratio χ×ω.

We combine two of Dennis Gabor’s inventions, holography and the Gabor transform, creating a novel method for numerical Fresnel diffraction. Different notions from signal processing were brought together, such as sparsity, phase space, chirplets, frames, and numerical Fresnel diffraction, resulting in a novel mathematical model for the propagation of coherent wavefronts of light. The algorithm is shown to be faster than FFT-based diffraction. Multiple experiments were performed, analyzing the PSF shape and the algorithm’s speed and accuracy.

In future work, we aim to extend the algorithm to handle propagation between arbitrarily oriented planes or surfaces. These can be achieved with affine transformations of phase space; as this requires non-uniform resampling of Fourier space [24], extensions to the current chirplet algorithm model will be needed. Furthermore, this could be generalized to represent arbitrary 3D surfaces by a non-planar collection of chirplets that can be propagated efficiently to the hologram plane. This could contribute to various applications utilizing CGH, such as diffractive simulations for waveguide-based displays and non-line-of-sight imaging.

As an alternative Fresnel diffraction operator with a flexible phase space mapping operator, the proposed method could be integrated with iterative CGH algorithms, using Gerchberg-Saxton or stochastic gradient descent, enabling the potential for more selective iterative phase space updates.

Future work may also include optimizing the Gabor frame parameters or utilizing more general frames and Gabor lattices. We plan to investigate how to integrate this model with different CGH algorithms, assess to what extent further optimizations are possible on GPU, and how suitable the algorithm is for hardware implementations on FPGA or ASIC.