## Abstract

Phase space tomography estimates correlation functions entirely from snapshots in the evolution of the wave function along a time or space variable. In contrast, traditional interferometric methods require measurement of multiple two–point correlations. However, as in every tomographic formulation, undersampling poses a severe limitation. Here we present the first, to our knowledge, experimental demonstration of compressive reconstruction of the classical optical correlation function, *i.e.* the mutual intensity function. Our compressive algorithm makes explicit use of the physically justifiable assumption of a low–entropy source (or state.) Since the source was directly accessible in our classical experiment, we were able to compare the compressive estimate of the mutual intensity to an independent ground–truth estimate from the van Cittert–Zernike theorem and verify substantial quantitative improvements in the reconstruction.

© 2012 Optical Society of America

## 1. Introduction

Correlation functions provide complete characterization of wave fields in several branches of physics, e.g. the mutual intensity of stationary quasi–monochromatic partially coherent light [1], and the density matrix of conservative quantum systems (*i.e.*, those with a time–independent Hamiltonian) [2]. Classical mutual intensity expresses the joint statistics between two points on a wavefront, and it is traditionally measured using interferometry: two sheared versions of a field are overlapped in a Young, Mach–Zehnder, or rotational shear [3,4] arrangement, and two–point ensemble statistics are estimated as time averages by a slow detector under the assumption of ergodicity [1, 5].

As an alternative to interferometry, phase space tomography (PST) is an elegant method to measure correlation functions. In classical optics, PST involves measuring the intensity under spatial propagation [6–8] or time evolution [9]. In quantum mechanics, analogous techniques apply [10–13]. However, the large dimensionality of the unknown state makes tomography difficult. In order to recover the correlation matrix corresponding to just *n* points in space, a standard implementation would require at least *n*^{2} data points.

Compressed sensing [14–16] exploits sparsity priors to recover missing data with high confidence from a few measurements derived from a linear operator. Here, sparsity means that the unknown vector contains only a small number of nonzero entries in some specified basis. Low–rank matrix recovery (LRMR) [17, 18] is a generalization of compressed sensing from vectors to matrices: one attempts to reconstruct a high–fidelity and low–rank description of the unknown matrix from very few linear measurements.

In this paper, we present the first, to our knowledge, experimental measurement and verification of the correlation function of a classical partially coherent field using LRMR. It is worth noting that LRMR came about in the context of compressive quantum state tomography (QST) [19], which utilizes different physics to attain the same end goal of reconstructing the quantum state. In PST, one performs tomographic projection measurements, rotating the Wigner space between successive projections by evolving the wave function [6, 7]. This is directly analogous to the classical optical experiment we are presenting here, where we perform intensity measurements (*i.e.*, tomographic projections in Wigner space) and utilize propagation along the optical axis to rotate the Wigner space between projections. The difference lies in the fact that in QST the state is recovered via successive applications of the Pauli dimensionality–reducing operator, and there is no need to evolve the state. Nevertheless, both approaches lead to the same Hermitian LRMR problem, as long as the assumption of a quasi–pure unknown state is satisfied. In [20], it was shown that estimation of a low–rank matrix of dimension *n* and rank *r* requires only *O*(*rn*ln*n*) to *O*(*rn*ln^{2} *n*) data points. A similar LRMR method was also used to recover the complex amplitude of an unknown object under known illumination [21–23]. Since the complex amplitude of the object is time–invariant, a rank–one solution was assumed in these works.

The low–rank assumption for classical partially coherent light anticipates a source composed of a small number of mutually incoherent effective sources, *i.e.* coherent modes [24], to describe measurements. This is essentially equivalent to the low entropy assumption [19], *i.e.* a nearly pure quantum state in the quantum analogue. This assumption is valid for lasers, synchrotron and table–top X–ray sources [25], and Köhler illumination in optical microscopes [1]. An additional requirement for LRMR to succeed is that measurements are “incoherent” with respect to the eigenvectors of the matrix, *i.e.* the measured energy is approximately evenly spread between modes [20, 21]. Diffraction certainly mixes the coherent modes of the source rapidly, so we expect LRMR to perform well for classical PST. The same expectation for QST has already been established [19].

## 2. Theory and method

The two–point correlation function of a stationary quasi–monochromatic partially spatially coherent field is the mutual intensity function [1]

where 〈·〉 denotes the expectation value over a statistical ensemble of realizations of the field*g*(

*x*).

The measurable quantity of the classical field, *i.e.* the intensity, after propagation by distance *z*, is [1]

*P*denotes the free–space propagation operator that combines both the quadratic phase and Fourier transform operations in Eq. (2), tr(·) computes the trace, and

