Rays, beams and diffraction in a discrete phase space: Wilson bases

John M. Arnold

doi:10.1364/OE.10.000716

1 Introduction

A recent development in the computational modelling of imaging systems has been the use of superpositions of gaussian beams as propagators for the forward propagation of some specified distribution of field u(x,0) to predict the field at another observation space u(x,z). The attraction of this method is that the gaussian beams are expected to propagate along paths predicted by geometrical optics, but the spectral content of the gaussian beams fills in data which is missing from a purely geometrical-optics description of the propagating wavefield, and hence yields a representation which is uniform at caustics, foci, diffraction shadow boundaries and other nonuniformities of classical geometrical optics. Some specific recent examples of gaussian beam tracking and superposition in practical imaging configurations can be found in Refs. [1, 2, 3, 4]. The use of these beam superposition schemes raises a fundamental question about how many beams are needed to accurately synthesise a certain wavefield when its initial specification is purely in terms of geometrical optics phase data, bearing in mind that eventually the beam superpositions have to be evaluated numerically and this requires discrete representations that are finite. Hence there is a question of how compact such a representation can be. This question is considered in this paper, using some simple examples of finite aperture distributions that are typical of diffraction theory.

The general configuration envisaged here is exemplified by the following situation. An initial scalar field u(x) is specified over a plane which we define to be the plane z = 0 in coordinates x = (x, y, z) (or sometimes (x _⊥, z), where x _⊥ represents transverse coordinates and z is the axial coordinate). The sources of the field are all assumed to lie in z < 0. The field is to be propagated forward from its initial plane z = 0 to some z ≥ 0 as a solution of the scalar wave equation ∇² u + k ² u = 0, with boundary conditions u(x _⊥,0) prescribed, and u an outgoing wave as ∥x∥ → ∞ in z > 0. For simplicity here we shall omit one transverse dimension (the y-coordinate) since the modifications of the theory to 3-dimensions are purely formal, and introduce no particular difficulties.

A typical situation of this kind occurs when there is an aperture $- \frac{1}{2} d < x < - \frac{1}{2} d$ located in a screen placed in the initial plane z = 0, illuminated from behind the aperture in z < 0; suppose also that the illuminating wave is focussed at some point in front of the screen at x = x ₀, z = f > 0. In this scenario, one wishes to compute the field in the vicinity of the focus, including the diffraction effects which take place at the edges of the aperture. This might typically be accomplished by the Fourier transform method [5]

u (x, z) = {(2 π)}^{- 1} \int_{- \infty}^{\infty} \tilde{u} (ξ) e^{i β z} e^{i ξ x} d ξ,

\tilde{u} (ξ) = \int_{- \infty}^{\infty} u (x, 0) e^{- i ξ x} d x,

with $β = \sqrt{k^{2} - ξ^{2}}$ . Computationally, this requires one FFT for the evaluation of the spectrum ũ(ξ), and one inverse FFT for each z over which the field is to be computed. The propagated field may itself be modified by another aperture at some subsequent plane z = z ₁ > 0, and the result may be again described by the Fourier method. In this way, complex optical imaging systems are modelled which consist of an arbitrary number of parallel planes at z = z_n, n ∊ Ƶ, containing apertures, and one seeks to compute the effect of the combined diffraction and focussing that takes place at each aperture plane, followed by free-space propagation between the component planes.

An ancient method of predicting the behaviour of optical systems of this type is that of ray tracing, in which the propagation of a disturbance from a source through various focussing elements is described along trajectories which emanate from the source and arrive at the observation point after deflections through the focussing elements. This method describes the phase evolution of a wave, and the ray trajectories are normals to surfaces of constant phase. A fundamental difficulty encountered by classical ray-tracing is that it does not account for diffraction that can be observed at the edges of sharp apertures. This defect is corrected by the geometrical theory of diffraction [6], which predicts the addition of new ray trajectories emanating from edges, which are also traced to the observation point. It is also now well-understood that if the phase of the wave can be predicted by the ray-tracing method, then the entire physical field can in principle be deduced from the phase in the form of an asymptotic expansion in inverse powers of the free-space wavenumber k, formally proposed as the Kline–Luneberg expansion, although the origins of this expansion are very old. Further difficulties are encountered in the transformation of geometrical phase information to field information in nonuniformities of the Kline–Luneberg representation at caustics, foci, diffraction shadow boundaries and other types of singularities, which require significant mathematical transformations to regularise the field representation. Nevertheless, all transformations for correcting these deficiencies can be expressed (or parameterised) explicitly in terms of information which is already known from the geometrical construction of phase information along ray trajectories [7].

It follows that the information required for computation of the physical imaged fields is all encoded in the phase of the wave, which can itself be determined by relatively straightforward geometrical computations which are vastly more efficient numerically than, for example, the Fourier Transform method. In a sense, the conventional methods of physical optics are highly redundant. Recent research on computational physical optics has sought to release this redundancy, and to synthesize new methods which approach the computational efficiency of geometrical ray tracing without the singularity-corrections required by the conventional methods.

