In this paper, we consider the continuous one-dimensional sparse phase retrieval problem to reconstruct a complex-valued signal from the modulus of its Fourier transform. Applications of this problem occur in electron microscopy, wave front sensing, laser optics [6, 7] as well as in X-ray crystallography and speckle imaging [8]. For the posed problem, we will show that for sparse signals the given Fourier intensities are already sufficient for an almost sure unique recovery, and we will give a construction algorithm to recover f.

We assume that the sparse signal is either of the form (1.1)f(t)=∑j=1Ncj(0)δ(t−Tj) or, for m > 0, (1.2)f(t)=∑j=1Ncj(m)Bj,m(t) with cj(m)∈ℂ, T_j ∈ ℝ for j = 1, …, N, where δ denotes the Delta distribution, and B_{j, m} is the B-spline of order m being determined by the (real) knots T_j < T_j+1 < … < T_j+m. We want to recover these signals from the Fourier intensities |f^(ω)|2 and will show that only O(N²) samples are needed to recover f, i.e., all coefficients cj(m), j = 1, …, N and knots T_j, j = 1, …, N + m, almost surely up to trivial ambiguities. The proposed procedure is constructive and consists in two steps. In a first step, we employ Prony's method to determine the coefficients and frequencies of the exponential sum |F[f](ω)|2. These frequencies are of the form T_j − T_k with j, k ∈ {1, …, N + m}. If these knot differences T_j − T_k are pairwise different for j ≠ k, then we can use this information in a second reconstruction step to compute the knots T_j and the coefficients c_j, and thus the desired signal.

While the general phase retrieval problem has been extensively studied for a long time, the special case of sparse phase retrieval grew to a strongly emerging field of research only recently, particularly often connected with ideas from compressed sensing. Most of the papers consider a discrete setting, where the N-dimensional real or complex k-sparse vector x has to be reconstructed from measurements of the more general form |〈aj,x〉|2 with vectors a_j forming the rows of a measurement matrix A ∈ ℂ^M×N. The needed number M of measurements depends on the sparsity k.

If A presents rows of a Fourier matrix, this setting is close to the sparse phase retrieval problem considered in optics, see e.g.,[9]. Here the problem is first rewritten as (non-convex) rank minimization problem, then a tight convex relaxation is applied and the optimization problem is solved by a re-weighted l₁-minimization method. The related approach in Eldar et al. [10] employs the magnitudes of the short-time Fourier transform and applies the occurring redundancy for unique recovery of the desired signal. A corresponding reconstruction algorithm is then based on an adaptation of the GESPAR algorithm in Shechtman et al. [11].

In Li and Voroninski [12], the measurement matrix A is taken with random rows and the PhaseLift approach [13] leads to a convex optimization problem that recovers the sparse solution with high probability. Employing a thresholded gradient descent algorithm to a non-convex empirical risk minimization problem that is derived from the phase retrieval problem, Cai et al. [14] have established the minimax optimal rates of convergence for noisy sparse phase retrieval under sub-exponential noise.

Other papers rely on the compressed sensing approach to construct special frame vectors a_j to ensure uniqueness of the phase retrieval problem with high probability, where the number of needed vectors is O(k), see e.g., [15–17].

We would like to emphasize that all approaches employing general or random measurement matrices in phase retrieval are quite different in nature from our phase retrieval problem based on Fourier intensity measurements. In this paper, we want to stick on considering Fourier intensity measurements because of their particular relevance in practice.

Early attempts to exploit sparsity of a discrete signal for unique recovery using Fourier intensities go back to unpublished manuscripts by Yagle [18, 19], where a variation of Prony's method is applied in a non-iterative algorithm to sparse signal and image reconstruction. Unfortunately, the algorithm proposed there not always determines the signal support correctly.

The continuous one-dimensional phase retrieval problem has been rarely discussed in the literature, see [5, 8, 20–22]. In the preprint [8], the authors also considered the recovery of sparse continuous signals of the form (1.1). However, in that paper the sparse phase retrieval problem is in turn transferred into a turnpike problem that is computationally expensive to solve. Moreover there exist cases, where a unique solution cannot be found, see [23]. Our method circumvents this problem by proposing an iterative procedure to fix the signal support (resp. the knots of the signal represented as a B-spline function) where the corresponding signal coefficients are evaluated simultaneously.