*x*

_{o}denotes the lateral coordinate at the observation plane. By changing variables

*x*= (

*x*

_{1}+

*x*

_{2})/2,

*x*′ =

*x*

_{1}–

*x*

_{2}and Fourier transforming the mutual intensity with respect to

*x*we obtain the Ambiguity Function (AF) [26–28]

*Ĩ*is the Fourier transform of the vector of measured intensities with respect to

*x*

_{o}. Thus, radial slices of the AF may be obtained from Fourier transforming the vectors of intensities measured at corresponding propagation distances, and from the AF the mutual intensity can be recovered by an additional inverse Fourier transform, subject to sufficient sampling.

To formulate a linear model for compressive PST, the measured intensity data is first arranged in Ambiguity space. The mutual intensity is defined as the “sparse” unknown to solve for. To relate the unknowns (mutual intensity) to the measurements (AF), the center–difference coordinate–transform is first applied, expressed as a linear transformation *T* upon the mutual intensity *J*, followed by Fourier transform , and adding measurement noise *e* as

The mutual intensity propagation operator is unitary and Hermitian, since it preserves energy. We use eigenvalue decomposition to determine the basis where the measurement is sparse. The resulting basis, *i.e.* the set of eigenvectors, is also known as coherent modes in optical coherence theory, whereas the whole process is known as coherent mode decomposition [24]. The goal of the LRMR method is to minimize the number of coherent modes to describe measurements. By doing LRMR, we impose two *physically* meaningful priors: *(1)* existence of the coherent modes [24], and *(2)* sparse representation of the partially coherent field in terms of coherent modes.

Mathematically, if we define all the eigenvalues *λ _{i}* and the estimated mutual intensity as

*J*̂, the method can be written as

*ℓ*

_{1}norm) of the matrix

*J*[17, 29]. The corresponding problem is stated as

*σ*= |

_{i}*λ*|, ||

_{i}*J*̂||

_{*}= ∑

_{i}*σ*. This problem is convex and a number of numerical solvers can be applied to solve it. In our implementation, we used the singular value thresholding (SVT) method [30]. The output estimate after each iteration of SVT typically has a sub-normalized total energy,

_{i}*i.e.*∑

_{i}*λ*< 1; we compensated for this by renormalizing at the end of each iteration [19].

_{i}## 3. Numerical simulations

First, we demonstrate the LRMR method with a numerical example using a 1D Gaussian-Schell model source (GSMS). Both the intensity distribution and the degree of coherence of GSMS follow a Gaussian distribution [31]

*σ*determines the spatial extent of the source, and

_{I}*σ*is proportional to the coherence length and determines the number of coherent modes in the input source. The eigenvalues of GSMS are never zero (analytical solution given in [31]). We defined the number of modes (rank of the source)

_{c}*r*as the first

*r*modes containing the 99% of the total energy.

One example is shown in Fig. 1(a). The parameters in this example are *σ _{I}* = 17 and

*σ*= 13 (rank

_{c}*r*= 6). Intensities are calculated at 40 different axial distances and the coverage in Ambiguity space is shown in Fig. 1(b). We simulate the case where data from both the near field and the far field are missing due to the finite range of camera scanning motion allowed in the actual experiment. The missing cone around the

*u*′-axis is due to missing data from near field, while the data missing from far field results in the missing cone around the

*x*′-axis. Both cones have an apex angle of 20 degrees.

For comparison, the data is first processed using the traditional filtered-backprojection (FBP) method [32]. Applying the Fourier-slice theorem to Eq. (5) implies that the 1D Fourier transform of a radial slice in the Ambiguity space (an intensity measurement) is related to a projection in the AF’s 2D Fourier space (the Wigner space [33,34]). The Wigner distribution function (WDF) is related to the mutual intensity by

*i.e.*, the intensity) when we use FBP for the reconstruction can be best explained with the help of Figure 1(b). Going from the Ambiguity space to the mutual intensity space involves Fourier transforming along horizontal lines, parallel to the

*u*′ axis. The diagonal in particular corresponds to the line

*x*′ = 0. It can be easily seen that, due to the missing cone, pretty much all the data are missing from that line, except near the origin; thus resulting in a low–pass filtering effect. The fact that the compressive reconstruction method manages to restore the physically correct values of the correlation along the diagonal corroborates that the missing cone is successfully retrieved in our LRMR reconstruction. The FBP reconstruction may also be compared quantitatively to the compressive reconstruction in terms of the global degree of coherence parameter $\overline{\mu}=\frac{\sqrt{{\sum}_{i}{\lambda}_{i}^{2}}}{{\sum}_{i}\left|{\lambda}_{i}\right|}$ [35, 36], which was found as 0.150 and 0.617, respectively; the true state has

