Theory and algorithm of the homeomorphic Fourier transform for optical simulations

Zongzhao Wang; Zongzhao Wang; Olga Baladron-Zorita; Olga Baladron-Zorita; Christian Hellmann; Frank Wyrowski

doi:10.1364/OE.388022

1. Introduction

The Fourier transform is one of the most important and widely used mathematical tools [1,2] in pretty much every technical discipline, including physical optics, where it can be most notably employed to analyze the spectral representation of an electromagnetic field (one-dimensional Fourier transform between the time domain, $t$, and the spectral frequency, $\omega$, domains), as well as to obtain its plane-wave composition (two-dimensional Fourier transform between the space domain, $x, y$, and the spatial frequency domain, $k_x, k_y$). Consequently, an efficient implementation of the Fourier transform can benefit a large number of applications [3–6]. In 1965, J. W. Cooley and J. W. Tukey published a paper [7] to introduce a most important technique: the Fast Fourier transform (FFT). The advent of the FFT represented a massive improvement on the computing efficiency since, while the standard Discrete Fourier Transform (DFT) incurs a numerical complexity of O($N^2$) for a function with $N$ sampling points, the FFT reduces that to $N \log N$ operations, which can be considered as approximately linear on the sample number of the function to be transformed.

These unique computing advantages make the FFT one of the most influential algorithms in the field of numerical simulations. However, even with the benefit of these advantages, in physical optics it is still often considered an insurmountable challenge to work with field components that present a strong wavefront phase. The periodic nature of the complex exponential function constrains all phase points to the range $(-\pi , \pi ]$, commonly referred to as “$2\pi$-modulo” phase or “$2\pi$-wrapped” phase [8,9]. The wrapped smooth phase becomes jaggedly discontinuous, and the steeper the slope of the original smooth phase function, the denser the teeth of its wrapped counterpart become. The most pronounced slope in the phase function restricts the period which must be used when sampling the field. In short, comparing the wrapped phase function to the smooth original, the wrapped phase leads to much larger sample numbers necessary for good resolution.

In order to overcome this sampling overkill that plagues the FFT algorithm, people look forward to an alternative approach to handle the smooth wavefront phase instead of the fully sampled wrapped phase. Several rigorous solutions to deal with the quadratic phase when performing the FFT have been proposed; for instance, the chirped-z transform [10–12] or the semi-analytical Fourier transform [13]. However, this kind of strategy can fall short in non-paraxial cases. This means we must explore a different line of thought. From the mathematical point of view, the Fourier-transform integrals are rapidly oscillating functions in the case of strong wavefront phases. There is a mathematical method, working on the basis of an asymptotic approximation to integrals of such rapidly varying functions, which occasionally comes up in wave theory, applied usually with the aim of obtaining an analytical solution to a specific problem (as is the case with the Debye integral [14,15]): the Method of Stationary Phase (MSP) [16,17]. Jakob J. Stamnes already revealed the connection between waves and rays with the application of the Method of Stationary Phase in his publication [18]. Despite the MSP being an algorithm with a much more general scope of application to the Fouier integral [19], its application is restricted usually only to the spherical phase or the quadratic phase [20].

In this work, we propose an approximated algorithm named “homeomorphic Fourier transform", whose aim is to efficiently compute the Fourier transform of a field that presents a strong smooth wavefront phase. The full derivation is presented in Section 2, taking as our starting point an alternative expression of the electromagnetic field in which a phase term, which is referred to as “wavefront phase", is extracted, and using the method of stationary phase to solve the integrand. After that, in Section 3 we consider the application of the algorithm to numerical simulations. In four groups of experiments, we select different wavefront phases and test their influence on the Fourier transform, from the most fundamental spherical phase, through the case where a single aberration term is present, to a general arbitrary numerical phase. Simulation results are presented alongside an accuracy and efficiency analysis to illustrate the advantages of the homeomorphic Fourier transform.

2. Theory

2.1 Field representation and hybrid sampling strategy

In physical optics, light is represented as an electromagnetic field, a vector field with six components, i.e., $\left \{ E_x, E_y, E_z, H_x, H_y, H_z \right \}$. In general, any electromagnetic field component is complex valued. To avoid inconveniently working with harmonic functions, we can write any field component in terms of its amplitude and phase,

(1)$$\begin{array}{lll} V\!\left( {\boldsymbol{\rho}} \right) & = & \left| V\!\left( {\boldsymbol{\rho}} \right) \right| \exp\left\{ \mathrm{i} \arg \left[ V\!\left( {\boldsymbol{\rho}} \right) \right] \right\} \\ & = & \left| V\!\left( {\boldsymbol{\rho}} \right) \right| \exp\left\{ \mathrm{i} \gamma \left({\boldsymbol{\rho}}\right) \right\} \end{array}$$

where the field component $V\!\left ( {\boldsymbol{\rho}} \right )$ is defined on a plane $\mathsf {P}$ which is assumed embedded in a homogeneous medium. In our notation, ${\boldsymbol {r}} = \left ( x, y, z \right )$ and ${\boldsymbol{\rho}} = {\boldsymbol {r}}_\perp = \left ( x, y \right )$ denote the three-dimensional position vector and its projection onto the transversal plane respectively. $\gamma \left ({\boldsymbol{\rho}}\right )$ indicates the phase of the field component. Besides, we introduce the symbol $N \left ( V \right ) = N_x \times N_y$ to indicate the number of sampling points for the given field $V\!\left ( {\boldsymbol{\rho}} \right )$.

In a numerical simulation, the FFT technique requires a strict equidistant Nyquist sampling, which applies not only to the amplitude but also to the “$2\pi$-modulo” phase $\exp \left \{ \mathrm {i} \gamma \left ({\boldsymbol{\rho}}\right ) \right \}$. It is very typical in optics to encounter fields whose sample number runs into several orders of magnitude, in a way that is caused by the phase part, i.e. $N (V) \sim N ( \exp \left \{ \mathrm {i} \gamma \left ({\boldsymbol{\rho}}\right ) \right \} )$. When this is the case, the computational advantage of the FFT can fall short. In what follows, we propose an alternative sampling approach: a hybrid sampling strategy.

First of all, we reformulate Eq. (1) to yield

(2)$$\begin{array}{lll} V\!\left( {\boldsymbol{\rho}} \right) & = & \left| V\!\left( {\boldsymbol{\rho}} \right) \right| \exp \left\{\mathrm{i} \left[ \phi\!\left( {\boldsymbol{\rho}} \right) + \psi\!\left( {\boldsymbol{\rho}} \right) \right] \right\} \\ & = & U\!\left( {\boldsymbol{\rho}} \right) \exp\left[ \mathrm{i}\psi\!\left( {\boldsymbol{\rho}} \right) \right] \textrm{,} \end{array}$$

where we have extracted the smooth wavefront phase $\psi \!\left ( {\boldsymbol{\rho}} \right )$ – the choice of which is arbitrary – and grouped the rest of the phase $\phi \!\left ( {\boldsymbol{\rho}} \right )$ alongside the amplitude $\left | V\!\left ({\boldsymbol{\rho}}\right )\right |$ into $U\!\left ({\boldsymbol{\rho}}\right )$. It should be noted that Eq. (2) entails no approximation whatsoever: it merely constitutes an alternative way to express $V\!\left ({\boldsymbol{\rho}}\right )$.

We have now rewritten the expression of our field component in a way which isolates the troublesome phase term $\exp \left [\mathrm {i}\psi \!\left ({\boldsymbol{\rho}}\right )\right ]$, so now, considering that the present discussion is not merely concerned with the underlying mathematics, but also with the implementation of the method for simulation purposes, in order to numerically store $V\!\left ({\boldsymbol{\rho}}\right )$ we could follow a hybrid sampling strategy, one which employs different sampling and interpolation techniques for the two terms involved (namely, $U\!\left ({\boldsymbol{\rho}}\right )$ and $\exp \left [\mathrm {i}\psi \!\left ({\boldsymbol{\rho}}\right )\right ]$): a strict equidistant Nyquist sampling applied to the residual complex amplitude $U\!\left ({\boldsymbol{\rho}}\right )$, and a non-equidistant sampling approach for the smooth wavefront phase $\psi \!\left ({\boldsymbol{\rho}}\right )$, with such varied options as spline or quadratic interpolation methods available for $\psi \!\left ({\boldsymbol{\rho}}\right )$. The application of this hybrid approach provides a workaround that allows us to avoid having to sample the exponential term associated with the wavefront phase, as $N \left ( \exp \left [ \mathrm {i} \gamma \! \left ({\boldsymbol{\rho}}\right ) \right ] \right ) \gg N \left ( \gamma \! \left ({\boldsymbol{\rho}}\right ) \right )$. Very often, in the absence of the wrapped wavefront phase, the sampling of $U\!\left ({\boldsymbol{\rho}}\right )$ does not pose much of a challenge, and $N \left (V \right ) \gg N \left ( U \right )$. Consequently, for strong wavefront phases, the hybrid sampling strategy can result in a dramatic decrease of the computational effort.

