EURASIP Journal on Applied Signal Processing 2005:10, 1520–1534 c ○ 2005 Hindawi Publishing Corporation Applications of the Wigner Distribution Function in Signal Processing

We present a review of the applications of the Wigner distribution function in various areas of signal processing: amplitude and phase retrieval, signal recognition, characterization of arbitrary signals, optical systems and devices, and coupling coefficient estimation in phase space. Although reference is made to specific signals and systems, the mathematical formulation is general and can be applied to either spatial, temporal, or spatio-temporal phase spaces, to coherent, partially coherent, or discrete signals. The universal and intuitive character of the Wigner distribution approach to signal characterization and processing and its simplicity in solving many issues are evidenced throughout the paper.


INTRODUCTION
Phase space methods become increasingly exploited in signal processing applications due to their intuitive character, universal validity, and last but not least simplicity in an increasing number of practical situations. These methods refer to any description of a signal or an optical system through a function that depends jointly on the canonical conjugate phase space variables, which are the transverse position vector r and the angular (spatial frequency) vector k for light beams, the time variable t, and frequency ω for optical pulses or, more generally, all four coordinates r, k, t, and ω. There are many phase space distribution functions that can be employed in signal processing applications; the spectrogram, the ambiguity function, or the Wigner distribution function (WDF) are just a few examples. These phase space distribution functions are actually closely related to one another, the use of one or another depend to a great extent on the specific application and/or on the type of experiments that are carried on.
This paper focuses on the applications of the WDF in signal processing; this phase space distribution function is often referred to as Wigner-Ville distribution, but the WDF term is used in a larger range of applications and therefore the term WDF will be used throughout this paper. The intention is not to present an exhaustive list of WDF properties (including its relation to other phase space distributions) or problems that can be tackled with this distribution function; there are already a large number of review papers and even books dedicated to these subjects [1,2,3,4,5,6,7,8,9]. Rather, it will be shown how to make use of the specific properties of the WDF and why it is desirable to use this phase space distribution in addressing some relevant issues in signal processing. Based on the available mathematical and experimental results, the paper offers a glimpse on the beauty, simplicity, and potentialities of the WDF approach to this research field; therefore, many details will be left aside and only the most significant results will be discussed in greater detail. The responsibility of doing justice to the so many briefly mentioned issues is entrusted to the references.

