Sampling Without Time:
Recovering Echoes of Light via Temporal Phase Retrieval

Ayush Bhandari†, Aurélien Bourquard ‡ and Ramesh Raskar†

† Media Laboratory and ‡ Research Laboratory of Electronics
Massachusetts Institute of Technology Cambridge, MA 02139–4307 USA.

<EMAIL_ADDRESS>$\bullet$<EMAIL_ADDRESS>$\bullet$<EMAIL_ADDRESS>
To Appear in the Proceedings of IEEE ICASSP, 2017111The $42^{\rm{nd}}$ IEEE
International Conference on Acoustics, Speech, and Signal Processing
(ICASSP)..
This work expands on the ideas discussed in [1] and [2].

Abstract

_This paper considers the problem of sampling and reconstruction of a
continuous-time sparse signal without assuming the knowledge of the sampling
instants or the sampling rate. This topic has its roots in the problem of
recovering multiple echoes of light from its low-pass filtered and auto-
correlated, time-domain measurements. Our work is closely related to the topic
of sparse phase retrieval and in this context, we discuss the advantage of
phase-free measurements. While this problem is ill-posed, cues based on
physical constraints allow for its appropriate regularization. We validate our
theory with experiments based on customized, optical time-of-flight imaging
sensors. What singles out our approach is that our sensing method allows for
temporal phase retrieval as opposed to the usual case of spatial phase
retrieval. Preliminary experiments and results demonstrate a compelling
capability of our phase-retrieval based imaging device._

###### Contents

  1. I Introduction
  2. II Image–Formation Model for ToF Sensors
  3. III Sampling Echoes of Light
    1. III-A Bandlimited Approximation of Probing Function
    2. III-B Sampling Theory Context
  4. IV Reconstruction via Phase Retrieval
  5. V Experimental Validation
    1. V-A Simulation
    2. V-B Practical Validation
  6. VI Conclusions

##  I Introduction

During a recent visit to meet a collaborator, the author (who happens to be an
avid photographer) saw a stark reflection of a local monument on the window
panes of the Harvard science center. This is shown in Fig. 1. This is a common
place phenomenon at a macro scale as well as a micro scale (for example
microscopy). At the heart of this problem is a fundamental limitation, that
is, all of the conventional imaging sensors are agnostic to the time
information. Alternatively stated, the image–formation process is insensitive
to the potentially distinct times that the photons can spend traveling between
their sources and the detector.

To elaborate, consider the _Gedankenexperiment_ version of Fig. 1 described in
Fig. 2 (a). Let us assume that a light source is co-located with the imaging
sensor (such as a digital camera). The reflection from the semi-reflective
sheet (such as a window pane) arrives at the sensor at time $t_{1}=2d_{1}/c$
while the mannequin is only observable on or after $t_{2}=2d_{2}/c$. Here, $c$
is the speed of light. While we typically interpret images as a two
dimensional, spatial representation of the scene, let us for now, consider the
photograph in Fig. 1 along the time-dimension. For the pixel
$\left(x_{0},y_{0}\right)$, this time-aware image is shown in Fig. 2 (b).
Mathematically, our time-aware image can be written as a $2$–sparse signal
(and as a $K$–sparse signal in general),

$m\left({\boldsymbol{r},t}\right)=\left({{\lambda_{1}}{\Gamma_{1}}}\right)\left(\boldsymbol{r}\right)\delta\left({t-{t_{1}{(\boldsymbol{r})}}}\right)+\left({{\lambda_{2}}{\Gamma_{2}}}\right)\left(\boldsymbol{r}\right)\delta\left({t-{t_{2}{(\boldsymbol{r})}}}\right)$
(1)

where $\left\\{{{\Gamma_{k}}\left(\boldsymbol{r}\right)}\right\\}$ are the
constituent image intensities,
$\left\\{{{\lambda_{k}}\left(\boldsymbol{r}\right)}\right\\}$ are the
reflection coefficients, $\boldsymbol{r}=\left(x,y\right)^{\top}$ is the 2D
spatial coordinate, and $\delta(\cdot)$ is the Dirac distribution.

Figure 1: Reflection from semi-reflective surfaces. The Memorial Church can be
seen imprinted on the glass facade of the Harvard University Science Center.

A conventional imaging sensor produces images by marginalizing time
information, resulting in the 2D photograph,

$\int_{0}^{\Delta\gg{t_{2}}}{m\left({\boldsymbol{r},t}\right)dt}=\left({{\lambda_{1}}{\Gamma_{1}}}\right)\left(\boldsymbol{r}\right)+\left({{\lambda_{2}}{\Gamma_{2}}}\right)\left(\boldsymbol{r}\right)\equiv
I\left(\boldsymbol{r}\right).$ (2)

Recovering $\left\\{{{\Gamma_{k}}}\right\\},k=1,2$ given $I$, more generally,
$K$ echoes of light,

