Compressive optical interferometry
Davood Mardani,1, ∗ H. Esat Kondakci,2, ∗ Lane Martin,2 Ayman F. Abouraddy,2 and George K. Atia1
1 Dept.

of Electrical and Computer Engineering, University of Central Florida, Orlando, FL 32816, USA
The College of Optics & Photonics, University of Central Florida, Orlando, Florida 32816, USA

2 CREOL,

arXiv:1706.06166v1 [] 19 Jun 2017

Compressive sensing (CS) combines data acquisition with compression coding to reduce the number of
measurements required to reconstruct a sparse signal. In optics, this usually takes the form of projecting
the field onto sequences of random spatial patterns that are selected from an appropriate random
ensemble. We show here that CS can be exploited in ‘native’ optics hardware without introducing
added components. Specifically, we show that random sub-Nyquist sampling of an interferogram helps
reconstruct the field modal structure. The distribution of reduced sensing matrices corresponding to
random measurements is provably incoherent and isotropic, which helps us carry out CS successfully.

Compressive sensing (CS) is a strategy for reducing the
number of measurements required to reconstruct a signal
by projecting it onto a basis of random functions, thereby
combining the two steps of data acquisition and compression
coding [1, 2]. This procedure is particularly effective when
the signal is sparse – i.e., the signal receives contributions
from a small number of basis functions in some representation [1]. Underlying this approach is an under-determined
linear transformation between the sparse representation
and the measurement basis called a sensing matrix that is
amenable to a stable pseudoinversion [1, 3]. CS has found
many uses in optics such as the ‘single-pixel’ camera [4, 5],
holography [6–9], optical imaging and tomography [10–13],
optical encryption [14, 15], spectral analysis [16], and efficient quantum state tomography [17–20]. In all these cases,
the native optical system is modified by incorporating a
new component that projects the field onto random patterns.
CS then reduces the number of data samples acquired by
the modified system that are needed to recover the target
sparse vector.
In this paper, we develop a compressive approach to interferometry for optical modal analysis in arbitrary bases
and realize it experimentally. We underscore four major differences from related work. First, unlike much of the prior
work which introduces random masks along the optical
path to obtain randomized measurements [6, 20], here we
exploit CS in the native optics hardware without modification to the underlying interferometer structure nor addition
of new hardware components. Second, this is the first work
to realize compression in the generalized interferometry framework introduced recently by some of the authors to carry
out modal analysis in arbitrary domains [21, 22]. Third,
we achieve compression gains both in the number of measurements used for reconstruction and the acquisition time.
This is in contrast to prior work [23] employing CS in optical interferometry, which solely focused on reducing the
number of measurements. Fourth, we establish analytical
guarantees for successful modal reconstruction under sensing constraints set by the limited degrees of freedom of an
optical interferometer.
We introduce CS into a prototypical optical system – a
two-path interferometer – and demonstrate experimentally

∗

These authors contributed equally.

that sampling the interferogram below the Nyquist rate still
enables spectral/modal recovery. We first cast the procedure of interferometric spectral/modal reconstruction as a
linear measurement problem. We find that the reduced sensing matrix associated with random sub-Nyquist sampling
satisfies sub-optimal conditions of isotropy and incoherence, thereby enabling compressive reconstruction from
fewer samples [24]. We denote this overall optical measurement scheme compressive interferometry since it combines
reduced data acquisition (reduced number of interferometric points sampled) with compressive coding (replacing
harmonic analysis with CS reconstruction).
For sake of generality, we employ a two-path interferometer in which the usual temporal delay is replaced by
a ‘generalized delay’ that operates in the Hilbert space defined by any modal basis of interest [21, 25]. Such a delay corresponds to an optical transformation that has this
modal basis as eigenfunctions, which reduces to the traditional temporal delay when spectral harmonics are of
interest. The ‘delays’ for Hermite-Gaussian (HG) or radial
Laguerre-Gaussian (LG) modal bases are optical realizations of the fractional Fourier transform (FrFT) [26, 27] or
fractional Hankel transform (FrHT) [28], respectively. In all
cases, the compressive interferometry strategy is to replace
harmonic analysis of the Nyquist-sampled interferogram
with random sub-Nyquist sampling without modifying the
native optical hardware. By constructing an intrinsically
stable, basis-neutral interferometer that includes realizations of fractional transforms [22], we confirm that CS helps
successfully reconstruct the modal weights of sparse optical beams in a basis of HG or LG modes from the random
sub-Nyquist samples. This may help in delay-intolerant applications that require real-time acquisition and processing,
such as real-time imaging and high-speed communications
exploiting spatial multiplexing.
A typical optical implementation of CS is depicted in
Fig. 1(a) [4, 17, 18]. Consider an optical field E( x ) =
∑n cn ψn ( x ) represented in an N-dimensional orthonormal
basis {ψn ( x )}, and E( x ) is normalized such ∑ |cn |2 = 1.
Here x may refer to one or more spatial dimensions, or even
time. The N × 1-vector x = {cn } of coefficients is s-sparse
if it receives contributions from ≤ s elements. In lieu of the
M ≥ 2N Nyquist-rate measurements normally required to
reconstruct x, CS leverages sparsity to reduce M by projecting the field onto a sequence of random masks selected from
a Gaussian ensemble to produce a new M × 1 measurement