THE WIGNER DISTRIBUTION FUNCTION AND SOME OF ITS PROPERTIES
The WDF has been originally introduced in quantum mechanics [10], but found itself rapidly in hosts of applications in the analysis and processing of both light beams and pulses [1,2,3,11,12]. The WDF of coherent signals that propagate along the z direction and that are characterized by a scalar field distribution ϕ(u, z) is defined at a z = const. plane as W ϕ (u, v; z) where * and T denote the operations of complex conjugation and transposition, respectively, and ϕ(v, z) = (2π) −n/2 ϕ(u, z) exp iuv T du (2) is the Fourier transform of ϕ(u, z). (For vectorial field distributions, a WDF matrix can be introduced, as shown in [13,14].) In the expressions above, u represents the transverse coordinate vector r = (x, y) in the spatial phase space appropriate for describing stationary and monochromatic light beams, symbolizes in the temporal phase space the time coordinate t in a reference frame that travels with the pulse, or designates the ensemble (r, t) for a general spatio-temporal phase space. The coordinate v is its canonically conjugate variable equal, respectively, to the projection k = (k x , k y ) of the wavevector on the transverse coordinates (known also as spatial frequency), to the frequency ω measured from the central frequency ω 0 of a light pulse, or to (k, ω) in general; n denotes the dimensionality of u and v (the phase space is 2n-dimensional). The parameter z in the WDF definition can be removed for simplicity if no ambiguity arises. It should be mentioned that in literature slightly different forms of the WDF are sometimes encountered, defined either with an opposite sign in the Fourier integral or even with scaled spatial frequency or frequency coordinates; due account of these differences must be taken when comparing the formulas in this paper with results derived with a different WDF definition.
In the paraxial approximation k = (kφ x , kφ y ), where k = 2π/λ is the wavenumber of the monochromatic field with wavelength λ, and φ x , φ y are the angles between the wavevector direction and the z-axis; thus, for monochromatic fields, k has the significance of a scaled angular coordinate. The WDF characterization of fields in the temporal phase space is sometimes referred to as the chronocyclic representation of light pulses [15]; it can adequately describe the propagation of signals when their spatial extent is of no concern. The use of spatio-temporal phase spaces is justified when optical systems, such as prisms or diffraction gratings couple the space and time coordinates through a frequencydependent refraction phenomenon, due to which different frequency components follow different paths; this situation is common for ultrashort light pulses characterized by a broad spectrum.
Extremely important for a large number of applications is that the WDF can be defined also for partially coherent light distributions such as [16] where the correlation function Γ(u 1 , u 2 ) = ϕ(u 1 )ϕ * (u 2 ) is the ensemble average of ϕ(u 1 )ϕ * (u 2 ); similarly to (1), the Fourier transform of the correlation function can also be used to define the WDF for partially coherent fields.
Digital signal processing applications can benefit from a discrete form of the WDF, which is defined in the temporal phase space, for example, as [17] where ϕ(t) is a complex-valued function that is sampled with a sampling interval τ, n, and l are discrete and ω is continuous. Note that the discrete-time WDF (4) is periodic in the frequency coordinate, with a period equal to one-half the sampling frequency ∆ω = 2π/τ, that is, W ϕ (nτ, ω) = W ϕ (nτ, ω + ∆ω/2). The practical considerations regarding the sampling rate that is needed to avoid aliasing are detailed in [17], which also investigates an efficient algorithm to implement in real time the temporal WDF. In an analogous manner a discrete-frequency WDF can be defined [18] through a formula similar to (4) but referring to the Fourier transform of the signal; a change in sign in the exponential term is also required. A WDF discrete in both n and θ variables, with properties similar to that of the continuous WDF in (1), has been introduced by Brenner [19]. The discretetime/frequency WDF is defined as where it is assumed that the signal can be approximated by a limited number N of samples, the values of the WDF being evaluated for t = nτ/2 and ω = m/Nτ (the product between the sample rates in time and frequency is equal to N −1 ). The discrete-time/frequency WDF is periodic in both variables with a period 2N, its properties and its relation to the continuous WDF being detailed in [18]. Irrespective of its definition, the WDF satisfies some properties that are tremendously important in practical applications. For the case of coherent light distributions these properties are as follows.
(i) The WDF is real, but can take negative values. Only Gaussian fields have WDFs that are everywhere positive. For a forward propagating field the negative WDF values were associated with small regions of backward flux, the local flux of the monochromatic field being defined as J(r) = k −1 W(r, k)k dk [20]. According to another interpretation [21], the negative values of the WDF originate from the phase space interference between neighboring minimumvalue phase space areas that correspond to the Gaussian beams into which an arbitrary field distribution can be decomposed.
(ii) The WDF is limited to the u and v intervals that limit the field distribution and its Fourier transform, respectively.
(iii) A shift of the field distribution in u or a shift of its Fourier transform in v leads to the same shift for the WDF.
(iv) The partial integral of the WDF over v is proportional to the field intensity: (v) The WDF integral over u is proportional to the intensity of the Fourier transform of the field distribution (the spectrum of the field): (vi) The field distribution and its Fourier transforms can be reconstructed from the WDF up to constant phase factors, according to with 0 the null vector, since the modulus of ϕ * (0) andφ * (0) can be found from (8) and (9), respectively, for u = 0 and v = 0.
(vii) The WDF satisfies the so-called overlap or Moyal relation For partially coherent sources with correlation functions Γ 1 (u, u ), Γ 2 (u, u ), the Moyal relation becomes which reduces to (10) The WDF is bilinear in the field distribution, which means that the WDF of a superposition of signals ϕ(u) = m a m ϕ m (u) contains not only the WDF of the individual signals (the auto-terms), W ϕ m (u, v) but also complex crossterms W ϕmϕn (u, v) that describe phase space interference effects: where W ϕmϕn (u, v) The interference terms in the WDF have an oscillatory behavior if the field distributions are coherent, and show no oscillations when the superposing beams are incoherent [22].
The oscillations indicate the eventuality of interference between the superposing states in the u and/or v domains, interferences in the u and v spaces, respectively, taking place only when the superposing field distributions have common projections along these axes. Note that the cross-terms in the WDF appear even when the individual WDFs do not overlap. The presence of the interference terms in the WDF of a multicomponent signal hampers in many situations the interpretation of the WDF and the reconstruction or synthesis of the field distribution. The geometry of the interference terms of the WDF and related signal representations in the temporal phase space is analyzed in [23]. Other properties of the WDF for coherent distributions can be found in [24], while the corresponding properties of the WDF for partially coherent distributions and for the discrete WDF are detailed in [16,19], respectively. An attractive property of the WDF is that its transformation law through first-order optical systems characterized in the geometrical optical approximation by a real symplectic matrix relating the ray vectors (uv) T at the output z = const. plane of an optical system to those at the input z = 0, is very simple [24,25,26]: A matrix is symplectic if it satisfies the relation S T JS = J with where 0 and I are the n-dimensional null and identity matrices, respectively; for n = 1 the symplecticity condition reduces to det S = 1. The WDF embodies the link between geometrical optics and wave optics: it remains constant along geometrical ray paths, but incorporates wave optical effects when considering diffraction and interference phenomena since the WDF is defined in terms of the field distribution. Note that (15) is valid not only for the better known case of optical systems for monochromatic beams but also for optical systems for light pulses; in particular the propagation of the WDF of pulses through temporal lenses [27] and dispersive media (the analogs of free spaces for monochromatic beams [28]) is also accounted for by (15). The influence of ray-pulse matrices [29] on the WDF is also described by (15). A similar law can be derived for each temporal frequency component of a polychromatic-pulsed paraxial beam propagating through nondispersive and achromatic optical systems [30] (see also [31] for the definition of the WDF for polychromatic paraxial fields), or for light propagation in misaligned optical systems [32]. Dragoman [33] has demonstrated that the ABCD transformation law for the WDF (15) holds also for particular optical systems characterized by complex S matrices when illuminated with partially coherent Gaussian-Schell beam; in this case, however, the A, B, C, D elements of the S matrix in (15) depend also on the incident field. The WDF transformation law in homogeneous, inhomogeneous, and weakly dispersive media has been derived by Bastiaans [34,35], the corresponding transformation in nonstationary and inhomogeneous optical media can be found in [36], while the evolution of the WDF through optical systems characterized by complex matrices is given in [37]. It is also possible to derive the propagation law of the temporal WDF of pulsed plane-wave signals in continuous random media [38]. All these propagation laws, although very useful in some applications, are overcome in popularity by the ABCD transformation law of the WDF expressed in (15). From (6), (7), and (15), it follows that the integral of the WDF over both u and v coordinates is proportional to the total energy content of the field distribution and that the energy is conserved at propagation through a first-order optical system characterized by a real and symplectic matrix. This conclusion no longer holds for propagation through optical systems characterized by a symplectic matrix with complex elements [39]; symplecticity of optical systems is not equivalent to energy conservation, with the notable exception of real matrices.
In the derivation of the elements of the symplectic matrix, S the reflection at various interfaces is neglected. An explicit accommodation of these reflections for layered structures is known to introduce a matrix that relates not the phase space coordinates but the field distribution and its derivative along the stratification direction [40]. A matrixlike relation exists even in this case between the vector WDF at the input and output plane, the elements of the WDF vector being the WDF of field, of its derivative, and the crossterms [41]; the WDF moments at the output plane are related to those at the input plane through a more complicated relation than (15).
The simplicity of the WDF transformation law through first-order systems reflects itself in the simplicity of the transformation laws of the WDF moments. For a two-dimensional stationary and monochromatic light distribution ϕ(x, y), for example, the moment of order l + m + n + p is defined as and transforms through first-order systems according to [42] In the expression above is the matrix that contains all moments of order j calculated at a z = const. plane, with ⊗ denoting the direct product of matrices and with the bar on top indicating that each matrix element is averaged in the sense of (17); M j (0) has the same significance but with the moments calculated at the input z = 0 plane of the optical system. Numerous moment invariants of several orders can be derived starting from (18) (see [6] and the references therein). In particular, (det M 2 ) 1/2 , which is invariant at propagation through first-order optical systems, can be identified with the phase space area occupied by the system. Its invariance can be viewed as a selection rule for various optical systems, such as the astigmatic mode converter consisting of a pair of properly oriented cylindrical lenses that transforms a Hermite-Gauss mode with rectangular symmetry into a Laguerre-Gauss mode with cylindrical symmetry [43]. Note that, due to (6), the moments of the spatial coordinates in (17) (those with n = p = 0), calculated with the WDF as the weighting function, are identical to the moments of spatial coordinates calculated with the field intensity as the weighting function, while the moments of the angular coordinates (with l = m = 0) in (17) are identical to the moments calculated with the field spectrum as weighting function, in agreement with (7). The moments defined in (17), with the WDF as weighting function, allow however a simplified calculation procedure for mixed, that is, spatial and angular, coordinates.
The propagation law (18) for WDF-based moments can be extended in certain cases even to optical systems that are not of the first-order type. For example, a modified definition of second-order moments for hard-edge diffracted fields leads for j = 2 to a propagation law similar to (18) even if hard-edge apertures are not first-order systems [44]. Analogously, it was demonstrated that pure phase transmittances can be represented as 4 × 4 matrices and a law similar to (18) has been derived for the second-order moments of the WDF [45].