$m\left({\boldsymbol{r},t}\right)=\sum\limits_{k=0}^{K-1}{\left({{\lambda_{k}}{\Gamma_{k}}}\right)\left(\boldsymbol{r}\right)\delta\left({t-{t_{k}{\color[rgb]{0,0,0}(\boldsymbol{r})}}}\right)}$
(3)

is an ill-posed problem. Each year, a number of papers [3, 4, 5, 6, 7, 8]
attempt to address this issue by using regularization and/or measurement-
acquisition diversity based on image statistics, polarization, shift, motion,
color, or scene features. Unlike previous works, here, we ask the question:
Can we directly estimate $\left\\{{{\Gamma_{k}}}\right\\}_{k=0}^{K-1}$ in (1)?
In practice, sampling (3) would require exorbitant sampling rates and this is
certainly not an option. Also, we are interested in
$\left\\{{{\Gamma_{k}}}\right\\}_{k=0}^{K-1}$ only, as opposed to
$\left\\{{{\Gamma_{k},t_{k}}}\right\\}_{k=0}^{K-1}$ [9] where $t_{k}$ is a
non-linear argument in (3). As a result, our goal is to recover the
intensities of a sparse signal. For this purpose, we explore the idea of
_sampling without time_ —a method for sampling a sparse signal which does not
assume the knowledge of sampling instants or the sampling rate. In this
context, our contributions are twofold:

  1. $1)$

For the general case of $K$-echoes of light, our work relies on estimating the
constituent images $\left\\{{{\Gamma_{k}}}\right\\}_{k=0}^{K-1}$ from the
filtered, auto-correlated, time-resolved measurements. This is the
distinguishing feature of our approach and is fundamentally different from
methods proposed in literature which are solely based on spatio-angular
information (cf. [3, 4, 5, 6, 7, 8] and references therein).

  2. $2)$

As will be apparent shortly, our work is intimately linked with the problem of
(sparse) phase retrieval (cf. [10, 11, 12, 13, 14] and references therein).
Our “sampling without time” architecture leads to an interesting measurement
system which is based on time-of-flight imaging sensors [1] and suggests that
$K^{2}-K+1$ measurements suffice for estimation of $K$ echoes of light. To the
best of our knowledge, neither such a measurement device nor bounds have been
studied in the context of image source separation [3, 4, 5, 6, 7].

Figure 2: Time-aware imaging. (a) Exemplary setting for imaging through semi-
reflective media. (b) Corresponding time-aware image at pixel
$\boldsymbol{r}_{0}=\left(x_{0},y_{0}\right)^{\top}$.

For the sake of simplicity, we will drop the dependence of $m$ and $\Gamma$ on
spatial coordinates $\boldsymbol{r}$. This is particularly well suited for our
case because we do not rely on any cross-spatial information or priors. Also,
scaling factors $\lambda_{k}$ are assumed to be absorbed in $\Gamma_{k}$.

Significance of Phase-free or Amplitude-only measurements: It is worth
emphasizing that resolving spikes from a superposition of low-pass filtered
echoes, is a problem that frequently occurs in many other disciplines. This is
a prototype model used in the study of multi-layered or multi-path models.
Some examples include seismic imaging [15], time-delay estimation [16],
channel estimation [17], optical tomography [18], ultrasound imaging [19],
computational imaging [1, 20, 21, 22, 23] and source localization [24]. Almost
all of these variations rely on amplitude and phase information or amplitude
and time-delay information. However, recording amplitude-only data can be
advantageous due to several reasons. Consider a common-place example based on
pulse-echo ranging. Let $p\left(t\right)=\sin\left({{\omega}t}\right)$ be the
emitted pulse. Then, the backscattered signal reads
$r\left(t\right)=\Gamma_{0}\sin\left({\omega t-{\theta_{0}}}\right)$. In this
setting, on-chip estimation of phase ($\theta_{0}$) or time delay
($t_{0}$)[22],

  1. $\bullet$

can either be computationally expensive or challenging since $t_{k}$ is a non-
linear parameter in (3), and hence, in $r\left(t\right)$.

  2. $\bullet$

requires more measurements. For instance, $2$ measurements suffice for
amplitude-only estimation ($\Gamma_{0}$) while amplitude-phase pair requires
$4$ measurements [2]. This aspect of phase estimation is an important
bottleneck as the the frame rate of an imaging device is inversely
proportional to number time-domain measurements.

  3. $\bullet$

is prone to errors. In many applications, multiple-frequency measurements,
$\omega=k{\omega_{0}}$, are acquired [22, 25, 2] assuming that phase and
frequency, ${\theta_{0}}=\omega{t_{0}}$ are linearly proportional. However,
this is not the case in practice. The usual practice is to oversample [26,
27].

In all such cases, one may adopt our methodology of working with intensity-
only measurements.

##  II Image–Formation Model for ToF Sensors