*μ*̄ = 0.618.

The coherent modes for a GSMS are Hermite–Gaussian sources [31]. The theoretical and LRMR estimated first nine coherent modes in this example are shown in Fig. 2(a) and 2(b), respectively. The theoretical eigenvalues are shown in Fig. 3(a). The FBP and LRMR estimated eigenvalues are compared in Fig. 3(b) and 3(c), respectively. The FBP estimates have several negative values, which does not satisfy the positive energy constraint. The absolute errors in LRMR estimates are plotted in Fig. 3(d).

Next, we study the noise performance of the LRMR method with a numerical example. In this simulation, the dimension of the input GSMS is 256×256 with parameters *σ _{I}* = 36 and

*σ*= 18 (rank

_{c}*r*= 9). We generate noisy data with different signal-to-noise ratio (SNR) from both an additive random Gaussian noise model and a Poisson noise model. However, we emphasize that the reconstruction algorithm does not make use of the noise statistics. For each SNR level, we repeat the simulation 100 times with different random noise terms, and then record the average relative mean-square-error (MSE) from the LRMR reconstruction. The ratio between the number of samples taken from the intensity measurements and the rank

*r*of the input mutual intensity matrix determines the oversampling rate [22]. This rate is plotted versus relative MSE for different SNR cases in Fig. 4. For good performance, the required oversampling rate is at least 5–6 (the theoretical oversampling rate is on the order of ln(256) = 5.5 according to [20]). Furthermore, the LRMR method is robust to noise in the sense that the reconstruction degrades gracefully as the SNR decreases.

## 4. Experimental result

The experimental arrangement is illustrated in Fig. 5. The illumination is generated by an LED with 620nm central wavelength and 20nm bandwidth. To generate partially coherent illumination, a single slit of width 355.6*μ*m (0.014″) is placed immediately after the LED and one focal length (75 mm) to the left of a cylindrical lens. One focal length to the right of the lens, we place the second single slit of width 457.2*μ*m (0.018″), which is used as a one–dimensional (1D) object.

The goal is to retrieve the mutual intensity immediately to the right of the object from a sequence of intensity measurements at varying *z*–distances downstream from the object, as described in the theory. We measured the intensities at 20 *z*–distances, ranging from 18.2mm to 467.2mm, to the right of the object. The data are given in Fig. 6. Each 1D intensity measurement consists of 512 samples, captured by a CMOS sensor with 12*μ*m pixel size. The dimension of the unknown mutual intensity matrix to be recovered is 512 × 512. Since only intensities at positive *z*, *i.e.* downstream from the object, are accessible, we can only fill up the top right and bottom left quadrants of Ambiguity space. The other two quadrants are filled symmetrically, *i.e.* assuming that if the field propagating to the right of the object were phase conjugated with respect to the axial variable *z*, it would yield the correct field to the left of the object, *i.e.* negative *z* [8, 13]. Under this assumption, a total of 40 radial slices are sampled in Ambiguity space, as shown in Fig. 7. The apex angle of the missing cone around the *u*′–axis is approximately 17.4 degrees, and the one around the *x*′–axis is approximately 28.6 degrees. The number of measurements is only 7.8% of the total number of entries in the unknown mutual intensity matrix.

The reconstructions from the FBP and LRMR methods are compared in Fig. 8(a) and 8(b), respectively. The FBP reconstruction suffers from the same artifacts detailed in the numerical simulations section. All these artifacts are greatly suppressed or completely removed in the LRMR reconstruction. In the real part of the reconstruction, the width of the square at the center is approximately 456*μ*m (38 pixels), which agrees with the actual width of the slit. The imaginary part is orders of magnitude smaller than the real part.

FBP estimated eigenvalues contain several negative values, and are shown in Fig. 9(a). This does not satisfy the positive energy constraint. LRMR estimated eigenvalues are compared in Fig. 9(b), and all eigenvalues are positive.

We further validated our compressive estimates by measuring the field intensity immediately to the right of the illumination slit [Fig. 10(a)]. Assuming that the illumination is spatially incoherent (a good assumption in the LED case), the mutual intensity of the field immediately to the left of the object is the Fourier transform of the measured intensity, according to the van Cittert–Zernike theorem [1, 5]. This calculated mutual intensity, based on the measurement of Fig. 10(a) and screened by the object slit, is shown in Fig. 10(b). The eigenvalues computed by coherent mode decomposition are shown in Fig. 10(c) and are in good agreement with the LRMR estimates, as compared in Fig. 10(d). It is seen that 99% of the energy is contained in the first 13 modes, which confirms our low–rank assumption. The FBP reconstruction may also be compared quantitatively to the compressive reconstruction in terms of the global degree of coherence parameters, which were experimentally found as 0.12 and 0.46, respectively; whereas the estimate yielded by the van Cittert–Zernike theorem is 0.49. The first nine eigenvectors of each individual mode are shown and compared in Fig. 11. Small errors in the compressive estimate are because the missing cone is still not perfectly compensated by the compressive approach, and because of other experimental imperfections.