In Section 2, we shortly recall the mathematical formulation of the considered sparse phase retrieval problem and the notion of trivial ambiguities of the phase retrieval problem that always occur.

Section 3 is devoted to the special case of phase retrieval for signals of the form (1.1). Using Prony's method, we give a constructive proof for the unique recovery of the N-sparse signal f up to trivial ambiguities using ³/₂N(N−1)+1 Fourier intensity measurements. Here we have to assume that the knot differences T_j − T_k are pairwise different.

In Section 4, the ansatz is generalized to the unique recovery of spline functions of the form (1.2) where we need to employ ³/₂(N + m)(N + m − 1) + 1 Fourier intensity measurements. In Section 5, we present an explicit algorithm for the considered sparse phase retrieval problem and illustrate it at different examples.

Let us assume that the knot differences T_j − T_k in (3.1) are pairwise different for j ≠ k. The squared Fourier intensity can then be written in the form

where we assume that the frequencies τ_ℓ, ℓ = −^N(N−1)/₂, …, ^N(N−1)/₂ are ordered by size. Obviously, the frequency differences τ_ℓ satisfy τ_−ℓ = −τ_ℓ. Further, each τ_ℓ > 0 corresponds to a difference T_j − T_k for some j > k, and the related coefficient γ_ℓ then equals to cj(0)c¯k(0). For the zero frequency τ₀ = 0, we have γ0:=∑j=1N|cj(0)|2.

In the first step, we want to recover all frequencies τ_ℓ and the corresponding coefficients γ_ℓ, ℓ = 0, …, ^N(N−1)/₂ of P(ω). However, at this stage, the bijective mapping between ℓ > 0 and (j, k) with j > k such that τ_ℓ = T_j − T_k will be still unknown and needs to be found in a second reconstruction step. In order to recover the frequency differences τ_ℓ and the unknown coefficients γ_ℓ from the exponential sum (3.2) we employ Prony's method [24, 25].

Let h > 0 be chosen such that hτ_ℓ < π for all ℓ = 1, …, ^N(N−1)/₂. Using the intensity values P(hk)=|F[f](hk)|2, k = 0, …, 2N(N − 1) + 1, the unknown parameters γ_ℓ and τ_ℓ in (3.2) can be determined by exploiting the algebraic Prony polynomial Λ(z) defined by

where λ_k denote the coefficients in the monomial representation of Λ(z). Obviously, Λ(z) is always a monic polynomial, which means that λ_N(N−1)+1 = 1.

Using the definition of the Prony polynomial Λ(z) in (3.3), we observe that

for m = 0, …, N(N − 1). Consequently, the vector of remaining coefficients λ:=(λ0,…,λN(N-1))T of the Prony polynomial Λ(z) can be determined by solving the system of linear equations

with H:=(P(h(k+m)))m,k=0N(N-1) and h:=(P(h(N(N-1)+1+m)))m=0N(N-1). Since the Hankel matrix H can be written as

with the Vandermonde matrix V :=(e−hkτℓ)ℓ=−N(N−1)/2,k=0N(N−1)/2,N(N−1)+1, the system of linear equations (3.4) possesses a unique solution if and only if the unimodular values e-ihτℓ differ pairwise for ℓ = −^N(N−1)/₂, …, ^N(N−1)/₂. This assumption has been ensured by choosing an h such that hτ_ℓ ∈ (−π, π), since the τ_ℓ had been supposed to be pairwise different.

Knowing the coefficients λ_k of Λ(z), we can determine the unknown frequencies τ_ℓ by evaluating the roots of the Prony polynomial (3.3). The coefficients γ_ℓ can now be computed by solving the over-determined equation system

with a Vandermonde-type system matrix.

The procedure summarized above is Prony's method, adapted to the non-negative exponential sum P(ω) in (3.2). In the numerical experiments in Section 5, we will apply the approximate Prony method (APM) in Potts and Tasche [26]. APM is based on the above considerations but it is numerically more stable and exploits the special properties γ-ℓ=γ¯ℓ and τ_−ℓ = −τ_ℓ for ℓ = 0, …, ^N(N−1)/₂.