APPLICATIONS TO SIGNAL RETRIEVAL
The complex field distribution can be determined, according to (8), from its WDF; the recovery of both amplitude and phase distributions of a signal is one of the most important tasks in signal processing since many detectors are only sensitive to the incident intensity. The WDF of onedimensional field distributions can be easily displayed on a two-dimensional screen using common optical elements. Optical set-ups that generate the spatial WDF of real signals, complex transparencies (holograms), or arbitrary complex one-dimensional fields are described in [46]; the product between the two displaced replicas of the field distribution that enter the definition of the WDF in (1) is implemented using rotated transparencies or holograms through which the light passes twice (but with a transverse coordinate inversion between passages) or using a rotated object illuminated by a tilted light beam, the subsequent Fourier transform being generated with the help of lenses and free space propagation regions. More complex set-ups involving joint transform correlator architectures [47] or incorporating a degenerate phase conjugation device for the implementation of the product between the field and its complex conjugate (and shifted) replica [48] can also be imagined. The latter set-up can be employed as well to generate samples of the WDF of two-dimensional complex signals. A white-light implementation of the WDF of one-dimensional field distributions containing an achromatic processor that consists of linear blazed zone plates and an achromatic cylindrical objective was proposed in [49], while acousto-optic processors for real-time implementation of temporal WDFs are described in [50]. Note that, since the WDF is real, a measurement of the intensity at the output plane of a WDF-generating setup yields the desired WDF values up to their sign, which can however be easily deduced from the (eventually) passages through zero of the intensity distribution if the WDF sign is known in one point; this is a trivial task for symmetric field distributions and a not-impossible task in general.
Except the rotationally symmetric beam for which the four-dimensional WDF can be recovered from the twodimensional WDF of the field distribution along a line that passes through the center [51], only samples of the fourdimensional WDF of an arbitrary two-dimensional light source can be produced. Several methods have been proposed to generate these samples, and they differ through the manner of implementing the product between the field distribution and its shifted and complex conjugate replica: a dual-channel display technique [52], two-tilted [52] or properly shifted and rotated transparencies [52,53], or a single transparency illuminated twice [53] can be used to generate desired samples of the WDF of real-valued fields, while a simultaneous optical production of several sections, W n (r n /2, k/2), of the WDF can be implemented using a fiber grating to shift multiple copies of the input object with r n [54].
Apart from the possibility of direct generation, the WDF function can be recovered from measurements of other distribution functions such as the spectrogram or the Radon or Hartley transforms (see the references in [6]). Alternatively, a powerful technique has been devised for the computation of the WDF, and hence for the retrieval of both amplitude and phase of the original field distribution, from measurements of light intensities at different planes or after passing through optical systems that rotate the original WDF in phase space. This method is known as tomography [55], and has been extensively used for the reconstruction of the WDF of one- [56] or two-dimensional optical beams [57], coherent and lowcoherent light sources [58,59], and optical pulses [60,61]. As expected, the two methods of complex field recovery, that is, the direct generation of the WDF and the tomography, provide comparable results [62].
For two-dimensional monochromatic point-symmetric fields, that is, for fields for which ϕ(−r) = ±ϕ(r), the WDF is the inverse Fourier transform of its scaled fractional Fourier transform spectrum, and hence can be determined from intensity measurements of the fractional Fourier transform of the incident field in a suitable number of points [63].
Closely related to the problem of signal retrieval is that of signal synthesis. A signal synthesis algorithm that involves the discrete-time WDF defined in (4) has been developed in [64]; it consists in finding the digital sequence ϕ(nτ) whose discrete-time WDF best approximates (in a least-squares sense) a desired function. A more detailed account on timevarying signal processing using WDF synthesis techniques can be found in [65].