2

Figure 1 | Concept of compressive optical interferometry.
(a) Schematic of the traditional CS scheme in optics, where the field E( x ) is subject to random projections before measurements. (b)
The mathematical concept used in traditional CS. The field is represented by a M × 1 vector x that is transformed by a M × N sensing
matrix Φ (with disordered entries corresponding to the random masks) to yield a M × 1 measurement vector y, from which x is
reconstructed. (c) The compressive interferometry scheme. The input beam E( x ) is directed to a Mach-Zhender interferferometer in
which a generalized delay α replaces the usual temporal delay. Two copies of E( x ) are created at a beam splitter and the modes
{ψn ( x )} underlying the beam acquire phase shifts of the form einα after passing through the generalized delay to yield a new beam
E( x; α). The original and delayed beams combine at another beam splitter, and their superposition is integrated by a bucket detector to
produce an interferogram. (d,e) Graphical depiction of the matrix form of a generalized interferogram. (d) With evenly spaced
Nyquist-rate sampling of the interferogram M ≥ 2N, the interferometric sensing matrix Φint has a well-defined deterministic structure.
(e) With sub-Nyquist randomly sampled points α, Φint now appears disordered and similar to Φ in (b).

vector y. Effectively, the field has undergone a linear transformation represented by a M × N sensing matrix Φ, y = Φx
[Fig. 1(b)]. CS exploits this under-determined linear system
to compressively recover x from M ≈ O(s log N )  N
linear measurements of y – provided some conditions on
Φ are satisfied. For example, Basis Pursuit (BP) recovers
x if Φ satisfies the Restricted Isometry Property (RIP) [2],
which requires that (1 − δ)kxk2 ≤ kΦxk2 ≤ (1 + δ)kxk2 ,
for any s-sparse x, where 0 < δ < 1 is a constant. Sensing matrices from Gaussian ensembles, for example, satisfy
the RIP with overwhelming probability [2]. The random
masks introduce new components into the optical system,
thus adding to its complexity while reducing the number
of measurements acquired.
The compressive interferometry scheme relies on a different strategy. Consider an interferogram P(α) traced by
scanning a ‘delay’ α placed in one arm of a two-path interferometer [Fig. 1(c)]. With respect to the modal basis
{ψn ( x )}, such a delay is an optical transformation having a continuous real control parameter α, Λ( x, x 0 ; α) =
∑n einα ψn ( x )ψn∗ ( x 0 ), which is a fractional transform hav-

ing the basis functions as eigenstates with eigenvalues einα ,
and Λ( x, x 0 ; 0) = δ( x − x 0 ) [21]. The parameter α acts as
a ‘delay’ in the Hilbert space spanned by {ψn ( x )} just as
a delay τ does in the time domain. That α is in fact a delay can be Rreadily seen by superposing the delayed field
E( x; α) = dxΛ( x, x 0 ; α) E( x 0 ) = ∑n cn einα ψn ( x ) and the
reference E( x; 0) to produce the interferogram
P(α) ∝

Z

dx | E( x; α) + E( x; 0)|2 = 1 +

N

∑ |cn |2 cos(nα).

n =1