Optical ToF sensors are active devices that capture 3D scene information. We
adopt the generic image–formation model used in [1] which is common to all ToF
modalities such as lidar/radar/sonar, optical coherence tomography, terahertz,
ultrasound, and seismic imaging. In its full generality, and dropping
dependency on the spatial coordinate for convenience, one can first formalize
this ToF acquisition model as:

$p\to\boxed{h}\to r\to\boxed{\phi}\to
y\xrightarrow{{{\textsf{Sampling}}}}\boldsymbol{y}\mapsto\underbrace{m=\mathcal{C}\left[\boldsymbol{y}\right]}_{{\textsf{Intensity}}}$
(4)

where,

  1. $\rhd$

$p\left(t\right)>0$ is a $T$-periodic probing function which is used to
illuminate the scene. This may be a sinusoidal waveform, or even a periodized
spline, Gaussian, or Gabor pulse, for instance.

  2. $\rhd$

$h\left(t,\tau\right)$ is the scene response function (SRF). This may be a
filter, a shift-invariant function $h_{\sf
SI}\left(t,\tau\right)=h\left(t-\tau\right)$, or a partial differential
equation modeling some physical phenomenon [25].

  3. $\rhd$

$r\left(t\right)=\int{p\left(\tau\right)h\left({t,\tau}\right)d\tau}$ is the
reflected signal resulting from the interaction between the probing function
and the SRF.

  4. $\rhd$

$y\left(t\right)=\int{r\left(\tau\right)\phi\left({t,\tau}\right)d\tau}$ is
the continuous-time signal resulting from the interaction between the
reflected function and the instrument response function (IRF), or $\phi$,
which characterizes the electro-optical transfer function of the sensor.

  5. $\rhd$

$\boldsymbol{y}$ is a set of discrete measurements of the form
${\left.{y\left(t\right)}\right|_{t=n\Delta}}$.

  6. $\rhd$

${\color[rgb]{0,0,0}\mathcal{C}\left[\boldsymbol{y}\right]=\boldsymbol{y}*\overline{\boldsymbol{y}}}$
is the cyclic auto-correlation of $\boldsymbol{y}$, where $*$ and
$\overline{\left(\cdot\right)}$ denote the convolution and time reversal
operators, respectively.

The interplay between $p$, $h$, and $\phi$ results in variations on the theme
of ToF imaging [1]. In this work, we will focus on an optical ToF setting.
Accordingly:

  1. $1)$

The probing function corresponds to a time-localized pulse.

  2. $2)$

The SRF, accounting for the echoes of light, is a $K$-sparse filter,

${h_{{\textsf{SI}}}}\left(t,\tau\right)\equiv
h_{K}\left(t-\tau\right)=\sum\limits_{k=0}^{K-1}{{\Gamma_{k}}\delta\left({t-\tau-{t_{k}}}\right)}.$
(5)

  3. $3)$

The IRF is fixed by design as
${\phi_{{\sf{SI}}}}\left(t,\tau\right)=p\left(t+\tau\right)$. This implements
the so-called homodyne, lock-in sensor [1].

Due to this specific shift-invariant structure of the SRF and IRF, we have,
$\boldsymbol{y}=p*\overline{p}*{\overline{h}_{K}}\equiv\mathcal{C}\left[p\right]*\overline{h}_{K}$.
Finally, the measurements read,

$m=\mathcal{C}\left[\boldsymbol{y}\right]=\left(\mathcal{C}\left[\varphi\right]*\mathcal{C}\left[{h_{K}}\right]\right),\quad\varphi=\mathcal{C}\left[p\right].$
(6)

##  III Sampling Echoes of Light

###  III-A Bandlimited Approximation of Probing Function

Due to physical constraints inherent to all electro-optical systems, it is
reasonable to approximate the probing signal as a bandlimited function [25,
9]. We use $L$–term Fourier series with basis functions
${u_{n}}\left({\color[rgb]{0,0,0}\omega_{0}t}\right)\stackrel{{\scriptstyle\mathrm{def}}}{{=}}{e^{\jmath{\omega_{0}}nt}}$,
$\omega_{0}T=2\pi$, to approximate $p$ with,

$p\left(t\right)\approx\widetilde{p}\left(t\right)=\sum\limits_{\left|\ell\right|\leqslant
L}{{{\widehat{p}}_{\ell}}{u_{\ell}}\left({\color[rgb]{0,0,0}\omega_{0}t}\right)},$
(7)

where the $\widehat{p}_{\ell}$ are the Fourier coefficients. Note that there
is no need to compute $p\left(t\right)$ explicitly as we only require the
knowledge of $\mathcal{C}\left[p\right]$ (cf. (6)). In Fig. 3, we plot the
calibrated $\mathcal{C}\left[p\right]$ and its Fourier coefficients. In this
case, $T=232.40$ ns. We also plot its approximation,
${\widetilde{\varphi}}=\mathcal{C}\left[\widetilde{p}\right]$ with $L=25$
together with the measured $\mathcal{C}\left[p\right]$ and $\varphi$.