APPLICATIONS TO IMAGE RECOGNITION
In applications involving pattern recognition and image analysis, the moments of the field distribution as well as the invariants that can be constructed with these moments are essential [66]. For example, the first-order moments of the spatial coordinates calculated with the monochromatic field intensity as the weighting function locate the centroid of the intensity distribution, while the second-order moments of the spatial coordinates characterize the orientation and the size of the image. Although the information contained in the original image can be recovered only if the infinite set of all moments is known, in most applications only a finite set of higher-order moments and moment invariants (in particular Zernike moment invariants [66]) are sufficient for image recognition and/or image classification purposes. Teague [67] showed that it is possible to calculate (with limited accuracy) image moments by optical means, that is, by using a Fourier transforming lens to perform the spatial integration. Why do we need then a WDF-based definition for moments? A powerful justification of the usefulness of WDF-based moments has been provided by Freeman and Saleh [68]. They showed that for many image recognition purposes it is better to use a set composed of moment invariants from both spatial and angular (spatial frequency) domains instead of the same number of moment invariants from the spatial or spatial frequency domain alone. This result, welcomed by experimentalists who know that noise hamper the accurate determination of moments to a degree that increases with the moment's order, validates the use of the mixed WDF-based moments in order to achieve the image recognition task with lower-order moments and moment invariants compared to the case when spatial or spatial frequency moments alone are used. Not to mention that the simple propagation laws of WDF-based moments through real and symplectic firstorder systems offer the easiest way to find moment invariants.
A last but not least reason for the employment of moments defined with the WDF as a weighting function is that they can, for arbitrary orders, be determined from the measurement of intensity moments of an appropriate number of fractional power spectra of the incident field distribution, obtained by measuring the intensity at the output plane of generally anamorphic fractional Fourier transform systems [69]. Other set-ups for measuring the ten second-order WDF-based moments of two-dimensional beams, which involve the measurements of second-order spatial moments at the output of astigmatic optical systems, have been put forward by Nemes and Siegman [70], and by Eppich et al. [71], whereas the phase space beam analyzer [72] measures simultaneously the beam radius and far-field divergence of onedimensional field distributions.

