Efficient modal analysis using compressive optical interferometry

Interferometry is routinely used for spectral or modal analysis of optical signals. By posing interferometric modal analysis as a sparse recovery problem, we show that compressive sampling helps exploit the sparsity of typical optical signals in modal space and reduces the number of required measurements. Instead of collecting evenly spaced interferometric samples at the Nyquist rate followed by a Fourier transform as is common practice, we show that random sampling at sub-Nyquist rates followed by a sparse reconstruction algorithm suffices. We demonstrate our approach, which we call compressive interferometry, numerically in the context of modal analysis of spatial beams using a generalized interferometric configuration. Compressive interferometry applies to widely used optical modal sets and is robust with respect to noise, thus holding promise to enhance real-time processing in optical imaging and communications. © 2015 Optical Society of America OCIS codes: (260.3160) Interference; (120.3180) Interferometry; (030.4070) Modes; (050.4865) Optical vortices. References and links 1. A. F. Abouraddy, T. M. Yarnall, and B. E. A. Saleh, “Generalized optical interferometry for modal analysis in arbitrary degrees of freedom,” Opt. Lett. 37, 2889–2891 (2012). 2. J. Wang, J. Yang, I. M. Fazal, N. Ahmed, Y. Yan, H. Huang, Y. Ren, Y. Yue, S. Dolinar, M. Tur, and A. E. Willner, “Terabit free-space data transmission employing orbital angular momentum multiplexing,” Nat. Photon. 6, 488–496 (2012). 3. N. Bozinovic, Y. Yue, Y. Ren, M. Tur, P. Kristensen, H. Huang, A. E. Willner, and S. Ramachandran, “Terabitscale orbital angular momentum mode division multiplexing in fibers,” Science 340, 1545–1548 (2013). 4. J. Demas and S. Ramachandran, “Sub-second mode measurement of fibers using C2 imaging,” Opt. Express 22, 23043–23056 (2014). 5. L. Allen, M. W. Beijersbergen, R. J. C. Spreeuw, and J. P. Woerdman, “Orbital angular momentum of light and the transformation of Laguerre-Gaussian laser modes,” Phys. Rev. A 45, 8185–8190 (1992). 6. A. F. Abouraddy, T. M. Yarnall, and B. E. A. Saleh, “Angular and radial mode analyzer for optical beams,” Opt. Lett. 36, 4683–4685 (2011). 7. V. Namias, “The fractional order Fourier transform and its applications to quantum mechanics,” IMA J. Appl. Math. 25, 241–265 (1980). 8. H. M. Ozaktas, Z. Zalevsky, and M. A. Kutay, The Fractional Fourier Transform (Wiley, 2001). 9. V. Namias, “Fractionalization of Hankel transforms,” IMA J. Appl. Math. 26, 187–197 (1980). 10. L. Yu, Y. Lu, X. Zeng, M. Huang, M. Chen, W. Huang, and Z. Zhu, “Deriving the integral representation of a fractional Hankel transform from a fractional Fourier transform,” Opt. Lett. 23, 1158–1160 (1998). 11. E. J. Candès, J. K. Romberg, and T. Tao, “Robust uncertainty principles: Exact signal reconstruction from highly incomplete frequency information,” IEEE Trans. Inf. Theory 52, 489–509 (2006). #242828 Received 8 Jul 2015; revised 15 Sep 2015; accepted 21 Sep 2015; published 22 Oct 2015 (C) 2015 OSA 2 Nov 2015 | Vol. 23, No. 22 | DOI:10.1364/OE.23.028449 | OPTICS EXPRESS 28449 12. E. J. Candès, “Compressed sensing,” in Proceedings of the International Congress of Mathematicians (2006), M. Sanz-Solé, J. Soria, J. L. Varona, and J. Verdera , eds. (Eur. Math. Soc., 2007), pp. 1433–1452. 13. E. J. Candès, J. Romberg, and T. Tao, “Stable signal recovery from incomplete and inaccurate measurements,” Comm. Pure Appl. Math. 59, 1207–1223 (2006). 14. D. L. Donoho, “Compressed sensing,” IEEE Trans. Inf. Theory 52, 1289–1306 (2006). 15. E. Candès and J. Romberg, “Sparsity and incoherence in compressive sampling,” Inverse Probl. 23, 969–985 (2007). 16. E. J. Candès, “The restricted isometry property and its implications for compressed sensing,” C. R. Acad. Sci. Paris 346, 589–592 (2008). 17. M. A. Davenport, P. T. Boufounos, M. B. Wakin, and R. G. Baraniuk, “Signal processing with compressive measurements,” IEEE J. Sel. Top. Signal Process. 4, 445–460 (2010). 18. M. F. Duarte, M. A. Davenport, D. Takhar, J. N. Laska, T. Sun, K. F. Kelly, and R. G. Baraniuk, “Single-pixel imaging via compressive sampling,” IEEE Signal Process. Mag. 25, 83–91 (2008). 19. G. A. Howland, J. Schneeloch, D. J. Lum, and J. C. Howell, “Simultaneous measurement of complementary observables with compressive sensing,” Phys. Rev. Lett. 112, 253602 (2014). 20. M. Mirhosseini, O. S. Magaña-Loaiza, S. M. H. Rafsanjani, and R. W. Boyd, “Compressive direct measurement of the quantum wave function,” Phys. Rev. Lett. 113, 090402 (2014). 21. S. S. Chen, D. L. Donoho, and M. A. Saunders, “Atomic decomposition by basis pursuit,” SIAM J. Sci. Comput. 20, 33–61 (1998). 22. M. Saunders, “PDCO: Primal-Dual Interior Method for Convex Objectives,” http://www.stanford.edu/group/SOL/software/pdco. 23. E. J. Candès, and T. Tao, “Decoding by Linear Programming,” IEEE Trans. Inf. Theory 51, 4203–4215 (2005). 24. J. Tropp, “Greed is good: Algorithmic results for sparseapproximation, IEEE Trans. Inf. Theory 50, 2231–2242 (2004). 25. W. Dai, O. Milenkovic, “Subspace pursuit for compressive sensing signal reconstruction, IEEE Trans. Inf. Theory 55, 2230–2249 (2009). 26. J. A. Rodrigo, T. Alieva, and M. L. Calvo, “Programmable two-dimensional optical fractional Fourier processor,” Opt. Express 17, 4976–4983 (2009).