###  III-B Sampling Theory Context

The shift-invariant characterization of the ToF image–formation model allows
to re-interpret the sampled version of (6) as the filtering of a sparse signal
$\mathcal{C}\left[h_{K}\right]$ with a low-pass kernel
$\psi=\mathcal{C}\left[\varphi\right]$ (cf. [28, 9]).

$\displaystyle m\left(t\right)$
$\displaystyle\stackrel{{\scriptstyle(\ref{acm})}}{{=}}\mathcal{C}\left[h_{K}\right]*\mathcal{C}\left[\varphi\right]\equiv\mathcal{C}\left[h_{K}\right]*\psi.$
(8)

Note that the properties of the auto-correlation operation imply that the
sparsity of $\mathcal{C}\left[h_{K}\right]$ is $K^{2}-K+1$, unlike $h_{K}$
that is $K$-sparse and is completely specified by $(K^{2}-K)/2+1$ due to
symmetry. Based on the approximation (7) and the properties of convolution and
complex exponentials, and defining
$\widehat{\psi}_{\ell}=|\widehat{p}_{\ell}|^{4}$, we rewrite $m(t)$ as,

${\color[rgb]{0,0,0}m\left(t\right)}=\sum\limits_{|\ell|\leq
L}\widehat{\psi}_{\ell}u_{\ell}\left({\color[rgb]{0,0,0}\omega_{0}}t\right)\underbrace{\int{\mathcal{C}\left[h_{K}\right]\left(\tau\right)u_{\ell}^{*}{\left({\color[rgb]{0,0,0}\omega_{0}}\tau\right)}}d\tau}_{\sf{Fourier\
Integral}},$ (9)

Finally, the properties of the Fourier transform imply that sampling the above
signal $m(t)$ at time instants $t=n\Delta$, results in discrete measurements
of the form
$\boldsymbol{m}=\mathbf{U}\mathbf{D}_{\hat{\psi}}\hat{\boldsymbol{s}}$, which
corresponds to the available acquired data in our acquisition setting.
Combining all the above definitions, it follows that:

  1. —$\bullet$

${\boldsymbol{m}}\in{\mathbb{R}^{N}}$ is a vector of filtered measurements
(cf. (8)).

  2. —$\bullet$

${\mathbf{U}}\in{\mathbb{C}^{N\times\left({2{L}+1}\right)}}$ is a DFT matrix
with elements
${\left[u_{n}\left({\color[rgb]{0,0,0}\omega_{0}}\ell\Delta\right)\right]_{n,\ell}}$.

  3. —$\bullet$

${{\mathbf{D}}_{\hat{\psi}}}\in{\mathbb{C}^{\left(2L+1\right)\times\left(2L+1\right)}}$
is a diagonal matrix with diagonal elements $\widehat{\psi}_{\ell}$. These are
the Fourier-series coefficients of $\mathcal{C}\left[\varphi\right]$.

  4. —$\bullet$

$\widehat{\boldsymbol{s}}\in\mathbb{R}^{2L+1}$ is a phase-less vector
containing the Fourier coefficients of $\mathcal{C}\left[h_{K}\right]$, which
is obtained by sampling the Fourier transform
$\widehat{s}\left({{\omega}}\right)$ of $\mathcal{C}\left[h_{K}\right]$ at
instants $\omega=\ell\omega_{0}$. The signal
$\widehat{s}\left({{\omega}}\right)$ directly depends on the quantities
$|\Gamma_{k}|$ of interest to be retrieved and is expressed as

$\displaystyle\widehat{s}\left({{\omega}}\right)$
$\displaystyle={{\left|{\sum\limits_{k=0}^{K-1}{{\Gamma_{k}}{\color[rgb]{0,0,0}u_{t_{k}}\left(\omega\right)}}}\right|^{2}}}\equiv|\widehat{h}_{K}\left(\omega\right)|^{2}\hfill$
(10)
$\displaystyle=\underbrace{\sum\limits_{k=0}^{K-1}{|\Gamma_{k}|^{2}}}_{a_{0}}+2\sum\limits_{k=0}^{K-1}{\sum\limits_{m=k+1}^{K-1}{\underbrace{\left|{{\Gamma_{k}}}\right|\left|{{\Gamma_{m}}}\right|}_{a_{k,m}}\cos\left({{\omega}{t_{k,m}}+\angle{\Gamma_{k,m}}}\right)}}$

where ${t_{k,m}}={t_{k}}-{t_{m}}$ and
$\angle{\Gamma_{k,m}}=\angle{\Gamma_{k}}-\angle{\Gamma_{m}}$.

It is instructive to note that the relevant unknowns
$\left\\{{{\Gamma_{k}}}\right\\}_{k=0}^{K-1}$ can be estimated from
$\widehat{\boldsymbol{s}}\in\mathbb{R}^{2L+1}$ in (10) which in turn depend
only on $L$ as opposed to sampling rate and sampling instants. This aptly
justifies our philosophy of _sampling without time_.