APPLICATIONS TO THE CHARACTERIZATION OF SIGNALS
The spatial and angular extent of optical field distributions is determined by its second-order moments. For stationary monochromatic and centered two-dimensional field distributions (for whichx =ȳ = k x = k y = 0) the squared width of the beam along x, y and the associated squared far-field divergence angles are determined, respectively, by x 2 , y 2 and k 2 x , k 2 y , ((x −x) 2 , (y −ȳ) 2 and (k x − k x ) 2 , (k y − k y ) 2 for noncentered light beams) while an overall characterization of the quality of the beam is provided by [73] The second-order moments and the beam quality parameter are tools that allow a quantitative comparison of the spatial and spatial frequency extent of beams that can otherwise differ substantially with respect to their shape or wavefront profile. Moreover, the propagation law (18) of the WDF-based moments through first-order systems is a useful tool for designing optical systems that can optimize the extent or the shape of the field distribution required in some applications. First-order optical systems leave invariant the beam quality parameter Q and modify the WDF-based moments according to (18). The beam quality factor Q takes a minimum value for a Gaussian field distribution. The square root of the ratio between the Q value of an arbitrary beam and that of a Gaussian field is known as the beam-propagation factor and is denoted by M 2 ; M 2 describes the far-field spreading of an arbitrary beam relative to a Gaussian field with the same waist. A high-quality beam, that is, a beam with a small Q value is simultaneously localized in spatial and angular coordinates. The minimum value of Q is 1/4 for one-dimensional monochromatic beams, and 1/2 for two-dimensional field distributions. The occurrence of a minimum value for Q corresponds to the existence of a classical counterpart of the uncertainty relations in quantum mechanics; uncertainty relations can also be defined for nonparaxial wave fields [74]. For multimode beams in an optical resonator, M 2 determines the number of transverse modes in one dimension [75].
An effective radius of curvature matrix R −1 for an arbitrary quasimonochromatic but centered optical beam can be introduced in terms of second-order WDF-based moments as R −1 = (1/k)u T ⊗ v/u T ⊗ u; R −1 has in the spatiotemporal phase space the meaning of a spatial curvature and a temporal chirp [30]. The matrix division in the definition of R −1 should be understood as multiplication of the matrix in the numerator with the inverse of the matrix in the denominator. A similar definition of an effective radius of curvature matrix can be introduced for each frequency component of a polychromatic pulse that propagates through nondispersive and achromatic optical systems [30]. A complex radius of curvature matrix q −1 can also be introduced for arbitrary beams as q −1 = R −1 − iλM 2 /(4πu T ⊗ u), the arbitrary beam satisfying the same propagation rule through first-order optical systems as a Gaussian beam with the same size if the wavelength λ is replaced with λM 2 [76,77,78]. In particular the second-order spatial moment (the beam width) obeys a parabolic law of propagation in free space, the minimum beam width occurring at the waist plane irrespective of the shape and wavefront profile of the beam. This remarkable result, which allows the characterization of arbitrary beams with concepts derived for Gaussian field distributions, is only possible when second-order moments, and in particular WDF-based second-order moments, are employed for beam characterization.
The twist parameter of partially coherent Gaussian field distributions [79], as well as the angular rotation of arbitrary fields under free space propagation can also be expressed in terms of second-order moments [80]. More precisely, for a two-dimensional monochromatic light beam, the angle between the principal axes, defined as those for which xy = 0, and the absolute axes, that is, those for which k x k y = 0, is given by cot 2θ = (x 2 − y 2 )/(2xy), the angular rotation offering a classification criterion of light beams. It is also worth mentioning that the Stokes matrices associated with partially polarized fields can be also expressed in terms of the secondorder moments of the WDF matrix of the vectorial electric field [14].
The symmetry of one-dimensional fields is described by the third-order moment x 3 , the moments x 2 k x and x 3 k x measure, respectively, the spatial range of the beam's symmetry and sharpness, x 2 k 2 x /(x 2 k 2 x ) expresses the degree of similarity of an arbitrary beam with a quasihomogeneous field distribution [81], while the kurtosis parameter K = x 4 /(x 2 ) 2 provides a quantitative measure for the classification of beams with respect to their sharpness [82,83]. K is useful to estimate the capabilities of laser beams in material processing applications.
Special relations between moments of the WDF occur for light beams with certain symmetries. For example, for onedimensional self-Fourier field distributions ϕ(x), which are characterized byφ(k x ) = a 1 ϕ(a 2 k x ) with constant a 1 , a 2 parameters such that |a 2 1 /a 2 | = 1, the WDF satisfies the symmetry condition W ϕ (x, k x ) = W ϕ (a 2 k x , −x/a 2 ), while its moments are related through x m k n x = (−1) m a m−n 2 x n k m x [84]. The symmetry of moments is thus an indication of the symmetry of the field. Field distributions with such a special symmetry can preserve it only at propagation through certain optical systems, characterized by a symplectic matrix that can be determined from (15). Bright solitons are examples of self-Fourier functions.
Not only the moments defined by (17), which are called global moments, have physical significance, but also local moments can be associated with parameters of the signal. For example, K(r) = kW(r, k)dk/ W(r, k)dk can be interpreted as the scaled (with k) average propagation angle at position r, R(k) = rW(r, k)dr/ W(r, k)dr can be identified with the average beam position at the transverse wavevector k, T(ω) = tW(t, ω)dt/ W(t, ω)dt can be seen as the group delay of the optical pulse, and Ω(t) = ωW(t, ω)dω/ W(t, ω)dω can be defined as the instantaneous frequency of a pulse [1].
Although the WDF is mainly used for characterizing optical beams and pulses that propagate through linear media, some studies have also focused on nonlinear propagation. A thorough analysis of the propagation of the WDF and of its moments through active media, for example, can be performed under the thin sheet condition, which restricts the diffraction along the amplifier and limits the magnitude and variation of the refractive index and gain coefficient of the active medium. Under these assumptions the field amplitude in the paraxial approximation can be expressed as a sum of two terms, one of them being approximately proportional to the z coordinate along the propagation direction. The WDF of this sum contains both auto-terms and a cross-term, each of them satisfying the free space propagation law in (15) with (in one dimension) A = D = 1, C = 0, and B = z [85]. As a result, the WDF does not propagate along straight lines in active media and the radiant intensity, that is, the average of the WDF along the spatial coordinate, as well as the beam quality parameter do not remain invariant at propagation. Moreover, an incident aligned field distribution, withx = k x = 0 at the input plane does not preserve this property at the output plane, in contrast to the behavior in homogeneous and passive media. Numerical simulations have shown that at propagation of Gaussian beams through active media under the thin sheet approximation, in both homogeneously and inhomogeneously broadened cases, the beam size increases compared to the propagation through passive media and the kurtosis parameter lowers (the beam becomes flatter) [86]. The form of the WDF also changes at propagation through nonlinear media. For example, at self-phase modulation the WDF of an incident Gaussian pulse develops an asymmetric pronounced dip and can even take negative values [87].
Studies performed using the WDF as a mathematical tool have shown that the degree of coherence of the beam has a determinant role in propagation. For example, partial incoherence weakens the nonlinearity, a higher intensity being required to achieve the same effect as that obtained with a coherent beam. In particular, partial incoherence tends to suppress the coherent instabilities, acting as a Landau-like damping effect on modulational instabilities, or tends to prevent the self-focusing collapse instability of bidimensional field distributions [88]. These effects have been demonstrated employing the WDF of the envelope function of partially coherent fields together with the Klimontovich statistical average. The evolution of the beam quality factor in a nonlinear Kerr medium with a quadratic refractive index profile has been derived in [89]; an invariant beam quality parameter can be introduced for nonlinear Kerr-like media for a properly defined effective radius of curvature. Moreover, Nasalski [90] has demonstrated that the first-order formalism can be preserved at propagation of a light pulse or beam through a nonlinear Kerr-like medium if the phase space coordinates are properly scaled: the low-power (linear) coordinates u, v should be replaced with U = u/w(ζ), V = vw(ζ), where w is the half-width pulse waist and ζ is the (temporal and spatial) chirp parameter, which vanishes at the actual waist plane. ζ replaces the z coordinate as propagation parameter and differs from it because of the self-shortening effect and because of the waist shift with respect to the linear propagation situation. Apart from these scaling operations of the phase space coordinates, an on-axis phase correction must be introduced in the expressions of both the field amplitude and its Fourier transform, and an additional scaling with w n/2 (ζ) and w −n/2 (ζ) of the field amplitude and its Fourier transform, respectively, must be taken into account, where n is the dimension of the nonlinear Schrödinger equation satisfied by the quasimonochromatic pulse.
The analytic expression of the WDF for bright solitons and for the two-soliton solution of the nonlinear Schrödinger equation, which is identical to the equation satisfied by optical solitons propagating in a dispersive and nonlinear medium, can be found in [91]. Bright solitons are characterized by a diamond-shaped WDF, while the WDF of the twosoliton solution evolves in time oscillating between profiles similar to the dark soliton and profiles with a large momentum distribution, which correspond to the most squeezed state in spatial coordinates. Most interesting, however, is the appearance of a distinct interference cross-term in the WDF of the two-soliton solution when the amplitudes of the two solitons are comparable. The WDF of black and gray solitons in the temporal phase space as well as the WDF moments of all envelope solitons (bright, black, gray) also have analytical expressions [92]. These expressions have provided test methods to decide if an optical pulse is or not a soliton; such phase space tests were applied to a bright soliton pulse source used for 270 terabit km/s transmission system [92] and for bright and dark magnetostatic solitons [93].
The transport equation of the WDF for field distributions that satisfy the nonlinear Schrödinger equation, the Korteweg-de Vries equation, as well as the Burgers equation was derived by Kamp [94], who separated the operator G(r, ∂/∂r, ϕ(r, t)) acting upon the field distribution into a linear and a nonlinear part: G(r, ∂/∂r, ϕ(r, t)) = F(r, ∂/∂r) + L[ϕ(r, t)]. Note that the WDF evolution law under a linear operation that acts on the field distribution according to is given by (See also [34,35,36,37].) The equation of motion of the WDF for an underlying governing equation that is a linear ordinary or partial differential equation has been derived in [95]. Starting from (22) the transformation laws for the WDF, the radiant intensity, the radiant emittance, the first-and second-order moments of the WDF, and the kurtosis parameter of a field that propagates through an inhomogeneous Kerr-type medium have been derived [96]. Quite surprisingly, it was found that the WDF transforms according to a law similar to (15), where the elements of the real symplectic matrix depend on the field distribution. This fact leads, however, to a variation of the beam quality factor with the propagation distance, in contrast to the case when the elements of the real symplectic matrix do not depend on the field distribution.
The soliton solutions of the nonlinear Schrödinger equation refer to envelope solitary waves in χ (3) materials, but stationary wave solutions of soliton type exist also in χ (2) materials. These are related to a variety of phenomena that include degenerate and nondegenerate parametric amplification and the coupling of intense electromagnetic fields to polaritons. The solitons in χ (2) media are solutions of coupled field equations in which the nonlinearity is second order in the field, and exist in both anomalously and normally dispersive regions. In absorptive media, the widths of χ (2) solitons are maintained over long distances while their amplitudes decay exponentially. The WDF of the bright and dark χ (2) solitons and their second-order moments are given in [97]; the obtained expressions are different from the corresponding ones in χ (3) solitons, providing a signature of the soliton type, and can be used to test different pulses.

