Recovery of Compactly Supported Functions from Spectrogram Measurements via Lifting

TL;DR

This paper introduces a new phase retrieval technique inspired by ptychography, enabling the approximate reconstruction of compactly supported functions from spectrogram data using a lifted formulation and numerical validation.

Contribution

It presents a novel phase retrieval method based on lifting for spectrogram measurements, extending phaseLift ideas to infinite-dimensional functions with practical truncation.

Findings

Numerical experiments show promising reconstruction results.

The method effectively handles infinite-dimensional spectrogram data.

The approach offers a new avenue for phase retrieval in imaging applications.

Abstract

A novel phase retrieval method, motivated by ptychographic imaging, is proposed for the approximate recovery of a compactly supported specimen function $f:\mathbb{R}\rightarrow\mathbb{C}$ from its continuous short time Fourier transform (STFT) spectrogram measurements. The method, partially inspired by the well known PhaseLift algorithm, is based on a lifted formulation of the infinite dimensional problem which is then later truncated for the sake of computation. Numerical experiments demonstrate the promise of the proposed approach.

Equations

$$b_{k,j}:={\left|{\int }_{-\infty }^{\infty }f\left(t\right)g\left(t-l_k\right)e^{-2\pi i{\omega }_{j}t}dt\right|}^2$$

$$\vec{b}=\left(\begin{array}[]{c}b_{1,1},\dots,b_{1,N},b_{2,1},\dots,b_{K,N}\end{array}\right)^{T}\in[0,\infty)^{NK},$$

$$\left|\mathcal{F}\left[f\cdot S_{l}g\right]\left(\omega \right)\right|=\frac{1}{2}\left|\sum _{m\in \mathbb{Z}}{e}^{-\pi ilm}\hat{f}\left(\frac{m}{2}\right)\hat{g}\left(\frac{m}{2}-\omega \right)\right|$$

$$\mathcal{F}\left[f\cdot S_{l}g\right]\left(\omega \right)={\int }_{-\infty }^{\infty }f\left(t\right)g\left(t-l\right)e^{-2\pi i\omega t}dt.$$

$${\left|\mathcal{F}\left[f\cdot S_{l}g\right]\left(\omega \right)\right|}^2={\left|{\int }_{-\infty }^{\infty }f\left(t\right)g\left(t-l\right)e^{-2\pi i\omega t}dt\right|}^2={\left|⟨f,h⟩\right|}^2$$

$$\hat{h}\left(k\right)$$

$${\left|⟨f,h⟩\right|}^2$$

$$\overline{\mathcal{F}\left[\overline{g\left(\cdot \right)}\right]\left(\eta \right)}=\hat{g}\left(-\eta \right).$$

$${\left|{\int }_{-\infty }^{\infty }\hat{f}\left(\omega -\eta \right)\hat{g}\left(-\eta \right)e^{-2\pi il\eta }d\eta \right|}^2$$

$$\mathcal{F}\left[\hat{g}\left(\cdot \right)e^{-2\pi il\left(\cdot \right)}\right]\left(\xi \right)=\hat{\hat{g}}\left(\xi +l\right)=g\left(-l-\xi \right),$$

$$\mathcal{F}\left[\frac{\sin \pi \left(m+2x\right)}{\pi \left(m+2x\right)}\right]\left(\xi \right)=\frac{e^{\pi im\xi }}{2}{\chi }_{\left(-1,1\right)}\left(\xi \right).$$

$${\left|\sum _{m\in \mathbb{Z}}\hat{f}\left(\frac{m}{2}\right){\mathcal{F}}^{-1}\left[g\left(-l-\left(\cdot \right)\right)\frac{e^{\pi im\left(\cdot \right)}}{2}{\chi }_{\left(-1,1\right)}\left(\cdot \right)\right]\left(-\omega \right)\right|}^2$$

$$\frac{1}{4}{\left|\sum _{m\in \mathbb{Z}}\hat{f}\left(\frac{m}{2}\right)e^{-\pi il\left(m-2\omega \right)}{\int }_{-a}^{a}g\left(u\right)e^{-2\pi iu\left(\frac{m}{2}-\omega \right)}du\right|}^2$$

$${\left|\mathcal{F}\left[f\cdot S_{l}g\right]\left(\omega \right)\right|}^2=\frac{1}{4}\sum _{k\in \mathbb{Z}}\sum _{j\in \mathbb{Z}}{A}_k\overline{A_j}$$