2 Geometrical optics, phase-space, and the Lagrange manifold

Consider a wavefront of a monochromatic wave arising from a phase function kσ(x, z) incident on the observation space z = 0, with free-space wavenumber k. The phase of the wave along the subspace z = 0 is described by the function kσ(x, 0), and the component of wavevector at x parallel to the space is k_x = k∂_xσ(x). The rays of geometrical optics are trajectories along wavefront normals; along rays the phase difference between two points on the same ray is ks, where s is distance along the trajectory between the two points. Since wavefronts are equiphase surfaces, the gradient of phase k ∇ σ is directed along rays with magnitude k. Hence the function k_x(x) defined above is just the x-component of this vector at each position along the space z = 0. In a 2-dimensional space with coordinates (x,ξ) the curve ξ = k_x(x) is the graph of the assignment of wavevector component k_x to each position x on the space z = 0. The space spanned by (x, ξ) is the phase-space, and the set of points in phase-space ξ = k_x(x) is the Lagrange submanifold Λ₀ for the particular wave being described. The importance of the Lagrange submanifold is that it carries all the geometrical information about the variation of phase of the wave on the subspace z = 0. The phase at any point in this subspace is kσ(x,0) = ∫^x k_x(x′)dx′. Each subspace at different values of z has a Lagrange submanifold Λ_z defined in the same way as Λ₀. The full Lagrange manifold Λ is the union of all submanifolds Λ_z as z varies.

Various imaging tranformations can be represented by symplectic (area-preserving) transformations in phase-space. The two most important are phase transformation in a lens (or a mirror), and propagation in free space. For a phase transformation, σ(x) → σ(x) + p(x) where p(x) is a polynomial in x; for an ideal thin lens, p(x) is quadratic: p(x) = -(x - a)²/2f with focal length f, and the centre of the lens is at x = a. Under the lens transformation, phase-space points transform by

S_{1} : x \mapsto x, ξ \mapsto ξ - k \partial_{x} p (x),

which is a shear transformation of the phase-space in the direction parallel to the ξ-axis. If p is of degree 2, this shear is linear. For transformation by propagation a distance z from the plane z = 0, the ray coordinates transform by

S_{2} : x \mapsto x - γ (ξ) z, ξ \mapsto ξ,

where $γ = \frac{ξ}{\sqrt{k^{2} - ξ^{2}}}$ , which is also a shear transformation, parallel to the x-axis. In the paraxial regime |ξ| ≪ k this shear is approximately linear. Under these symplectic transformations, which are defined on the whole phase-space, the Lagrange submanifold transforms to a new Lagrange submanifold, which is the locus in phase-space of the transformed ray system. In the case of propagation, the new Lagrange submanifold is just the Lagrange submanifold Λ_z for the rays at the observation space z. The use of symplectic imaging transformations for geometrical optics in phase-space is described by Dragt, Forest and Wolf [8].

A particular occurence in imaging systems that cannot be described by these simple transformations is that of aperturing, where an incident field illuminates a slit or hole in a screen. It is known from the Geometrical Theory of Diffraction (GTD) [6] that the proper ray description of this phenomenon involves the addition of a new set of rays, diffracted rays, emanating from apparent sources on the edge of the aperture. In the 2-D case that we are considering here, the phase-space picture of a set of rays all originating at a single spatial point x = d is a straight line parallel to the ξ-axis passing through the point (d, 0) on the x-axis where the rays originate. If the aperture is a slit of width d, there are two edges located at $x = \pm \frac{1}{2} d$ , and hence there are two parallel diffraction lines in phase-space. In addition, the incident field is truncated beyond the edges of the aperture in real-space, so the full phase-space picture of aperture diffraction involves only that part of the Lagrange manifold of the incident ray system lying between the two diffraction lines corresponding to the two edges, and terminated at these two lines. This entire phase-space system is propagated forward to some z > 0 by the shear transformation (4); there, it may be again subjected to lensing and apertured diffraction following the same rules, then propagated forward to the next component plane, and so on.

3 The windowed Fourier Transform

The Windowed Fourier Transform (WFT) is obtained from the standard Fourier Transform (1) by the introduction of a movable window, described by a window function g(x):

u (x) = {(2 π)}^{- 1} \int_{- \infty}^{\infty} \int_{- \infty}^{\infty} U (x_{0}, ξ_{0}) g (x - x_{0}) e^{i ξ_{0} x} d x_{0} d ξ_{0},

\tilde{u} (ξ) = {(2 π)}^{- 1} \int_{- \infty}^{\infty} \int_{- \infty}^{\infty} U (x_{0}, ξ_{0}) \tilde{g} (ξ - ξ_{0}) e^{- i (ξ - ξ_{0}) x_{0}} d x_{0} d ξ_{0}