Introduction
Optical beams offer the potential for carrying a high-information content by exploiting the large-dimensional space spanned by its physical degrees of freedom (DoFs). The spatial DoF has attracted particular interest with recent advances in the synthesis and analysis of beams having complex spatial profiles [1]. Indeed, spatial multiplexing for high-speed communications in free space [2] and in multimode fibers [3] has brought to the fore the importance of accurate and rapid modal analysis in a desired basis [4], such as the that of orbital angular momentum (OAM) modes [5].
We recently proposed an interferometric procedure that allows -in principle -for an optical beam to be analyzed in terms of a complete and orthogonal -but otherwise arbitrary -modal basis [1,6]. We call this approach henceforth 'generalized interferometry'. Such a strategy exploits a two-path interferometer (such as a Mach-Zhender interferometer, MZI), but replaces the usual optical delay with a 'generalized phase operator' (GPO) -a unitary spatial transformation parameterized by a continuous real number that plays the role of a 'generalized delay' in modal space. The GPO is in fact an optical transformation whose eigenfunctions are the functional elements of the modal basis of interest. Indeed, the GPO generalizes to an arbitrary basis the notion of temporal delay. The GPO delay parameter in the generalized interferometer is swept, an interferogram is recorded, and its Fourier transform reveals the beam's modal content. In the case of discrete modal bases, the GPO is a fractional optical transform; e.g., the fractional Fourier transform (fFT) [7,8] or fractional Hankel transform (fHT) [9,10] for Hermite-Gaussian (HG) or radial Laguerre-Gaussian (LG) modes, respectively. In practice, measurements are acquired by sampling the GPO delay -the order of the associated fractional transform -at the Nyquist rate to avoid aliasing in modal analysis. This requires collecting a large number of samples and implies more latency, which may be intolerable for delay-sensitive applications.
In this paper, we introduce a novel strategy we call compressive interferometry that exploits the modal sparsity of typical optical beams -the small number of modes carrying significant energy out of the potentially accessible modes -to reduce the number of measurements required to determine the modal weights. Recent work has made use of compressive sensing (CS) approaches [11][12][13][14][15][16][17] for optical field estimation and imaging -and not modal analysisby introducing a sequence of random masks in the field path [18][19][20]. Therein, compressive sensing is achieved by applying M random masks (each with N random pixel patterns) to emulate the random ensembles that are typically used in standard CS. In contrast, we exploit CS here in the context of 'native' optics hardware, namely two-path interferometry, with no need for extra hardware -for the purpose of modal analysis. Instead of recording an interferogram sampled uniformly at the Nyquist rate, the number of required samples is significantly reduced by sub-Nyquist sampling of the interferogram at randomly chosen values of the GPO delay parameter and exploiting the linear CS model to harness the sparse representation of the beam in modal space. These randomly collected measurements are viewed as the result of operating with a random transformation -or restricted sensing matrix -on the sparse vector of modal weights. Unlike the sensing matrices typically studied in CS that are drawn from random ensembles [11], the sensing matrix here has only a few degrees of freedom constrained by the interferometer itself. Nevertheless, we show that the modal weights are recovered via a reconstruction algorithm such as Basis Pursuit (BP) [21].
The paper is organized as follows. Generalized interferometry for modal analysis in arbitrary basis is presented in Section 2. Modal analysis is then posed as a sparse recovery problem and our compressive interferometry approach is presented in Section 3. Section 4 provides representative numerical examples for the reconstruction of 1D and 2D optical beams in HG, LG, and OAM modal bases using the proposed approach.