In this classical experiment, we have the benefit that direct observation of the one–dimensional object is available; thus, we were able to carry out quantitative analysis of the accuracy of the compressive estimate. In the quantum analogue of measuring a complete quantum state, direct observation would have of course not been possible, but the accuracy attained through the compressive estimate should be comparable, provided the low entropy assumption holds [19].

## 5. Discussion

In conclusion, we experimentally demonstrated compressive reconstruction of the mutual intensity function of a classical partially coherent source using phase space tomography. By exploiting the physically justifiable assumption of a quasi–pure source, both measurement and post–processing dimensionality are greatly reduced. We used the van Cittert–Zernike theorem to estimate the true mutual intensity function as a way to cross–validate the compressive reconstruction, and found indeed good agreement.

Here we followed a much simplified version of the approach described in [21] which showed that the complex operators describing the measurements should be uniformly distributed in the *n*–dimensional unit sphere, whereas we simply utilized free space propagation. The phase masks described in [21] to implement optimal sampling are outside the scope of the present work.

## Acknowledgments

The authors thank Baile Zhang, Jonathan Petruccelli, and Yi Liu for helpful discussions. Financial support was provided by Singapore’s National Research Foundation through the Center for Environmental Sensing and Modeling (CENSAM) and BioSystems and bioMechanics (BioSyM) independent research groups of the Singapore-MIT Alliance for Research and Technology (SMART) Centre, by the US National Institutes of Health and by the Chevron–MIT University Partnership Program.

## References and links

**1. **L. Mandel and E. Wolf, *Optical Coherence and Quantum Optics* (Cambridge University Press, 1995).

**2. **K. Blum, *Density Matrix Theory and Applications* (Plenum Press, 1981).

**3. **K. Itoh and Y. Ohtsuka, “Fourier-transform spectral imaging: retrieval of source information from three-dimensional spatial coherence,” J. Opt. Soc. Am. A **3**, 94–100 (1986). [CrossRef]

**4. **D. L. Marks, R. A. Stack, and D. J. Brady, “Three-dimensional coherence imaging in the Fresnel domain,” Appl. Opt. **38**, 1332–1342 (1999). [CrossRef]

**5. **J. W. Goodman, *Statistical Optics* (Wiley-Interscience, 2000).

**6. **K. A. Nugent, “Wave field determination using three-dimensional intensity information,” Phys. Rev. Lett. **68**, 2261–2264 (1992). [CrossRef] [PubMed]

**7. **M. G. Raymer, M. Beck, and D. McAlister, “Complex wave-field reconstruction using phase-space tomography,” Phys. Rev. Lett. **72**, 1137–1140 (1994). [CrossRef] [PubMed]

**8. **C. Q. Tran, A. G. Peele, A. Roberts, K. A. Nugent, D. Paterson, and I. McNulty, “X-ray imaging: a generalized approach using phase-space tomography,” J. Opt. Soc. Am. A **22**, 1691–1700 (2005). [CrossRef]

**9. **M. Beck, M. G. Raymer, I. A. Walmsley, and V. Wong, “Chronocyclic tomography for measuring the amplitude and phase structure of optical pulses,” Opt. Lett. **18**, 2041–2043 (1993). [CrossRef] [PubMed]

**10. **K. Vogel and H. Risken, “Determination of quasiprobability distributions in terms of probability distributions for the rotated quadrature phase,” Phys. Rev. A **40**, 2847–2849 (1989). [CrossRef] [PubMed]

**11. **D. T. Smithey, M. Beck, M. G. Raymer, and A. Faridani, “Measurement of the Wigner distribution and the density matrix of a light mode using optical homodyne tomography: application to squeezed states and the vacuum,” Phys. Rev. Lett. **70**, 1244–1247 (1993). [CrossRef] [PubMed]

**12. **U. Leonhardt, “Quantum–state tomography and discrete Wigner function,” Phys. Rev. Lett. **74**, 4101–4105 (1995). [CrossRef] [PubMed]