Figure 3: Bandlimited approximation
$\widetilde{\varphi}\approx\mathcal{C}\left[p\right]$ with $L=25$ (red curves)
and its cyclic auto-correlation $\mathcal{C}\left[\widetilde{\varphi}\right]$
(blue curve). Shown are the spatial-domain signals (top) and the corresponding
Fourier-domain magnitudes in log scale (bottom).

##  IV Reconstruction via Phase Retrieval

Given data model
$\boldsymbol{m}=\mathbf{U}\mathbf{D}_{\hat{\psi}}\hat{\boldsymbol{s}}$ (9), we
aim to recover the image parameters $|\Gamma_{k}|$ at each sensor pixel. For
this purpose, we first estimate samples of $\hat{s}(\omega)$ as
$\widehat{\boldsymbol{s}}=\mathbf{D}_{\psi}^{-1}\mathbf{U}^{+}\boldsymbol{m}$,
where $\mathbf{U}^{+}$ is the matrix pseudo-inverse of $\mathbf{U}$. This is
akin to performing a weighted deconvolution, knowing that $\mathbf{U}$ is a
DFT matrix. Next, we solve the problem of estimating $|\Gamma_{k}|$ in two
steps (1) First we estimate $a_{k,m}$, and then, (2) based on the estimated
values, we resolve ambiguities due to $|\cdot|$.

Parameter Identification via Spectral Estimation: Based on the coefficients
$\widehat{\boldsymbol{s}}$, one can then retrieve the amplitude and frequency
parameters that are associated with the oscillatory terms as well as the
constant value in (10). The oscillatory-term and constant-term parameters
correspond to $\\{a_{k,m},t_{k,m}\\}$ and $a_{0}$, respectively. All parameter
values are retrievable from $\widehat{\boldsymbol{s}}$ through spectral
estimation [29]; details are provided in [2] for the interested reader.

Note that, given the form of (10) and our acquisition model, the sparsity
level of the sequence $\widehat{\boldsymbol{s}}$—corresponding to the total
amount of oscillatory and constant terms—is $(K^{2}-K)/2+1$. The magnitude
values $a_{k,m}$ and $a_{0}$ can thus be retrieved if at least
$L\geq(K^{2}-K)+1$ autocorrelation measurements are performed [30, 2], which
is the case in the experiments described in Section V.

Resolving Ambiguities: Based on the aforementioned retrieved parameters, one
wishes to then deduce $\\{|\Gamma_{k}|\\}_{k=0}^{K-1}$. The estimated cross
terms $\\{a_{k,m},t_{k,m}\\}$ allow to retrieve the values of $|\Gamma_{k}|$
through simple pointwise operations. Here we will focus on the case, $K=2$.
Due to space limitations, we refer the interested readers to our companion
paper [2] for details on the general case ($K>2$). The case when $K=2$,

$\displaystyle\widehat{s}\left(\omega\right)\stackrel{{\scriptstyle(\ref{mag})}}{{=}}$
$\displaystyle{\left|{{\Gamma_{0}}u_{{t_{0}}}^{*}\left(\omega\right)+{\Gamma_{1}}u_{{t_{1}}}^{*}\left(\omega\right)}\right|^{2}}+\underbrace{{\varepsilon_{\omega}}{{\left({{\Gamma_{k}},{t_{k}}}\right)}_{k>2}}}_{\approx
0}$ $\displaystyle=$
$\displaystyle\underbrace{{{\left|{{\Gamma_{0}}}\right|}^{2}}+{{\left|{{\Gamma_{1}}}\right|}^{2}}}_{{a_{0}}}+\underbrace{2\left|{{\Gamma_{0}}}\right|\left|{{\Gamma_{1}}}\right|}_{{a_{01}}}\cos\left({\omega{t_{0,1}}+\angle{\Gamma_{0,1}}}\right),$

may be a result of two distinct echoes (Fig. 2) or due to approximation of
higher order echoes by two dominant echoes (due to inverse-square law). In
this case, we effectively estimate two terms: $a_{0}$ and $a_{01}$ (10). The
set of retrieved magnitude values then amount to solving[2]:

$\left|{\left|{{\Gamma_{0}}}\right|\pm\left|{{\Gamma_{1}}}\right|}\right|=\sqrt{{a_{0}}\pm{a_{01}}},\qquad{a_{0}},{a_{01}}>0.$
(11)

Thanks to the _isoperimetric property_ for rectangles: $a_{0}\geq a_{01}$, the
r.h.s above is always positive unless there is an estimation error in which
case, an exchange of variables leads to the unique estimates,

${\left\\{{|{{\widetilde{\Gamma}}_{k}}|}\right\\}_{k=1,2}}=\frac{{\sqrt{{{\widetilde{a}}_{0}}+{{\widetilde{a}}_{01}}}\pm\sqrt{{{\widetilde{a}}_{0}}-{{\widetilde{a}}_{01}}}}}{2}.$