Generalized interferometry
To set the stage for developing the concept of compressive interferometry, we first describe the reconstruction of a beam's modal content using the modified MZI shown in Fig. 1(a) [1,6]. In traditional two-path interferometry, an optical delay is swept to produce an interferogram whose Fourier transform yields the spectrum. In other words, a traditional interferometer is capable of analyzing an optical beam in terms of the temporal harmonics e −iωt . We pose the following question: does the same general strategy apply to other modal bases, such as spatial modes? In that case, what optical device should replace the usual optical delay in order to implement a 'delay' in the space spanned by the modal basis of interest?
Consider a paraxial optical beam E(x, y) in a given plane written as a linear superposition of the elements of orthonormal modal bases {ψ n (x)} n and {φ m (y)} m with coefficients c nm , Here, the upper limits N 1 and N 2 are set by the optical system; e.g., the numerical aperture. Our goal is to analyze the input beam into the two modal bases {ψ n (x)} n and {φ m (y)} m and identify the modal weights c nm . To achieve this goal, we search for optical transformations that generalize the traditional role of the optical delay in a typical two-path interferometer. It has been shown [1,6] that this goal can be achieved in a generalized interferometer by replacing the optical delay with GPOs corresponding to the basis sets {ψ n (x)} n and {φ m (y)} m having delay parameters α and β , respectively. The GPO associated with {ψ n (x)} n , for example, is the optical system implementing the unitary spatial transformation having {ψ n (x)} n as its eigenfunctions, which adds a phase nα to the n th mode; a second transformation is similarly constructed for {φ m (y)} m . The impact of such a GPO on the associated modal basis is analogous to that of a traditional temporal delay on a time harmonic. A harmonic e −iωt emerges from the delay without change -except for picking up a phase of the form e iωτ , where τ is the delay value. Likewise, a mode ψ m (x) emerges from the GPO described by Eq. (2) with no change other than picking up a phase of the form e inα with α now playing the role of the 'delay'. The specific optical arrangement to implement a GPO is determined by the modal basis of interest. Several examples will help clarify this point. For example, if {ψ n (x)} n corresponds to HG modes along a Cartesian coordinate, then the associated GPO is the fractional Fourier transform (fFT). This can be seen immediately by realizing that HG modes of order m are eigenfunctions of the fFT of order α (with α ranging from 0 to 4) with eigenvalue e i π 2 αm [1]. Alternatively, for radial LG modes in polar coordinates, the GPO in Eq. (2) is identified as the order-α fHT [6]. Furthermore, when considering OAM states of order , the associated GPO corresponds to an angular rotation α, which leaves an OAM mode unchanged except for a phase of the form e i θ . In all these cases, the real, continuous parameter α -the order of the fractional transform -now plays the role of the temporal delay in a traditional interferometer. The interferogram is given by and we recover the modal coefficients by taking a 2D Fourier transform (FT) of P(α, β ). Note that this equation applies to all the different modal bases described above, and hence this strat-egy of generalized interferometry can be considered 'basis neutral'. To cast this generalized interferometry scheme into the form of a linear sensing model [12], we write the interferogram in matrix form. Assume, for simplicity, a 1D modal basis with maximum dimension N (known as the ambient dimension) and that M measurement points are sampled along the interferogram, then the discretized components of Eq. (3) may be written in matrix form, ⎡ The M×1 vector y is called the measurement vector, the M×N matrixΦ the sensing matrix, and the N×1 vector x the modal-coefficient vector. Usually M > N; e.g., M ≥ 2N to achieve the Nyquist rate. Each sample P(α i ) can be written as an inner product, where Φ i is a row vector whose n-th element is cos nα i , n = 1,...,N. The mathematical model in Eq. (4) is depicted graphically in Fig. 1(b). The interferogram samples for two modal bases (2D beam profile) are where α i and β j are the GPO parameters corresponding to the two spatial DoFs, the measurements are assembled in an M×1 vector y (M = M 1 ×M 2 ), x is an N-dimensional sparse modalcoefficient vector (N =N 1 ×N 2 ), and Φ i, j is a 1×N row vector with elements cos(nα i +mβ j ), n= 1, ..., N 1 ; m = 1, ..., N 2 , in the M×N sensing matrixΦ. It is natural to ask whether one may recover the modal weights with fewer measurements when the beam is sparse in modal space [ Fig. 1(c)].