Then, from a practical point of view, let us analyze the practicalities of the hybrid sampling strategy. The most important two questions are how to obtain the smooth wavefront phase and how to deal with it in actual simulations of optical systems. The complete solution to obtain the optimal smooth wavefront phase from arbitrary complex field data is to combine the phase unwrapping algorithm and an appropriate interpolation algorithm. However, it is common for the wavefront phase to be directly established with the model of the source, e.g., an ideal spherical wave. And then, for most optical components, the physical response of the element on the wavefront phase is known. For instance, in the case of lenses and interfaces. However, for Fourier transform operations, there is no straightforward solution: the Fast Fourier Transform (FFT) algorithm cannot handle a field sampled in a hybrid manner. With an eye to taking the numerical advantage of the hybrid sampling strategy into the Fourier transform, we have developed the following algorithm, which we refer to as “homeomorphic Fourier transform”.

2.2 Derivation of the homeomorphic Fourier transform

Let us next substitute Eq. (2) in the definition of the Fourier-transform integral,

(3)$$\begin{array}{lll} \tilde{V}\! \left( {\boldsymbol{\kappa}} \right) & = & \frac{1}{2\pi} \iint^{\infty}_{-\infty} V\! \left( {\boldsymbol{\rho}} \right) \exp\left( -\mathrm{i} {\boldsymbol{\rho}} \cdot {\boldsymbol{\kappa}} \right) \mathrm{d}^2\rho = \\ & = & \frac{1}{2\pi} \iint^{\infty}_{-\infty} U\!\left( {\boldsymbol{\rho}} \right) \exp\left[ \mathrm{i}\psi\!\left( {\boldsymbol{\rho}} \right) - \mathrm{i} {\boldsymbol{\rho}} \cdot {\boldsymbol{\kappa}} \right] \mathrm{d}^2\rho \textrm{,} \end{array}$$

where ${\boldsymbol {\kappa }} = \left (k_x, k_y \right )$ is the transversal spatial frequency vector and the tilde in $\tilde {V}$ makes explicit reference to a function defined in the spatial frequency domain, spanned by ${\boldsymbol {\kappa }}$, as opposed to a function defined in the space domain – spanned by ${\boldsymbol{\rho}}$, as already introduced above – which would have no tilde.

It is relevant at this point to highlight the similarities between Eq. (3) and the stationary phase approximation formula, the working assumptions for and derivation of which can be found in the work by Leonard Mandel and Emil Wolf [16]. It is directly from said cited work that we extract the following two preconditions for the application of the stationary phase method:

• The term $U\!\left ({\boldsymbol{\rho}}\right )$ must be slowly varying.
• In the domain where the transversal position vector ${\boldsymbol{\rho}}$ is defined, there is one and only one critical point of the first kind defined by the formula ${\boldsymbol {\nabla }}\left [ \psi \!\left ( {\boldsymbol{\rho}} \right ) - {\boldsymbol{\rho}} \cdot {\boldsymbol {\kappa }} \right ] = 0$.

Mathematically, the first condition means that the integrand in Eq. (3) is dominated by the exponential term. The second condition implies a bijective mapping from the spatial domain to the spatial frequency domain, i.e. a homeomorphism. In other words, the conditions require that the field be in a region of “normal congruence” and that caustics be avoided. From the physical point of view, the gradient of the phase term has the same dimension as the spatial frequency. When these two conditions are satisfied, we obtain a bijective mapping relation between the local gradient and the spatial frequency. It is referred to as “local spatial frequency” in [2].

The application of the stationary phase method to Eq. (3) produces the following result:

(4)$$\tilde{V}\left( {\boldsymbol{\kappa}} \right) \approx a\!\left[ {\boldsymbol{\rho}} \left( {\boldsymbol{\kappa}} \right) \right] U\!\left[ {\boldsymbol{\rho}} \left( {\boldsymbol{\kappa}} \right) \right] \exp\left\{ \mathrm{i} \psi\! \left[ {\boldsymbol{\rho}} \left( {\boldsymbol{\kappa}} \right) \right] - \mathrm{i} {\boldsymbol{\kappa}} \cdot {\boldsymbol{\rho}}\!\left( {\boldsymbol{\kappa}} \right) \right\} \textrm{,}$$

with the bijective mapping relation ${\boldsymbol{\rho}} \rightarrow {\boldsymbol {\kappa }}$

(5)$${\boldsymbol{\nabla}}\psi\!\left({\boldsymbol{\rho}}\right) = {\boldsymbol{\kappa}} \textrm{,}$$

and with