U (x_{0}, ξ_{0}) = \int_{- \infty}^{\infty} u (x) g^{*} (x - x_{0}) e^{- i ξ_{0} x} d x

= {(2 π)}^{- 1} \int_{- \infty}^{\infty} \tilde{u} (ξ) {\tilde{g}}^{*} (ξ - ξ_{0}) e^{i ξ x_{0}} d ξ

with the condition

\int_{- \infty}^{\infty} {∣ g (x) ∣}^{2} d x = 1 .

Suppose that we have given a monchromatic wavefunction of geometrical optics type

u (x) ~ a (x) \exp [i k σ (x)],

with functions a(x) and σ(x) having variations on a typical length scale l, and suppose that kl is very large, where k is the free-space wavenumber. Then as position x varies, the oscillations of the function exp[ikσ(x)] go through very many cycles in the typical length scale l; one period of oscillation is a length λ(x) ~ 2π(k∂_xσ(x))^-1 which is generally much less than l. If the window scale parameter L is chosen much smaller than the length scale l but much larger than the local wavelength λ, then the amplitude a(x) and the phase σ(x) are approximately constant over the extent of the window function g(x - x ₀) but the function exp[ikσ(x)] goes through many oscillations. Near a fixed point x = x ₀ the phase function kσ(x) has a Taylor expansion

k σ (x) = k {σ (x_{0}) + (x - x_{0}) σ' (x_{0}) + \frac{1}{2} {(x - x_{0})}^{2} σ ″ (x_{0}) + . . .} .

Over the window’s extent, the function (10) can barely be distinguished from the Fresnel wave function

A \exp [i k_{x} (x_{0}) (x - x_{0}) + \frac{1}{2} i k σ ″ (x_{0}) {(x - x_{0})}^{2}]

where A = u(x ₀) and k_x = k∂_xσ(x)|_x=x₀, and therefore the WFT coefficient (7) is approximately

U (x_{0}, ξ_{0}) ~ u (x_{0}) \tilde{G} (x_{0}, ξ_{0} - k_{x} (x_{0})) .

with

\tilde{G} (x_{0}, ξ) = \int_{- \infty}^{\infty} g (x) \exp (\frac{1}{2} i k σ ″ (x_{0}) x^{2})) \exp (- i ξ x) d x .

If the window function g(x) is appropriately chosen, then this distribution is peaked sharply at the position (x ₀, k_x(x ₀)) in phase-space, which is located on the Lagrange submanifold determined by geometrical optics. The decay of the function G̃(x ₀,ξ ₀ - k_x(x ₀)) as ξ ₀ varies away from the point on the Lagrange manifold above the spatial point x ₀ is typically over the order of Δξ ₀ ~ L ^-1. The principle condition for this localisation behaviour to occur is that λ ≪ L ≪ l; a suitable choice for L to satisfy this criterion would be $L = \sqrt{l λ}$ . To summarise, under these conditions, the WFT coefficients of a function u(x) of geometrical optics type exhibit two fundamental properties: (i) asymptotic (λ → 0) localisation around the GO Lagrange submanifold in phase-space; and, (ii) determination by a local Fresnel wave (quadratic phase) approximation to the wavefunction u(x) near each window parameter value x ₀ such that the amplitude, phase and derivatives of phase of u are sampled at x = x ₀.

A peculiarity of the WFT representation for functions u(x) is that there are infinitely many coefficient functions U(x ₀, ξ ₀) other than (7) that will reconstruct any given function u(x). This arises from the fact that the basis functions for the representation are not generally linearly independent:

g (x - x_{1}) e^{i ξ_{1} x} = {(2 π)}^{- 1} \int_{- \infty}^{\infty} \int_{- \infty}^{\infty} K (x_{0}, ξ_{0}; x_{1}, ξ_{1}) g (x - x_{0}) e^{i ξ_{0} x} d x_{0} d ξ_{0},

K (x_{0}, ξ_{0}; x_{1}, ξ_{1}) = \int_{- \infty}^{\infty} g (x - x_{0}) g (x - x_{1}) e^{i (ξ_{0} - ξ_{1}) x} d x .

It follows that if the basis functions g(x) exp(iξx) in (5) are replaced by (14), then the coefficient functions transform to

U (x_{0}, ξ_{0}) \to {(2 π)}^{- 1} \int_{- \infty}^{\infty} \int_{- \infty}^{\infty} K (x_{0}, ξ_{0}; x_{1} {, ξ}_{1}) U (x_{1}, ξ_{1}) d x_{0} d ξ_{0} .

The transformed coefficients still represent the same function u(x). This transformation may be represented in operator form as

U \to U' = KU,

where K represents the operator on phase-space functions whose kernel is (15). It is not difficult to prove, by direct computation using (15), that K is a projection operator, satisfying

K^{2} = K .