Compressive sensing model for generalized interferometry
Sparse signals are often vectors in a high-dimensional space with only few non-zero elements -known as the support of the signal. A vector with at most s non-zero elements is called 's-sparse'. The linear CS model shows that we can recover this sparse signal in a space of ambient dimension N from only M N linear measurements -that is, at sub-Nyquist sampling rate [12,14]-provided that the sensing matrix satisfies some sufficient conditions. Generalized interferometry [Eq. (3)] is posed here as a linear measurement problem, y =Φx [Eq. (4)]. When the number of interferometric measurements M is smaller than the ambient dimension N, i.e., M N, modal analysis is equivalent to solving an under-determined system of linear equations, which is generally an ill-posed problem since it has less equations than unknowns. However, searching for a solution is feasible if the beam is sparse in modal space. In particular, we can search for the most sparse solution that minimizes the 0 -norm, x 0 , of the modal-coefficient vector x, subject to the data constraint, i.e., min x 0 subject to y =Φx, where x 0 is equal to the number of non-zero entries of x [11]. Nonetheless, the problem in Eq. (7) is generally NP-hard as it involves a combinatorial search over all s-sparse vectors. The complexity can be significantly reduced using a convex relaxation of Eq. (7) min x 1 subject to y =Φx, where x 1 = ∑ N n=1 |x n | is the 1 -norm of x, and x n = |c n | 2 , n = 1, ..., N are the energies of the different modes. This 1 -minimization, known as the Basis Pursuit (BP), can be reduced to a simpler Linear Programming (LP) problem, then solved using an appropriate technique such as the Primal-Dual Interior-Point Method for Convex Objectives (PDCO) [22] adopted in this paper. Furthermore, the solution to the minimization in Eq. (8) recovers x exactly provided that x is sufficiently sparse, and that the sensing matrixΦ satisfies some sufficient conditions such as the Restricted Isometry Property (RIP) [23].Φ is said to satisfy the RIP with isometry constant for all s-sparse vectors x. This definition requires a low-rank matrixΦ to behave as a unitary matrix with respect to all s-sparse vectors. Although it is generally hard to verify if a given matrix satisfies the RIP, random matrices drawn from Gaussian ensembles satisfy the RIP with overwhelming probability when M ≈ O(s log N) [16]. Following standard notation, we write f (x) = O(g(x)) as x → ∞, if and only if there is a positive constant C and a number x 0 such that Despite the fact that we cannot guarantee that the restricted sensing matrixΦ in Eq. (4) inherits the RIP, we conjecture that if x admits an s-sparse representation in an N-dimensional modal space, it may still be possible to reconstruct x from M N interferogram samples. This conjecture is reasonable becauseΦ is generated through random GPO settings, and therefore is rich enough to capture the corresponding behavior of random ensembles. That is, by viewing modal analysis as a sparse recovery problem, sub-Nyquist sampling may provide enough information to successfully solve an under-determined system of linear equations and recover the modal content. We have validated our conjecture with the simulations provided below. We also note that, while BP is our algorithm of choice in this paper, other reconstruction algorithms such as Matching Pursuit (MP) [24] and Subspace Pursuit [25] could also be applied to recover the active modes using the proposed compressive interferometry approach.
To assess the robustness of the modal reconstruction and the impact of noise, we rewrite the compressive interferometry measurements as y =Φx + n, where n is an additive noise vector modeled as white Gaussian noise N(0, σ 2 n I), σ 2 n = (||y|| 2 2 /M)/10 SNR 10 , and SNR denotes the signal-to-noise ratio in dB. In our case, it is unclear a priori whetherΦ will be favorable for reconstruction, e.g., satisfy the RIP, since some of the rows may not be linearly independent. In the presence of noise, this could degrade the performance of existing reconstruction algorithms that typically involve inversion of submatrices. To recover the independent rows ofΦ, we use a singular value decomposition (SVD)Φ =ÛΣV T , whereÛ andV are unitary bases for the column and row spaces ofΦ, respectively, andΣ is the diagonal matrix of singular values. Zero or small singular values indicate thatΦ has dependent or nearly-dependent rows. Hence, we discard those measurements that correspond to the small singular values by thresholding the ratios of the singular values to the largest one. The threshold is chosen based on the noise level and the spread of the eigenvalues of the sensing matrix. Therefore, the compressive interferometry problem can be remodeled using a truncated SVD as, y =Û T o y =Σ oV T x + n , whereΣ o , and U o denote the truncated versions ofΣ,Û, respectively, the superscript 'T' denotes transposition, and n is additive Gaussian noise with the same distribution as n. Our numerical results confirm that the errors in modal reconstruction are significantly reduced using this truncated SVD.

