Open Access
Add to BibSonomyAbstract
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 n2 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(rnlnn) to O(rnln2 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]
This can be expressed in operator form as where 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 xo denotes the lateral coordinate at the observation plane. By changing variables x = (x1 + x2)/2, x′ = x1 – x2 and Fourier transforming the mutual intensity with respect to x we obtain the Ambiguity Function (AF) [26–28] Equation (2) can be written as [6–8, 26, 28], where Ĩ is the Fourier transform of the vector of measured intensities with respect to xo. 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 Ĵ, the method can be written as
Direct rank minimization is NP–hard; however, it can be accomplished by solving instead a proxy problem: convex minimization of the “nuclear norm” (ℓ1 norm) of the matrix J [17, 29]. The corresponding problem is stated as where the nuclear norm is the sum of the singular values σi = |λi|, ||Ĵ||* = ∑i σ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.e. ∑i λi < 1; we compensated for this by renormalizing at the end of each iteration [19].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]
and where σI determines the spatial extent of the source, and σc 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) 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 σc = 13 (rank 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.
Fig. 1 (a) Input mutual intensity of a GSMS with paramters σI = 17 and σc = 13, (b) data point locations in the Ambiguity space, mutual intensities estimated by (c) FBP and (d) LRMR methods.
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
To implement the FBP method, each intensity projection is first filtered by a Ram–Lak kernel apodized by a Hamming window; the estimated WDF is obtained by back-projecting all the filtered intensities, and then an inverse Fourier transform is applied to produce the mutual intensity. Figure 1(c) shows the reconstructed mutual intensity following this procedure. Three types of artifacts can be seen in this reconstruction. First, the reconstructed mutual intensity has lower values along the diagonal of the matrix due to the missing cones. However, this is un-physical because a correlation function should always have maximum value at zero separation. Second, the estimated degree of coherence is lower than the original field. The third artifact is the high frequency noise around the diagonal of the matrix, which is due to undersampling between the radial slices. All these artifacts have been greatly suppressed or completely removed by LRMR, whose reconstruction result is shown in Fig. 1(d). The disappearance of the correlation peak along the diagonal (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 [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).
Fig. 2 The first nine coherent modes of the mutual intensity in Fig. 1(a). (a) Theoretical modes, and (b) LRMR estimates.
Fig. 3 Eigenvalues of the mutual intensity in Fig. 1(a). (a) Theoretical values, (b) FBP estimates, (c) LRMR estimates, and (d) absolute errors in the LRMR estimates versus mode index.
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 σc = 18 (rank 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.
Fig. 4 Oversampling rate versus relative MSE of LRMR estimates. The input field is a GSMS with parameters σI = 36 and σc = 18. The noisy data is generated with different SNR from (a) an additive random Gaussian noise model, and (b) a Poisson noise model.
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.
Fig. 7 (a) Real and (b) imaginary parts of the radial slices in Ambiguity space from Fourier transforming the vectors of intensities measured at corresponding propagation distances.
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.
Fig. 10 (a) Intensity measured immediately to the right of the illumination slit; (b) real part of van Cittert–Zernike theorem estimated mutual intensity immediately to the right of the object slit; (c) eigenvalues of the mutual intensity in (b); (d) absolute error between the eigenvalues in Fig. 9(b) and 10(c) versus mode index.
Fig. 11 (a) LRMR estimated coherent modes of the mutual intensities in Fig. 8(b), and (b) coherent modes of the mutual intensities in Fig. 10(b), calculated via use of the van Cittert–Zernike theorem, and assumption of incoherent illumination.
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]