##  V Experimental Validation

###  V-A Simulation

Noting that measurements
$\boldsymbol{m}=\mathbf{U}\mathbf{D}_{\hat{\psi}}\hat{\boldsymbol{s}}$ and
$\hat{\boldsymbol{s}}$ are linked by an invertible, linear system of
equations, knowing $\boldsymbol{m}$ amounts to knowing $\hat{\boldsymbol{s}}$.
We have presented detailed comparison of simulation results in [2] for the
case where one directly measures $\hat{\boldsymbol{s}}$. In this paper, we
will focus on practical setting where the starting point of our algorithm is
(9).

Figure 4: Experimental results on acquisition and reconstruction of echoes of
light in a scene. (a) Single pixel, time-stamped samples
$y=\mathcal{C}\left[p\right]*\overline{h}_{2}$ serve as our ground truth
because the phase information is known. We also plot the estimated, 2-sparse
SRF (5) in red ink. (b) Fourier domain frequency samples of the SRF,
$\widehat{h}_{2}\left(\ell\omega_{0}\right)$ with $L=[-20,20]$ and
$\omega_{0}=2\pi/\Delta,\Delta=70$ ps. We plot the measured data in green and
the fitted data in black. The frequencies are estimated using Cadzow’s method.
(c) Same pixel, time-stampless, low-pass filtered, auto-correlated samples (8)
together with the estimated, auto-correlated SRF
$\mathcal{C}\left[h_{2}\right]$ in red ink. (d) Fourier domain samples
$\widehat{s}\left(\ell\omega_{0}\right)$ (10) and its fitted version (black
ink). (e) Ground truth images for the experiment. (f) Estimated images using
temporal phase retrieval.

###  V-B Practical Validation

Our experimental setup mimics the setting of Fig. 2. A placard that reads
“Time–of–Flight” is covered by a semi-transparent sheet, hence $K=2$. The
sampling rate is $\Delta\approx 70\times 10^{-12}$ seconds using a our custom
designed ToF imaging sensor [1]. Overall, the goal of this experiment is to
recover the magnitudes ${\left\\{{|{\Gamma_{k}}|}\right\\}_{k=0}^{k=1}}$ given
auto-correlated measurements
$m\left(n\Delta\right)\stackrel{{\scriptstyle(\ref{acm})}}{{=}}\mathcal{C}\left[y\right]\left(n\Delta\right),\mathbf{m}\in\mathbb{R}^{2795}$.
To be able to compare with a “ground truth”, we acquire time-domain
measurements $y\left(n\Delta\right)$ before autocorrelation, whose Fourier-
domain phases are intact. In Fig. 4(a), we plot the non-autocorrelated
measurements $\mathbf{y}$ while phase-less measurements $\mathbf{m}$ are shown
in Fig. 4(c) from which we note that the samples are symmetric in time domain
due to $m=\mathcal{C}\left[y\right]$ (cf. 6). In the first case (cf. Fig.
4(a)) where $y=\mathcal{C}\left[p\right]*\overline{h}_{K}$, we have [1],

$\widehat{y}\left(\omega\right)={\left|{\widehat{p}\left(\omega\right)}\right|^{2}}h_{K}^{*}\left(\omega\right),\quad{h_{K}}\left(\omega\right)=\sum\limits_{k=0}^{K-1}{{\Gamma_{k}}{u_{{t_{k}}}}\left(\omega\right)}.$

Similar to $\mathbf{m}$ in (9), in this case, the measurements read
$\mathbf{y}=\mathbf{U}\mathbf{D}_{\widehat{\varphi}}\widehat{\mathbf{h}}$ and
we can estimate the complex-valued vector $\widehat{\mathbf{h}}$. The phase
information in $\widehat{\mathbf{h}}$ allows for the precise computation of
$\\{\Gamma_{k},t_{k}\\}_{k=0}^{K-1}$ [1, 9]. These “intermediate” measurements
serve as our ground truth. The spikes corresponding to the SRF (5) are also
marked in Fig. 4(a). Fourier-domain measurements $\widehat{h}_{2}^{*}$ linked
with $y\left(n\Delta\right)$ are shown in Fig. 4(b). Accordingly,
${\widehat{h}^{*}_{2}}\left(\omega\right),\
\omega=\ell\omega_{0},l=\\{-L,\ldots,L\\}$ where $\omega_{0}=2\pi/\Delta$ and
$L=20$. We estimate the frequency parameters
$\\{\Gamma_{k},t_{k}\\}_{k=0}^{K-1}$ using Cadzow’s method [31]. The resulting
fit is plotted in Fig. 4(b).