(1)
The modal coefficients x = {|cn |2 } are revealed by taking
a Fourier transform (FT) with respect to α. This requires
sampling α at the Nyquist rate dictated by the highest-order
mode N – even if only a few modes contribute significantly
[Fig.1 (d)].
To show that CS algorithms may recover the sparse vector
x from M  N measurements, we first cast P(α) from
Eq. 1 with selected values of α into a linear model. A M × N
interferometric sensing matrix Φint relates the N × 1 vector
x of contributions from each Hilbert-space basis-element to

3
a typical realization of the sensing matrix for a subset of
an ensemble of 106 sparse-vector realizations x with s = 4,
where maxx |η (x)| ≤ δ. A histogram constructed from the
full ensemble is plotted in Fig. 2(b). However, the sensing
matrix Φint with randomly selected α indeed belongs to a
random ensemble satisfying weaker incoherence and isotropy
properties despite its limited degrees of freedom imposed
by the interferometric process. Provably, such matrices also
yield exact reconstruction with high probability via BP [29].

Figure 2 | Restricted Isometry Property.
√
(a) Calculated values of η (x) = k( 2/M)Φint xk2 /kxk2 − 1 for
s-sparse (s = 4) vectors x. (b) Histogram of η (x) for an ensemble
of 106 realizations of s-sparse vectors. (c) The sensing matrix
Φint from Eq. 2 with randomly sampled α. (d) Φint is incoherent
with incoherence parameter µ = 1 and satisfies the isotropy
property E{φ† φ} = 0.5 I.

To clarify, let φ denote a row of the sensing matrix selected from a random ensemble F . The incoherence parameter µ(F ) is defined as the smallest number for which
maxn |φn |2 ≤ µ(F ), n = 1, 2, . . . , N, where φn are the entries of φ. The distribution F is said to obey the isotropy
property if E[φ† φ] = aI, where E[.] is the expected value, I
the identity matrix, a some constant, and φ† the conjugate
transpose of φ. Φint satisfies the conditions of incoherence
and isotropy for α randomly selected from a uniform distribution over [0, 2π ] with µ = 1, as confirmed in Fig. 2(c)
that shows Φint consisting of M rows selected from the full
matrix in Fig. 1(d). The matrix appears random in spite
of being sampled from a structured sensing system. The
matrix Φint is shown to satisfy the isotropy property in
Fig. 2(d).

We now move on to implementing the compressive interferometric scheme experimentally. The generalized delay
transformation Λ( x, x 0 ; α) has the modal basis {ψn ( x )} as
the M × 1 measurement-vector y that contains the sampled
eigenfunctions and its order parameter α represents the
interferogram [Fig. 1(e)],

 
 delay in the Hilbert space spanned by this basis. We carry

| c1 |2
out our experiments using two modal bases, 1D HG funcP ( α1 ) − 1
cos α1 cos 2α1 . . . cos Nα1
 P(α2 ) − 1   cos α2 cos 2α2 . . . cos Nα2  |c2 |2  tions and radial LG functions. The delay transformations

 



=
 ..  . Λ( x, x 0 ; α) in these two modal bases correspond to the FrFT
..
..
..
..
..





.
.
.
.
.
.  [26, 27] and the FrHT [28], respectively [21, 25]. The usual
P(α M ) − 1
cos α M cos 2α M . . . cos Nα M
| c N |2
|
{z
} |
{z
}| {z } temporal delay is replaced with an optical realization of the
appropriate fractional transform, whose order α is swept
y
x
Φint
in the interval [0, 2π ] to produce an interferogram P(α).
(2)
We construct these fractional transforms using spatial light
Here Φint corresponds to any interferometric measurement,
modulators (SLMs) that realize generalized cylindrical or
temporal or otherwise. This matrix differs fundamentally
spherical lenses to produce the FrFT or FrHT, respectively,
from those typically employed in CS applications in optics
of desired order [30]. We then exploit the polarization[Fig. 1(a) and 1(b)] where controllably random transformaselectivity of liquid-crystal-based SLMs [31] to create an
tions Φ are realized by an array of random patterns that
intrinsically stable common-path interferometer. One polarare judiciously selected such that Φ has the RIP, for examization – that undergoes the transformation implemented
ple. In interferometry, this freedom in designing Φint is
by the SLMs – is subjected to the generalized ‘delay’ while
not available because Φint has a fixed structure with few
the orthogonal polarization component is left invariant afcontrollable degrees of freedom (the values of the delay
ter the SLMs and thus corresponds to a reference [22]. A
parameter α) imposed by the constrained sensing structure
polarization projection at the output can then superpose
(the interferometer itself).
the reference and the delayed beams, thereby yielding the
Given the constrained form of Φint in Eq. 2, it can be
target interferogram P(α) as the delay parameter α is swept
shown that selecting random values of α from a uniform
evenly at the Nyquist rate.
distribution at a sub-Nyquist rate yields Φint that satisfies
the RIP with high probability when M and N are sufficiently
In the compressive interferometry approach, α is ranlarge. For smaller values of M and N, however, the RIP condomly sampled at a sub-Nyquist rate. We study the perstant δ is large. The sufficient conditions for all existing reformance of this strategy by comparing the modal reconconstruction algorithms require smaller values of δ to guarstruction in the HG and radial LG bases to that obtained
antee successful reconstruction.
via the FT. M measurements are collected randomly by
√ For example, a sufficient
condition for BP requires δ < 2 − 1. Having the large δ asselecting generalized phases α j , j = 1, 2, . . . , M from a unisociated with Φint therefore does not guarantee reconstrucform random distribution [0, 2π ], and we use the BP altion on the basis of these requirements.
As an illustration, in
gorithm to reconstruct the vector x̃ of modal coefficients.
√
Fig. 2(a) we plot η (x) = (k( 2/M )Φint xk2 /kxk2 ) − 1 with
The reconstructed vector is then compared to the vector xFT

