# Inverting Spectrogram Measurements via Aliased Wigner Distribution Deconvolution and Angular Synchronization

We propose a two-step approach for reconstructing a signal x∈C^d from subsampled short-time Fourier transform magnitude (spectogram) measurements: First, we use an aliased Wigner distribution deconvolution approach to solve for a portion of the rank-one matrix x x^*. Second, we use angular syncrhonization to solve for x (and then for x by Fourier inversion). Using this method, we produce two new efficient phase retrieval algorithms that perform well numerically in comparison to standard approaches and also prove two theorems, one which guarantees the recovery of discrete, bandlimited signals x∈C^d from fewer than d STFT magnitude measurements and another which establishes a new class of deterministic coded diffraction pattern measurements which are guaranteed to allow efficient and noise robust recovery.

## Authors

• 13 publications
• 1 publication
• 2 publications
• 9 publications
06/04/2021

### Phase Retrieval for L^2([-π,π]) via the Provably Accurate and Noise Robust Numerical Inversion of Spectrogram Measurements

In this paper, we focus on the approximation of smooth functions f: [-π,...
08/22/2018

### Blind Phaseless Short-Time Fourier Transform Recovery

The problem of recovering a pair of signals from their blind phaseless s...
05/16/2022

### Phase retrieval of analytic signals from short-time Fourier transform measurements

Analytic signals belong to a widely applied class of signals especially ...
11/26/2021

### Non-Convex Recovery from Phaseless Low-Resolution Blind Deconvolution Measurements using Noisy Masked Patterns

This paper addresses recovery of a kernel h∈ℂ^n and a signal x∈ℂ^n from ...
02/03/2015

### Recovery of Piecewise Smooth Images from Few Fourier Samples

We introduce a Prony-like method to recover a continuous domain 2-D piec...
12/20/2021

### Toward Fast and Provably Accurate Near-field Ptychographic Phase Retrieval

Ptychography is an imaging technique which involves a sample being illum...
10/04/2019

### Admissible Measurements and Robust Algorithms for Ptychography

We study an approach to solving the phase retrieval problem as it arises...
##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 1. Introduction

The phase retrieval problem, i.e., reconstructing a signal from phaseless measurements, is at the core of many scientific breakthroughs related to the imaging of cells [46], viruses [44], and nanocrystals [15], and also advances in crystallographic imaging [28], optics [48], astronomy [18], quantum mechanics [16], and speech signal processing [2, 24]. As a result, many sophisticated algorithms, which achieve great empircal success, have been developed for solving this problem in applications throughout science and engineering (see [19, 21, 24] for widely used examples). Motivated by the success of these methods, the mathematical community has recently began to study the challenging problem of designing measurement masks and corresponding reconstruction algorithms with rigorous convergence guarantees and noise robustness properties (see, e.g., the work of Balan, Candès, Strohmer, and others [1, 2, 13, 25]). In this paper, we aim to extend the mathematical analysis of phaseless measurement maps and noise-robust reconstruction algorithms to include a broad class of phaseless Gabor measurements such as those that are utilized in, e.g., ptychographic imaging [14, 17, 40, 41].

Specifically, we will develop and analyze several algorithms for recovering (up to a global phase) a signal from the magnitudes of its inner products with shifts of masks that are locally supported in either physical space or Fourier space. The local support of these masks in physical space corresponds to the use of concentrated beams in ptychographic imaging to measure small portions of a large sample, whereas the local support of these masks in Fourier space simulates the recovery of samples belonging to a special class of deterministic coded diffraction patterns (CDP).

Following [33, 34], we will assume that we have a family of measurement masks, or windows, such that for all the nonzero entries of either or are contained in the set for some fixed where for any integer we let

 [n]0={0,1,…,n−1}

denote the set of the first nonnegative integers. Let be an integer which divides and let be the matrix-valued measurement map defined by its coordinate functions

 Y′k,ℓ\coloneqqY′k,ℓ(x)\coloneqq|⟨Sℓamk,x⟩|2+N′k,ℓ, (1.1)

for and where is the circular shift operator on defined for and by

 (Sℓy)j\coloneqqy((j+ℓ)modd), (1.2)

and represents an arbitrary perturbation due to, e.g., measurement noise or imperfect knowledge of the masks

Our goal is to reconstruct from these measurements. It is clear that for all so at best we can hope to reconstruct up to a global phase, i.e., up to the equivalence relation

 x∼x′ if x=eiϕx′ for some ϕ∈R.

Algorithms 1 and 2, presented in Section 4, will accomplish this goal in the special case where the masks are obtained by modulating a single mask Specifically, we let and for we let

 mk=Wkm (1.3)

where is the modulation operator given by

 (Wkm)j\coloneqqe2πijkdmj. (1.4)

As we will see, assuming that our masks have this form will allow us to recover even when the shift size is strictly greater than one. Towards this end, we let be the matrix-valued measurement map defined by its coordinate functions

 Yk,ℓ\coloneqqYk,ℓ(x)\coloneqq|⟨SℓWkm,x⟩|2+Nk,ℓ, (1.5)

where analogously to (1.1), represents an arbirtrary perturbation. We note that is the special case of where and the masks have the form (1.3). For positve integers, and which divide we let

 dK[K]0\coloneqq{0,dK,2dK,…,d−dK},and dL[L]0\coloneqq{0,dL,2dL,…,d−dL},

and we let be the partial measurement matrix obtained by restricting to rows in and columns in so that the -th entry of is given by

 (YK,L)k,ℓ=YkdK,ℓdL. (1.6)

Similarly, we let be the matrix obtained by restricting to rows and columns in and

Letting so that we note that

 (YK,L)k,ℓ=∣∣⟨x,SℓaWωkm⟩∣∣2+Nωk,ℓa=∣∣∣⟨x,e2πiℓaωkdWωkSℓam⟩∣∣∣2+Nωk,ℓa=∣∣⟨x,WωkSℓam⟩∣∣2+Nωk,ℓa. (1.7)

Therefore, forms a matrix of STFT magnitude measurements. Furthermore, when and (1.5) also encompasses a large class of masked Fourier magnitude measurements (i.e., CDP measurements) of the form

 (YK,L)k,ℓ=∣∣⟨x,WωkSℓam⟩∣∣2+Nk,ℓ=|(Fd Diag(m′ℓ) x)k|2+(NK,L)k,ℓ, (1.8)

where and is the discrete Fourier transform matrix whose entries are defined by

 (Fd)j,k\coloneqqe−2πijkd. (1.9)

Measurements similar to (1.8) are considered in, e.g., the recent works by Candès and others [3, 10, 11, 25]. However, their masks are usually generated randomly, whereas we will consider deterministically designed mask constructed as shifts of a single base mask

Our method for recovering is based on a two-step approach. Following the example of, e.g., [2, 13], we can lift the nonlinear, phaseless measurements (1.1) to linear measurements of the Hermitian rank-one matrix . Specifically, it can be shown that

 Y′k,ℓ(x) =⟨xx∗,Sℓamkmk∗S∗ℓa⟩+N′k,ℓ,

where the inner product above is the Hilbert-Schmidt inner product. Restricting, for the moment, to the case

(i.e. ) and assuming that the nonzero entries of are contained in the set for all , one can see that every matrix will have all of its nonzero entries concentrated near the main diagonal. Specifically, we have unless either or Therefore, letting be the restriction operator given by

 Tδ(G)ij={Gijif |i−j|<δ or |i−j|>d−δ0otherwise,

we see

 Y′k,ℓ(x)=⟨xx∗,Sℓmkmk∗S∗ℓ⟩+N′k,ℓ=⟨Tδ(xx∗),Sℓmkmk∗S∗ℓ⟩+N′k,ℓ,(k,ℓ)∈[K]0×[d]0.

Our lifted, linearized measurements are therefore given by where is defined by

 (A(X))(k,ℓ)=⟨X,Sℓmkmk∗S∗ℓ⟩for(k,ℓ)∈[K]0×[d]0andX∈Tδ(Cd×d). (1.10)

As a result, one can approximately solve for up to a global phase factor by evaluating on in order to recover a Hermitian approximation to , and then applying a noise robust angular synchronization method (e.g., see [45, 47]

) to obtain an estimate of

from . See [33, 34] for further details.

The following theorem summarizes previous work using this two-stage approach for the case where the nonzero entries of the masks are contained in the set for all .

###### Theorem 1 (See [33, 34]).

For , let and set and so that in (1.1). There exists a practical nonlinear reconstruction algorithm that takes in measurements for all and outputs an estimate that always satisfies

 minϕ∈[0,2π]∥∥x−eiϕxe∥∥2≤C(∥x∥∞min|x|2)(dδ)2κ∥N′∥F+Cd14√κ∥N′∥F. (1.11)

Here is the condition number of the linear map in (1.10) and is an absolute universal constant.

Furthermore, it is possible to choose masks such that (see [33]), and it is also possible to construct a single mask such that if then for the measurements that appear in (1.5) (see [34]). Also, if is sufficiently smalll, the algorithm mentioned above is guaranteed to require just total flops to achieve (1.11) up to machine precision.

Note that (1.11) guarantees that the algorithm in [33] referred to by Theorem 1 exactly inverts (up to a global phase) the measurement map in (1.1) for all nonvanishing in the noiseless setting (i.e., when ). Furthermore, the error between the recovered and original signal degrades gracefully with small amounts of arbitrary additive noise, and when the algorithm runs in essentially FFT-time. Indeed, a thorough numerical evaluation of this method has demonstrated it to be significantly more computationally efficient than competing techniques when, e.g., (see [33]). While the the first term of the error bound obtained in Theorem 1 exhibits quadratic dependence in we note that the main results of [32]

imply, at least heuristically, that polynomial dependencies on

are actually unavoidable in any upper bound like (1.11) when the masks are locally supported. As a result, both (1.11) as well as the new error bounds developed below generally must exhibit such polynomial dependences on .

### 1.1. Main Results

One of the main drawbacks of Theorem 1 is that it only holds for shifts of size i.e., when . In real ptychographic imaging applications, however, the equivalent of our parameter will in fact often at least with being moderately large. Therefore, we will consider recovery scenarios where both and are strictly less than in (1.6), and consider classes of and for which we can still guarantee noise-robust recovery results. This motivates our first new result, which allows us to recover bandlimited signals.

###### Theorem 2 (Convergence Gaurantees for Algorithm 2).

Let with and , and let

 μ2=min|p|≤γ−1,|q|≤δ−1∣∣∣Fd(ˆm∘Sp¯¯¯¯¯¯ˆm)q∣∣∣ > 0. (1.12)

Assume that , and also that and divide . Furthermore, suppose that the phaseless measurements (1.6) have noise dominated by the norm of so that

 ∥NK,L∥F≤β∥x∥22 (1.13)

for some . Then Algorithm 2 in Section 4 outputs an estimate to with relative error

 (1.14)

where is the partial Fourier matrix with entries and singular value . Furthermore, if is sufficently small, then Algorithm 2 is always guaranteed to require at most total flops to achieve (1.14) up to machine precision.

As mentioned earlier, Theorem 2 allows us to recover even when and are both strictly less than . Indeed, the total number of measurements, is independent of the sample size (though it does exhibit dependence on the parameters, and and it also requires that ). Nonetheless, the fact that Algorithm 2 exhibits robustness to arbitrary noise indicates that one can use it to quickly obtain a low-pass approximation to a sufficiently smooth using fewer than STFT magnitude measurements. We also note that Proposition 2, stated in Section 4, shows that locally supported masks with are relatively simple to construct.

Our second result utilizes the connection between CPD measurements and STFT magnitude measurements (see (1.7) and (1.8)) to provide a new class of deterministic CPD measurement constructions along with an associated noise-robust recovery algorithm. Unlike previously existing deterministic constructions (see, e.g., Theorem 3.1 in [10]) the following result presents a general means of constructing deterministic CDP masks using shifts of a single bandlimited mask .

###### Theorem 3 (Convergence Gaurantees for Algorithm 1).

Let with for some Let and let

 μ1\coloneqqmin|p|≤γ−1|q|≤ρ−1∣∣∣Fd(ˆm∘Sp¯¯¯¯¯¯ˆm)q∣∣∣ > 0. (1.15)

Fix an integer and assume that divides . Then, when Algorithm 1 in Section 4 will output an estimate of such that

 minϕ∈[0,2π]∥∥x−eiϕxe∥∥2≤Cd7/2∥ˆx∥∞∥Nd,L∥FL12μ1κ52⋅min|ˆx|2+C′d32L14√∥Nd,L∥Fμ1 (1.16)

for some absolute constants . Furthermore, if is sufficently small, then Algorithm 1 is always guaranteed to require just total flops to achieve (1.16) up to machine precision.

When Theorem 3 guarantees that CDP measurements suffice in order to recover any signal with a nonvanishing discrete Fourier transform in the noiseless setting as long as . Analogously to Proposition 2, Proposition 1, also stated in Section 4, provides a straightforward way to construct masks with For general and , Theorem 3 shows that one can reconstruct signals using CPD measurements based on windows with Fourier support in just -time when is sufficiently small. We also note that in addition to the theoretical guarantees provided by Theorems 2 and 3, Section 5 demonstrates that both Algorithm 1 and 2 are fast, accurate, and robust to noise in practice as well.

### 1.2. Related Work

The connections between theoretical time-frequency analysis and phaseless imaging (e.g., ptychography) have been touched on in the physics community many times over the past several decades. As noted well over two decades ago in [41] and later in [14], continuous spectrogram measurements can be written as the convolution of the Wigner distribution functions of the specimen and the probe Furthermore, [17] has pointed out that this allows one to recover the specimen of interest if enough samples are drawn so that the Heisenberg boxes sufficiently cover the time-frequency plane. In this work, we use similar ideas formulated in the discrete setting to efficiently invert the types of structured lifted linear maps as per (1.10) that appear in [33, 34], and use angular synchronization approaches to recover the signal up to a global phase. Specifically, we produce two new, efficient algorithms for inverting discrete spectrogram measurements that are provably accurate and robust to arbitrary additive measurement errors.

In [7], Bendory and Eldar prove results similar to some of those summarized in Theorem 1 in the case there and and they also demonstrate numerically that their algorithms are robust to noise. In this paper, we prove noise-robust recovery results, where we allow However, we make additional assumptions about either the support of or the supports of and We also note the very recent and excellent work of Rayan Saab and Brian Preskitt [39] as well as that of Melnyk, Filbir, and Krahmer [35] which both prove results similar to Theorem 2. As in Theorem 2, the results of [35, 39] can guarantee recovery with shift sizes . Their results primarily differ from Theorem 2 in that they don’t used Wigner Distribution Deconvolution (WDD) based methods. As a result, they consider different classes of masks and signals than we do.

Other related work includes that of Salanevich and Pfander [38, 42] which builds upon the work of Alexeev et al. [1] to establish noise robust recovery results for Gabor frame-based measurements. Their noise robust approach has similar characteristics to the approach taken here with the primary differences being that they require additional measurements beyond those provided by shifts and modulations of a single mask (see, e.g., equation (8) in [38]), and in some sense utilize the reverse of the approach taken here: Instead of first solving a linear system to obtain an approximation of (a portion of) , and then using angular synchronization to obtain an approximation to , the methods of [38, 42] instead first use angular synchronization methods to obtain frame coefficients of , and then reconstruct using the recovered frame coefficients.

The rest of the paper is organized as follows. In Section 2, we establish necessary notation and state a number of preliminary lemmas. Then, in Section 3, we establish several discrete and aliased variants of WDD, some of which can be used when the mask is locally supported in physical space, and others for when is locally supported in Fourier space. In Section 4, we prove Theorems 2 and 3 which provide recovery guarantees for our proposed methods and also state propositions which describe ways to design masks so that the assumptions of these theorems are valid. Finally, in Section 5, we evaluate our algorithms numerically and show that they are fast and robust to additive measurement noise.

## 2. Notation and Preliminary Results

For we let

 supp(x)\coloneqq{n∈[d]0:xn≠0}

denote the support of where, as in Section 1, We let denote the reversal of about its first entry, i.e.,

 (Rx)n=˜xn\coloneqqx−nmodd% for 0≤n≤d−1,

and we recall from (1.2) and (1.4) the circular shift and modulation operators given by and In order to avoid cumbersome notation, if is not an element of we will write in place of For , we define the Fourier transform of by

 ˆxk\coloneqq(Fdx)k=d−1∑n=0xne−2πinkd,

where as in (1.9), denotes the discrete Fourier transform matrix with entries for For and , we define circular convolution and Hadamard (pointwise) multiplication by

 (x∗dy)ℓ \coloneqqd−1∑n=0xnyℓ−n, and(x∘y)ℓ\coloneqqxℓyℓ,

and we define their componentwise quotient and componentwise absolute value by

 (xy)n=xnyn% and|x|n=|xn|.

For a matrix we let denote its -th column, and let denote its Frobenius norm. When proving the convergence of our algorithms, we will use the fact that, up to a reorganization of the terms, a banded matrix, whose nonzero entries are contained within entries of the main diagonal is equivalent to a matrix whose columns are the diagonal bands of the square, banded matrix. Towards this end, if and is a matrix with columns indexed from to so that column zero is the middle column, we let be the banded matrix with entries given by

 (C2κ−1(M))j,k={Mj,k−jif% |j−k|<κ or |j−k|>d−κ0otherwise (2.1)

for By construction, the columns of are the diagonal bands of with the middle column lying on the main diagonal. For example, in the case where

Below, we will state a number of lemmas, some of which are well known, which we will use in the proofs of our main results. Proofs are provided in the appendix. Our first lemma summarizes a number of properties of the discrete Fourier transform and the operators above.

###### Lemma 1.

For all and ,

The following lemma is the discrete analogue of the convolution theorem.

###### Lemma 2.

(Convolution Theorem) For all

 F−1d(ˆx∘ˆy)=x∗dy,

and

 (Fdx)∗d(Fdy)=dFd(x∘y).

In much of our analysis, we will have to consider the Hadamard product of a vector with a shifted copy of itself. The next three lemmas will be useful when we need to manipulate terms of that form.

###### Lemma 3.

Let and let . Then,

 (Fd(x∘Sω¯¯¯x))α=1de2πiωαd(Fd(ˆx∘S−α¯¯¯ˆx))ω.
###### Lemma 4.

Let and let Then,

###### Lemma 5.

Let and let . Then,

 ((x∘S−ℓy)∗d(¯¯¯˜x∘Sℓ¯¯¯˜y))k=((x∘S−k¯¯¯x)∗d(˜y∘Sk¯¯¯˜y))ℓ.

For a positive integer which divides , we introduce the subsampling operator

 Zs:Cd→Cds,

defined by

 (Zsx)n\coloneqqxns for n∈[ds]0.

The following lemma shows that taking the Fourier transform of a subsampled vector produses an aliasing effect.

###### Lemma 6.

(Aliasing) Let be a positive integer which divides Then for and ,

 (Fds(Zsx))ω=1ss−1∑r=0ˆxω−rds.

## 3. Aliased Wigner Distribution Deconvolution for Fast Phase Retrieval

As in Section 1, we let denote an unknown quantity of interest and let denote a known measurement mask, and consider measurements of the form (1.5). By (1.7), we see we may write as a noisy windowed Fourier magnitude measurement of the form

 Yk,ℓ=∣∣ ∣∣d−1∑n=0xnmn−ℓe−2πinkd∣∣ ∣∣2+Nk,ℓ,for 0≤k,ℓ≤d−1. (3.1)

Let and denote the -th columns of the measurement matrix and the noise matrix respectively, and, as in Section 1, let be the partial measurement matrix obtained by restricting to rows in and columns in so the the entries of are given by (1.6), and let be the analogous matrix obtained by restricting to rows and columns in and

Our goal is to recover (up to a global phase) from these measurements with an error that may be bounded in terms of the magnitude of the noise Our method will be based on the following result that is an aliased and discrete variant of the Wigner Distribution Deconvolution (WDD) approach presented in the continuous setting by Chapman in [14]. Together with Lemmas 9, 10, and 11, it will allow us to recover portions of the rank one matrices and

###### Theorem 4.

Let be the partial measurement matrix defined in (1.6), and let be the corresponding partial noise matrix. Let and be the matrices defined by

 ˜Y\coloneqqFLYTK,LFTKand˜N\coloneqqFLNTK,LFTK.

Then for any and ,

 ˜Yα,ω =KLd3dK−1∑r=0dL−1∑ℓ=0(Fd(ˆx∘SℓL−α¯¯¯ˆx))ω−rK(Fd(ˆm∘Sα−ℓL¯¯¯¯¯¯ˆm))ω−rK+˜Nα,ω (3.2) =KLd2dK−1∑r=0dL−1∑ℓ=0e−2πi(ℓL−α)(ω−rK)/d(Fd(ˆx∘SℓL−α¯¯¯ˆx))ω−rK(Fd(m∘Sω−rK¯m))ℓL−α+˜Nα,ω (3.3) =KLd2dK−1∑r=0dL−1∑ℓ=0e2πi(ℓL−α)(ω−rK)/d(Fd(x∘Sω−rK¯¯¯x))α−ℓL(Fd(ˆm∘Sα−ℓL¯¯¯¯¯¯ˆm))ω−rK+˜Nα,ω (3.4) =KLddK−1∑r=0dL−1∑ℓ=0(Fd(x∘Sω−rK¯¯¯x))α−ℓL(Fd(m∘Sω−rK¯¯¯¯¯m))ℓL−α+˜Nα,ω,. (3.5)

To aid in the readers understanding, before proving Theorem 4, we will first give a short proof of the following lemma which is the special case of (3.5) where . It is the direct analogue of Chapman’s WDD approach as formulated in the continuous setting in [14].

###### Lemma 7.

Let be the measurement matrix with entries defined as in (3.1) and let be the corresponding noise matrix. Then, the -th column of is given by

 ˜Yω=d⋅Fd(x∘Sω¯¯¯x)∘R(Fd(m∘Sω¯¯¯¯¯m))+˜Nω, (3.6)

where

###### The Proof of Lemma 7.

As noted in (3.1), we may write as the STFT of with window Therefore, by Lemma 1, parts 5 and 8, we see that for any ,

 yℓ =|Fd(x∘S−ℓm)|2+ηℓ =Fd((x∘S−ℓm)∗d(¯¯¯˜x∘Sℓ¯¯¯¯¯¯˜m))+ηℓ. (3.7)

Thus, taking a Fourier transform of and applying Lemma 1, part 1, yields

 (Fdyℓ)ω=d((x∘S−ℓm)∗d(¯¯¯˜x∘Sℓ¯¯¯¯¯¯˜m))−ω+(Fdηℓ)ω,

and so, by Lemma 5,

 (Fdyℓ)ω =d((x∘Sω¯¯¯x)∗d(˜m∘S−ω¯¯¯¯¯¯˜m))ℓ+(Fdηℓ)ω. (3.8)

Since taking the transpose of the above equation implies

Therefore, the -th columns of and satisfy

 (YTFTd)ω=d(x∘Sω¯¯¯x)∗d(˜m∘S−ω¯¯¯¯¯¯˜m)+(NTFTd)ω,

so, taking the Fourier transform of both sides and applying Lemmas 2 and 4 yields

 (FdYTFTd)ω =dFd(x∘Sω¯¯¯x)∘Fd(˜m∘S−ω¯¯¯¯¯¯˜m)+(FdNTFTd)ω

Recalling that and completes the proof. ∎

The following lemma applies analysis similar to the previous lemma to subsampled column vectors using Lemma 6.

###### Lemma 8.

For and ,

 (FKZdK(yℓ))ω=KdK−1∑r=0((x∘Sω−rK¯¯¯x)∗d(˜m∘