The action achieved by the projector K is to project the equivalence class of all coefficient functions representing a given function u(x) onto the single representative coefficient function given by (7). Hence it follows that, if there be given a coefficient function computed by any method (not necessarily (7)) which represents a given function u(x), then that coefficient function can always be reduced to the canonical representative (7) by application of the projection operator K to it.

4 Examples of windowed Fourier Transform representations

In the following examples it is very convenient to adopt a length scale normalised to the width of the window in the windowed Fourier Transform. Suppose for example that the window function is chosen as the Gaussian function

g (x) = 2^{\frac{1}{4}} L^{- \frac{1}{2}} \exp (\frac{- π x^{2}}{L^{2}})

with some length scaling parameter L for the width of the window. Setting L = 1 is equivalent to a choice of length units for x in which x = 1 represents a real length unit of L. The unscaled variables can always be recovered from the scaled variables by the transformations x → x/L, g(x) → L ^-1/2 g(x/L), etc. With this scaling, the centres of the translated window functions are at positions x ₀/L = m. It is further convenient to scale wavenumbers to normalised spatial frequencies by ξ ₀ = 2πn/L. Then area in phase space scales by dx ₀dξ ₀ = 2πdmdn. These normalisations permit, by a slight abuse of notation, the writing of the WFT coefficient as U(m,n). Noncanonical WFT coefficients will be written with a prime, U′(m,n) and the canonical versions will be unprimed, U(m,n).

4.1 Fresnel wave Consider the Fresnel wave

Consider the Fresnel wave

u (x) = \exp (\frac{- i k x^{2}}{2 f}) .

with f < 0 for a diverging wave and f > 0 for a converging wave. The windowed Fourier Transform coefficients for this function are

U (m, n) = 2^{\frac{1}{4}} {(1 - i κ)}^{\frac{1}{2}} \exp [i κπ m^{2} - 2 i πmn - \frac{π}{1 - i κ} {(n - κ m)}^{2}]

with m = x ₀/L, n = ξ ₀ L/2π, κ = -kL ²/2πf. The magnitude of this function is very sharply concentrated around the line n = κm, or ξ ₀ = -kx ₀/f, in the phase-space. This line also defines local wavenumbers in the geometrical optics sense: the local spatial frequency at x = x ₀ of the phase variation of the Fresnel wave (20) is

ξ_{0} = \frac{d}{d x} (- \frac{k x^{2}}{2 f}) |_{x = x_{0}} = - \frac{k x_{0}}{f} .

The application of the projection operator can also be illustrated by this example. It can be shown by direct computation of eq. (7) that an alternative coefficient function for the Fresnel wave is

U' (m, n) = 2^{- \frac{1}{4}} {(1 + i κ)}^{\frac{1}{2}} \exp [i κπ m^{2} - 2 i πmn] δ (n - κ m)

where δ(x) is the Dirac delta function. The coefficient (23) obviously differs from the canonical coefficient function (21) by the contraction of the gaussian factor into a delta function, in addition to some rescaling of constant multipliers. It is a straightforward direct computation to show that the action of the projection operator K on the non-canonical coefficient function U′ in (23) is to reduce it to the canonical coefficient function U in (21). This example is very instructive, because it shows that a noncanonical coefficient function for the Fresnel wave can be located in phase-space entirely over the line n = κm, with zero values off this line. The position of this line is predicted solely by geometrical optics. In addition, the factor exp [iκπm ²] in (23) samples the phase of the Fresnel wave at the spatial point x = m, which is the position of the peak of the window function g(x - m) appearing in the reconstruction formula (5). Since all wavefunctions of geometrical-optics type are local Fresnel waves, the localisation of the WFT coefficients to the Lagrange submanifold is a generic property of this type. To summarise, in this example of the Fresnel wave it is clear that a noncanonical WFT coefficient function can be constructed entirely from geometrical optics, the reduction to the canonical coefficient function being carried out by the projection operator K.

4.2 Edge-diffracted waves

Consider an incident wave u ₀(x, z) generated by sources in z ≤ 0 illuminating a screen at z = 0 containing an aperture A, and seek the form of the total field in z ≥ 0 as predicted by the method of physical optics (PO). In the PO method, the aperture field is approximated by u(x,0) ~ u ₀(x, 0)χ_A(x), where