Let us now investigate the question, how many intensity values are at least necessary for the recovery of P(ω) in (3.2). Counting the number of unknowns of P(ω) in (3.2), we only need to recover the ³/₂N(N − 1) + 1 real values γ₀ and Re γ_ℓ, Im γ_ℓ, τ_ℓ, for ℓ = 1, …^N(N−1)/₂. We will show now that using the special structure of the real polynomial P(ω) in (0.2) and of the Prony polynomial Λ(z) in (0.3), we indeed need only ³/₂N(N − 1) + 1 exact equidistant real measurements P(kh), k = 0, …, ³/₂N(N − 1) to recover all parameters determining P(ω). This can be seen as follows.

Reconsidering Λ(z) in (3.3) with τ₀ = 0 and τ_ℓ = −τ_−ℓ, we obtain

where all occurring coefficients λ_k are real. Moreover, since

is antisymmetric, it follows that

and particularly λ_N(N−1)+1 = −λ₀ = 1. In order to determine the unknown coefficients λ_k, k = 1, …, ^N(N−1)/₂ of

we employ (3.2) and observe that for m = 0, …, ^N(N−1)/₂ − 1,

Therefore, the vector of unknown coefficients λ :=(λ1,…,λN(N−1)/2)T can be evaluated from the system ∑k=1N(N−1)/2λk[P(h(k+m))−P(h(N(N−1)+1+m−k))]=[P(hm)−P(h(N(N−1)+1+m))]      (m=0,…,N(N−1)/2−1). The frequency differences τ_ℓ are then extracted from the zeros of Λ(z), and the coefficients γ_ℓ, ℓ = 0, …, ^N(N−1)/_2, are computed as in (3.5) but with k = 0, …, ³/₂N(N − 1).

Having determined the frequency differences τ_ℓ as well as the corresponding coefficients γ_ℓ of (3.2), we want to reconstruct the parameters T_j and cj(0), j = 1, …, N, of f in (1.1) in a second step.

Theorem 3.1. Let f be a signal of the form (1.1), whose knot differences T_j − T_k differ pairwise for j, k ∈ {1, …, N} with j ≠ k, and whose coefficients satisfy |c1(0)|≠|cN(0)|. Further, let h be a step size such that h(T_j − T_k) ∈ (−π, π) for all j, k. Then f can be uniquely recovered from its Fourier intensities |F[f](hℓ)| with ℓ = 0, …, ³/₂N(N − 1) up to trivial ambiguities.

Proof. Applying Prony's method to the given data |F[f](hℓ)|, we can compute the frequency differences τ_ℓ and the related coefficients γ_ℓ of the squared Fourier intensity (3.2). We denote by T :={τℓ}ℓ=1N(N−1)/2 the list of obtained positive frequencies ordered by size. Now, we need to recover the mapping ℓ → (j, k) such that τ_ℓ = T_j − T_k, where we can assume that j > k for ℓ > 0.

Obviously, the maximal distance τ_^N(N−1)/2 is now equal to the length T_N − T₁ of the unknown f in (1.1). Due to the trivial shift ambiguity, we can assume without loss of generality that T₁ = 0 and T_N = τ_N(N−1)/2. Further, the second largest distance τ_{(N(N−1)/2)−1} corresponds either to T_N−1 − T₁ or to T_N − T₂. Due to the trivial reflection and conjugation ambiguity, we can assume that T_N−1 = T_N−1 − T₁ = τ_{(N(N−1)/2)−1}. By definition, there exists a τℓ*>0 in our sequence of parameters T such that τℓ*+τ(N(N−1)/2)−1=τN(N−1)/2, and τℓ* is hence equal to the knot difference T_N − T_N−1. Thus, we obtain cN(0)c¯1(0)=γN(N−1)/2,  cN−1(0)c¯1(0)=γ(N(N−1)/2)−1, andcN(0)c¯N−1(0)=γℓ*. These equations lead us to cN(0)=γN(N−1)/2c¯1(0),  cN−1(0)=γ(N(N−1)/2)−1c¯1(0), and thus to |c1(0)|2=γN(N−1)/2γ(N(N−1)/2)−1¯γℓ*. Since f can only be recovered up to a global rotation, we can assume that c1(0) is real and non-negative, which allows us to determine the coefficients c1(0), cN(0), and cN-1(0) in a unique way.