**13. **C. Kurtsiefer, T. Pfau, and J. Mlynek, “Measurement of the Wigner function of an ensemble of Helium atoms,” Nature (London) **386**, 150–153 (1997). [CrossRef]

**14. **E. Candès, J. Romberg, and T. Tao, “Robust uncertainty principles: exact signal reconstruction from highly incomplete frequency information,” IEEE Trans. Inform. Theory **52**, 489–509 (2006). [CrossRef]

**15. **E. Candès, J. Romberg, and T. Tao, “Stable signal recovery from incomplete and inaccurate measurements,” Comm. Pure Appl. Math. **59**, 1207–1223 (2006). [CrossRef]

**16. **D. L. Donoho, “Compressed sensing,” IEEE Trans. Inform. Theory **52**, 1289–1306 (2006). [CrossRef]

**17. **E. J. Candès and B. Recht, “Exact matrix completion via convex optimization,” Found. Comput. Math. **9**, 717–772 (2009). [CrossRef]

**18. **E. J. Candès and T. Tao, “The power of convex relaxation: near-optimal matrix completion,” IEEE Trans. Inform. Theory **56**, 2053–2080 (2010). [CrossRef]

**19. **D. Gross, Y.-K. Liu, S. T. Flammia, S. Becker, and J. Eisert, “Quantum state tomography via compressed sensing,” Phys. Rev. Lett. **105**, 150401 (2010). [CrossRef]

**20. **D. Gross, “Recovering low-rank matrices from few coefficients in any basis,” IEEE Trans. Inf. Theory **57**, 1548–1566 (2011). [CrossRef]

**21. **E. J. Candès, T. Strohmer, and V. Voroninski, “Phaselift: exact and stable signal recovery from magnitude measurements via convex programming,” ArXiv: 1109.4499v1 (2011).

**22. **E. J. Candès, Y. Eldar, T. Strohmer, and V. Voroninski, “Phase retrieval via matrix completion,” ArXiv: 1109.0573 (2011).

**23. **Y. Shechtman, Y. C. Eldar, A. Szameit, and M. Segev, “Sparsity based sub-wavelength imaging with partially incoherent light via quadratic compressed sensing,” Opt. Express **19**, 14807–14822 (2011). [CrossRef] [PubMed]

**24. **E. Wolf, “New theory of partial coherence in the space-frequency domain. Part I: spectra and cross spectra of steady-state sources,” J. Opt. Soc. Am. **72**, 343–351 (1982). [CrossRef]

**25. **D. Pelliccia, A. Y. Nikulin, H. O. Moser, and K. A. Nugent, “Experimental characterization of the coherence properties of hard x-ray sources,” Opt. Express **19**, 8073–8078 (2011). [CrossRef] [PubMed]

**26. **K.-H. Brenner, A. Lohmann, and J. Ojeda-Castañeda, “The ambiguity function as a polar display of the OTF,” Opt. Commun. **44**, 323–326 (1983). [CrossRef]

**27. **K.-H. Brenner and J. Ojeda-Castañeda, “Ambiguity function and Wigner distribution function applied to partially coherent imagery,” Opt. Acta. **31**, 213–223 (1984). [CrossRef]

**28. **J. Tu, “Wave field determination using tomography of the ambiguity function,” Phys. Rev. E **55**, 1946–1949 (1997). [CrossRef]

**29. **E. J. Candès and Y. Plan, “Matrix completion with noise,” ArXiv: 0903.3131 (2009).

**30. **J.-F. Cai, E. J. Candès, and Z. Shen, “A singular value thresholding algorithm for matrix completion,” ArXiv: 0810.3286 (2008).

**31. **A. Starikov and E. Wolf, “Coherent-mode representation of Gaussian Schell-model sources and of their radiation fields,” J. Opt. Soc. Am. **72**, 923–928 (1982). [CrossRef]

**32. **A. C. Kak and M. Slaney, *Principle of Computerized Tomographic Imaging* (Society for Industrial and Applied Mathematics, 2001). [CrossRef]

**33. **M. J. Bastiaans, “The Wigner distribution function applied to optical signals and systems,” Opt. Commun. **25**, 26–30 (1978). [CrossRef]

**34. **M. J. Bastiaans, “Application of the Wigner distribution function to partially coherent light,” J. Opt. Soc. Am. A **3**, 1227–1238 (1986). [CrossRef]

**35. **A. Starikov, “Effective number of degrees of freedom of partially coherent sources,” J. Opt. Soc. Am. **72**, 1538–1544 (1982). [CrossRef]

**36. **M. J. Bastiaans, “New class of uncertainty relations for partially coherent light,” J. Opt. Soc. Am. A **1**, 711–715 (1984). [CrossRef]