χ_{A} (x) = {\begin{matrix} 1, x \in A \\ 0, x \notin A \end{matrix}

is the characteristic function of the aperture.

Suppose that a gaussian beam u ₀(x, 0) = 2^1/4 exp[-π(x - a)² + 2πib(x - a)] is incident on a plane at z = 0, such that the waist of the beam is centred on the point x = a and the tilt of the beam is specified by the parameter b; in addition suppose that there is a hard screen placed in z = 0 so that the edge of the screen is at x = 0, and the screen lies in the half-space x ≤ 0. For this problem, it is known that a noncanonical WFT coefficient for the aperture field u(x, 0) in the physical optics approximation is [13]

U' (m, n) = \frac{1}{2 π i} \frac{1}{(n - b)} δ (m - a) \exp (- 2 π i ab)

(It is here assumed that b has a small positive imaginary part for the purposes of computing integrals over the spectral parameter n with respect to the pole in U(m, n) at n = b.) The corresponding result for diffraction of the same beam by a slit in $- \frac{1}{2} d \leq x \leq \frac{1}{2} d$ is

U' (m, n) = \frac{1}{2 π i} \frac{1}{(n - b)} δ (m - a) \exp [π i (n - b) d] \exp (- 2 π i ab)

More complex diffraction results can be built up from these elements; for example, the result for plane-wave diffraction through a slit is obtained by integrating (25) over the parameter a, after removing the factor exp(-2πiab), to obtain

U' (m, n) = \frac{1}{2 π i} \frac{1}{n - b} {\exp [π i (n - b) d] - \exp [- π i (n - b) d]}

Finally, the noncanonical WFT coefficient for the Fresnel wave diffracted through a slit is

U' (m, n) = \frac{1}{2 π i} \frac{1}{n - κ m} \exp [- i κπ m^{2}] \exp [π i (n - κ m) d]

A characteristic feature of these noncanonical WFT coefficients is that they are expressed in closed form, requiring no integrals to be carried out.

The canonical WFT coefficients corresponding to each of these noncanonical WFT functions are easily obtained by computing their projections with the operator K. (In fact, for all of these examples the canonical coefficient can be easily computed directly from (7) with u(x) = u ₀(x,0)χ_A(x).)

The recovery of the canonical coefficient is particularly instructive for the diffracted plane wave. To simplify, we consider just the case of a single edge at x = 0, in which case we retain only the first term of (26) with d = 0, and remove the second term. Thus we have a noncanonical coefficient

U' (m, n) = \frac{1}{2 π i} \frac{1}{n - b}

for the aperture field of a plane wave propagating in direction b diffracted at a single edge at x = 0. After applying the projection operator K, the canonical WFT coefficient for this case is

U (m, n) = 2^{\frac{1}{4}} \exp [- π {(n - b)}^{2} + 2 π i m (n - b)] \int_{- m - i (n - b)}^{\infty} \exp (- π s^{2}) d s .

Consider now this expression asymptotically in a number of limits.

m → ∞, |n - b| < m.
Then the integral in (29) tends to 1, and we have asymptotically
$U (m, n) ~ 2^{\frac{1}{4}} \exp [- π {(n - b)}^{2} + 2 π i m (n - b)], m \to \infty$
This function is strongly localised around the half-line {m > 0, n = b}, which is also part of the Lagrange submanifold for the incident plane wave.
m → -∞, |n - b| < |m|.
Then the integral in (29) tends to 0, and we have asymptotically
$U (m, n) ~ 0, m \to - \infty$
n → ±∞, |m| < |n - b|.
Then the integral in (29) tends to
$\int_{- m - i (n - b)}^{\infty} \exp (- π s^{2}) d s ~ - {(2 π (m + i (n - b)))}^{- 1} \exp [- π {(m + i (n - b))}^{2}],$
and we have asymptotically
$U (m, n) ~ - 2^{\frac{1}{4}} {(2 π i (n - b))}^{- 1} \exp [- π m^{2}], n \to \pm \infty .$
This function is strongly localised around the line m = 0, which is the Lagrange sub-manifold for the edge-diffracted (Keller) wave. The variation of the function in the factor (n - b)^-1 along the line m = 0 is essentially the Keller PO diffraction coefficient.
The function U(m, n) tends to zero with gaussian decay along any direction in the (m, n)-plane other than along the lines in the cases (i)–(iii) above. The four regions of the (m, n)-plane defined in cases (i)–(iii) are essentially the Stokes sectors of the analytic function erfc $(- \sqrt{2 π} ζ)$ with respect to the complex variable ζ, = m + i(n - b).

4.3 Geometrical construction of the canonical WFT coefficient

The results of the computation of the canonical WFT coefficients of the examples in this section lead to the following interpretation: the Lagrange manifold is a union of line segments (more generally, arcs of curves) in phase-space derived from the geometry of phases of the constituent geometrical optics components of the wave; the canonical WFT coefficient is approximated locally near each segment of the Lagrange manifold by functions each of which has a maximum on that segment of the Lagrange submanifold, with gaussian decay on either side of the Lagrange submanifold; if the tangent to the Lagrange submanifold at point (m ₀,n ₀) ∊ Λ₀ is the line n - n ₀ = κ(m - m ₀) then the decay of the canonical coefficient is greatest in the direction perpendicular to the tangent, and follows the factor exp[-πρ ²] where ρ = (1 + κ ²)^-1/2{(n -n ₀) - κ(m - m ₀)}. The variation along the Lagrange submanifold is determined by the phase and amplitude variations of the GO component described by each segment of the Lagrange submanifold. Points in phase-space where different segments of Λ₀ meet are transition points, where the phase-space variation of the WFT coefficient must be filled in by a special transition function. A similar construction in phase-space has also been proposed by Alonso and Forbes [9].

Using these rules, we now revisit the problem of aperture diffraction, this time to construct the canonical WFT coefficient for the problem of a converging (Fresnel) wave u ₀(x) = exp(iκπx ²) with κ < 0 in a finite aperture $A = {x : - \frac{1}{2} d \leq x \leq \frac{1}{2} d}$ . The Lagrange submanifold for the incident wave in this problem consists of the single line n = κm, from Example 4.1 above. The Lagrange submanifold for the diffracted aperture field u(x) = u ₀(x)χ_A(x) consists of the part of the line n = κm in the interval $- \frac{1}{2} d \leq m \leq \frac{1}{2} d$ , augmented by the diffraction lines $m = \pm \frac{1}{2} d$ . The segment of the line n = κm meets the two diffraction lines in the two points $n = \pm \frac{1}{2} κ d$ , but does not extend outside the diffraction lines. The incident field u ₀ has the WFT coefficient (21), so this is the variation that must be attached along the line segment n = κm in the interval $- \frac{1}{2} d \leq m \leq \frac{1}{2} d$ :

U (m, n) ~ 2^{\frac{1}{4}} {(1 - i κ)}^{- \frac{1}{2}} \exp [i κπ m^{2} - 2 i πmn - \frac{π}{1 - i κ} {(n - κm)}^{2}]

(- \frac{1}{2} d \leq m \leq \frac{1}{2} d) .

Along each diffraction line of the Lagrange submanifold $m = \pm \frac{1}{2} d$ there is attached a function similar to (32):

U (m, n) ~ - 2^{\frac{1}{4}} {(2 π i (n \pm \frac{1}{2} κ d))}^{- 1} \exp [- π {(m \mp \frac{1}{2} d)}^{2}], (m ~ \pm \frac{1}{2} d) .

It is assumed that the canonical WFT coefficient has negligible values other than these. There is also required a transition function which covers the region around each of the two points $(\pm \frac{1}{2} d, \mp \frac{1}{2} κ d)$ where different segments of the Lagrange manifold meet; these points are the origin in phase-space of the shadow boundaries that are well-known in real-space diffraction theory.

5 Wilson basis

Even the canonical WFT coefficients infinitely oversample the phase-space. For the numerical implementation of these WFT representations, the phase-space must be discretised in some manner. This is achieved optimally by the Wilson basis [11], which we shall now show samples the canonical WFT coefficients at the optimum rate of 1 coefficient per unit area of phase space. The Wilson basis is an orthogonal basis for L ²(R) whose coefficients are related to sampled canonical WFT coefficients. For an arbitrary function u(x) the Wilson expansion is

u (x) = \sum_{m \in Z} \sum_{n \in Z} A_{mn} w_{mn} (x)

with

w_{mn} (x) = {\begin{matrix} ϕ (x - m), n = 0 \\ \sqrt{2} ϕ (x - \frac{1}{2} m) \cos (2 πnx), m + n = 0 (\mod 2), n \neq 0 \\ \sqrt{2} ϕ (x - \frac{1}{2} m) \sin (2 πnx), m + n = 1 (\mod 2) . n \neq 0 \end{matrix}

Under certain conditions on the wavelet function ϕ(x) the function set {w_mn(x)} can be made to be orthonormal. The conditions for orthonormality are

\sum_{m \in Ƶ} ϕ^{*} (x - k - \frac{1}{2} m) ϕ (x - \frac{1}{2} m) = 2 δ_{0 k}, \forall k \in Ƶ, \forall x \in ℜ,

\sum_{m \in Ƶ} {(- 1)}^{m} ϕ^{*} (x - k - \frac{1}{2} m - \frac{1}{2}) ϕ (- x - \frac{1}{2} m) = 0, \forall k \in Ƶ, \forall x \in ℜ,

The second of these conditions, (38), can always be satisfied if ϕ(x) is a real even function, and we assume this from now onwards.

When the Wilson basis is orthonormal, there is a simple formula for the computation of the coefficients,

A_{mn} = \int_{- \infty}^{\infty} w_{mn}^{*} (x) u (x) d x .

After substituting the forms (36) for the Wilson basis functions, this reduces to

A_{mn} = {\begin{matrix} \int_{- \infty}^{\infty} ϕ^{*} (x - m) u (x) d x, n = 0 \\ \sqrt{2} \int_{- \infty}^{\infty} ϕ^{*} (x - \frac{1}{2} m) \cos (2 πnx) u (x) d x, m + n = 0 (\mod 2), n \neq 0 \\ \sqrt{2} \int_{- \infty}^{\infty} ϕ^{*} (x - \frac{1}{2} m) \sin (2 πnx) u (x) d x, m + n = 1 (\mod 2) . n \neq 0 \end{matrix}

These expressions are now clearly related to samples of the WFT coefficients when the window function for the WFT is ϕ ^*(x) (actually, ϕ is real-valued).

The method of constructing a Wilson basis [12] seeks to satisfy conditions (37)–(38), by a superposition of the form

ϕ (x) = \sum_{(r, s) \in Ƶ^{2}} C_{r s} g (x - \frac{1}{2} r) \exp (2 π i sx)

with a decaying function g(x) and suitable choices for the coefficients C_rs. This expansion is a frame expansion for ϕ, and is oversampled by a factor of 2 (there are twice as many coefficients required than are strictly necessary by sampling theorems). With the choice g(x) = (2ν)^1/4 exp(-νπx ²), it has proved possible to compute the coefficients C_rs numerically for selected values of ν [12]; they show good decay as the index parameters r and s change away from r = 0, s = 0. An expansion of the form (41) for ϕ implies an expansion of the form

u (x) = \sum B_{r s} g (x - \frac{1}{2} r) \exp (2 π i sx)

for the function u(x), which has the form of a 2-fold oversampled Gabor series, or frame expansion. When this form is substituted in the expressions for the Wilson basis (36) and the Wilson coefficients (40), it is deduced that

w_{mn} (x) = {\begin{matrix} \sum_{r, s} C_{r - 2 m s} g_{r s} (x), n = 0 \\ 2^{- \frac{1}{2}} \sum_{r, s} {C_{r - m s - n} + C_{r - m s + n}} g_{rs} (x), n \neq 0, m + n = 0 (\mod 2) \\ 2^{- \frac{1}{2}} \sum_{r, s} {C_{r - m s - n} - C_{r - m s + n}} g_{rs} (x), n \neq 0, m + n = 1 (\mod 2) \end{matrix}

u (x) = \sum_{r, s, r', s'} U_{r s} T_{r s r' s'} g_{r' s'} (x)

with

g_{r s} (x) = g (x - \frac{1}{2} r) \exp (2 π i sx)

U_{r s} = \int_{- \infty}^{\infty} u (x) g_{r s}^{*} (x) d x

T_{r s r' s'} = \sum_{m = - \infty}^{\infty} \sum_{n = 0}^{\infty} D_{mn r s}^{*} D_{mn r' s'}

D_{mn r s} = {\begin{matrix} C_{r - 2 m s,} n = 0 \\ 2^{- \frac{1}{2}} (C_{r - m s - n} + {(- 1)}^{m + n} C_{r - m s + n}), n \neq 0 \end{matrix}

The coefficients U_rs in (46) are sampled values of the canonical WFT coefficients of the function u(x), with the samples being taken at the phase-space points $(\frac{1}{2} r, s)$ with integer values of r and s. It follows that the coefficients of the orthogonal Wilson expansion are linearly related to samples of WFT coefficients on a rectangular mesh in phase-space.

6 Conclusions

The construction of canonical Windowed Fourier Transform (WFT) coefficients has been demonstrated for several simple diffraction problems. The canonical coefficients exhibit strong localisation near the Lagrange submanifold, which is the locus in phase-space of elementary geometrical optics, and which encodes all the information about the phase of a wave on the reference space (z = 0 in this case). Canonical coefficients exhibit gaussian decay of the coefficients for parameter values (m, n) away from the Lagrange manifold. In addition, a variety of noncanonical coefficients can be obtained which have simpler analytical form than the canonical coefficients, and which can be transformed to the canonical coefficients by the application of a projection operator; the canonical coefficients of a distinct wavefunction u(x) are unique.

The WFT coefficients oversample the phase-space, since they are continuous distributions on the phase-space, but a minimally specified set of data should have 1 coefficient per unit area in phase-space. The minimal sampling of the wavefield in phase-space is achieved in a Wilson basis, which is an orthogonal basis whose coefficients are linearly related to samples of the canonical WFT coefficients on a rectangular grid which has twice the density of the critical sampling grid.

The orthogonality of the Wilson basis is a very strong property. Under forward propagation the Wilson basis elements change, since their component gaussian beams evolve, but they remain orthogonal at any z. Parseval’s Theorem applies, so that the power transported by any finite subset of the basis elements is just the sum of the square magnitudes of the coefficients of the subset. Orthogonal expansions do not exhibit the instabilities typical of, for example, critically sampled Gabor expansions [2, 10, 14]. Their coefficients are determined by simple projection integrals (39), in contrast to the highly complicated bi-orthogonal projections associated with critically sampled Gabor expansions. In short, Wilson bases provide the minimal discrete basis parameterised in phase-space for expansion of wavefunctions of geometrical optics type, and their orthogonality guarantees the minimal mean-square error when the basis set is truncated to a finite one consisting of a few elements whose ϕ-functions in (36) are parameterised by positions and frequencies near the Lagrange manifold for the function being expanded.

This paper has concentrated only on the representation of a given initial field at z = 0 by a superposition of orthogonal functions that are discretely parameterised in phase-space. The field is propagated forwards to another observation space by propagation of the individual constituent window functions ϕ(x - m) exp[±2πinx]. This in turn can be accomplished by the forward propagation of the constituent gaussian windows g_rs(x) in the frame superposition (44); since propagation between two parallel planes in free space is a unitary operation, the forward propagated Wilson basis is again an orthogonal Wilson basis at the new observation plane, because

w_{mn} (x, 0) \to w_{mn} (x, z)

where the arrow represents forward propagation by a distance z and the tilde represents the spectral domain conjugate to the spatial domain of x. The propagation of these window functions occurs essentially along the ray trajectories of geometrical optics. Thus, a new Wilson basis is constructed at the forward observation space; the coefficients of the wavefield u(x, z) at the new observation space in the new Wilson basis w_mn(x, z) are identical to those of the initial field in the initial Wilson basis at z = 0.

A fundamental advantage of the proposed method over conventional ray asymptotics, apart from its potential numerical efficiency, is that no special transition functions are needed at nonuniformities of the GO representation of the field, such as caustics, when these develop along ray trajectories: the Wilson representation is itself the transition function. This feature is more significant in 3-dimensional field tracking, where the catalogue of canonical caustics (catastrophes) is quite large [15], and each nonuniformity of the field has to be classified as one of the canonical catastrophes in order to efficiently compute an asymptotic transition function; this is obviated in the localised discrete-basis approach described here.

References and links

1. H. T. Chou, P. H. Pathak, and R. J. Burkholder, “Novel guassian beam method for the rapid analysis of large reflector antennas,” IEEE Trans. Ant. Prop. 49, 880–893 (2001). [CrossRef]

2. D. Lugara and C. Letrou, “Alternative to Gabor’s representation of plane aperture radiation,” Elect. Lett. 34, 2286–2287 (1998). [CrossRef]

3. V. Galdi, L. B. Felsen, and D. Castanon, “Quasi-ray Gaussian beam algorithm for time-harmonic two-dimensional scattering by moderately rough interfaces,” IEEE Trans. Ant. Prop. 49, 1305–1314 (2001). [CrossRef]

4. C. Rieckman, M. R. Rayner, and C. G. Parini, “Diffracted Gaussian beam analysis of quasioptical multi-reflector systems,” Elect. Lett. 36, 1600–1601 (2000). [CrossRef]

5. J. W. Goodman, Introduction to Fourier optics (McGraw-Hill, New York, 1968).

6. G. L. James, The geometrical theory of diffraction for electromagnetic waves (IEEE Electromagnetic Waves Series, vol. 1, 3d ed., Peter Peregrinus Ltd., 1986).

7. J. M. Arnold, “Geometrical methods in the theory of wave propagation: a contemporary review,” IEE Proc. J 133, 165–188 (1986).

8. A. J. Dragt, E. Forest, and K. B. Wolf, “Foundations of a Lie-algebraic theory of geometrical optics,” Chap. 4 in Lie methods in optics, ed. J. Sanchez Mondragon and K. B. Wolf (Springer Lecture Notes in Physics, vol. 250, Springer-Verlag, Berlin, 1986). [CrossRef]

9. M. A. Alonso and G. W. Forbes, “Phase space distributions for high-frequency fields,” J. Opt. Soc. Am. A 17, 2288–2300 (2000). [CrossRef]

10. D. Gabor, “Theory of communication,” Jour. IEE (III) 93, 429–457 (1946).

11. K. G. Wilson, “Generalised Wannier functions,” Cornell University preprint (1987).

12. I. Daubechies, S. Jaffard, and J. L. Journe, “A simple Wilson orthonormal basis with exponential decay,” SIAM Jour. Math Anal. 22, 554–572 (1991). [CrossRef]

13. J. M. Arnold, “Phase space localisation in high-frequency scattering using elementary wavelets,” URSI General Assembly, Kyoto, Japan, Paper B8-2, p60, (1993).

14. J. M. Arnold, “Phase-space localisation and discrete representations of wave fields,” Jour. Opt. Soc. Am. A 12, 111–123 (1995). [CrossRef]

15. J. F. Nye, Natural focusing and fine structure of light (IOP Publishing, Ltd., Bristol, 1999).

Rays, beams and diffraction in a discrete phase space: Wilson bases

Abstract

1 Introduction

2 Geometrical optics, phase-space, and the Lagrange manifold

3 The windowed Fourier Transform

4 Examples of windowed Fourier Transform representations

4.1 Fresnel wave Consider the Fresnel wave

4.2 Edge-diffracted waves

4.3 Geometrical construction of the canonical WFT coefficient

5 Wilson basis

6 Conclusions

References and links

Cited By

Equations (56)

Optics Express