In parallel to Figs. 4(a),(b), normalized, autocorrelated data
$\mathbf{m}=\mathcal{C}\left[\mathbf{y}\right]$ are marked in Fig. 4(c). The
signal $\mathbf{y}$ is also shown as reference in green dashed line. We also
plot the estimated $\mathcal{C}\left[h_{2}\right]$ (autocorrelated SRF). In
Fig. 4(d), we plot measured and deconvolved vector
$\widehat{\boldsymbol{s}}=\mathbf{D}_{\psi}^{-1}\mathbf{U}^{+}\boldsymbol{m}$
from which we estimate $a_{0},a_{k,m}$ (10). The result of fitting using
Cadzow’s method with $41$ samples is shown in Fig. 4(d).

The reconstructed images, $|\Gamma_{0}|$ and $|\Gamma_{1}|$, due to amplitude-
phase measurements (our ground truth, Fig. 4(a)) are shown in Fig. 4(e).
Alternatively, the reconstructed images with auto-correlated/intensity-only
information are shown in 4(f). One can observe the great similarity between
the images obtained in both cases, where only a few outliers appear in the
phase-less setting. The PSNR values between the maps reconstructed with and
without the phase information are of $30.25$ dB for $|\Gamma_{1}|$ and $48.88$
dB for $|\Gamma_{0}|$. These numerical results indicate that, overall, the
phase loss that occurs in our autocorrelated measurements still allows for
accurate reconstruction of the field of view of the scene.

Finally, in order to determine the consistency of our reconstruction approach,
the PSNR between the original available measurements and their re-synthesized
versions—as obtained when reintroducing our reconstructed profiles
$|\Gamma_{k}|$ into our forward model—are also provided. As shown in Figs.
4(a) and (c), the PSNR values in the oracle and phase-less settings correspond
to $43.68$ dB and $42.41$ dB, respectively, which confirms that our
reconstruction approach accurately takes the parameters and structure of the
acquisition model into account.

##  VI Conclusions

In this paper, we have introduced a method to satisfactorily recover the
intensities of superimposed echoes of light, using a customized ToF imaging
sensor for acquisition and temporal phase retrieval for reconstruction. Up to
our knowledge, this is the first method that performs time-stampless sampling
of a sparse ToF signal, in the sense that we only measure the amplitudes
sampled at uniform instants without caring about the particular sampling times
or sampling rates. This innovation can potentially lead to alternative
hardware designs and mathematical simplicity as phases need not be estimated
in hardware.