APPLICATIONS TO THE CHARACTERIZATION OF OPTICAL SYSTEMS AND DEVICES
A measurement of the WDF after it passes through an optical system is a valuable means to characterize the system if the incident field distribution is known. Due to the ABCD transformation law and due to the fact that the WDF can be defined for either coherent or partially coherent field distributions, the WDF incorporates information about the coherence of illumination and the optical system through which it passes. Therefore it can be used to evaluate the performance and tolerances of optical processors [98] and, in particular, to estimate the parameters of various optical systems. For example, according to (15), the WDF of a collimated, monochromatic, and one-dimensional field, W(x, k x ) = δ(k x ), becomes W(x, k x ) = δ(k x + k x x/ f ) after passing through a thin lens characterized by a transmission function exp(ikx 2 /2 f ), so that the tilt angle of the WDF after the lens is a measure of its focal length. (The free space produces a shear transformation of the WDF along x.) Analogously, the passage of the same field distribution through a pure phase object with a transmission function exp[iφ(x)] that cannot be in general characterized through a real symplectic matrix yields an output WDF of the form W(x, k x ) = δ(k x + dφ/dx), from which the profile of the pure phase object can be recovered. An example of the phase profile recovery from the WDF measurement can be found in [99]. One example of special interest is the case of spherically aberrated lenses. The type of aberration, as well as the value of the aberration coefficient can be determined from the form of the WDF of the output field distribution [100]. The presence of aberrations manifests itself in the change of shape of the output WDF, in addition to the coordinate transformation predicted by the ABCD law in (15). As demonstrated in [101], the output WDF of an aberrated optical system is obtained by applying an exponential differential operator to the WDF of the same optical system in the absence of aberrations, the first term of the operator predicting the coordinate transformation of the unaberrated WDF while the second term predicts the shape distortion.
In many applications the overall effect of aberrations can be quantified by the changes in Q. For a Gaussian incident field, for example, the effect of spherical aberrations can be calculated analytically [102]; in the two-dimensional case, spherical aberrations always deteriorate the beam quality if the incident field distribution is real. On the contrary, a Gaussian aperture always improves the beam quality parameter for one-dimensional partially coherent Gaussian-Schell beams since it reduces the beamwidth (the far-field divergence can be enlarged or decreased depending on the position of the aperture with respect to the waist plane) [103].
The WDF can also be used to express the intensity distribution along different paths in the image space of an optical system [104], and can hence be used for an efficient analysis of the performance of the imaging system in the presence of spherical aberrations. Moreover, the WDF is also useful for the design of pupil filters that would generate specific axial responses in the presence of spherical aberrations [105]. The WDF of the desired pupil filters is obtained through tomography from the desired axial irradiance distributions, and the pupil filter is subsequently determined from its WDF.
Due to the simplicity of the ABCD transformation law of the WDF through first-order optical system, it is easier to study optical systems using the WDF approach than to compute the field evolution through such systems. And this approach can also be used for polychromatic light beams, as long as each frequency component is treated separately. This is the reason why it is desirable to employ the WDF approach to follow the propagation of polychromatic coherent fields through a system designed to achromatize Fresnel diffraction patterns, for example, [106]. The residual chromatic aberrations can be obtained as well using this approach to achromatic white-light self-imaging.
The propagation of partially coherent fields through apertures is another subject that can be linked to the WDF. More precisely, the angular spectrum of the cross-spectral density of the partially coherent field at an observation plane is a WDF that depends on the path vectors between the aperture and the observation plane [107]. The complex degree of spatial coherence acts as a weighting factor for the contribution of aperture radiator pairs to the WDF at each position in the observation plane.
When pulses are launched into single-mode fibers,t and ω represent the mean arrival time and the mean frequency of the pulse, respectively, the squared temporal and spectral pulse widths being determined by (t −t) 2 and (ω −ω) 2 , respectively, [108]. Moreover, the transmission capacity of a fiber at a z = const. plane can be defined as the inverse of the pulse width at that plane, that is, as At propagation through dispersive fibers with a secondorder dispersion coefficient, the beam quality parameter Q = (t −t) 2 (ω −ω) 2 − ((t −t)(ω −ω)) 2 remains constant and the moments of the time and frequency coordinates transform according to (18), whereas for fibers with third-order dispersion terms, Q is a second-order polynomial in the fiber's length z; the evolution of the transmission capacities of the fibers in the two cases can be calculated analytically [108].
The study of active devices, in particular lasers, can also benefit from the computation of moments of the field distribution and of the beam quality parameter. For example, the beam quality parameter and the second-order moments for the fields emitted by heterojunction monomode diode lasers with very thin active layers can be calculated analytically and can be related to physical constants that characterize the laser device [109]. These parameters of the emitted field can then be used to design first-order optical systems that can transform the field emitted by monomode diode lasers into a circular field distribution at a desired plane. (This transformation can only be achieved at some planes since the emitted field has different degrees of global coherence on each transverse axis.) An efficient coupling of such laser beams into circularly symmetric optical fibers, for example, require thus the introduction of precisely designed first-order optical systems between the laser and the fiber's end.
The WDF was also applied by Gase [110] to the design of multimode laser resonators consisting of two mirrors and a laser rod between them, the laser rod acting as a convergence lens due to the temperature gradient. The ABCD transformation law of the WDF and of its moments was used to match the radius of the laser beam, considered as an isotropic Gaussian-Schell beam, to the free-rod diameter, this condition describing the stability of the resonator; note that the radius of the laser beam depends on the beam quality parameter. The focusing performances of multimode laser beams as well as the disturbance effect of an aperture diaphragm were also treated with the WDF formalism.
Holograms can also be related to the WDF. A hologram is an interference pattern recorded by the superposition of the signal ϕ s (r) and the reference beam ϕ r (r), the signal wavefront being recovered by reading the hologram, that is, by illuminating it with a replica of the reference signal. It can be demonstrated that the interference pattern of a hologram is identical to the partial integral over k of the interference term (or cross-term) W ϕsϕr (r, k) of the field distribution ϕ s (r) + ϕ r (r) obtained by superposing the signal and reference beams [111]. This cross-term in the WDF has a strongly oscillating behavior and therefore has been called "the smile function" of the superposition state. Note that once a hologram is recorded it is possible to read the information carried by it by measuring the WDF of the hologram without the need to employ a replica of the reference beam in order to recover the signal wavefront. This fact has been used successfully to extract the location of three-dimensional objects from inline holograms without reconstructing the three-dimensional field [112]. The procedure works if the concentration of particles in inline holograms is not too high, since otherwise the cross-terms in the WDF may interfere with the auto-terms and the interpretation of the WDF becomes extremely difficult.
Another effect that can be analyzed with the WDF is the moiré pattern, which is an undesired feature of pixelated displays that originates from undersampling and which is encountered in image processing, topography, Talbot interferometry, and submicrometer alignment methods. The multiplicative moiré patterns with well-defined local frequencies, resulting from the superposition of nonperiodic masks, for example, can be analyzed disregarding the cross-terms in the WDF of the superposed masks [113], but the cross-terms must be accounted for in interferometric applications. The WDF interpretation of moiré patterns provides an intuitive insight into the information stored in the phase modulation of nonperiodic moiré fringes.
Optical set-ups that shape ultrashort pulses, for example, can only be characterized by a WDF defined on a spatiotemporal phase space. A typical set-up of this kind consists of a grating that spatially separates the different frequency components of an ultrashort pulse, followed by a collimating lens, a mask that modifies in a desired manner each spectral component, and another lens that focuses the collimated beam onto a subsequent grating, which superimposes spatially the diverse spectral components. Each of these devices, with the exception of the mask, can be described by a real symplectic S matrix such that the evolution of the WDF through the set-up can be described in a much simpler way than the evolution of the spatially-and frequency-dependent field distribution. As mentioned before, a spatio-temporal phase space is needed since the grating couples the spatial and temporal features of the pulse. The effect of the mask with a transmission function t(u) on the incident field distribution ϕ in (u), described as ϕ out (u) = t(u)ϕ in (u), corresponds to a change of the WDF given by (see also [12]) The convolution in the v coordinate describes the action of any filtering device in u, a symmetrical effect, that is, a convolution of the WDF in u, occurring for any filtering action in the v domain. The WDF treatment of the pulse shaper is detailed in [114].
Another optical system that has been studied with the spatio-temporal WDF is a Lukosz multiplexer with moving gratings, which sends a one-dimensional signal through a finite-width aperture by segmenting it into parts separated by small differences in wavelengths, and then reconstructs the signal [115]. More precisely, the multiplexer consists of a moving grating, a Fourier transforming system, the finitewidth aperture, another Fourier transforming system, and a counterpropagating grating. Since a Fourier-transforming system formed from a lens and free spaces is described by an ABCD matrix and the action of the finite-width aperture can be seen as that of a filtering device, a WDF approach to the analysis of the multiplexer would be justified if the moving grating could be described easily in phase space. (Note that, unlike in the beam shaper the grating is illuminated with an extended object.) Fortunately, such a simple description exists: the WDF W in (x, k x ; λ) of a one-dimensional input signal of wavelength λ, ϕ(x, λ) is multiplied by the grating and transformed into where Λ is the grating period and W n (x) is the intensity of each copy, which is directly related to the shape of the grating. A moving grating viewed with optical sensors that integrate over time should be represented in phase space by the same expression as in (24), but the coefficients W n (x) become now time dependent and hence should be averaged over time.
Other optical set-ups analyzed with spatio-temporal WDFs include a white-light spatial Fourier transformer, a rotated zone plate [116], and a grating compressor for optical pulses consisting of two parallel diffraction gratings separated by a free space region [117]. One of the most recently emerged applications of temporal phase space distribution functions is the synthesis [118] or reconstruction [119] of periodic grating profiles, in particular of fiber Bragg gratings. Such devices can function, for example, as wavelength-selective components, dispersion compensators, or mode converters. The problem is to reconstruct or synthesize the refractive index profile from the measured or desired values of the complex field reflection coefficient r(ω) or the reflection impulse response r(t); these two functions are related by a Fourier transform. Once one of these is known, it is easy to compute the WDF but the next step, that is, the finding of the refractive index profile, is generally difficult due to the presence of cross-terms in the WDF of a multicomponent signal. These cross-terms, which reflect the correlations between pairs of auto-terms, have larger amplitudes than the auto-terms and are situated in the middle (in both temporal and frequency coordinates) between the corresponding auto-terms. It is even possible that cross-terms that connect pairs of auto-terms overlap with other autoterms, the resulting WDF becoming indeed very difficult to interpret. So, for fiber grating reconstruction or synthesis, the WDF is not the best choice: the spectrogram, for example, or the Choi-Williams distribution [4] are better suited since they are not plagued by the presence of cross-terms. On the other hand, the resolution limitation of the spectrogram technique imposed by the choice of an appropriate window function is always worse than that achieved with the WDF approach; a suitable compromise is then to use the crossterm-free spectrogram as a mask over the WDF and so to reduce the unwanted cross-terms of the latter. Our point here is that the WDF is not always the best suited phase space distribution, although it satisfies a number of attractive properties that are not shared by other distribution functions (see, e.g., [1]).
Not only the refractive index profile of fiber gratings can be recovered using the WDF, but also this approach works well for the recovery of spatial nonperiodic profiles of the refractive index of single-or multimoded guided structures. The refractive index profile reconstruction can be achieved by solving a linear system of equations that involve moments of the gradients along r and k of the measured or desired WDF. The procedure works for either longitudinally constant [120] or longitudinally variant but continuous refractive index profiles [121], and for different propagating modes in the optical structure. If the WDF gradients ∇ r W and ∇ k W are approximated with [W(r+dr, k)−W(r, k)]/dr and [W(r, k + dk) − W(r, k)]/dk, respectively, the recovery precision depends on the steps in r and k of the measured or simulated WDF.
It is important to note that the usefulness of the WDF for the characterization of optical systems is not limited to the case when these are of the first-order type, with an eventual inclusion of aberration effects. The effect of complex integrated structures on light propagation is usually estimated with the beam-propagation method, which implies the division of the guiding structure into a number of steps such that in each step the refractive index profile can be decomposed in longitudinally-and transverse-varying parts. The computation of the output field is carried out by applying in each step different operators to the field distribution and to its Fourier transform, which implies that at each step a computation of the direct and the inverse Fourier transform of the field is required [122]. The succession of these operations can be replaced with a much simpler WDF formulation of the beam-propagation method, which has the advantage that it depends on both spatial and spatial frequency coordinates and thus the operators at each step can be applied directly on the WDF, the successive time-consuming direct and indirect Fourier transforms being no longer needed [123]. The operators that must be applied on the WDF are simple differential operators; it is even possible to obtain in the first-order approximation an analytic formula for the WDF transformation through some integrated structures, such as a linear taper with a transverse parabolic refractive index profile.