4

Figure 3 | Reconstructed modal distributions using the FT and CS approaches.
(a) The modal weights |cn |2 calculated for fields in the HG-basis by applying the FT to Nyquist-rate evenly sampled interferograms (blue)
and CS to sub-Nyquist randomly sampled (yellow) interferograms. (b) Same as (a) for LG modes (LG0 and LG1 ) and (c) for
superposition of HG modes (HG0 +HG1 and HG1 + iHG2 ). For FT, we use 128 evenly spaced values of α from 0 to 2π. For CS, we use
M = 30 randomly selected α ∼ U [0, 2π ]. The insets show the ideal (exact) and approximated modes that are implemented
experimentally. For HG0 and LG0 , the modes are produced without approximation.

Figure 4 | Reconstruction error.

Reconstruction error versus number of measurements M. Each
curve results from averaging over 100 runs of the experiment. The
average curve (solid black) is the mean and the shaded area
designate one standard deviation spread on either side of the
required M for 1000 randomly generated examples of sparse
vectors with support size s ≤ 4.

In conclusion, we have devised a compressive approach
to interferometry whereby the modal content of light beams
can be recovered by sampling intereferograms randomly
at sub-Nyquist rate using native optics hardware, which
was demonstrated experimentally. Two fundamental differences distinguish our strategy from prior work employing
CS in optics. First, the interferometric sensing matrix is
imposed through the structure of the interferometer. As
such, compression has to be carried out under sensing constraints enforced by the limited degrees of freedom of the
sensing system. Second, previous work in the context of
interferometry focused on reducing the number of measurements used for reconstruction, but not on compressive
data acquisition. Leveraging the notion of generalized delay,
our approach enables compressive modal analysis in arbitrary modal basis. Furthermore, we established provable
guarantees for successful reconstruction despite the sensing constraints set by the limited degrees of freedom of the
interferometer.
Funding. Office of Naval Research (ONR) contracts
N00014-14-1-0260 and N00014-17-1-2458.

obtained from the FT of 128 uniformly sampled interferometric measurements by computing the scaled error metric,
||xFT − x̃||2 /||xFT ||2 . For each example, the average reconstruction error is evaluated by averaging over 100 runs. The
number of potential modes is N = 64. Figures 3(a) and 3(b)
show the reconstructed modal coefficients using FT and CS
for input beams formed from HG or LG modes. In Fig. 3(c),
we depict the reconstructed modal coefficients for two different superpositions of HG modes. The performance of the
BP algorithm is comparable to that of FT using M = 30 compressive measurements (≈ 25% of the measurements in FT),
leading to substantial savings in acquisition time without
any added complexity to the system. The reconstruction
error is plotted in Fig. 4 as a function of M. The average
error approaches zero when M ≈ 25. Statistical analysis
of the observed error reveals that the required number of
measurements lies between M = 15 to M = 30 to attain
nearly error-free reconstruction.

5

[1]
[2]
[3]
[4]
[5]
[6]
[7]
[8]
[9]
[10]
[11]
[12]
[13]
[14]
[15]
[16]
[17]