Having fixed T_N = T_N − T₁ = τ_N(N−1)/2 and T_N−1 = T_N−1−T₁ = τ_{(N(N−1)/2)−1} we notice that the third largest distance τ_{(N(N−1)/2)−2} is either equal to T_N−T₂ or to T_N−2−T₁ = T_N−2. As before, there exists a frequency τℓ* such that τ(N(N-1)/2)-2+τℓ*=τN(N-1)/2.

Case 1: If τ_{(N(N−1)/2)−2} = T_N−T₂, then we have τℓ∗=τN(N−1)/2−τ(N(N−1)/2)−2=(TN−T1)−(TN−T2)    =T2−T1 with the related coefficient γℓ*=c2(0)c¯1(0). Moreover, we have γ(N(N-1)/2)-2=cN(0)c¯2(0) such that (3.6)c2(0)=γℓ*c¯1(0)=γ¯(N(N−1)/2)−2c¯N(0). Case 2: If τ_{(N(N−1)/2)−2} = T_N−2−T₁, then we have τℓ∗=τN(N−1)/2−τ(N(N−1)/2)−2 = (TN−T1)−(TN−2−T1)    =TN−TN−2 with the related coefficient γℓ*=cN(0)c¯N-2(0) and γ(N(N-1)/2)-2=cN-2(0)c¯1(0). Thus, (3.7)cN−2(0)=γ¯ℓ∗c¯N(0)=γ(N(N−1)/2)−2c¯1(0). However, only one of the two equalities in (3.6) and (3.7) can be true, since if both were true then γℓ*c¯N(0)=c¯1(0)γ¯(N(N-1)/2)-2 and c¯1(0)γ¯ℓ*=c¯N(0)γ(N(N-1)/2)-2 lead to |cN(0)c1(0)|=|γ(N(N−1)/2)−2γℓ*|=|c1(0)cN(0)| contradicting the assumption that |cN0|≠|c10|. Consequently, either the equation in (3.6) or the equation in (3.7) holds true and we can either determine T₂ with c2(0) or T_N−2 with cN-2(0). Removing all frequency differences τ_ℓ from the sequence of distances T that correspond to the differences T_j − T_k of the recovered knots, we can repeat this approach to find the remaining coefficients and knots of f inductively.                                                                                                          □

If we identify the space of complex-valued signals of the form (1.1) with the real space ℝ^3N, the condition that two knot differences T_j1−T_k1 and T_j2−T_k2 are equal for fixed indices j₁, j₂, k₁, and k₂ defines a hyperplane with Lebesgue measure zero. An analogous observation follows for the condition |c1(0)|=|cN(0)|. The signals excluded in Theorem 3.1 hence form a negligible null set.

Corollary 3.2. Almost all signals f in (1.1) can be uniquely recovered from their Fourier intensities |F[f]| up to trivial ambiguities.

Remark 3.3. 1. Since the proof of Theorem 3.1 is constructive, it can be used to recover an unknown signal (1.1) analytically and numerically. If the number N of spikes is known beforehand then the assumption of Theorem 3.1 can be simply checked during the computation. If the assumption regarding pairwise different distances T_j − T_k is not satisfied, then the application of Prony's method in the first step yields less than N(N − 1) + 1 pairwise distinct frequency differences τ_ℓ. The second assumption |cN0|≠|c10| can be verified in the second step, where c1(0), cN-1(0), and cN(0) are evaluated.

2. A similar phase retrieval problem had been transferred to a turnpike problem in Ranieri et al. [8]. The turnpike problem deals with the recovery of the knots T_j from an unlabeled set of distances. Although this problem is solvable under certain conditions, a backtracing algorithm can have exponential complexity in the worst case, see [27].

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.