APPLICATIONS TO THE CALCULATION OF THE COUPLING COEFFICIENT
The overlap or Moyal relation of the WDF (property (vii) of the WDF) can be used to express the coupling efficiency of light distributions to a region of space with different refractive or dispersive properties. The optimization of the coupling efficiency between different optical systems is an issue of enormous practical importance in any signal processing and signal transmission application.
A definition of the coupling coefficient between a light source and a waveguide, for example, in terms of spatial phase space distribution functions is more appropriate for the analysis of the coupling problem than a definition in terms of spatial or spatial frequency distributions since the numerical aperture (an angular parameter) is as important as the dimensions of the waveguide core (a spatial parameter) in the characterization of waveguide field propagation. The expression of the coupling efficiency between a coherent source, with field distribution ϕ s (r) and WDF W s (r, k), and a waveguide mode ϕ w (r) with a WDF W w (r, k) is given by [124] η = (2π) n W s (r, k)W w (r, k)dr dk W s (r, k)dr dk W w (r, k)dr dk , where the integrals extend over the illumination plane. The behavior of the coupling coefficient in (25), which is a real parameter, with spatial and angular misalignments between the source and the waveguide mode can be intuitively grasped from the graphical representations of W s (r, k) and W w (r, k) (see [124,125] for more details). Note that in phase space, the spatial and angular misalignments correspond to displacements of the respective WDF along the r and k axes, respectively. The connection between the value of the coupling coefficient and the graphic representations of the WDF of a quantum-well modulator and the first four modes of a planar optical waveguide, for example, has been discussed in [126]. According to (25), coupling is achieved if the WDFs (or the fields) of the source and waveguide overlap. For partially coherent light distributions, the expression of the coupling coefficient in (25) becomes [127] η = W sw (r, k)dr dk 2 W s (r, k)dr dk W w (r, k)dr dk , where W sw (r, k) = (2π) −n/2 ϕ s r + r 2 ϕ * w r − r 2 exp ir k T dr = (2π) −n/2 Γ sw r + r 2 , r − r 2 exp ir k T dr (27) while the coupling coefficient between a coherent multimode light source described by the field distribution ϕ s (r) = n a n ϕ n s (r) and a superposition ϕ w (r) = m b m ϕ m w (r) of waveguides modes ϕ m w (r) can be expressed as [127] η = a l a * k b n b * m η lknm .
Here η lknm = (2π) n W lk s (r, k)W nm w (r, k)dr dk W s (r, k)dr dk W w (r, k)dr dk (29) with W lk i (r, k) = (2π) −n/2 ϕ l i r + i = s, w, W s (r, k) and W w (r, k) having the same significance as in (25).