(6)$$a\!\left({\boldsymbol{\rho}}\right) = \left\{ \begin{array}{ccl} \sqrt{\frac{\mathrm{i}}{\psi_{xx}\left({\boldsymbol{\rho}}\right)}} \sqrt{-\frac{\mathrm{i}\psi_{xx}\left({\boldsymbol{\rho}}\right)} {\psi^2_{xy}\left({\boldsymbol{\rho}}\right)-\psi_{xx}\left({\boldsymbol{\rho}}\right)\,\psi_{yy}\left({\boldsymbol{\rho}}\right)}} & , & \psi_{xx}\left( {\boldsymbol{\rho}} \right) \neq 0 \\ \frac{1}{\left| \psi_{xy}\left({\boldsymbol{\rho}}\right) \right|} & , & \psi_{xx}\left( {\boldsymbol{\rho}} \right) = 0 \end{array} \right. \textrm{,}$$

where $\psi _{x_ix_j} \overset {\textrm {def}}{=} \frac {\partial \psi }{\partial x_i \partial x_j}$.

The sign (in a broad sense of the word) of the amplitude scaling factor $a$ defined in Eq. (6), $\operatorname {sgn}\left [ a\!\left ( {\boldsymbol{\rho}} \right )\right ]$, can take one of three different values: $\mathrm {i}$, $-\mathrm {i}$ and $1$. Each of the three options has different physical implications; the first two have already been analysed in literature [18] and correspond, respectively, to an overall convergent or an overall divergent field. The third case, $\mathrm {sign}\left [ a \!\left ( {\boldsymbol{\rho}} \right ) \right ] = 1$, denotes a situation in which the field is convergent in one dimension and divergent in the other.

It should be noted that the derivation presented above incurs no approximation whatsoever up to and including Eq. (3). It is only when the two preconditions for the validity of the stationary phase method are employed in the simplification of the Fourier integral to produce Eq. (4) that we finally abandon the path of full rigor. The resulting expression (namely, as already stated, Eq. (4)) constitutes an approximation, albeit a good and accurate one as long as the two aforementioned preconditions hold.

2.3 Inverse homeomorphic Fourier transform

The derivation of the inverse operation can be performed in an analogous manner. Let us then start from the spectrum of the field, $\tilde {V}\!\left ({\boldsymbol {\kappa }}\right )$, and rewrite it by extracting a smooth wavefront phase term, like we did for the direct homeomorphic Fourier transform in Eq. (2):

(7)$$\begin{array}{lll} \tilde{V}\!\left( {\boldsymbol{\kappa}} \right) & = & \left| \tilde{V}\!\left( {\boldsymbol{\kappa}} \right) \right| \exp\left\{ \mathrm{i} \arg \left[ \tilde{V}\!\left( {\boldsymbol{\kappa}} \right)\right] \right\} \\ & = & \tilde{A}\!\left( {\boldsymbol{\kappa}} \right) \exp\left[ \mathrm{i}\tilde{\psi}\! \left( {\boldsymbol{\kappa}} \right) \right] \end{array}$$

where the term $\tilde {\psi }\! \left ( {\boldsymbol {\kappa }} \right )$ is the smooth wavefront phase in the frequency domain. Identifying this expression with the direct homeomorphic Fourier transform presented in Eq. (4) it is possible to write $\tilde {A}\!\left ({\boldsymbol {\kappa }}\right )$ and $\tilde {\psi }\!\left ({\boldsymbol {\kappa }}\right )$ as follows:

(8)$$\left\{ \begin{array}{lll} \tilde{A}\!\left({\boldsymbol{\kappa}}\right) & = & a\!\left[ {\boldsymbol{\rho}} \!\left( {\boldsymbol{\kappa}} \right) \right] U\!\left[ {\boldsymbol{\rho}} \!\left( {\boldsymbol{\kappa}} \right) \right] \\ \tilde{\psi}\!\left( {\boldsymbol{\kappa}} \right) & = & \psi\!\left[ {\boldsymbol{\rho}} \!\left( {\boldsymbol{\kappa}} \right) \right] - {\boldsymbol{\kappa}} \cdot {\boldsymbol{\rho}} \!\left( {\boldsymbol{\kappa}} \right) \end{array} \right. \textrm{.}$$

We can then use Eq. (7) to substitute for $\tilde {V}$ in the definition of the inverse Fourier transform integral

(9)$$\begin{array}{lll} V\!\left({\boldsymbol{\rho}}\right) & = & \frac{1}{2\pi} \iint^{\infty}_{-\infty} \tilde{V}\!\left( {\boldsymbol{\kappa}} \right) \exp\left( \mathrm{i}{\boldsymbol{\rho}} \cdot {\boldsymbol{\kappa}} \right) \mathrm{d}^2\kappa = \\ & = & \frac{1}{2\pi} \iint^{\infty}_{-\infty} \tilde{A}\!\left({\boldsymbol{\kappa}}\right) \exp\left[ \mathrm{i}\tilde{\psi}\!\left( {\boldsymbol{\kappa}} \right) + \mathrm{i} {\boldsymbol{\rho}} \cdot {\boldsymbol{\kappa}} \right] \mathrm{d}^2\kappa \textrm{.} \end{array}$$

Again, carrying out a similar process as that for the direct transform, the smooth phase term $\tilde {\psi }\!\left ({\boldsymbol {\kappa }}\right )$ will provide the inverse bijective mapping relation ${\boldsymbol {\kappa }} \rightarrow {\boldsymbol{\rho}}$,

(10)$$\tilde{{\boldsymbol{\nabla}}}\tilde{\psi}\!\left({\boldsymbol{\kappa}}\right) ={-}{\boldsymbol{\rho}} \textrm{.}$$

Applying the stationary phase method, the integral in Eq. (9) becomes

(11)$$V\!\left({\boldsymbol{\rho}}\right) = \tilde{a}\!\left[ {\boldsymbol{\kappa}}\!\left( {\boldsymbol{\rho}} \right) \right] \tilde{A}\!\left[ {\boldsymbol{\kappa}}\!\left( {\boldsymbol{\rho}} \right) \right] \exp\left\{ \mathrm{i}\tilde{\psi}\!\left[ {\boldsymbol{\kappa}}\!\left( {\boldsymbol{\rho}} \right)\right] + \mathrm{i}{\boldsymbol{\rho}}\cdot{\boldsymbol{\kappa}}\!\left({\boldsymbol{\rho}}\right)\right\} \textrm{,}$$

with

(12)$$\tilde{a}\!\left({\boldsymbol{\kappa}}\right) = \left\{ \begin{array}{ccl} \sqrt{\frac{\mathrm{i}}{\tilde{\psi}_{k_xk_x}\left({\boldsymbol{\kappa}}\right)}} \sqrt{-\frac{\mathrm{i}\tilde{\psi}_{k_xk_x}\left({\boldsymbol{\kappa}}\right)}{\tilde{\psi}^2_{k_xk_y}\left({\boldsymbol{\kappa}}\right)-\tilde{\psi}_{k_xk_x}\left({\boldsymbol{\kappa}}\right)\,\tilde{\psi}_{k_yk_y}\left({\boldsymbol{\kappa}}\right)}} & , & \tilde{\psi}_{k_xk_x}\left({\boldsymbol{\kappa}}\right) \neq 0 \\ \frac{1}{\left| \tilde{\psi}_{k_xk_x}\left({\boldsymbol{\kappa}}\right) \right|} & , & \tilde{\psi}_{k_xk_x}\left({\boldsymbol{\kappa}}\right) = 0 \end{array} \right. \textrm{.}$$

It is then possible to identify the terms

(13)$$\left\{ \begin{array}{lll} U\!\left({\boldsymbol{\rho}}\right) & = & \tilde{a}\!\left[{\boldsymbol{\kappa}}\!\left({\boldsymbol{\rho}}\right)\right] \tilde{A}\!\left[ {\boldsymbol{\kappa}}\!\left({\boldsymbol{\rho}}\right) \right] \\ \psi\!\left({\boldsymbol{\rho}}\right) & = & \tilde{\psi}\!\left[{\boldsymbol{\kappa}}\!\left({\boldsymbol{\rho}}\right)\right] + {\boldsymbol{\rho}}\cdot{\boldsymbol{\kappa}}\!\left({\boldsymbol{\rho}}\right) \end{array} \right. \textrm{,}$$

in Eq. (11), so that after subsequent application of first the direct and then the inverse homeomorphic Fourier transforms, the original expression of the field $V\!\left ({\boldsymbol{\rho}}\right )$ is retrieved, thus proving the consistency of the method.

2.4 Pointwise Fourier transform in the case of a non-bijective phase

As already covered in the previous section, one of the preconditions for the validity of both the direct and the inverse homeomorphic Fourier transform (HFT and IHFT respectively, for short) is that the smooth wavefront phase be capable of generating a bijective mapping relation between the Fourier domains, ${\boldsymbol{\rho}}$ and ${\boldsymbol {\kappa }}$. The characteristics of the bijection and how they relate to the generation of a homeomorphism are shown in the Appendix.

However, it is possible in an optics context to encounter situations where this pre-requisite for bijection is not fulfilled, a phenomenon which often happens in the presence of strong aberration. In such kind of situation, different points in the starting domain might map to the same point in the other domain, namely a surjective mapping. Let us recall that, when performing the decomposition of the field outlined in Eq. (2), we have total freedom when deciding exactly what part of the total phase is included in the wavefront term, and what part remains in $U$. This choice will evidently affect the eventual accuracy of the HFT, so that, ideally, the aim is for as much of the smooth phase as possible to be contained in $\psi \!\left ({\boldsymbol{\rho}}\right )$. But we must bear in mind that the condition for bijectivity must prevail, as on it hinges the consistency of the definition of the HFT/IHFT method. Therefore, should we come across a situation in which the initially given $\psi$ term violates the bijectivity condition, we need a strategy to eliminate the source of that transgression while leaving as much of the $\psi$ as possible intact. This will go a long way in the task of developing a robust algorithm and extending the validity of the HFT.

Starting, as before, from the field $V\!\left ({\boldsymbol{\rho}}\right )$, but already provided in the decomposed form $V\!\left ({\boldsymbol{\rho}}\right ) = U\!\left ({\boldsymbol{\rho}}\right ) \exp \left [\mathrm {i} \psi \!\left ({\boldsymbol{\rho}}\right )\right ]$, we aim to extract the bijective part $\psi ^{\textrm {map}}\!\left ({\boldsymbol{\rho}}\right )$ from the full smooth phase term $\psi \!\left ({\boldsymbol{\rho}}\right )$, thus:

(14)$$\begin{array}{lll} V\!\left({\boldsymbol{\rho}}\right) & = & U\!\left({\boldsymbol{\rho}}\right) \exp\left[\mathrm{i} \psi\!\left({\boldsymbol{\rho}}\right)\right] = \\ & = & U\!\left( {\boldsymbol{\rho}} \right) \exp\left[ \mathrm{i} \psi^{\textrm{res}}\!\left( {\boldsymbol{\rho}} \right) \right] \exp\left[ \mathrm{i} \psi^{\textrm{map}}\!\left( {\boldsymbol{\rho}} \right) \right] \textrm{.} \end{array}$$

For the numerical implementation of this strategy we propose the following roadmap:

1. Perform phase fitting on the initial wavefront phase. As a first step, we sequentially extract, using a fitting technique, the linear phase, spherical phase and Zernike terms from the initial smooth phase $\psi \!\left ({\boldsymbol{\rho}}\right )$, resulting in $\psi \!\left ({\boldsymbol{\rho}}\right ) = \psi ^{\textrm {fit}}\!\left ({\boldsymbol{\rho}}\right ) + \Delta \psi \!\left ({\boldsymbol{\rho}}\right )$, where the difference $\Delta \psi$ between the original phase $\psi$ and the fitted result $\psi ^{\textrm {fit}}$ should be negligible.
2. Elementary optimization. Since the analytical expression of the fitting model is known, we can now examine whether the result of the fitting $\psi ^{\textrm {fit}}$ is bijective, according to the homeomorphism criterion (as presented in the Appendix). Numerically, we compute the second derivatives of the fitted phase $\psi ^{\textrm {fit}}$ at the position of each sampling point, and check whether they have the same sign. Should this check reveal a non-bijective phase, an iterative approach is introduced to omit the highest-order term in the list of Zernike polynomials and afterwards re-examine the homeomorphism criterion. This step is repeated until a bijective phase $\psi ^{\textrm {map}}$ is obtained.
3. Advanced optimization. In the elementary optimization we roughly omit the highest-order Zernike term, but it is important to note that Zernike polynomials are global functions which operate on the entirety of the definition domain, and many commonly encountered cases of non-bijective phases are caused by local phase defects. It would then make sense, taking the result of the elementary optimization as the starting point, to perform a more thorough investigation of the possible sources of the violation of bijectivity, with more of a view to identifying the more local ones and extracting them from the $\psi$ in a more pinpointed manner.

We would now, after the application of the steps described above, have succeeded in the extraction of the bijective phase part $\psi ^{\textrm {map}}\!\left ({\boldsymbol{\rho}}\right )$, which provides the mapping relation ${\boldsymbol{\rho}} \rightarrow {\boldsymbol {\kappa }}$,

(15)$${\boldsymbol{\kappa}}^{\textrm{map}} = {\boldsymbol{\nabla}} \psi^{\textrm{map}}\!\left({\boldsymbol{\rho}}\right) \textrm{.}$$

We could now apply the stationary phase method, but this time using the mapping relation from Eq. (15). The resulting expression is the same as the standard homeomorphic Fourier transform. It reveals that the only difference is that the general mapping relation has been replaced by the bijective mapping. The approximation intrinsic to the forced regularization of the mapping relation will certainly lead to a dislocation in the spectrum position and a resulting, inevitable, error. This error can be numerically investigated and compared against a given accuracy threshold, on the basis of which the final decision could be made whether to go ahead with the application of the homeomorphic Fourier transform in that specific case, or to discard it in favor of the rigorous variant.

3. Application

So far in this paper, we have focused on the theoretical derivation of the homeomorphic Fourier transform. However, one of the most promising aspects of this concept is its application in computer simulations of optical systems. This new technique has been implemented in the physical optics modeling and design software VirtualLab Fusion [21]. In this coming section, therefore, we will present several illustrative examples of the application of this technique.

3.1 Fourier transform of a field with spherical wavefront phase

Let us start from a simple optical setup: an ideal plane wave, with a wavelength of $532$ nm, illuminates a house-shaped mask behind which a divergent spherical phase $\psi ^{\textrm {sph}}\!\left ({\boldsymbol{\rho}}\right ) = \operatorname {sgn}\left (R\right )k_0\check {n}\sqrt {{\boldsymbol{\rho}}^2+R^2}$ is superimposed onto the field. Here, $\check {n}$ indicates the complex refractive index of the surrounding medium (standard air, $\check {n} = 1.0003$) and the radius of curvature $R$ of the spherical wavefront is assumed to take on a positive value, i.e. $R>0$.

For this first set of experiments the amplitude distribution is fixed, as per Fig. 1(a), which shows the $E_x$ component of the field behind the mask. The radius of curvature $R$ of the spherical phase which is then superimposed on this amplitude is employed as the variable. The result of the fast Fourier transform corresponding to different values of $R$ is then shown in Fig. 1(b)–(g). Mere visual investigation reveals a clear tendency: in the case, pictured in panel (b), of $R=+\infty$ (or, in other words, in the limit case of a plane wavefront) the original pattern and its Fourier transform bear no resemblance whatsoever to each other, virtue of the integral nature of the Fourier-transform operation; however, as the magnitude of $R$ progressively decreases, and the importance of the spherical phase increases as a direct result thereof, a one-to-one mapping behavior seems to take over and one can ever more easily discern a house-shaped distribution also in the Fourier transform. For the specific configuration studied here, panel (e) happens to be a particularly evident instance on that score – the mapping behavior is crystal clear by this point already, but the distortion is still temperate. Still, even for stronger distortion, as is the case for the results of panels (f) and (g), the mapping behavior is incontestable, since distortion does not inherently violate the homeomorphism: the definition of a homeomorphism hinges on the fact that the coordinates may be changed by the transformation, but the neighbors of the points remain the same. In addition, from panel (e) to (g), the distortion effect becomes progressively evident. It is because that the bijective mapping from ${\boldsymbol{\rho}}$ to ${\boldsymbol {\kappa }}$ follows the spherical phase function which cannot be approximated to a linear mapping in the high NA nonparaxial regime.Serving as a contrast, we compute the homeomorphic Fourier transform of the same house-shaped fields, shown in Fig. 2. Since a comparison of both operations only makes sense when the homeomorphic Fourier transform constitutes a good approximation of the rigorous integral, we have performed an analysis of the accuracy of the HFT using the FFT as reference. Here, according to Wyrowski and Kuhn [22], the deviation between the FFT result and the HFT result is evaluated by

(16)$$\sigma := \frac{\sum_{x,y} \left| \tilde{V}^{\textrm{FFT}} \left({\boldsymbol{\kappa}}\right) - \tilde{V}^{\textrm{HFT}} \left({\boldsymbol{\kappa}} \right) \right|^2}{\sum_{x,y} \left|\tilde{V}^{\textrm{FFT}} \left({\boldsymbol{\kappa}}\right) \right|^2} \,$$

where the reference $\tilde {V}^{\textrm {FFT}} \left ({\boldsymbol {\kappa }}\right )$ is computed by the rigorous fast Fourier transform technique and $\tilde {V}^{\textrm {HFT}} \left ({\boldsymbol {\kappa }}\right )$ denotes the result obtained by the homeomorphic Fourier transform. For each case, the deviation between the two approaches is presented in Tab. 1. The deviation exhibits initially high values at the low-NA which rapidly drop as the NA increases, with an asymptotic tendency as the NA continues to increase even further. That is, when the smooth wavefront phase term is strong enough to dominate the Fourier transform, the HFT will yield a result that is in practice indistinguishable from that of the rigorous FFT integral.

Fig. 1. Fast Fourier transform of a field with house-shaped amplitude and different values of the spherical-phase radius. Panel (a) shows the amplitude distribution of the $E_x$ component in the spatial domain, common starting point to all results. A spherical phase with radius of curvature $R$ is then superimposed on this amplitude, and its Fourier transform computed. Panels (b) to (g) present the $\tilde {E}_x$ component (i.e. the result of said Fourier-transform operation) for different values of $R$, as indicated below each of the individual panels. The required numbers of sampling points of fast Fourier transform for each case are given in Tab. 1.

Download Full Size | PDF

Fig. 2. Homeomorphic Fourier transform of the same house-shaped amplitude with a superimposed spherical phase in comparison to Fig. 1. Panel (a) shows the amplitude distribution of the $E_x$ component in the spatial domain. Panels (b) to (g) present the result of Homeomorphic Fourier-transform operation for different values of $R$. Please note the window size of panel (b) and the scale of panels (b) to (d) are different from Fig. 1. The required numbers of sampling points of homeomorphic Fourier transform for each case are given in Tab. 1.

Download Full Size | PDF

Table 1. Comparison of the homeomorphic Fourier transform and the fast Fourier transform: number of sampling points and deviation. The numerical aperture (NA) is defined as $n\sin \theta = \frac {na}{\sqrt {a^2+R^2}}$. Here, $a$ is the radius of the aperture and $n$ is the refractive index of the surrounding medium. Corresponding field distribution results are respectively shown in Fig. 1 and Fig. 2.

View Table | View all tables in this article

The investigation of the behavior of the fast Fourier transform and the homeomorphic Fourier transform with respect to the significance of the smooth wavefront (of which a spherically shaped one is just one example) for a given field distribution illustrates, is the effect of the stationary phase condition on the Fourier transform. In other words, when the aforementioned smooth wavefront is strong enough, the stationary phase assumption will be fulfilled and the Fourier transform will cease to exhibit an integral behavior, and effect a mere homeomorphism between the space and the spatial frequency domains: when this is true, the substitution of the (albeit strictly speaking approximated) homeomorphic Fourier transform – as presented in Section 2. – for the rigorous Fourier integral will be justified.

Let us at this point retrieve another dimension of the discussion which was already mentioned in the previous section: numerical effort. When, in Section 2., we presented the theoretical derivation of the homeomorphic Fourier transform, we indicated how using a hybrid sampling strategy for a field given in the form of Eq. (2) could effectively massively reduce the sampling effort compared against an indiscriminate Nyquist sampling of the entire field, including the $2\pi$-wrapped wavefront phase. In the following experiment we seek to compare the numerical effort of computing the Fast Fourier Transform (FFT) against that required to calculate the corresponding homeomorphic Fourier transform (HFT) of the same function. Even though the number of sampling points has been given in Tab. 1, for better visualization of the results and gauging of conclusions, we have plotted in Fig. 3 the numerical effort of the FFT and of the HFT in terms of the number of sampling points as function of the numerical aperture (NA). The graphs in Fig. 3 and Tab. 1 reveal that, as the numerical aperture increases as $R$ decreases, there is a steep rise in the number of sampling points required by the sampling conditions of the FFT operation. In the case of the HFT, however, the number of necessary sampling points remains constant with increasing NA. This is due to the fact that the residual field $U$ and the smooth phase $\psi$ are sampled separately, following the aforementioned hybrid sampling strategy. How many samples it is necessary to handle in order to compute the HFT of a given function depends on which of the two terms demands a larger sampling number; in the particular case considered here, taking into account the fact that the smooth phase is a simple spherical phase, and sampling the residual field $U\!\left ({\boldsymbol{\rho}}\right )$ necessitates resolving all the details in the house-shaped aperture, the sampling heft of the latter clearly trumps that of the former. Finally, to clarify the performance of the HFT in practical tasks, we record the CPU time at two specific positions for both FFT and HFT, shown in Fig. 3. It reveals that the small number of sampling points translates into a short computational time. HFT does contribute to reduce the numerical complexity and could have broad application prospects.

Fig. 3. Comparison of the numerical effort of the homeomorphic and the Fast Fourier transforms (HFT and FFT respectively, for short) for the case of a house-shaped field with a spherical phase, for different values of the numerical aperture (NA) of said spherical phase. The abscissa records increasing values of the NA, which correspond to decreasing values of the radius of curvature $R$ of the spherical phase. The ordinate corresponds to a logarithmic scale of the number of sampling points. The dashed line for the number of sampling points used to sample the function in the case of the HFT operation, and, the solid line for the number of sampling points necessitated by the computation of the FFT. Finally, the CPU time at two specific positions, where $R = 10$ mm and $R = 5$ mm, are pointed out.

Download Full Size | PDF

3.2 Fourier transform of a field with a simple aberrant phase

The concept of aberrations has a long and prominent history in the field of optical design, and is used, among other things, to characterize the shape of wavefront phases, usually as a means of identifying and quantifying in what way and how much that phase differs from a certain reference. Although there are alternatives, one of the most often encountered representations is that provided by the set of Zernike polynomials.

In this second set of experiments, we intend to study the influence of aberrant phases on the Fourier transform. The selected setup is similar to the one previously used in this work: an input plane wave is truncated by an aperture shaped as a sequence of letters spelling “Light”, after which aperture an aberrant phase is superimposed on the field. In this case, as before, the amplitude distribution will remain unchanged throughout the sequence of experiments (as shown in Fig. 4); it is the Zernike polynomial terms imposed on the field what we will vary to study the effect of different types of aberrations on the Fourier transform. We can then mathematically describe the smooth wavefront phase by means of its Zernike decomposition:

(17)$$\psi\!\left({\boldsymbol{\rho}}\right) = \psi^{\textrm{Zer}}\!\left({\boldsymbol{\rho}}\right) = k \sum_{m=0}^{M}\sum_{n=0}^{N}c^m_n Z^m_n\!\left(r, \theta\right) \textrm{,}$$

where $k=\frac {2\pi }{\lambda }\check {n}$, $\lambda$ being the wavelength, $r=\frac {\left |{\boldsymbol{\rho}}\right |}{\left |{\boldsymbol{\rho}}_\textrm {max}\right |}$ and $\theta = \arctan \left (\frac {y}{x}\right )$; $c^m_n$ denotes the coefficients of the corresponding Zernike term. For the experiment we have selected six types of lower-order Zernike phases, which we enumerate and describe in Table 2. They all have the mathematical characteristic of being quadratic polynomial terms so that the resulting mapping relation ${\boldsymbol{\rho}}\rightarrow {\boldsymbol {\kappa }}$ is an exactly linear transformation.

Fig. 4. Homeomorphic Fourier transform (FFT results are not shown, since deviation between HFT and FFT just 0.2 %) of a field with amplitude shaped in the form of a letter sequence spelling “Light”, in the cases of different aberrant phases. Panel (a) shows the amplitude distribution of the $E_x$ component in the space domain, common starting point to all results. Different aberrations are then superimposed on this amplitude, and its Fourier transform computed. The remaining panels present the $\tilde {E}_x$ component (i.e. the result of said Fourier-transform operation) for (b) positive defocus, (c) negative defocus, (d) positive $x$ astigmatism, (e) negative $x$ astigmatism, (f) positive $y$ astigmatism and (g) negative $y$ astigmatism.

Download Full Size | PDF

Table 2. Simulation parameters of the aberrant phase for the example presented in Section 3.2.

View Table | View all tables in this article

In Fig. 4, the result of the Fourier transform for different aberrations is presented. One can see that, regardless of the scaling, the Fourier-transform operation delivers either a rotation or a mirroring of the original pattern. It should be noted that, for all the cases gathered in Fig. 4, the homeomorphic Fourier transform would produce virtually the same results as the rigorous FFT, with $\sim {0.2}{\% }$ deviation when compared against the FFT reference.

Additionally, we discussed the appearance of distortion in the mapping in the case of the spherical phase, but it is worth mentioning that no such effect arises for the aberrations considered in this experiment. This is related to the fact, already pointed out above, that all the aberration terms we have selected produce a linear mapping relation between ${\boldsymbol{\rho}}$ and ${\boldsymbol {\kappa }}$, and the fact that the scaling factor of Zernike polynomials is symmetric for both the $x$ and $y$ dimensions.

3.3 Fourier transform of a field with a general, bijective, wavefront phase

We have used the previous two examples to demonstrate, respectively and in isolation, the application of the homeomorphic Fourier transform to fields exhibiting spherical phases or simple aberrant phases. In what follows we shall extend the case study to a more general scene: a spherical wave is truncated by a circular aperture with radius $a=12$ mm; at the plane where the aperture is located the spherical phase has a radius of curvature of $R = 6$ mm or, in other words, the distance from the point source to the plane of the aperture is 6 mm. The amplitude of the field behind this mask is presented in Fig. 5(a). The numerical aperture can be computed by the formula $\textrm {NA} = \frac {an}{\sqrt {a^2+R^2}} = 0.3$, a value for which the HFT proves sufficiently accurate. The Fourier transform of the initial field – i.e. with only spherical phase – is given in Fig. 5(b) which, in accordance with what we already know about the behavior of fields for which the HFT is an accurate choice, exhibits a very similar distribution to the original pattern of the field in the space domain. This case will serve as a reference for the following tests, where we will add different Zernike terms onto the initial spherical phase and observe their influence on the Fourier transform. We have selected six types of complex aberrant phases, so that one or a combination of several of them can be added to the starting-point case of a purely spherical phase.

Fig. 5. Homeomorphic Fourier transform (FFT results are not shown, since deviation between HFT and FFT just 0.2 %) of a field with a general bijective wavefront phase. Panel (a) should the amplitude distribution of the $E_x$ component in the space domain, common starting point to all results. A spherical wavefront with different aberrations is then superimposed on this amplitude, and the Fourier transform computed. The remaining panels present the $\tilde {E}_x$ component (i.e. the result of said Fourier-transform operation) for (b) a purely spherical phase, (c) spherical phase and trefoil $x$, (d) spherical phase and coma $x$, (e) spherical phase, trefoil $x$ and coma $x$, (f) spherical phase and tetrafoil $y$, (g) spherical phase and tertiary astigmatism $y$, (h) spherical phase, tetrafoil $y$ and tertiary astigmatism $y$, (i) spherical phase and pentafoil $x$, (j) spherical phase and tertiary spherical, and (k) spherical phase, pentafoil $x$ and tertiary spherical. The corresponding Zernike terms are indicated under each panel.

Download Full Size | PDF

We show the simulation results of this experiment in Fig. 5. We have organized them as follows: the second and third columns gather the Fourier-transform results for the cases where a single higher-order Zernike phase term has been added to the initial, purely spherical, phase. The fourth column corresponds to the Fourier transform of the same initial field, but this time the superimposed aberration term is a combination of the higher-order terms employed in the same row in columns two and three, so that the effect of both can be observed at the same time and in a combined manner. The specific simulation parameters can be consulted in Table 3. Figure 5 illustrates how the inclusion of different aberration terms and combinations thereof modify and reshape the pattern and amplitude distribution of the spectrum.

Table 3. Simulation parameters of the aberrant phase for the example in Section 3.3.

View Table | View all tables in this article

It is important to note that for all the examples presented in this set of experiments the result of performing the Fourier transform with the homeomorphic operation was always compared against that obtained with the rigorous FFT algorithm. The deviation remained bounded throughout at under 0.2 %. A direct implication of this is that the homeomorphic Fourier transform algorithm can be applied to fields with any wavefront phase, as long as it sustains the bijectivity condition. In other words, there is no additional restriction on the type or shape of the wavefront phase, apart from the requirement for a bijective mapping.

3.4 Fourier transform of a field with a non-bijective wavefront phase

Up to this point we have only discussed those cases where the wavefront phase of the test field produced a bijective mapping when subjected to the algorithm of the HFT. Last but not least, we would like to delve into the pragmatically quite interesting case of how to extend the application of the HFT to situations in which the given wavefront phase violates the bijectivity condition. The basic workflow of how to go about this task numerically was already outlined in Section 2.4. Here let us focus on a specific, tangible, instance of its application.

We take as the starting point for this experiment a field linearly polarized along $E_x$ with a slowly varying amplitude, as shown in Fig. 6(a). Its smooth wavefront phase is shown, numerically sampled, in Fig. 7.

Fig. 6. Fourier transform of a field with a non-bijective wavefront phase. Panel (a) shows the amplitude of the $E_x$ component in the space domain, panel (b) shows the result of the Fourier transform as obtained via the rigorous Fast Fourier Transform, and panel (c) shows the result of the Fourier transform as per the approximation provided by the homeomorphic Fourier transform with bijective regularization of the mapping relation.

Download Full Size | PDF

Fig. 7. Analysis of the wavefront phase for the simulation task: Homeomorphic Fourier transform of a field with a non-bijective wavefront phase. Panel (a) shows the initial smooth phase, panel (b) the residual phase after performing the fitting, panel (c) the phase from which the bijective mapping will be obtained, panel (e) serves to examine the homeomorphism condition for the fitted phase (the violation of the condition is evident from the change in sign of the second derivatives) and panel (f), analogously to (e), serves to examine how the same condition is fulfilled in the case of the final regularized phase employed to compute the bijective mapping.

Download Full Size | PDF

The next step, rather than straight away attempting to apply the HFT, is to examine whether the homeomorphism condition is fulfilled by the wavefront phase of the field in question. The mathematical criterion to establish or disprove this conformity rests on the sign of the second derivatives of the wavefront in the entire definition domain of the field: if the sign does not change, then the condition holds, the mapping relation can be obtained, and the HFT computed. If, on the other hand, the sign changes from positive to negative (or vice versa) at any point in the definition domain, the homeomorphism condition is violated and further measures need to be taken before the HFT algorithm can be applied.

This takes us to the practicalities of calculating, numerically, the point-wise second derivatives of the wavefront. There are several methods to attain this end, which can depend on the initial parametrization of the three-dimensional surface in question. For instance, a spline interpolation can provide the local gradient information. In this example, however, we choose to follow the proposed method of first performing a fitting for the linear phase, spherical phase, and any additional Zernike terms, in that order, starting from the initial numerical phase data. We consequently obtain an analytical expression, $\psi ^{\textrm {fit}}\!\left ({\boldsymbol{\rho}}\right )$, for the phase. For the specific case at hand, the result of this fit (expression and coefficients for the different terms) is presented in Table 4. An analysis of the quality of the fit can be extracted from Fig. 7(b), where we show the subtraction of the smooth, analytic function produced by the fit, $\psi ^{\textrm {fit}}\!\left ({\boldsymbol{\rho}}\right )$, from the initial, discretely sampled phase $\psi \!\left ({\boldsymbol{\rho}}\right )$. Direct comparison of the scale used here ($\sim 3 \times 10^{-6}$rad) and that used for the initial phase in panel (a), $\sim 970$rad, confirms the accuracy of the fit.

Table 4. Simulation parameters of the aberrant phase for the example in Section 3.4.

View Table | View all tables in this article

Let us now calculate the second derivatives of $\psi ^{\textrm {fit}}\!\left ({\boldsymbol{\rho}}\right )$ in a point-wise manner; the sign of the result is graphically shown, color-coded, in Fig. 7(e): while over most of the definition domain the second derivatives exhibit a positive sign, a small area in the bottom right is negative-valued. The homeomorphism condition is thus infringed, and the wavefront phase must be regularized before moving on with the computation of the homeomorphic Fourier transform.

We follow the elementary optimization approach presented in Section 2.4, whereby the highest Zernike term present in the fit will be successively removed from the wavefront part of the phase and the homeomorphism condition re-examined until it is fulfilled. In the specific example studied here, only one iteration is necessary, and it suffices to remove the secondary astigmatism in $x$, leaving us with the final expression for the wavefront to be employed in the computation of the HFT, $\psi ^{\textrm {map}}$, which fulfills the homeomorphism condition, as shown in panels (c) and (f) of Fig. 7.

We are finally in a situation to compute the now bijective mapping relation ${\boldsymbol{\rho}}\rightarrow {\boldsymbol {\kappa }}^{\textrm {map}}$ from the regularized wavefront phase $\psi ^{\textrm {map}}$, according to Eq. (5). Please note that the remainder of the phase after the part which can provide us with a bijective mapping has been extracted is not lost in the calculation of the HFT: the expression of the total field $V\!\left ({\boldsymbol{\rho}}\right )$ will be reorganized, according to Eq. (2), so that this remainder of the phase will be henceforth included in the $U\!\left ({\boldsymbol{\rho}}\right )$ term. The only effect of this regularization of the wavefront phase is on the mapping relation.

Once the mapping is obtained, we can finally perform the HFT operation on the test field. Panels (b) and (c) of Fig. 6 show the results of performing the FFT and HFT respectively on said test field. Visually, the similarity between them is evident, with some small differences appearing in the bottom right (tellingly, the same area where the violation of the homeomorphism condition arose). The deviation of the HFT with respect to the FFT reference is about 1 %. The closer look at the area of interest which is also presented in the figure reveals that the differences between the HFT and the FFT results stem from what could be described as a “diffraction-like” pattern in the FFT which is absent in the HFT case. This observation in the purely mathematical perspective we have maintained throughout this work could be obviated; however, when looking to apply the mathematical tool that is the HFT in an optics context this detail will prove immensely relevant: one of the main effects (if not the main effect) of employing the HFT in the simulation of an optical system is that it will neglect any diffraction effects which might be relevant at the plane where it is performed.

4. Conclusions

Though the numerical effort of the FFT is approximately linear on the sample number of the field component $V\!\left ({\boldsymbol{\rho}}\right )$ to be transformed, due to the “$2\pi$-wrapped” phase, the sample number $N\!\left ( V \right )$ can explode when the field presents a strong wavefront phase. We use a hybrid sampling strategy to deal with the smooth wavefront phase $\psi \! \left ({\boldsymbol{\rho}}\right )$ and the remnant of the field, $U\!\left ({\boldsymbol{\rho}}\right )$, separately. Since we avoid sampling the wrapped phase $\exp \left [ \textrm {i}\psi \left ({\boldsymbol{\rho}}\right )\right ]$, the computational effort is dramatically reduced, $N\!\left ( U \right ) \ll N\!\left ( V \right )$. To allow for the application of the hybrid sampling approach to the Fourier transform operation, we present the theoretical derivation and the validity conditions of the homeomorphic Fourier transform: an approximated algorithm best suited to fields with a strong smooth wavefront phase. The fundamental basis of the proposed algorithm is the method of stationary phase. To develop a robust algorithm for more practical tasks, we extend the homeomorphic Fourier transform to cover fields with a non-bijective wavefront phase. Finally, several numerical examples and the analysis of simulation results are shown to illustrate the promising potential of this approach.

5. Outlook

In Section 2.4, we proposed a robust algorithm and extended the validity of the HFT. However, as shown in the last example, the regularization of the mapping relation results in a dislocation in the spectrum position and leads to an inevitable error. So far, we have investigated the error numerically and compared it against a given accuracy threshold. But there is no alternative approach to deal with the Fourier transform of a field with a strong non-bijective wavefront phase efficiently. Next, we would like to pursue not only efficient but also rigorous solutions for that specific case.

A. Homeomorphism condition for the smooth phase function

Let us consider two subsets $X \subseteq \mathbb {R}^2$ and $K \subseteq \mathbb {R}^2$. A mapping function between two such subsets can mathematically described as $f: {\boldsymbol{\rho}} \in \textrm{X} \Rightarrow f \left ( {\boldsymbol{\rho}}\right ) = {\boldsymbol {\kappa }} \in K$, where ${\boldsymbol{\rho}} = \left (x, y\right )$ and ${\boldsymbol {\kappa }}=\left (k_x, k_y\right )$. One such mapping function is precisely the vector mapping relation effected by the homeomorphic Fourier transform, which establishes a link between the two domains via the gradient of the wavefront phase function,

(A1)$${\boldsymbol{\kappa}} = f\!\left({\boldsymbol{\rho}}\right) = {\boldsymbol{\nabla}}\psi\!\left({\boldsymbol{\rho}}\right)\,\textrm{.}$$

It is made evident in the body of the work above that to identify whether this vector mapping constitutes a homeomorphism or not is key when it comes to applying the homeomorphic Fourier transform. In this section of the appendix we shall analyze this condition in a more in-depth manner from the strictly mathematical point of view, with a view to ascertaining a valid, easy-to-check criterion.

Let us then begin from the topological definition of a homeomorphism: a continuous bijection whose inverse is also continuous. A bijection describes a function which pairs each element contained in the domain $X$ with exactly one element in the domain $K$ and which leaves no element unpaired. For one-dimensional cases, it has been proved by Kenneth George that a strictly monotone real function is bijective [23]. Mathematically, the strict monotonicity of a one-dimensional function entails its first derivatives being all non-negative or non-positive.

Working from the one-dimensional result, we investigate the extrapolation of the criteria to the two-dimensional case. First, let us reformulate the vector mapping functions in Eq. (A1) in terms of their scalar components, thus:

(A2)$${\boldsymbol{\kappa}} = \left(k_x, k_y\right) :=\left\{ \begin{array}{l} k_x = \frac{\partial}{\partial x} \psi\!\left({\boldsymbol{\rho}}\right) = \psi_x\!\left({\boldsymbol{\rho}}\right)\,\textrm{,} \\ k_y = \frac{\partial}{\partial y} \psi\!\left({\boldsymbol{\rho}}\right) = \psi_y\!\left({\boldsymbol{\rho}}\right)\,\textrm{.} \end{array} \right.$$

In what follows the derivation process is demonstrated for the example of a spherical function. As shown in Fig. 8, each of the two components of the mapping vector $\left (k_x, k_y\right )$ corresponds, respectively, to a parametric surface. It is evident that, in order to prove the bijection, we must ensure that for two different positions ${\boldsymbol{\rho}}_1 \neq {\boldsymbol{\rho}}_2$ in the definition domain $X$, the corresponding mapping vectors are also different, ${\boldsymbol {\kappa }}\!\left ({\boldsymbol{\rho}}_1\right ) \neq \kappa \!\left ({\boldsymbol{\rho}}_2\right )$.

Fig. 8. Spherical phase function $\psi \! \left ( {\boldsymbol{\rho}}\right )$ and its partial derivative functions.

Download Full Size | PDF

Let us then pick out a random, arbitrary position ${\boldsymbol{\rho}}_i$ from the domain $X$. The corresponding mapping vector ${\boldsymbol {\kappa }}_i = \left (k_{xi}, k_{yi}\right )$ can be obtained from the application of Eq. (A2). Next, fixing the $k_x$ component so that henceforth it remains constant, we obtain a 3D parametric curve which corresponds to the intersection of $\psi _x\!\left ({\boldsymbol{\rho}}\right )$ and the plane $\psi _x\!\left ({\boldsymbol{\rho}}\right ) = k_{xi}$. As can be seen in Fig. 9, it is possible to introduce a new variable $t_\textrm {min} < t < t_\textrm {max}$ to describe this parametric curve in the form ${\boldsymbol{\rho}}\!\left (t\right ) = \left [ x\!\left (t\right ), y\!\left (t\right )\right ] = \psi _x^{-1}\!\left (k_{xi}\right )$ where ${\boldsymbol{\rho}}\!\left (t_\textrm {min}\right )$ and ${\boldsymbol{\rho}}\!\left (t_\textrm {max}\right )$ are the two end points of the curve.

Fig. 9. (a) The blue curve corresponds to the intersection of the first derivative of the spherical phase function $\psi _x\!\left ({\boldsymbol{\rho}}\right )$ with the constant plane given by $\psi _x\!\left ({\boldsymbol{\rho}}\right ) = k_{xi}$. (b) Along the path ${\boldsymbol{\rho}}\!\left (t\right ) = \left [v\!\left (t\right ), y\!\left (t\right )\right ]$, we can find another parametric curve on the surface $\psi _y\!\left ({\boldsymbol{\rho}}\right )$.

Download Full Size | PDF

The curve ${\boldsymbol{\rho}}\!\left (t\right )$ is obviously a subset of the definition domain $X$, whose mapping vectors all have the same $k_x$ component. Therefore, if we are to have a bijection, it is inescapable that the $k_y$ component corresponding to the points of this curve be all different from each other, that is, that along the path ${\boldsymbol{\rho}}\!\left (t\right )=\left [x\!\left (t\right ), y\!\left (t\right )\right ]$ no two $k_y = \psi _y\!\left [{\boldsymbol{\rho}}\!\left (t\right )\right ]$ have the same value. Or, in other words, the one-dimensional function $\psi _y\!\left [{\boldsymbol{\rho}}\!\left (t\right )\right ]$ is strictly monotonic. This same test can then be applied to all the dense set of ${\boldsymbol{\rho}}$ paths which arise from considering all possible different fixed $k_{xi}$. If the condition holds for all points in $X$, then it can be concluded that the mapping in question is indeed a homeomorphism.

Before we move on to the rigorous demonstration, let us list some important arguments and summarize the most fundamental concepts discussed up to this point:

(i) A one-dimensional function is bijective if it is strictly monotonic, and that is true if and only iff the first derivatives are all non-negative or non-positive.
(ii) We say that $\psi \!\left ({\boldsymbol{\rho}}\right )$ is a smooth function when the derivatives (of all orders) are continuous and differentiable, i.e. when $\psi \!\left ({\boldsymbol{\rho}}\right )\in C^\infty$.
(iii) The condition $\psi ^2_{xy}\!\left ({\boldsymbol{\rho}}\right ) - \psi _{xx}\!\left ({\boldsymbol{\rho}}\right )\psi _{yy}\!\left ({\boldsymbol{\rho}}\right ) \neq 0$ implies that $\psi _{xy}\!\left ({\boldsymbol{\rho}}\right )$ and $\psi _{xx}\!\left ({\boldsymbol{\rho}}\right )\psi _{yy}\!\left ({\boldsymbol{\rho}}\right )$ cannot be 0 simultaneously.
(iv) The mapping from the scalar variable $t$ to the vector ${\boldsymbol{\rho}}$, $t\rightarrow {\boldsymbol{\rho}}\!\left (t\right )$, is bijective.

The $\psi _x$ component is a constant along the intersection parametric curve. Therefore, we have

(A3)$$\frac{\mathrm{d}\psi_x\!\left[x\!\left(t\right), y\!\left(t\right)\right]}{\mathrm{d}t} = \psi_{xx}\frac{\mathrm{d}x}{\mathrm{d}t} + \psi_{xy}\frac{\mathrm{d}y}{\mathrm{d}t} = 0\,\textrm{.}$$

Additionally, we can write

(A4)$$\frac{\mathrm{d}\psi_y\!\left[x\!\left(t\right), y\!\left(t\right)\right]}{\mathrm{d}t} = \psi_{xy}\frac{\mathrm{d}x}{\mathrm{d}t} + \psi_{yy}\frac{\mathrm{d}y}{\mathrm{d}t}\,\textrm{,}$$

so that we can substitute Eq. (A3) into Eq. (A4), thus:

(A5)$$\frac{\mathrm{d}\psi_y\!\left[x\!\left(t\right), y\!\left(t\right)\right]}{\mathrm{d}t} = \left\{ \begin{array}{cll} \frac{\psi^2_{xy}-\psi_{xx}\psi_{yy}}{-\psi_{xx}}\frac{\mathrm{d}y}{\mathrm{d}t} & , & \psi_{xx} \neq 0 \\ \psi_{xy}\frac{\mathrm{d}x}{\mathrm{d}t} & , & \psi_{xx} = 0 \end{array} \right.$$

As already established, bijection requires that $\frac {\mathrm {d}}{\mathrm {d}t}\psi _y$ be non-negative or non-positive throughout. To investigate how this condition applies to Eq. (A5), we consider three types of cases:

1. $\psi _{xx} = \forall {\boldsymbol{\rho}}\!\left (t\right )$
- • From Eq. (A3), we get $\frac {\mathrm {d}y}{\mathrm {d}t} = 0$. From this we can conclude that $y\!\left (t\right ) = \textrm {constant}$.
- • Considering condition (iv), the mapping $t\rightarrow x\!\left (t\right )$ is bijective.
- • From condition (i), $\frac {}{}$ is non-negative or non-positive in its entirety.
- • If $\psi _{xy}$ is negative or positive throughout, then $\frac {\mathrm {d}}{\mathrm {d}t}\psi _y$ is non-negative or non-positive throughtout.
2. $\psi _{xx}\neq 0$ and $\psi _{xy}\neq 0 \forall {\boldsymbol{\rho}}\!\left (t\right )$
- • From condition (ii), $\psi _{xx}$ and $\psi _{xy}$ are each either positive or negative throughout.
- • From Eq. (A3) we deduce that $\frac {\mathrm {d}y}{\mathrm {d}t} = \frac {\mathrm {d}x}{\mathrm {d}t} = 0$ or that $\frac {\mathrm {d}y}{\mathrm {d}t}$ and $\frac {\mathrm {d}x}{\mathrm {d}t}$ are positive or negative throughout.
- • Condition (iv) preempts the possibility of $\frac {\mathrm {d}y}{\mathrm {d}t} = \frac {\mathrm {d}x}{\mathrm {d}t} = 0$.
- • $\frac {\mathrm {d}y}{\mathrm {d}t}$ is positive or negative throughout ${\boldsymbol{\rho}}\!\left (t\right )$.
- • If $\psi _{xy}^2 - \psi _{xx}\psi _{yy}$ is either negative or positive throughout, $\frac {\mathrm {d}}{\mathrm {d}t} \psi _y$ is negative or positive throughout.
3. At some singular points $\psi _{xx}=0$ or $\psi _{xy}=0$, e.g. curves in Fig. 10
- • From condition (iii), we conclude that $\psi _{xy}\!\left ({\boldsymbol{\rho}}\right )$ and $\psi _{xx}\!\left ({\boldsymbol{\rho}}\right )$ cannot be zero simultaneously.
- • We can divide the curve ${\boldsymbol{\rho}}\!\left (t\right )$ into a series of fragments separated by the zero-valued points.
- • For each fragment, the situation is analogous to that analyzed in the second kind (point 2. above), so each fragment can be analyzed independently according to those criteria.
- • Because of the continuity – condition (i) – we can connect the result of any two neighboring fragments. We find that the term $\frac {1}{\psi _{xx}}\frac {\mathrm {d}y}{\mathrm {d}t}$ is non-negative or non-positive throughout.
- • If $\psi _{xy}^2-\psi _{xx}\psi _{yy}$ is either negative or positive throughout, then $\frac {\mathrm {d}}{\mathrm {d}t}\psi _y$ is either non-negative or non-positive throughout.

Fig. 10. Example of the third kind of smooth function phase. At $P_1$, $P_2$ and $P_3$ the second derivative of $\psi \!\left ({\boldsymbol{\rho}}\right )$ is equal to zero. Based on these three points, the definition domain is divided into four segments. Comparing the variation of the derivatives between two neighboring segments, we can analyze the continuity of the zero position.

Download Full Size | PDF

From the above analysis we can conclude that, for all three different kinds of situations, the homeomorphism condition is contained in the requirement that $\psi _{xy}^2-\psi _{xx}\psi _{yy}$ be either negative or positive throughout.

Funding

European Social Fund (2017 SDP 0049).

Acknowledgments

We gratefully acknowledge financial support by the LightTrans GmbH and the funding from the European Social Fund (ESF) (Thüringen-Stipendium Plus: 2017 SDP 0049).

Disclosures

The authors declare no conflicts of interest.

References

1. R. N. Bracewell, The Fourier transform and its applications (McGraw-Hill, 1986).

2. J. W. Goodman, Introduction to Fourier optics (Roberts and Company Publishers, 2005).

3. D. Asoubar, S. Zhang, F. Wyrowski, and M. Kuhn, “Efficient semi-analytical propagation techniques for electromagnetic fields,” J. Opt. Soc. Am. A 31(3), 591–602 (2014). [CrossRef]

4. S. Zhang, D. Asoubar, C. Hellmann, and F. Wyrowski, “Propagation of electromagnetic fields between non-parallel planes: a fully vectorial formulation and an efficient implementation,” Appl. Opt. 55(3), 529–538 (2016). [CrossRef]

5. A. V. Pfeil and F. Wyrowski, “Wave-optical structure design with the local plane-interface approximation,” J. Mod. Opt. 47(13), 2335–2350 (2000). [CrossRef]

6. M. Kuhn, F. Wyrowski, and C. Hellmann, Non-sequential optical field tracing, Advanced Finite Element Methods and Applications (Springer, 2013257–273

7. J. W. Cooley and J. W. Tukey, “An algorithm for the machine calculation of complex fourier series,” Math. Comp. 19(90), 297 (1965). [CrossRef]

8. K. Itoh, “Analysis of the phase unwrapping algorithm,” Appl. Opt. 21(14), 2470 (1982). [CrossRef]

9. D. J. Bone, “Fourier fringe analysis: the two-dimensional phase unwrapping problem,” Appl. Opt. 30(25), 3627–3632 (1991). [CrossRef]

10. L. Rabiner, R. W. Schafer, and C. Rader, “The chirp z-transform algorithm,” IEEE Trans. Audio Electroacoust. 17(2), 86–92 (1969). [CrossRef]

11. A. Andrianov, A. Szabo, A. Sergeev, A. Kim, V. Chvykov, and M. Kalashnikov, “Computationally efficient method for fourier transform of highly chirped pulses for laser and parametric amplifier modeling,” Opt. Express 24(23), 25974–25982 (2016). [CrossRef]

12. R. Tong and R. W. Cox, “Rotation of nmr images using the 2d chirp-z transform,” Magn. Reson. Med. 41(2), 253–256 (1999). [CrossRef]

13. Z. Wang, S. Zhang, O. Baladron-Zorita, C. Hellmann, and F. Wyrowski, “Application of the semi-analytical fourier transform to electromagnetic modeling,” Opt. Express 27(11), 15335–15350 (2019). [CrossRef]

14. P. Debye, “Das verhalten von lichtwellen in der nähe eines brennpunktes oder einer brennlinie,” Ann. Phys. 335(14), 755–776 (1909). [CrossRef]

15. G. C. Sherman and W. C. Chew, “Aperture and far-field distributions expressed by the debye integral representation of focused fields,” J. Opt. Soc. Am. 72(8), 1076–1083 (1982). [CrossRef]

16. L. Mandel and E. Wolf, Optical coherence and quantum optics (Cambridge university press, 1995).

17. J. Brüning and V. W. Guillemin, Mathematics Past and Present: Fourier Integral Operators (Springer, 1994).

18. J. J. Stamnes, “Waves, rays, and the method of stationary phase,” Opt. Express 10(16), 740–751 (2002). [CrossRef]

19. K. G. Larkin, “Natural demodulation of two-dimensional fringe patterns. ii. stationary phase analysis of the spiral phase quadrature transform,” J. Opt. Soc. Am. A 18(8), 1871–1881 (2001). [CrossRef]

20. R. Rolleston and N. George, “Stationary phase approximations in fresnel-zone magnitude-only reconstructions,” J. Opt. Soc. Am. A 4(1), 148–153 (1987). [CrossRef]

21. Physical optics simulation and design software “Wyrowski VirtualLab Fusion,” developed by Wyrowski Photonics UG, distributed by LightTrans GmbH, Jena, Germany.

22. F. Wyrowski and M. Kuhn, “Introduction to field tracing,” J. Mod. Opt. 58(5-6), 449–466 (2011). [CrossRef]

23. K. G. Binmore and K. G. Binmore, Mathematical Analysis: a straightforward approach (Cambridge University Press, 1982).

spherical radius	$R = \infty$	$R = 200$ mm	$R = 50$ mm	$R = 10$ mm	$R = 2.5$ mm	$R = 1$ mm
numerical aperture	NA $= 0$	NA $= 0.006$	NA $= 0.026$	NA $= 0.13$	NA $= 0.45$	NA $= 0.71$
sampling of FFT	$(59 \times 55)$	$(97 \times 95)$	$(281 \times 23)$	$(1299 \times 915)$	$(4557 \times 3323)$	$(7753 \times 6091)$
sampling of HFT	$(59 \times 55)$	$(59 \times 55)$	$(59 \times 55)$	$(59 \times 55)$	$(59 \times 55)$	$(59 \times 55)$
deviation	$99.99 %$	$17.17 %$	$4.1 %$	$0.28 %$	$0.067 %$	$0.021 %$

No.	Name & Type	Expression (Zernike)	Value
Fig. 4 (b)	defocus	$Z_{2}^{0} = \sqrt{3} (2 r^{2} - 1)$	$c_{2}^{0} = 100 λ$
Fig. 4 (c)	defocus	$Z_{2}^{0} = \sqrt{3} (2 r^{2} - 1)$	$c_{2}^{0} = - 100 λ$
Fig. 4 (d)	vertical astigmatism	$Z_{2}^{- 2} = \sqrt{6} r^{2} \cos 2 θ$	$c_{2}^{- 2} = 200 λ$
Fig. 4 (e)	vertical astigmatism	$Z_{2}^{- 2} = \sqrt{6} r^{2} \cos 2 θ$	$c_{2}^{- 2} = - 200 λ$
Fig. 4 (f)	oblique astigmatism	$Z_{2}^{2} = \sqrt{6} r^{2} \sin 2 θ$	$c_{2}^{2} = 200 λ$
Fig. 4 (g)	oblique astigmatism	$Z_{2}^{2} = \sqrt{6} r^{2} \sin 2 θ$	$c_{2}^{2} = - 200 λ$

No.	Name & Type	Expression (Zernike)	Value
Fig. 5 (c)	trefoil $x$	$Z_{3}^{3} = \sqrt{8} r^{3} \cos 3 θ$	$c_{3}^{3} = 300 λ$
Fig. 5 (d)	coma $x$	$Z_{3}^{1} = \sqrt{8} (3 r^{3} - 2 r) \cos θ$	$c_{3}^{1} = 100 λ$
Fig. 5 (f)	tetrafoil $y$	$Z_{4}^{- 4} = \sqrt{10} r^{4} \sin 4 θ$	$c_{4}^{- 4} = 150 λ$
Fig. 5 (g)	tertiary astigmatism $y$	$Z_{6}^{- 2} = \sqrt{14} (15 r^{6} - 20 r^{4} + 6 r^{2}) \sin 2 θ$	$c_{2}^{- 2} = 40 λ$
Fig. 5 (i)	pentafoil $x$	$Z_{5}^{5} = \sqrt{12} r^{5} \cos 5 θ$	$c_{5}^{5} = 150 λ$
Fig. 5 (j)	tertiary spherical	$Z_{8}^{0} = \sqrt{9} (70 r^{8} - 140 r^{6} + 90 r^{4} - 20 r^{2} + 1)$	$c_{8}^{0} = 30 λ$

name & type	fitting expression	fitting coefficients	$ψ^{map} (ρ)$
spherical phase	$ψ^{sph} = sign (R) k_{0} \overset{ˇ}{n} \sqrt{ρ^{2} + R^{2}}$	$R = 9.1$ mm	✓
piston	$Z_{0}^{0} = 1$	$c_{0}^{0} = 0.88 λ$	✓
tilt $y$	$Z_{1}^{- 1} = 2 r \sin θ$	$c_{1}^{- 1} = - 2.52 λ$	✓
tilt $x$	$Z_{1}^{1} = 2 r \cos θ$	$c_{1}^{1} = - 1.17 λ$	✓
defocus	$Z_{2}^{0} = \sqrt{3} (2 ρ^{2} - 1)$	$c_{2}^{0} = - 1.97 λ$	✓
astigmatism $x$	$Z_{2}^{2} = \sqrt{6} r^{2} \cos 2 θ$	$c_{2}^{2} = 27.6 λ$	✓
trefoil $y$	$Z_{3}^{- 3} = \sqrt{8} r^{3} \sin 3 θ$	$c_{3}^{- 3} = 6.14 λ$	✓
coma $y$	$Z_{3}^{- 1} = \sqrt{8} (3 r^{3} - 2 r) \sin θ$	$c_{3}^{- 1} = 4.5 λ$	✓
coma $x$	$Z_{3}^{1} = \sqrt{8} (3 r^{3} - 2 r) \cos θ$	$c_{3}^{1} = 6 λ$	✓
trefoil $x$	$Z_{3}^{3} = \sqrt{8} r^{3} \cos 3 θ$	$c_{3}^{3} = 3.84 λ$	✓
tetrafoil $y$	$Z_{4}^{- 4} = \sqrt{10} r^{4} \sin 4 θ$	$c_{4}^{- 4} = - 3.53 λ$	✓
spherical aberration	$Z_{4}^{0} = \sqrt{5} (6 r^{4} - 6 r^{2} + 1)$	$c_{4}^{0} = 0.001 λ$	✓
secondary astigmatism $x$	$Z_{4}^{2} = \sqrt{10} (4 r^{4} - 4 r^{2}) \cos 2 θ$	$c_{4}^{2} = 0.15 λ$	✗

spherical radius	$R = \infty$	$R = 200$ mm	$R = 50$ mm	$R = 10$ mm	$R = 2.5$ mm	$R = 1$ mm
numerical aperture	NA $= 0$	NA $= 0.006$	NA $= 0.026$	NA $= 0.13$	NA $= 0.45$	NA $= 0.71$
sampling of FFT	$(59 \times 55)$	$(97 \times 95)$	$(281 \times 23)$	$(1299 \times 915)$	$(4557 \times 3323)$	$(7753 \times 6091)$
sampling of HFT	$(59 \times 55)$	$(59 \times 55)$	$(59 \times 55)$	$(59 \times 55)$	$(59 \times 55)$	$(59 \times 55)$
deviation	$99.99 %$	$17.17 %$	$4.1 %$	$0.28 %$	$0.067 %$	$0.021 %$

Theory and algorithm of the homeomorphic Fourier transform for optical simulations

Abstract

1. Introduction

2. Theory

2.1 Field representation and hybrid sampling strategy

2.2 Derivation of the homeomorphic Fourier transform

2.3 Inverse homeomorphic Fourier transform

2.4 Pointwise Fourier transform in the case of a non-bijective phase

3. Application

3.1 Fourier transform of a field with spherical wavefront phase

3.2 Fourier transform of a field with a simple aberrant phase

3.3 Fourier transform of a field with a general, bijective, wavefront phase

3.4 Fourier transform of a field with a non-bijective wavefront phase

4. Conclusions

5. Outlook

A. Homeomorphism condition for the smooth phase function

Funding

Acknowledgments

Disclosures

References

Cited By

Figures (10)

Tables (4)

Equations (22)

Optics Express