## References

  * [1] A. Bhandari and R. Raskar, “Signal processing for time-of-flight imaging sensors,” _IEEE Signal Processing Magazine_ , vol. 33, no. 4, pp. 45–58, Sep. 2016.
  * [2] A. Bhandari, A. Bourquard, S. Izadi, and R. Raskar, “Time-resolved image demixing,” in _Proc. of IEEE Intl. Conf. on Acoustics, Speech and Sig. Proc. (ICASSP)_ , Mar. 2016, pp. 4483–4487.
  * [3] N. Ohnishi, K. Kumaki, T. Yamamura, and T. Tanaka, “Separating real and virtual objects from their overlapping images,” in _Lecture Notes in Computer Science_. Springer Science $+$ Business Media, 1996, pp. 636–646.
  * [4] A. Levin, A. Zomet, and Y. Weiss, “Separating reflections from a single image using local features,” in _Proc. of IEEE Comp. Vis. and Patt. Recognition_ , vol. 1, 2004, pp. I–306.
  * [5] A. M. Bronstein, M. M. Bronstein, M. Zibulevsky, and Y. Y. Zeevi, “Sparse ICA for blind separation of transmitted and reflected images,” _Int. J. Imaging Syst. Technol._ , vol. 15, no. 1, pp. 84–91, 2005.
  * [6] N. Kong, Y.-W. Tai, and S. Y. Shin, “High-quality reflection separation using polarized images,” vol. 20, no. 12, pp. 3393–3405, 2011.
  * [7] Y. Li and M. S. Brown, “Single image layer separation using relative smoothness,” in _Proc. of IEEE Comp. Vis. and Patt. Recognition_ , 2014, pp. 2752–2759.
  * [8] P. Chandramouli, M. Noroozi, and P. Favaro, “ConvNet-based depth estimation, reflection separation and deblurring of plenoptic images,” University of Bern, submitted, 2016.
  * [9] A. Bhandari, A. M. Wallace, and R. Raskar, “Super-resolved time-of-flight sensing via FRI sampling theory,” in _Proc. of IEEE Intl. Conf. on Acoustics, Speech and Sig. Proc. (ICASSP)_ , Mar. 2016, pp. 4009–4013.
  * [10] Y. C. Eldar, N. Hammen, and D. G. Mixon, “Recent advances in phase retrieval [Lecture Notes],” vol. 33, no. 5, pp. 158–162, Sep. 2016.
  * [11] H. Qiao and P. Pal, “Sparse phase retrieval with near minimal measurements: A structured sampling based approach,” in _Proc. of IEEE Intl. Conf. on Acoustics, Speech and Sig. Proc. (ICASSP)_ , Mar. 2016, pp. 4722–4726.
  * [12] B. Rajaei, E. W. Tramel, S. Gigan, F. Krzakala, and L. Daudet, “Intensity-only optical compressive imaging using a multiply scattering material and a double phase retrieval approach,” in _Proc. of IEEE Intl. Conf. on Acoustics, Speech and Sig. Proc. (ICASSP)_ , Mar. 2016.
  * [13] K. Huang, Y. C. Eldar, and N. D. Sidiropoulos, “Phase retrieval from 1D Fourier measurements: Convexity, uniqueness, and algorithms,” vol. 64, no. 23, pp. 6105–6117, Dec. 2016.
  * [14] Y. M. Lu and M. Vetterli, “Sparse spectral factorization: Unicity and reconstruction algorithms,” in _Proc. of IEEE Intl. Conf. on Acoustics, Speech and Sig. Proc. (ICASSP)_ , May 2011.
  * [15] S. Levy and P. K. Fullagar, “Reconstruction of a sparse spike train from a portion of its spectrum and application to high-resolution deconvolution,” _Geophysics_ , vol. 46, no. 9, pp. 1235–1243, Sep. 1981.
  * [16] J. J. Fuchs, “Multipath time-delay detection and estimation,” vol. 47, no. 1, pp. 237–243, Jan. 1999.
  * [17] Y. Barbotin, A. Hormati, S. Rangan, and M. Vetterli, “Estimation of sparse MIMO channels with common support,” vol. 60, no. 12, pp. 3705–3716, Dec. 2012\.
  * [18] C. S. Seelamantula and S. Mulleti, “Super-resolution reconstruction in frequency-domain optical-coherence tomography using the finite-rate-of-innovation principle,” vol. 62, no. 19, pp. 5020–5029, Oct. 2014\.
  * [19] P. Boufounos, “Compressive sensing for over-the-air ultrasound,” in _Proc. of IEEE Intl. Conf. on Acoustics, Speech and Sig. Proc. (ICASSP)_ , May 2011.
  * [20] A. Velten, R. Raskar, D. Wu, B. Masia, A. Jarabo, C. Barsi, C. Joshi, E. Lawson, M. Bawendi, and D. Gutierrez, “Imaging the propagation of light through scenes at picosecond resolution,” _Commun. ACM_ , vol. 59, no. 9, pp. 79–86, Aug. 2016.
  * [21] D. Shin, F. Xu, F. N. C. Wong, J. H. Shapiro, and V. K. Goyal, “Computational multi-depth single-photon imaging,” _Optics Express_ , vol. 24, no. 3, p. 1873, Jan. 2016.
  * [22] A. Bhandari, A. Kadambi, R. Whyte, C. Barsi, M. Feigin, A. Dorrington, and R. Raskar, “Resolving multipath interference in time-of-flight imaging via modulation frequency diversity and sparse regularization,” _Optics Letters_ , vol. 39, no. 6, pp. 1705–1708, Mar. 2014.
  * [23] A. Kadambi, R. Whyte, A. Bhandari, L. Streeter, C. Barsi, A. Dorrington, and R. Raskar, “Coded time of flight cameras: sparse deconvolution to address multipath interference and recover time profiles,” _ACM Trans. Graphics_ , vol. 32, no. 6, pp. 1–10, Nov. 2013.
  * [24] J. J. Fuchs and H. Chuberre, “A deconvolution approach to source localization,” vol. 42, no. 6, pp. 1462–1470, Jun. 1994.
  * [25] A. Bhandari, C. Barsi, and R. Raskar, “Blind and reference-free fluorescence lifetime estimation via consumer time-of-flight sensors,” _Optica_ , vol. 2, no. 11, pp. 965–973, Nov. 2015.
  * [26] K. Hibino, B. F. Oreb, D. I. Farrant, and K. G. Larkin, “Phase-shifting algorithms for nonlinear and spatially nonuniform phase shifts,” _Jour.l of the Opt. Soc. of America A_ , vol. 14, no. 4, p. 918, Apr. 1997\.
  * [27] P. Hariharan, B. F. Oreb, and T. Eiju, “Digital phase-shifting interferometry: a simple error-compensating phase calculation algorithm,” _Appl. Opt._ , vol. 26, no. 13, p. 2504, Jul. 1987.
  * [28] H. Pan, T. Blu, and M. Vetterli, “Towards generalized FRI sampling with an application to source resolution in radioastronomy,” vol. 65, no. 4, pp. 821–835, Feb. 2017.
  * [29] P. Stoica and R. L. Moses, _Introduction to spectral analysis_. Prentice Hall, Upper Saddle River, 1997, vol. 1.
  * [30] D. Potts and M. Tasche, “Parameter estimation for exponential sums by approximate prony method,” _Signal Processing_ , vol. 90, no. 5, pp. 1631–1642, May 2010.
  * [31] J. A. Cadzow, “Signal enhancement–A composite property mapping algorithm,” vol. 36, no. 1, pp. 49–62, 1988.