CONCLUSIONS
The paper introduces the WDF and its most important properties as a mathematical tool in several areas of signal processing that include signal retrieval, signal recognition, characterization of signals and optical systems, and coupling coefficient estimation in phase space. The mathematical formalism can be applied to either spatial, temporal, or spatiotemporal phase spaces, to coherent, partially coherent or digital signals, offering a unified view for the analysis of field propagation through various optical systems. In many cases the WDF approach to the study of field propagation and coupling is very much simplified due to its desirable properties and due to the extremely simple propagation law (see (15)) through first-order optical systems. This propagation law can be employed to characterize optical systems described by symplectic matrices, active devices such as lasers, or to estimate the value of aberration effects. The moments of the WDF, subject also to a very simple evolution law expressed by (18), are related to spatial, temporal, angular, and/or frequency characteristics of the field distribution and can be either directly measurable or can be calculated from the experimentally obtained WDF. These moments and/or their invariant combinations (such as the beam quality factor) are quantitative parameters that allow the comparison of beams with totally different shapes or wavefronts. The WDF is thus a universal tool, in the sense that it can be applied to (almost) any field distribution and an intuitive and simple representation of field evolution (due to its ABCD law of transformation) and coupling (because it characterizes simultaneously the fields' overlap in both spatial and spatial frequency domains). The continuous growth of the number of applications of the WDF is hence predictable, desirable, and enjoyable.