Sparse modal reconstruction
We now illustrate the performance of compressive interferometry in different examples. The modal-recovery error is defined as a scaled distance between the true modal-coefficient vector x and the recovered vectorx, ||x−x|| 2 /||x|| 2 . We reconstruct sparse beams from a reduced number of measurements using the BP algorithm [21] with the GPO delay samples (α and β ) randomly chosen from a uniform probability distribution U[−π, π]. BP uses an 1 -relaxation of the 0 -norm to promote sparse solutions.
The numerical results show that when the GPO delays are chosen randomly,Φ preserves most of the information required to recover x when M N and a sparse reconstruction algorithm replaces the FT.
First consider the 1D beams E 1 (x) and E 2 (x) formed of superpositions of s = 4 HG modes selected out of N = 100 modes [Figs. 2(a)-2(d)]. Since the modal basis of interest is that of HG modes, the associated GPO is the fFT, and the order of the fFT plays the role of the GPO delay parameter. Using compressive interferometry in presence of noise (SNR=10 dB), the recovery error is seen to drop sharply when M > 35 samples are randomly selected [ Fig. 2(e)] whether the beam comprises low-or high-order modes. This reconstruction compares favorably with that obtained from the FT of an interferogram produced by the GPO-based setup in Fig. 1(a) sampled at the Nyquist rate (200 evenly spaced samples) as shown in Fig. 2(f) over a wide range of SNR values. Compressive interferometry therefore substantially reduces M while remaining robust to additive noise. Figure 3 shows the performance of the proposed compressive approach for a setting with unequal modal weights. Using about M = 35 measurements is sufficient for modal reconstruction at SNR=10 dB and SNR=20 dB as shown in Fig. 3 (c). Also, SVD is shown to improve the performance of modal reconstruction.
We next evaluate the performance of the proposed compressive interferometry approach with 2D beams. We use two different 2D modal bases: in Figs. 4(a) and 4(b) we examine beams E 3 (x, y) and E 4 (x, y) formed of superpositions of HG modes along Cartesian coordinates x and y, while in Figs. 4(c) and 4(d) we consider beams E 5 (r, θ ) and E 6 (r, θ ) formed of superpositions of LG-OAM modes along the radial r and azimuthal θ DoFs in a polar coordinate system. In the case of 2D beams, two separate GPOs are required, one for each dimension. For 2D HG modes, the two GPOs are two fFT systems, one for the x coordinate (with delay parameter α x ) and the other for the y coordinate (with delay parameter α y ). The interferogram P(α x , α y ) is recorded by sweeping the values of α x and α y -at the Nyquist rate -and a 2D Fourier transform reveals the modal weights |c nm | 2 . For the case of radial LG and OAM modes in polar coordinates, the GPOs consist of a fHT for the radial coordinate and a angular rotation for the azimuthal coordinate, respectively. We proceed to examine the impact of the compressive sensing model on reducing the number of measurements required. In each example, we consider ambient dimensions N 1 =10 and N 2 = 10 for each DoF. In each coordinate system, we consider a beam formed of low-order modes and one formed of high-order modes. In all cases (with SNR = 10 dB), we find that M ≈ 35 measurements are sufficient to reduce the reconstruction error. We consider the tradeoff between the number of measurements M, the SNR, and the error performance in Figs. 4(e) and 4(f). Several general features emerge from our numerical results. First, in Fig. 4(e) we see that when a sufficient number of measurements are selected, the reconstruction error drops with increasing SNR (M =35, 40). If the number of measurements is too small, however, then the reconstruction error is large and independent of SNR (M =10, 15). Second, we observe a phase transition in the dependence of M on SNR for a fixed reconstruction error [ Fig. 4(f)]. That is, below a threshold value of M, the reconstruction is not satisfactory for any value of SNR. This is a feature of all CS approaches when an insufficient number of measurements are collected [13]. Finally, note that fractional transforms, such as the fFT, may be implemented using spatial light modulators [26] such that random sampling of the GPO delay is performed electronically and requires no moving parts.

Conclusion
In conclusion, we have presented a general procedure, compressive interferometry, based on the CS model applied to a specific interferometric arrangement, but formulated in such a way that it may be readily adapted to a variety of settings. We leverage the sparse representation of optical beams in modal space by posing spatial-mode analysis as a sparse recovery problem. This strategy requires sub-Nyquist-rate randomized measurements, is resistant to noise, and is applicable to a variety of modal sets -and thus holds potential to enhance real-time processing in optical imaging and communications.