D. L. Donoho, IEEE Trans. Inf. Theory 52, 1289 (2006).
E. J. Candes, C. R. Acad. Sci. 346, 589 (2008).
E. J. Candes, Proc. Int. Congress of Mathematicians (2006).
M. F. Duarte, M. A. Davenport, D. Takhar, J. N. Laska, T. Sun, K. F.
Kelly, and R. G. Baraniuk, IEEE Signal Process. Mag. 25, 83 (2008).
F. Magalhaes, F. M. Araujo, M. V. Correia, M. Abolbashari, and
F. Farahi, Appl. Opt. 50, 405 (2011).
L. Martinez-Leon, P. Clemente, Y. Mori, V. Climent, J. Lancis, and
E. Tajahuerce, Opt. Express 25, 4975 (2017).
P. Clemente, V. Durán, E. Tajahuerce, P. Andrés, V. Climent, and
J. Lancis, Opt. Lett. 38, 2524 (2013).
B. Deepan, C. Quan, and C. J. Tay, Proc. SPIE 9234, 923419 (2014).
Y. Rivenson, A. Stern, and J. Rosen, Opt. Lett. 38, 2509 (2013).
A. Liutkus, D. Martina, S. Popoff, G. Chardon, O. Katz, G. Lerosey,
S. Gigan, L. Daudet, and I. Carron, Sci. Rep. 4, 5552 (2014).
C. G. Graff and E. Y. Sidky, Appl. Opt. 54, C23 (2015).
M. Süzen, A. Giannoula, and T. Durduran, Opt. Express 18, 23676
(2010).
R. Yao, Q. Pian, and X. Intes, Biomed. Opt. Express 6, 4887 (2015).
J. Li, J. S. Li, Y. Y. Pan, and R. Li, Sci. Rep. 5, 10374 (2015).
N. Rawat, B. Kim, I. Muniraj, G. Situ, and B.-G. Lee, Appl. Opt. 54,
1782 (2015).
Z. Wang and Z. Yu, Opt. Express 22, 25608 (2014).
G. A. Howland, J. Schneeloch, D. J. Lum, and J. C. Howell, Phys. Rev.
Lett. 112, 253602 (2014).

[18] M. Mirhosseini, O. S. Magaña-Loaiza, S. M. H. Rafsanjani, and R. W.
Boyd, Phys. Rev. Lett. 113, 090402 (2014).
[19] A. Kalev, R. L. Kosut, and I. H. Deutsch, npj Quantum Inf. 1, 15018
(2015).
[20] G. A. Howland, S. H. Knarr, J. Schneeloch, D. J. Lum, and J. C. Howell,
Phys. Rev. X 6, 021018 (2016).
[21] A. F. Abouraddy, T. M. Yarnall, and B. E. A. Saleh, Opt. Lett. 37, 2889
(2012).
[22] L. Martin, D. Mardani, H. E. Kondakci, W. D. Larson, S. Shabahang,
A. K. Jahromi, T. Malhotra, A. N. Vamivakas, G. K. Atia, and A. F.
Abouraddy, Sci. Rep. 7, 44995 (2017).
[23] N. Mohan, I. Stojanovic, W. C. Karl, B. E. A. Saleh, and M. C. Teich,
Proc. SPIE 7570, 75700L (2010).
[24] D. Mardani, A. F. Abouraddy, and G. K. Atia, Opt. Express 23, 28449
(2015).
[25] A. F. Abouraddy, T. M. Yarnall, and B. E. A. Saleh, Opt. Lett. 36, 4683
(2011).
[26] V. Namias, IMA J. Appl. Math. 25, 241 (1980).
[27] H. M. Ozaktas, Z. Zalevsky, and M. A. Kutay, The Fractional Fourier
Transform with Applications in Optics and Signal Processing (Wiley, 2001).
[28] V. Namias, IMA J. Appl. Math. 26, 187 (1980).
[29] E. J. Candes and Y. Plan, IEEE Trans. Inf. Theory 57, 7235 (2011).
[30] J. A. Rodrigo, T. Alieva, and M. L. Calvo, Opt. Express 17, 4976 (2009).
[31] K. H. Kagalwala, G. D. Giuseppe, A. F. Abouraddy, and B. E. A. Saleh,
Nat. Photon. 7, 72 (2013).

