Inverting Spectrogram Measurements via Aliased Wigner Distribution Deconvolution and Angular Synchronization
Michael Perlmutter, Sami Merhi, Aditya Viswanathan, Mark Iwen
TL;DR
This paper introduces a novel two-step method combining aliased Wigner distribution deconvolution and angular synchronization for efficient phase retrieval from subsampled spectrograms, with theoretical guarantees for signal recovery.
Contribution
It presents two new algorithms for phase retrieval from spectrogram measurements and proves guarantees for signal recovery with fewer measurements and deterministic coded diffraction patterns.
Findings
Successful numerical performance of the proposed algorithms.
Guarantees for recovering signals from fewer than d measurements.
Introduction of deterministic coded diffraction patterns with robust recovery.
Abstract
We propose a two-step approach for reconstructing a signal 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 Second, we use angular syncrhonization to solve for (and then for 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 from fewer than STFT magnitude measurements and another which establishes a new class of deterministic coded diffraction pattern measurements which are guaranteed to allow…
Peer Reviews
No public reviews on file for this paper yet. If you reviewed it on a platform where reviews are public (OpenReview, ICLR, NeurIPS, ICML), you can paste yours below so the community can read it here.
Videos
No videos yet. Explain this paper in a talk, walkthrough, or lecture? Add one.
Inverting Spectrogram Measurements via Aliased Wigner Distribution Deconvolution and Angular Synchronization
Michael Perlmutter, Sami Merhi, Aditya Viswanathan, Mark Iwen
Abstract.
We propose a two-step approach for reconstructing a signal 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 Second, we use angular syncrhonization to solve for (and then for 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 from fewer than 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.
Key words and phrases:
Phase Retrieval, Spectrogram Measurements, Short-Time Fourier Transform (STFT), Wigner Distribution Deconvolution, Angular Synchronization, Ptychography.
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, \mathbf{m_{0}},\mathbf{m_{1}},\ldots,$$\mathbf{m_{K-1}}\in\mathbbm{C}^{d} such that for all the nonzero entries of either or are contained in the set for some fixed where for any integer we let
[TABLE]
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
[TABLE]
for and where is the circular shift operator on defined for and by
[TABLE]
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
[TABLE]
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
[TABLE]
where is the modulation operator given by
[TABLE]
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
[TABLE]
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
[TABLE]
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
[TABLE]
Similarly, we let be the matrix obtained by restricting to rows and columns in and
Letting so that we note that
[TABLE]
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
[TABLE]
where and is the discrete Fourier transform matrix whose entries are defined by
[TABLE]
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
[TABLE]
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
[TABLE]
we see
[TABLE]
Our lifted, linearized measurements are therefore given by Y^{\prime}_{k,\ell}(\mathbf{x})=\mathcal{A}\big{(}T_{\delta}(\mathbf{x}\mathbf{x}^{*})\big{)}_{(k,\ell)}+N^{\prime}_{k,\ell} where is defined by
[TABLE]
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
[TABLE]
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
[TABLE]
Assume that , and also that and divide . Furthermore, suppose that the phaseless measurements (1.6) have noise dominated by the norm of so that
[TABLE]
for some . Then Algorithm 2 in Section 4 outputs an estimate to with relative error
[TABLE]
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
[TABLE]
Fix an integer and assume that divides . Then, when Algorithm 1 in Section 4 will output an estimate of such that
[TABLE]
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
[TABLE]
denote the support of where, as in Section 1, We let denote the reversal of about its first entry, i.e.,
[TABLE]
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
[TABLE]
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
[TABLE]
and we define their componentwise quotient and componentwise absolute value by
[TABLE]
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
[TABLE]
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
[TABLE]
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 ,
- (1)
** 2. (2)
** 3. (3)
** 4. (4)
** 5. (5)
** 6. (6)
** 7. (7)
** 8. (8)
**
The following lemma is the discrete analogue of the convolution theorem.
Lemma 2**.**
(Convolution Theorem) For all
[TABLE]
and
[TABLE]
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,
[TABLE]
Lemma 4**.**
Let and let Then,
[TABLE]
Lemma 5**.**
Let and let . Then,
[TABLE]
For a positive integer which divides , we introduce the subsampling operator
[TABLE]
defined by
[TABLE]
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 ,
[TABLE]
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
[TABLE]
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
[TABLE]
Then for any and ,
[TABLE]
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
[TABLE]
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 ,
[TABLE]
Thus, taking a Fourier transform of and applying Lemma 1, part 1, yields
[TABLE]
and so, by Lemma 5,
[TABLE]
Since taking the transpose of the above equation implies
[TABLE]
Therefore, the -th columns of and satisfy
[TABLE]
so, taking the Fourier transform of both sides and applying Lemmas 2 and 4 yields
[TABLE]
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 ,
[TABLE]
Proof.
As in the proof of Lemma 7, the columns of the and satisfy
[TABLE]
Therefore, subtracting from both sides, taking the Fourier transform, applying Lemma 6, and then using (3.8) we see
[TABLE]
The lemma follows from the linearity of the Fourier transform and of the subsampling operator ∎
Now we shall prove Theorem 4.
The Proof of Theorem 4.
Noting that is obtained by subsampling the rows and columns of we see that the -th column of is given by
[TABLE]
Therefore, applying Lemma 8 we see
[TABLE]
Thus, the -th column of is given by
[TABLE]
Taking the Fourier transform of both sides and applying Lemmas 6, 2, and 4 we see that
[TABLE]
for all
Using the linearity of the Fourier transform and the definitions of and completes the proof of (3.5). (3.2), (3.3), and (3.4) follow by using Lemma 3 to see that
[TABLE]
and
[TABLE]
∎
3.1. Solving for Diagonal Bands of the Rank-One Matrices
We wish to use Theorem 4 to solve for diagonal bands of the rank-one matrix In the case where one can use (3.6) to see that for
[TABLE]
However, in general, the right-hand side of (3.2)-(3.5) are linear combinations of multiple terms and therefore, it is not as straightforward to solve for these diagonal bands. In this subsection, we present several lemmas which make different assumptions on the spatial and frequency supports of and and identify special cases where these sums reduce to a single nonzero term. In these cases, we will then be able to solve for diagonal bands of either or by formulas similar to the one above. We will use Lemmas 9 and 10 in the proofs of Theorems 2 and 3. We state Lemma 11 in order to demonstrate that Wigner deconvolution approach can also be applied to the setting considered in [33]. We will provide the proof of Lemma 10. The proofs of Lemmas 9 and 11 are nearly identical.
The first lemma in this section assumes that is bandlimited and the spatial support of is contained in an interval of length It allows us to recover diagonal bands of the rank-one matrix
Lemma 9**.**
Let with and . Let and divide and let be the partial measurement matrix defined as in (1.6) and let be the corresponding subsampled noise matrix. As in the statement of Theorem 4, let
[TABLE]
Then for any and ,
[TABLE]
Moreover, if for some and for some , and if or and or the sum above collapses to only one term, so that
[TABLE]
The next lemma is similar to the previous one, but replaces the assumptions that has compact spatial support and that is bandlimited and with the assumption that is bandlimited. It allows us to recover diagonals of
Lemma 10**.**
Let and assume . Let divide let be the partial measurement matrix defined as in (1.6), and let be the corresponding partial noise matrix. As in the statement of Theorem 4, let
[TABLE]
Then for any and ,
[TABLE]
Moreover, if for some , then for all and all such that either or the sum above reduces to a single term and
[TABLE]
Proof.
(3.9) follows from Theorem 4 by setting in (3.2). To prove the second claim, we note that by the assumption that unless
[TABLE]
If and , this can only occur if Indeed, if then
[TABLE]
and if then
[TABLE]
Therefore, all other terms in the above sum are zero, and the right-hand side of (3.9) reduces to the desired result. Likewise, if then (3.10) can only hold when ∎
As in Lemma 9, the following lemma assumes that the spatial support of is contained in an interval of length and allows us to recover diagonals of However, it differs in that it assumes that but does not assume that that is -bandlimited.
Lemma 11**.**
Let with . Let divide let be the partial measurement matrix defined as in (1.6), and let be the corresponding subsampled noise matrix. As in the statement of Theorem 4, let
[TABLE]
Then for any and ,
[TABLE]
Moreover, if for some , and if either or then for all , the sum above reduces to only one term and
[TABLE]
Remark 1**.**
For convenience, in Lemmas 9, 10, and 11 we have assumed that the support of or were contained in the first or entries. However, inspecting the proofs we see these results remain valid if these intervals are replaced with any other intervals of the same length.
4. Recovery Guarantees
In this section, we will present two algorithms which allow us to reconstruct from our matrix of noisy measurements and prove Theorems 2 and 3, presented in the introduction, which guarantee that these algorithms converge. Before providing the proofs of these theorems, we will first state two propositions which show that it is possible to design masks in such a way that the mask dependent constants and are nonzero. For proofs of these propositions, please see the appendix.
Proposition 1**.**
Let be bandlimited with , so that its Fourier transform may be written as
[TABLE]
for some real numbers . As in (1.15), let
[TABLE]
for some . If
[TABLE]
and
[TABLE]
then
Proposition 2**.**
Let be a compactly supported mask with , given by
[TABLE]
for some real numbers . As in (1.12), let
[TABLE]
for some . If
[TABLE]
and
[TABLE]
then
We will now prove Theorem 3, which we restate below for the convenience of the reader. See 3
The Proof of Theorem 3.
Let and let Then by Lemma 10, if then
[TABLE]
and therefore
[TABLE]
Substituting we see
[TABLE]
for all Likewise, for
[TABLE]
so, since substituting implies
[TABLE]
for all
In order to write the equations above in a compact form, we will construct three matrices, , and As in Section 2, for notatational convenience, we will index the columns of these matrices from to so that column zero is the middle column. For we let the -th column of be the diagonal band of which is terms off of the main diagonal, i.e.
[TABLE]
and we define the columns of and by
[TABLE]
and
[TABLE]
By construction, (4.5) and (4.6) imply
[TABLE]
Using the fact is unitary, we see
[TABLE]
where is as in (1.15). Let be the Hermitianizing operator
[TABLE]
and note that Since operator defined in (2.1) is linear and is Hermitian, (4.8) implies
[TABLE]
Let be the signum function,
[TABLE]
and let and , be the (componentwise) normalized versions of and , respectively, i.e.,
[TABLE]
and note that
[TABLE]
For all and we have that
[TABLE]
Therefore, we can apply (4.11) and the fact that for all
[TABLE]
to see
[TABLE]
for all Thus, by (4.9) and the fact that we see that
[TABLE]
Therefore, by Corollary 2 of [33], we have
[TABLE]
where is the vector of true phases of , and is the lead eigenvector of .
As in Algorithm 1, define by
[TABLE]
for Lemma 3 of [34] implies that
[TABLE]
and so
[TABLE]
Therefore, by (4.9), we see
[TABLE]
Together with (4.13) this yields
[TABLE]
(1.16) now follows from the fact that for all ∎
Theorem 2, restated below, provides recovery guarantees for Algorithm 2 under the assumptions that is compactly supported in space and that is bandlimited. The proof is somewhat similar to the proof of Theorem 3 but uses Lemma 9 in place of Lemma 10 and uses the Lemma 8 of [33] during the angular synchronization step.
See 2
The Proof of Theorem 2.
Analogously to the proof of Theorem 3, we apply Lemma 9 with and to the cases where and and then substitute, and to see that
- (1)
if , then
[TABLE] 2. (2)
if , then
[TABLE] 3. (3)
if , then
[TABLE] 4. (4)
if , then
[TABLE]
We will write the above equations in matrix form, with rows indexed from to and columns indexed from to , as
[TABLE]
where is matrix with entries defined by
[TABLE]
and the entries of and are given by
[TABLE]
with
[TABLE]
for Let be as in (2). Then, by the same reasoning as in (4.9), we see
[TABLE]
Furthermore, for all and such that and
[TABLE]
Therefore, we see where is the Vandermone matrix are given below by
[TABLE]
and is the partial autocorrelation type matrix
[TABLE]
Therefore, by (4.14) we have
[TABLE]
has full rank since it is a Vandermonde matrix with distinct nodes, and therefore since it has a left inverse given by Thus,
[TABLE]
and so, by (4.16)
[TABLE]
where is the minimal singular value of
Let be a reshaping operator, such that if is a matrix,
[TABLE]
for so that
[TABLE]
and let where is the Hermitianizing operator defined in (4.10). (1.13) and the fact that for all imply that so the fact that is Hermitian, together with (4.25) implies
[TABLE]
By Lemma 8 of [33], if is the lead eigenvalue of and is an associated normalized eigenvector, then
[TABLE]
(1.14) follows by taking the inverse Fourier transform of both sides. ∎
5. Numerical Experiments
We now present numerical experiments which demonstrate the robustness and efficiency of the proposed algorithms and provide comparisons to existing phase retrieval methods. These results were generated using the open source BlockPR MATLAB software package (freely available at [31]) on a desktop computer (iMac, 2017) with an Intel® Core™i7-7700 (7th generation, quad core) processor, 16GB RAM, and running macOS High Sierra and MATLAB R2018b. In all of our plots, each data point was obtained by averaging the results of trials.
Unless otherwise stated, we used i.i.d. mean zero complex Gaussian random test signals with measurement errors modeled using a (real) i.i.d. Gaussian noise model. We will report both the signal to noise ratio (SNR) and reconstruction error in decibels (dB) with
[TABLE]
where and denote the true signal, recovered signal, (Gaussian) noise variance, and number of measurements respectively.
We will present selected results comparing the proposed formulation against other popular phase retrieval algorithms such as PhaseLift [13] (implemented as a trace-regularized least-squares problem using the first order convex optimization package TFOCS [5, 6], Hybrid Input-Output/Error Reduction (HIO+ER) alternating projection algorithm [4, 20], and Wirtinger Flow [12]. We note that more accurate results using PhaseLift may be obtained using other solvers and software packages (such as CVX [22, 23]), albeit at a prohibitively expensive computational cost. For the HIO+ER algorithm, the following two projections were utilized: (i) projection onto the measured magnitudes, and (ii) projection onto the span of these measurement vectors. The initial guess was set to be the zero vector, although use of a random starting guess did not change the qualitative nature of the results. As is common practice, (see, for example, [20]) we implemented the HIO+ER algorithm in blocks of twenty-five HIO iterations followed by five ER iterations in order to accelerate the convergence of the algorithm. To minimize computational cost while ensuring convergence, the total number of HIO+ER iterations was limited to (see Figure 1).
5.1. Empirical Validation of Algorithm 1
In Algorithm 1, whose convergence is guaranteed by Theorem 3, we assume that our measurements are obtained using a bandlimited mask with To demonstrate the effictiveness of this algorithm, we performed numerical experiments on the following two types of masks:
[TABLE]
where denotes an i.i.d uniform random distribution on the interval , and
[TABLE]
The exponential mask in (5.2) is closely related to the deterministic masks first introduced in [34]. The mask-dependent constant (see (1.15) in Theorem 3) for the random mask, with and (and averaged over trials), was . The behavior for other choices of and was similar. For the exponential mask, this constant was . The qualitative and quantitative performance of the algorithm was similar with both families of masks.
We performed experiments with both Algorithm 1, as presented in Section 4, and also with a modified version which uses a post-processing procedure to obtain improved accuracy. The modified algorithm replaces Steps (5) and (7) of Algorithm 1 (referred to as Diag. Mag. Est. and Norm. Ang. Sync. in Figure 2(a)) with the eigenvector based magnitude estimation procedure (Eig. Mag. Est.) in Section 6.1 of [33], and the graph Laplacian based angular synchronization method (Graph Ang. Sync.) described in Algorithm 3 of B. Preskitt’s dissertation [39]. As seen in Figure 2(a), which plots the reconstruction error at various noise levels with , and a random mask constructed as in (5.1) with these changes offered improved reconstruction accuracy.
Figure 2(b) demonstrates the importance of the number of shifts As expected, the reconstructions using larger (which entails using more measurements, each corresponding to greater overlap between successive masked regions of the specimen) offered improved accuracy. In order to ensure that divides we varied the value of slightly for different values of As in Figure 2(a), we used random masks, constructed as in (5.1), with We observe that for larger values of , performance improved by about dB. We also note that, in practice, a suitable value of can be chosen depending on whether the proposed method is used as a reconstruction procedure or as an initializer for another algorithm.
In Figure 3(a), we compare the performance of the proposed method to other popular phase retrieval methods. Reconstruction errors for recovering a signal of length using shifts and a random mask with are plotted for different levels of noise. We see that the proposed method performs well in comparison to the other algorithms, and even nearly matches the significantly more expensive algorithms such as PhaseLift which are based on semidefinite programming (SDP). We note that the Wirtinger Flow method is sensitive to the choice of parameters and iteration counts. We used fewer total iterations (150 at 10dB SNR) at higher noise levels and more iterations (4500 at 60dB) at lower levels in order to ensure that the algorithm converged to the level of noise. We are not aware of any methodical procedure for setting the various algorithmic parameters when utilizing the (local) measurement constructions considered in this paper. We also note that, of the algorithms considered, Algorithm 1 is the only one that has a theoretical convergence guarantee which applies to this class of spectrogram-type measurements.
Figure 3(b) plots the corresponding execution time for the various algorithms as a function of the problem size . In this case, random masks were chosen with along with shifts. The figure confirms the essentially FFT–time computational cost of Algorithm 1. Furthermore, it also shows that while the post-processing procedure of modifying steps (5) and (7) does increase the computational cost of the algorithm, it does not increase it drastically.222The modified Step (7) uses MATLAB’s eigs command which can be computationally inefficient for this problem for large ; we defer a more detailed analysis and more efficient implementations to future work In particular, even with these modifications, the proposed method provides best–in–class computational efficiency, and is significantly faster than the HIO+ER, Wirtinger Flow, and PhaseLift algorithms.
5.2. Empirical Validation of a Lemma 11 Based Approach
We next provide numerical results validating an approach based on Lemma 11 that applies the Wigner deconvolution method to the setting considered in [33]. As in Theorem 2, we assume that and we also add the assumption that In this setting, we may apply Lemma 11 and then solve for diagonal bands of in a manner analogous to Algorithms 1 and 2. We then can recover by applying the same angular synchronization procedure as in Algorithm 1. We note that, because we are using Lemma 11 rather than Lemma 9, we do not need to assume that is bandlimited as we do in Theorem 2. As in Section 5.1, we conducted experiments with both deterministicly constructed and randomly constructed masks, and found that we obtained similar results for both families of masks. The figures below use the exponential mask construction first introduced in [33],
[TABLE]
and therefore allow us to directly compare the performance of the proposed method with the algorithm introduced in [33]. The mask-dependent constant (see (1.12 in Theorem 2) for this mask, with and , was , with similar behavior for different choices of and . Figure 4(a) plots the reconstruction error with , and as in (5.3). Results with and without the post-processing modifications described in Section 5.1 are provided, along with results from [33] with and without the modified (see §6.1 in [33]) magnitude estimation and HIO+ER post-processing ( iterations).
As can be seen in Figure 5.4, the post-processing procedure yields a small improvement of about -dB in the reconstruction error, especially at low noise levels. We observe that the Wigner deconvolution based approach yields numerical performance which is comparable to [33] in the settings where the theoretical guarantees of [33] are applicable, while also adding the additional flexibility of allowing shifts of length under certain assumptions on either or as discussed in Theorems 2 and 3.
Next, we investigate the reconstruction accuracy as a function of , the number of Fourier modes. Figure 4(b) plots reconstruction error in recovering a test signal for and respectively, with the exponential masks defined as in (5.3) with As in Figure 2(b), we vary the signal length slighlty, in order to ensure that divides As expected, the plot shows that reconstruction accuracy improves when increases, i.e., when more measurements are acquired.
For completeness, we include noise robustness and execution time plots comparing the performance of the proposed method to the HIO+ER, PhaseLift, and Wirtinger Flow algorithms in Figures 5(a) and 5(b) respectively. From Figure 5(a), we see that the proposed method (both with and without the modified magnitude estimation/angular synchronization procedures) performs well in comparison to HIO+ER and the other algorithms across a wide range of SNRs. Furthermore, Figure 5(b) demonstrates the essentially FFT–time computational cost of the method as well as the best-in-class computational efficiency when compared to other competing algorithms.
5.3. Empirical Validation of Algorithm 2
We now provide numerical results validating Algorithm 2, whose convergence is guaranteed by Theorem 2. We begin by noting that the Vandermonde matrix defined in (4.24) often has a large condition number which poses a challenge in the accurate evaluation of Step (4) in Algorithm 2. One possible solution is to utilize the Tikhonov regularized solution (see [27] for example), where the regularization parameter is chosen using a procedure such as the L-curve method [26]. However, empirical simulations suggest that this procedure is not sufficiently robust to achieve reconstruction accuracy up to the level of added noise. Therefore, we replace Steps (4)–(6) in Algorithm 2 by a modified non-stationary iterated Tikhonov method inspired by the work of Buccini et al. in [9], which we detail in Algorithm 3. This procedure works by iteratively computing a Tikhonov regularized solution to the equation in Step (4) of Algorithm 2; however, at each step, the solution is applied to the residual of , with a geometrically decreasing regularization parameter. Buccini et al. showed that a similar iterative procedure has benefits over traditional Tikhonov regularization for more standard linear systems. While our problem setting is different, our empirical results suggest a similar benefit. We defer a more detailed theoretical analysis to future work.
Figure 5.6 presents empirical evaluation of the noise robustness and computational efficiency of Algorithm 2 with the Modified Iterated Tikhonov Method of Algorithm 3. Figure 6(a) plots the reconstruction error with signals of length with frequency support of length using (complex random) masks with spatial support of length We used Fourier modes and shifts, and utilized the following iterated Tikhonov parameters: , , and chosen using the L-curve method. We note that using standard Tikhonov regularization (Alg. 2 in Figure 6(a)) yields rather poor results. An aggressive regularization parameter has to be chosen to surmount the ill-conditioning effects in Step (4) of Algorithm 2. Consequently, even a few iterations of the HIO+ER algorithm performs better than Algorithm 2. However, using the modified iterated Tikhonov procedure (Alg. 2 (w/ Alg. 3) in Figure 6(a)) yields significantly improved results, with a clear improvement in noise robustness over even the HIO+ER algorithm. Furthermore, Figure 6(b) plots the execution time as a function of the problem size for both Algorithms 2 and 3. The plot confirms that the modified iterative Tikhonov procedure of Algorithm 3 does not impose a significant computational burden.333We note that Step (1) of Algorithm 3 is computationally tractable since is typically small, and that the matrix in Step 2(c) can be pre-computed. Indeed, both Algorithms 2 and Algorithm 2 with the Modified Iterated Tikhonov Method of Algorithm 3 are faster than the HIO+ER algorithm. We note that more efficient implementations (involving fast computations of Vandermonde systems) of all the algorithms in Figure 6(b) may be possible; we defer this to future research.
6. Future Work
In future work, one might develop variants of the algorithms presented here for two-dimensional problems along the lines of [29]. Additionally, one might also develop variations of these algorithms for recovering compactly supported functions from sampled spectrogram measurements (see [36]) in the continuous setting. Furthermore, another, perhaps less direct, extension of these works would be to attempt to apply the Wigner distribution methods used here to the sparse phase retrieval problem. In, e.g., [30] it was shown that sparse vectors with can be recovered up to a global phase from only magnitude measurements of the form . Thus, somewhat surprisingly, sparse phase retrieval problems generally do not require significantly more measurements to solve than compressive sensing problems. One may be able to generate new sparse phase retrieval methods for STFT magnitude measurements of the type considered here by replacing the standard Fourier techniques used in the methods above with sparse Fourier transform methods [8, 37, 43]. It has been shown that sparse phase retrieval problems can be solved in sublinear-time [47]. The further development of sublinear-time methods for solving sparse phase retrieval problems involving STFT magnitude measurements could prove valuable in the future for use in extremely large imaging scenarios.
Acknowledgements
Mark Iwen was supported in part by NSF DMS-1912706 and NSF CCF-1615489. Sami Merhi was supported in part by NSF CCF-1615489.
Appendix
In this section, we will prove the lemmas from Section 2 as well as Propositions 1 and 2.
The Proof of Lemma 1.
Let and let
Part 1:
[TABLE]
Part 2:
[TABLE]
Part 3:
[TABLE]
Part 4:
[TABLE]
Part 5:
[TABLE]
Part 6:
[TABLE]
Part 7:
[TABLE]
Part 8: For all
[TABLE]
∎
The Proof of Lemma 2.
For ,
[TABLE]
Therefore, so multiplying by proves the first claim. To verify the second claim, note that by Lemma 1 part 1,
[TABLE]
∎
The Proof of Lemma 3.
Let and let Observe that
[TABLE]
∎
The Proof of Lemma 4.
For any and any it is straightforward to check that
[TABLE]
Therefore,
[TABLE]
∎
The Proof of Lemma 5.
Let and let Then,
[TABLE]
∎
The Proof of Lemma 6.
For , the Fourier inversion formula states that
[TABLE]
Therefore, for all
[TABLE]
∎
The Proof of Proposition 1.
Let be a bandlimited mask, whose Fourier transform may be written as
[TABLE]
for some real numbers which satisfy (4.1) and (4.2). Let and recall that is defined by
[TABLE]
For , we have
[TABLE]
and for ,
[TABLE]
Therefore, for any and any
[TABLE]
where is some real number depending on and Using the assumptions (4.1) and (4.2) we see that
[TABLE]
With this,
[TABLE]
where the last inequality follows by 6.2. Therefore, is nonzero for all and and so
∎
The Proof of Proposition 2.
Let
[TABLE]
be a compactly supported mask, where are real numbers which satisfy (4.3) and (4.4). Let and recalll that is defined by
[TABLE]
By Lemma 3, it suffices to show that
[TABLE]
for all and all If , then
[TABLE]
and if , then
[TABLE]
Therefore, for all and all
[TABLE]
where is some real number depending on and By the same reasoning as in the proof of Proposition 1, this combined with (4.3) and (4.4) implies that for all and all
∎
The reference list from the paper itself. Each links out to its DOI / PubMed record.
- 1[1] B. Alexeev, A. S. Bandeira, M. Fickus, and D. G. Mixon. Phase Retrieval with Polarization. SIAM Journal on Imaging Sciences , 7(1):35–66, 2014.
- 2[2] R. Balan, P. Casazza, and D. Edidin. On signal reconstruction without phase. Applied and Computational Harmonic Analysis , 20(3):345–356, 2006.
- 3[3] A. S. Bandeira, Y. Chen, and D. G. Mixon. Phase retrieval from power spectra of masked signals. Information and Inference: a Journal of the IMA , 3(2):83–102, 2014.
- 4[4] H. H. Bauschke, P. L. Combettes, and D. R. Luke. Phase retrieval, error reduction algorithm, and Fienup variants: A view from convex optimization. Journal of the Optical Society of America. A, Optics, Image science, and Vision , 19(7):1334–1345, 2002.
- 5[5] S. Becker, E. J. Candès, and M. Grant. Templates for convex cone problems with applications to sparse signal recovery. Mathematical Programming Computation , 3(3):165–218, Aug. 2011.
- 6[6] S. Becker, E. J. Candes, and M. Grant. TFOCS: Templates for first-order conic solvers, version 1.3.1. http://cvxr.com/tfocs , Sep. 2014.
- 7[7] T. Bendory, Y. C. Eldar, and N. Boumal. Non-convex phase retrieval from stft measurements. IEEE Transactions on Information Theory , 64(1):467–484, 2017.
- 8[8] S. Bittens, R. Zhang, and M. A. Iwen. A deterministic sparse fft for functions with structured fourier sparsity. Advances in Computational Mathematics , 45(2):519–561, Apr 2019.