$$S_{\omega }=\left\{\left(j,k\right)\in \mathbb{Z}\times \mathbb{Z}|\left|k-2\omega \right|\le 2\delta \text{ and }\left|j-2\omega \right|\le 2\delta \right\},$$

$${\left|\mathcal{F}\left[f\cdot S_{l}g\right]\left(\omega \right)\right|}^2=\frac{1}{4}\sum _{(j,k)\in S_{\omega }}{A}_k\overline{A_j}+error.$$

$$\sum _{|j-2\omega |\le 2\delta }{e}^{\pi ilj}\overline{\hat{f}\left(\frac{j}{2}\right)}\overline{\hat{g}\left(\frac{j}{2}-\omega \right)}=e^{2\pi il\omega }{\vec{X}}_l^{\ast }\overrightarrow{Y_{\omega }}$$

$$\vec{X}_{l}=\left(\!\!\!\begin{array}[]{c}e^{\pi il\left(2\delta\right)}\hat{g}\left(-\delta\right)\\ e^{\pi il\left(2\delta-1\right)}\hat{g}\left(\frac{1}{2}-\delta\right)\\ \vdots\\ e^{\pi il\cdot 0}\hat{g}\left(0\right)\\ \vdots\\ e^{\pi il\left(1-2\delta\right)}\hat{g}\left(\delta-\frac{1}{2}\right)\\ e^{\pi il\left(-2\delta\right)}\hat{g}\left(\delta\right)\end{array}\!\!\!\right),\,\vec{Y}_{\omega}=\left(\!\!\!\begin{array}[]{c}\overline{\hat{f}\left(\omega-\delta\right)}\\ \overline{\hat{f}\left(\omega-\delta+\frac{1}{2}\right)}\\ \vdots\\ \overline{\hat{f}\left(\omega\right)}\\ \vdots\\ \overline{\hat{f}\left(\omega+\delta-\frac{1}{2}\right)}\\ \overline{\hat{f}\left(\omega+\delta\right)}\end{array}\!\!\!\right).$$

$${\left|\mathcal{F}\left[f\cdot S_{l}g\right]\left(\omega \right)\right|}^2$$

$$\!\left[\!\!\!\!\begin{array}[]{ccccc}\left|\hat{f}\!\left(\omega-\delta\right)\right|^{2}\!\!\!\!&\!\!\cdots\!\!\!&\overline{\hat{f}\!\left(\omega-\delta\right)}\hat{f}\!\left(\omega\right)\!\!\!&\!\!\cdots\!\!\!&\overline{\hat{f}\!\left(\omega-\delta\right)}\hat{f}\!\left(\omega+\delta\right)\!\!\\ \vdots\!\!\!&\!\!\ddots\!\!\!&\vdots\!\!\!&\!\!\vdots\!\!\!&\vdots\!\!\\ \overline{\hat{f}\!\left(\omega\right)}\hat{f}\!\left(\omega-\delta\right)\!\!\!\!&\!\!\cdots\!\!\!&\left|\hat{f}\!\left(\omega\right)\right|^{2}\!\!\!&\!\!\cdots\!\!\!&\overline{\hat{f}\!\left(\omega\right)}\hat{f}\!\left(\omega+\delta\right)\!\!\\ \vdots\!\!\!&\!\!\vdots\!\!\!&\vdots\!\!\!&\!\!\ddots\!\!\!&\vdots\!\!\\ \overline{\hat{f}\!\left(\omega+\delta\right)}\hat{f}\!\left(\omega-\delta\right)\!\!\!&\!\!\cdots\!\!\!&\overline{\hat{f}\!\left(\omega+\delta\right)}\hat{f}\!\left(\omega\right)\!\!\!&\!\!\cdots\!\!\!&\left|\hat{f}\!\left(\omega+\delta\right)\right|^{2}\!\!\end{array}\!\!\right]\!\!\!.$$

$$\vec{X}_{l}=\left(\begin{array}[]{ccccc}m_{-\delta}^{l},&m_{-\delta+\frac{1}{2}}^{l},&\dots,&m_{\delta-\frac{1}{2}}^{l},&m_{\delta}^{l}\end{array}\right)^{T}$$

$$\left[\begin{array}[]{cccccccc}m_{0}^{l}&m_{\frac{1}{2}}^{l}&\cdots&m_{\delta}^{l}&0&0&\cdots&0\\ m_{-\frac{1}{2}}^{l}&m_{0}^{l}&\cdots&m_{\delta-\frac{1}{2}}^{l}&m_{\delta}^{l}&0&\cdots&0\\ \vdots&\vdots&\vdots&\vdots&\vdots&\vdots&\vdots&\vdots\\ 0&0&\cdots&0&m_{-\delta}^{l}&m_{-\delta+\frac{1}{2}}^{l}&\cdots&m_{\frac{1}{2}}^{l}\\ 0&0&\cdots&0&0&m_{-\delta}^{l}&\cdots&m_{0}^{l}\end{array}\right]$$

$$\mathbf{G}=\left(\begin{array}[]{c}\mathbf{X}_{l_{1}}\\ \mathbf{X}_{l_{2}}\\ \vdots\\ \mathbf{X}_{l_{K}}\end{array}\right)$$

$${\mathbf{F}}_{i,j}=\{\begin{array}{ll}\overline{\hat{f}\left(\frac{i-2n-1}{2}\right)}\hat{f}\left(\frac{j-2n-1}{2}\right),& \text{if }\left|i-j\right|\le 2\delta ,\\ 0,& \text{otherwise,}\end{array}$$

$$\vec{b}\approx \text{mbox}diag(\mathbf{GF}{\mathbf{G}}^{\ast }),$$

$$\vec{b}\approx \mathbf{M}\vec{F},$$

Full text

I Introduction

The problem of signal recovery (up to a global phase) from phaseless STFT measurements appears in many audio engineering and imaging applications. Our principal motivation here, however, is ptychographic imaging (see, e.g., [14, 11]) in the 1-D setting where a compactly supported specimen, $f:\mathbb{R}\rightarrow\mathbb{C}$, is scanned by a focused illuminating beam $g:\mathbb{R}\rightarrow\mathbb{C}$ which translates across the specimen in fixed overlapping shifts $l_{1},\dots,l_{K}\in\mathbb{R}$. At each such shift of the beam (or, equivalently, the specimen) a phaseless diffraction image is then sampled in bulk by a detector. Due to the underlying physics the collected measurements are then modeled as sampled STFT magnitude measurements of $f$ of the form

[TABLE]

for a finite set of $KN$ shift and frequency pairs $(l_{k},\omega_{j})\in\{l_{1},\dots,l_{K}\}\times\{\omega_{1},\dots,\omega_{N}\}$. Our objective is to approximate $f$ (up to a global phase) using these $b_{k,j}$ measurements.

There has been a good deal of work on signal recovery from phaseless STFT measurements in the last couple of years in the discrete setting, where $f$ and $g$ are modeled as vectors ab initio, and then recovered from discrete STFT magnitude measurements. In this setting many related recovery techniques have been considered including iterative methods along the lines of Griffin and Lim [12, 18] and alternating projections [11], graph theoretic methods for Gabor frames based on polarization [15, 13], and semidefinite relaxation-based methods [8], among others [5, 2, 7, 6].

Herein we will instead consider the approximate recovery of $f$ (as a compactly supported function) from samples of its continuous STFT magnitude measurements $b_{k,j}$ as per (I.1). Besides perhaps better matching the continuous models considered in some applications such as ptychography, and allowing one to more naturally consider approaches that utilize, e.g., irregular sampling, we also take recent work on phase retrieval in infinite dimensional Hilbert spaces [19, 3, 1] as motivation for exploring numerical methods to solve this problem.

In particular, the recent work of Daubechies and her collaborators implies that the stability of such continuous phase retrieval problems is generally less well behaved than their discrete counterparts [3, 1]. Specifically, [1] characterizes a class of functions for which infinite dimensional phase retrieval (up to a single global phase) from Gabor measurements is unstable, and then proposes the reconstruction of these worst-case functions up to several local phase multiples as a stable alternative. We take this initial work on stable infinite dimensional phase retrieval from Gabor measurements as a further motivation to explore new fast numerical techniques for the robust recovery of compactly supported functions from their continuous spectrogram measurements.

I-A The Problem Statement and Specifications

Given a vector of stacked spectrogram samples from (I.1),

[TABLE]

our goal is to approximately recover a piecewise smooth and compactly supported function $f:\mathbb{R}\rightarrow\mathbb{C}$. Of course $f$ can only be recovered up to certain ambiguities (such as up to a global phase, etc.) which depend not only on $f$, but also the window function $g$ (see, e.g., [1]). Without loss of generality, we will assume that the support of $f$ is contained in $[-1,1]$. Given our motivation from ptychographic imaging we will, herein at least, primarily consider the unshifted beam function $g$ to also be (approximately) compactly supported within a smaller subset $[-a,a]\subset[-1,1]$. Furthermore, we will also assume that $g$ is smooth enough that its Fourier transform decays relatively rapidly in frequency space compared to $\hat{f}$. Examples of such $g$ include both suitably scaled Gaussians, as well as compactly supported $C^{\infty}$ bump functions [9].

I-B The Proposed Numerical Approach

The proposed method aims to recover samples from the Fourier transform of $f$ at frequencies in $\Omega=\{\omega_{1},\dots,\omega_{N}\}$, giving $\vec{f}\in\mathbb{C}^{N}$ with $f_{j}=\widehat{f}(\omega_{j})$, from which $\widehat{f}$ can then be approximately recovered via standard sampling theorems (see, e.g., [17]). The inverse Fourier transform of this approximation of $\widehat{f}$ then provides our approximation of $f$.

Recovery of the samples from $\widehat{f}$, $\vec{f}\in\mathbb{C}^{N}$, is performed in two steps using techniques from [7, 6] adapted to this continuous setting: first, a truncated lifted linear system is inverted in order to learn a portion of the rank-one matrix $\vec{f}\vec{f}^{*}$ from a finite set of STFT spectrogram samples, then, an eigenvector based angular synchronization method is used in order to recover $\vec{f}$ from the portion of $\vec{f}\vec{f}^{*}$ computed in the first step. Note that this truncated lifted linear system is both banded and Toeplitz, with band size determined by the decay of $\widehat{g}$. If $g$ is effectively bandlimited to $[-\delta,\delta]$ the proposed lifting-based algorithm can be implemented to run in $\mathcal{O}\left(\delta N(\log N+\delta^{2})\right)$-time, which is essentially FFT-time in $N$ for small $\delta$.

II Our Lifted Formulation

The following theorem forms the basis of our lifted setup.

Theorem 1**.**

Suppose $f:\mathbb{R}\rightarrow\mathbb{C}$ is piecewise smooth and compactly supported in $\left[-1,1\right]$. Let $g\in L^{2}\left(\left[-a,a\right]\right)$ be supported in $\left[-a,a\right]\subset[-1,1]$ for some $a<1$, with $\left\|g\right\|_{L^{2}}=1$. Then for all $\omega\in\mathbb{R}$,

[TABLE]

for all shifts $l\in[a-1,1-a]$.

Proof.

Denote by $S_{l}g$ the right shift of $g$ by $l$. The short-time Fourier transform (STFT) [10] of $f$ given $g$, at a shift $l$ and frequency $\omega$, is defined by

[TABLE]

The squared magnitude of the Fourier transform above is called a spectrogram measurement:

[TABLE]

where $h\left(t\right)=\overline{g\left(t-l\right)}e^{2\pi i\omega t}$. We calculate

[TABLE]

By Plancherel’s theorem, we have

[TABLE]

where in the last equality we have used

[TABLE]

And so, by Shannon’s Sampling theorem [16], applied to $\hat{f}$, we see that $\left|\mathcal{F}\left[f\cdot S_{l}g\right]\left(\omega\right)\right|^{2}$ is equal to

[TABLE]

where $\star$ denotes convolution.

Recall that $\mathcal{F}\left[f\star g\right]=\hat{f}\hat{g}$ so that $f\star g=\mathcal{F}^{-1}\left[\hat{f}\hat{g}\right].$ We calculate the Fourier transform

[TABLE]

and the Fourier transform $\mathcal{F}\left[\mbox{sinc}\pi\left(m+2\left(\cdot\right)\right)\right]\left(\xi\right)$ as

[TABLE]

With this, the spectrogram measurements $\left|\mathcal{F}\left[f\cdot S_{l}g\right]\left(\omega\right)\right|^{2}$ are given by

[TABLE]

Since $l$ is such that $\left[-l-1,-l+1\right]\cap\left[-a,a\right]=\left[-a,a\right]$, we have that $\left|\mathcal{F}\left[f\cdot S_{l}g\right]\left(\omega\right)\right|^{2}$ equals

[TABLE]

We have now proven the theorem. ∎

Using Theorem 1 we may now write

[TABLE]

where $A_{n}:=e^{-\pi iln}\hat{f}\left(\frac{n}{2}\right)\hat{g}\left(\frac{n}{2}-\omega\right).$

II-A Obtaining a Truncated, Finite Lifted Linear System

If $\hat{g}$ decays quickly we may truncate the sums above for a given frequency $\omega$ with minimal error. To that end, we pick the indices $j$ and $k$ so that $\left|\frac{k}{2}-\omega\right|\leq\delta$ and $\left|\frac{j}{2}-\omega\right|\leq\delta$ for some fixed $\delta\in\mathbb{N}$. If we denote

[TABLE]

then

[TABLE]

We may write

[TABLE]

where $\vec{X}_{l}\in\mathbb{C}^{4\delta+1}$ and $\vec{Y}_{\omega}\in\mathbb{C}^{4\delta+1}$ are the vectors

[TABLE]

This notation allows us to write our measurements in a lifted form

[TABLE]

Here, $\vec{Y}_{\omega}\vec{Y}_{\omega}^{*}$ is the rank-one matrix

[TABLE]

For each $\vec{X}_{l}\in\mathbb{C}^{4\delta+1}$, rewrite it as

[TABLE]

so that $m_{k}^{l}=e^{-\pi il\left(2k\right)}\hat{g}\left(k\right)$. Then construct the Toeplitz matrix $\mathbf{X}_{l}\in\mathbb{C}^{N\times N}$ as

[TABLE]

where $N$ is the number of frequencies $\omega$ being considered. Then we construct the block matrix ${\bf G}\in\mathbb{C}^{NK\times N}$ as

[TABLE]

where $K$ is the number of shifts of the window $g$.

Let $\mathbf{F}\in\mathbb{C}^{N\times N}$ be defined as

[TABLE]

where $n=\frac{N-1}{4}$. Note that ${\bf F}$ is composed of overlapping segments of the rank-1 matrices $\vec{Y}_{\omega}\vec{Y}_{\omega}^{*}$ for $\omega\in\{-n,\dots,n\}$. Thus, our measurements can be written as

[TABLE]

where $\vec{b}$ is defined in (I.2). By consistently vectorizing (II.2), we can obtain a simple linear system which can be inverted to learn $\vec{F}$, a vectorized version of ${\bf F}$. In particular, we have

[TABLE]

where the matrix $\mathbf{M}\in\mathbb{C}^{NK\times N^{2}}$ can be computed by, e.g., passing the canonical basis elements for $\mathbb{C}^{N\times N}$, ${\bf E}_{ij}$, through (II.2).

We solve the linear system (II.3) as a least squares problem; experiments have shown that $\mathbf{M}$ is of rank $NK$. The process of recovering the Fourier coefficients of $f$ from $\vec{F}$ is known as angular synchronization, and is described in detail in [6].

III Numerical Results

We test the Phase Retrieval algorithm above for two different choices of signal $f$. The first is a Gaussian signal $f\left(x\right)=2^{\frac{1}{4}}e^{-25\left(\frac{4x}{3}\right)^{2}}\chi_{\left[-1,1\right]}$, and the second is a modified Gaussian $f\left(x\right)=2^{\frac{1}{4}}e^{-8\pi x^{2}}\cos\left(24x\right)\chi_{\left[-1,1\right]}$. In both cases, the window used is the Gaussian $g\left(x\right)=c\cdot 2^{\frac{1}{4}}e^{-16\pi x^{2}}\chi_{\left[-\frac{1}{2},\frac{1}{2}\right]}$ where $c$ is a constant chosen so that $\left\|g\right\|_{L^{2}}=1$.

We use a total of 11 shifts of $g$ in each experiment. Since $g$ is supported on $\left[-\frac{1}{2},\frac{1}{2}\right]$, any two consecutive shifts are separated by $\frac{0.5}{11}$ (see Figure  III.1). We choose 61 values of $\omega$ from $\left[-15,15\right]$ sampled in half-steps, and set $\delta=7$.

The reconstructions in physical space are shown at selected grid points in Figures  III.2 and  III.3. The relative $\ell^{2}$ error in physical space is $1.47\times 10^{-3}$ for the first experiment and $1.872\times 10^{-2}$ for the second.

IV Future Work

While this paper addresses the 1D problem, the extension of this method to the 2D setting is an appealing avenue for future research. Indeed, preliminary results indicate that the underlying discrete method that forms the basis for this paper extends to two dimensions without too much difficulty. Furthermore, empirical results suggest that the method proposed here demonstrates robustness to noise, although we defer a detailed analysis (and derivation of an associated robust recovery guarantee) to future work.

Acknowledgement

This work was supported in part by the National Science Foundation grant NSF DMS-1416752.