Stable Gabor phase retrieval in Gaussian shift-invariant spaces via biorthogonality

We study the phase reconstruction of signals $f$ belonging to complex Gaussian shift-invariant spaces $V^\infty(\varphi)$ from spectrogram measurements $|\mathcal{G} f(X)|$ where $\mathcal{G}$ is the Gabor transform and $X \subseteq \mathbb{R}^2$. An explicit reconstruction formula will demonstrate that such signals can be recovered from measurements located on parallel lines in the time-frequency plane by means of a Riesz basis expansion. Moreover, connectedness assumptions on $|f|$ result in stability estimates in the situation where one aims to reconstruct $f$ on compact intervals. Driven by a recent observation that signals in Gaussian shift-invariant spaces are determined by lattice measurements [Grohs, P., Liehr, L., Injectivity of Gabor phase retrieval from lattice measurements, Appl. Comput. Harmon. Anal. 62 (2023), pp. 173-193] we prove a sampling result on the stable approximation from finitely many spectrogram samples. The resulting algorithm provides a provably stable and convergent approximation technique. In addition, it constitutes a method of approximating signals in function spaces beyond $V^\infty(\varphi)$, such as Paley-Wiener spaces.


Introduction
Phaseless signal reconstruction deals with the problem of recovering a function f from intensity measurements of the form {|Af (x)| : x ∈ X} where A is a linear transformation and X is a subset of the domain of Af .This is known as the phase retrieval problem, a non-linear inverse problem arising in numerous applications in science and engineering, such as coherent diffraction imaging [13,14,31], speech recognition [5] and quantum mechanics [12].In the situation where f belongs to an infinite-dimensional Banach or Hilbert space, phase retrieval is known to be never uniformly stable [2,6], hence constituting a challenging problem for numerical reconstruction approaches.See the surveys [15,20] for a current overview on uniqueness, stability and algorithms for phase retrieval.Recent attention was drawn to signals having a shift-invariant structure, that is, (1) f where ϕ ∈ L p (R) is the so-called generator of V p β (ϕ) and p ∈ [1, ∞].The constant β > 0 is a regularity parameter, commonly referred to as the stepsize.In a series of papers, Chen et.al. studied the phase retrievability of a real-valued map f from its modulus |f |, A = Id, assuming that the generator has compact support [7,8,9].Shenoy et.el.investigated the situation where A = F is the Fourier transform and the defining sequence (c n ) n of f obeys special properties [32].In the situation where ϕ = sinc is the cardinal sine, i.e. the considered signals lie in the classical Paley-Wiener space, Thakur demonstrated how to recover a real-valued bandlimited function from unsigned samples lying on a suitable dense sampling set [33].Similar results for real-valued maps were derived for Gaussian generators as well as totally positive generators of Gaussian type and are due to Gröchenig [19] and Romero [30].Besides, it was shown that if ϕ = φ σ is the Gaussian with variance σ 2 then, in general, complex-valued maps are not determined from samples of their absolute value.Hence, additional information is needed to achieve injectivity for complex signals from measurements |Af | with a suitably chosen signal-transformation A. The authors of the present article proved in [22] that every complex-valued map in V 1 β (φ) is uniquely determined up to a global phase factor from measurements of the form |V g f (X)| where V g is the short-time Fourier transform with Gaussian window g = φ σ , (3) and X ⊂ R 2 is a separated set of sampling points in the time-frequency plane.The resulting phase retrieval problem with transformation A = V g is known as the Gabor phase retrieval problem which plays an important role in applications such as ptychography, a popular and highly successful approach in coherent diffraction imaging.The previous observations serve as the starting point for a more penetrating investigation of the Gabor phase retrieval problem in a Gaussian shift-invariant setting and motivate the work in hand.In particular, we derive an explicit inversion formula and stability estimates as well as propose a new reconstruction algorithm which approximates complex-valued signals in V ∞ β (φ) from finitely many (deterministic) samples of |V g f | in a provable and stable manner.
Classical algorithmic approaches on solving the phase retrieval problem are based on fixed point iterations, iterative projection methods and gradient descent methods.Examples include the alternating projection method, the Gerchberg-Saxton algorithm [18], error reduction [17], the Douglas-Rachford algorithm [16], Hybrid Input-Output Algorithm [17] and the Averaged Alternating Reflection Algorithm [29,28].The difficulty of analyzing convergence and robustness of these methods lies in the non-convexity of the phase retrieval problem.Heuristically, these methods perform well if one chooses a right initialization which is close to a solution of the problem.However, it is in general unclear how to perform the initialization step and provide convergent results in a noise regime.Moreover, the discretization of these methods in an infinite-dimensional setting constitutes a difficult task.In contrast, the derived algorithm of the present paper overcomes these difficulties by providing a provably convergent approximation routine from noisy samples.Moreover, we quantify exactly how many samples of |V g f | need to be taken in order to guarantee a given reconstruction error.
1.1.Notation.Throughout this article, the shift-invariant space V p β (ϕ), the Gaussian φ σ and the short-time Fourier transform are defined as in ( 1), ( 2) and (3), respectively.The cross-ambiguity function of f and g is the map A(f, g)(x, ω) = e πixω V g f (x, ω).The spectrogram of f w.r.t. to a window g is given by the map (x, ω) → |V g f (x, ω)| 2 .If g = φ σ is a Gaussian then we set G := V φ := V φ σ and call G the Gabor transform of f .For clarity of the exposition and simplicity of notation, we will drop the dependence on σ.For a function f : R → C and an ω ∈ R we define its tensor product f ω via (4) where T ω f = f (• − ω) is the shift-operator.Moreover, the modulation operator M ω is defined by M ω f (t) = e 2πiωt f (t).In Section 2 we will exhibit that under suitable assumptions on ϕ the system of translates (T βn ϕ) n∈Z forms a Riesz basis for the space V 2 β (ϕ).Denoting by S the frame operator associated with (T βn ϕ) n∈Z and making use of the property that the operators S and T βn commute for every n ∈ Z implies that there exists a map φ such that (T βn φ) n∈Z is the dual frame of (T βn ϕ) n∈Z .The function φ given by (5) φ := S −1 ϕ is called the dual generator of ϕ.If (T βn ϕ) n∈Z forms a Riesz basis for V 2 β (ϕ) and f = n c n T βn ϕ then we call c = (c n ) n ⊂ C the defining sequence of f .In order to study stability properties of the Gabor phase retrieval problem in a Gaussian shift-invariant regime we define a mixed norm as follow.Let α > 0, p ∈ [1, ∞] and suppose that F : αZ × R → C is measurable.The mixed norm ∥ • ∥ α,p is defined by (6) ∥F ∥ α,p := This is the ℓ p -norm of the sequence (∥F (αn, •)∥ L 1 (R) ) n∈Z .If F is defined on R 2 we set ∥F ∥ α,p := ∥F | αZ×R ∥ α,p where F | αZ×R is the restriction of F to αZ × R. If not specified otherwise, all indices will run over the integers Z. Finally, we say that two functions f, h : R → C equal up to a global phase (or: up to a unimodular constant) if there exists a τ ∈ T := {z ∈ C : |z| = 1} such that f = τ h.

Main results.
The first main result, Theorem 1.1, provides an explicit inversion formula of a function f ∈ V ∞ β (φ) from its spectrogram |Gf (X)| where X = β 2 Z × R is a set of vertical lines in the time-frequency plane.Theorem 1.1.Let f ∈ V ∞ β (φ) and let p ∈ R such that f (p) ̸ = 0. Then there exists a unimodular constant τ ∈ T such that and φ ω is the dual generator of the Gaussian tensor product φ ω .
In the previous theorem the maps { φ ω : ω ∈ R} constitute a parametrized system of dual generators.In particular, the statement shows that two functions f, h ∈ V ∞ β (φ) agree up a global phase provided that their spectrograms agree on X = β 2 Z × R. Note, that this is a uniqueness statement within the signal class V ∞ β (φ).A similar statement does not hold if one replaces V ∞ β (φ) by L 2 (R), as we have shown in [21].
Next, we study stability properties of the above inversion.Heuristically, Theorem 1.2 assumes that the map f has the property that |f | is not too small on large intervals.More precisely, it will be assumed that there exists J ∈ N, p 1 , . . ., p J ∈ R and γ > 0 such that

Condition (P) excludes examples of the form f
where it is known that exponential instabilities arise for growing s [3].These examples further show that an assumption on |f | of the above form cannot be dropped.
and the implicit constant depends only on the step-size β and the standarddeviation σ of the Gaussian φ = φ σ .
We seek to discretize the previous regime in the sense that it yields an algorithmic approximation of functions in V ∞ β (φ) from finitely many measurements located on a grid of the form )×(2H+1) denotes a noise matrix then the given samples S take the form The contribution of the noise η will be quantified in terms of the ℓ ∞ -operator norm of η, i.e. the maximum absolute row sum of η.We denote this norm by ∥η∥ ∞ .The so-called numerical approximation routine R depends only on S and is defined as follows: for and phases ν 0 , . . ., ν J−1 ∈ T by Assuming that L j and ν j are well-defined, i.e. c j > 0 and L j (p j+1 − p j ) ̸ = 0 we define R : . This function has the property that it approximates f from finitely many noisy samples.Figure 2 illustrates the reconstruction performance of the numerical reconstruction routine R applied to finitely many noisy spectrogram samples of a complex-valued signal.The following theorem proves stability and convergence of this algorithm.
and that there exist and assume that the parameters h, H, N of the grid X have the property that The numerical approximation routine R combined with Theorem 1.3 results in a provably stable and convergent approximation algorithm as soon as condition (P) is satisfied.Note that this is an assumption on |f | which is a-priori unknown if only samples of |Gf | are available.However, this information is encoded in the spectrogram of f : we will show that the modulus |f | of any f ∈ V ∞ β (φ) can be approximated in a uniformly and globally stable way without imposing any additional assumption on f , see the discussion following Corollary 3.7, as well as Lemma 3.8.Based on a reconstruction of |f | the resulting algorithm presented in Section 4 first detects points p 1 , . . ., p J such that condition (P) is satisfied and finally uses the numerical approximation routine R to provide a reconstruction of f .In this way, the information on the location of points p j can be omitted.
Finally, we study the growth behavior of the number of spectrogram samples needed to approximate a function on an interval [−s, s].Assume that f has the property that there exists a partition for some γ > 0 and all j ∈ Z. Suppose we aim to approximate f on the interval [−s, s] using the numerical approximation routine R = R(s).How does the minimal number of spectrogram samples N (s) grow so that R(s) achieves the bound (9) for a fixed ε > 0?
for all j ∈ Z. Assume that the spectrogram samples S are noiseless, η = 0, and ε is given as in Theorem 1.3.Then the minimal number of spectrogram samples N (s) needed to achieve the bound (9) on the interval [−s, s] using the numerical approximation routine R satisfies Theorems 1.1, 1.2, 1.3 and 1.4 are a direct consequence of Theorems 3.6, 3.10, 4.6 and Corollary 4.7 that will be proved in Sections 3.1, 3.2, 4.2 and 4.2, respectively.1.3.Outline.The article is structured as follows: in Section 2 we present the necessary general theory on shift-invariant spaces which will be used throughout the article.Section 3 starts with an abstract discussion on the reconstruction of arbitrary maps from their tensor product.These ideas will be applied to Gaussian shift-invariant spaces and lead to an explicit inversion formula from spectrogram measurements which makes use of a biorthogonal expansion.The section ends with a stability analysis of the derived reconstruction formula.In the final section, Section 4, we present and analyze an algorithmic way of approximating functions in Gaussian shift-invariant spaces from finitely many spectrogram samples.Moreover, we elaborate that the obtained algorithm approximates signals in Paley-Wiener spaces.

Preliminaries on shift-invariant spaces
2.1.General theory.For β > 0, p ∈ [1, ∞] and ϕ ∈ L p (R) let the shiftinvariant space V p β (ϕ) be defined as in (1).Recall that a sequence (f n ) n ⊂ H in a separable Hilbert space H is said to be a Riesz basis for H if it is complete in H and there exist positive constants 0 < A ≤ B such that for arbitrary n ∈ N and arbitrary c 1 , . . ., c n ∈ C one has (10) A The optimal constants A and B are called the lower and upper bound of (f n ) n .Every Riesz basis with lower and upper bounds A and B is a frame with frame constants A and B. Associate to a Riesz basis (f n ) n its frame operator S : H → H, defined by Sf = n∈Z ⟨f, f n ⟩f n .Then S is a positive, invertible operator and ( Now let H = V 2 β (ϕ) be a shift-invariant space with a square-integrable generator ϕ.Define the 1-periodic map Φ β : R → R via (11) Φ where φ = Fϕ denotes the Fourier transform of ϕ which is defined on and extends to a unitary operator on L 2 (R).The map Φ β characterizes the property of (T βn ϕ) n being a Riesz basis for V If ϕ satisfies condition (12) and if S : is the frame operator associated with (T βn ϕ) n then S commutes with βZ-shifts, implying that (S −1 T βn ϕ) n = (T βn φ) n with φ = S −1 ϕ [10, Lemma 9.4.1].Accordingly, the canonical dual frame (S −1 T βn ϕ) n has the same structure as (T βn ϕ) n .Definition 2.2.Let β > 0 and ϕ ∈ L 2 (R) such that (T βn ϕ) n is a Riesz basis for V 2 β (ϕ).If S is the frame operator associated with (T βn ϕ) n then ϕ := S −1 ϕ is called the dual generator of ϕ.
Clearly, the dual generator generates the same shift-invariant space, i.e.
In the situation of Theorem 2.1, the map ϕ can be derived explicitly.To do so, define in a similar fashion as above the The dual generator is then given as a formula of φ and Ψ.
Proof.This follows from a generalization of [10,Proposition 9.4.2].For the convenience of the reader, we provide a proof in the Appendix 5.1.□ Note that the structure of V 2 β (ϕ) can be enriched by imposing further properties on the generator ϕ.For instance, if ϕ is continuous and belongs to the Wiener Amalgam space then V 2 β (ϕ) is a reproducing kernel Hilbert space (RKHS).We refer to [4,10,11] for results in this direction.

Gaussian generators.
We now specialize ϕ to be the Gaussian φ σ with variance σ 2 as defined in Equation (2).For the sake of simplicity we write φ = φ σ with the implicit understanding that φ has variance σ 2 .Properties of V 2 β (φ) can be expressed in terms of the Jacobi theta function of third kind.Recall that this is the map ϑ 3 : C × (0, 1) → C, For every fixed c ∈ (0, 1), ϑ 3 is an entire, π-periodic map in z.Using the results of the previous subsection gives Corollary 2.4.Let φ be a centred Gaussian with variance σ 2 .Then the 1-periodization of φ is given by )).In particular, the system (T βn φ) n constitutes a Riesz basis for V 2 β (φ).Moreover, if the function Λ : R → C is defined as Proof.The identity where the second equality follows from Poisson's summation formula.In a similar fashion, we obtain )), thereby proving the statement using Theorem 2. with c = φ( β √ 2 ) then the identity F(s → e −2π 2 σ 2 s 2 )(t) = (σ √ 2π) −1 φ(t) implies that φ can be written as Hence, up to a multiplicative constant, the defining sequence of the dual generator is given by the Fourier coefficients of the reciprocal theta function and therefore the decay of φ matches with the decay of (a n ) n .In fact, the Fourier coefficients of 1 ϑ 3 satisfy precisely an exponential decay (see a derivation by Janssen [27, p. 178] based on results of Whittaker and Watson [35, p. 489]): if the constant ξ = ξ(c) is defined by 2 ) .
This explicit formula for the Fourier coefficients will be used in the upcoming sections to derive decay estimates on the dual generator leading to stability and approximation results of the Gabor phase retrieval problem.We conclude the present section with the observation that if f has a Gaussian shift-invariant structure and is merely assumed to be bounded then the inner products ⟨f, T βn φ⟩ := R f (t)T βn φ(t) dt remain well-defined and the map where the series on the right converges uniformly on compact intervals of the real line.
Proof.See the Appendix 5.

□
We call n ⟨f, T βn φ⟩T βn φ the biorthogonal expansion of f ∈ V ∞ β (φ).Note that V ∞ β (φ) consists of smooth functions and is the largest among the spaces V p β (φ) due to the inclusion

Uniqueness and explicit reconstruction
Gf the Gabor transform of f .We aim to recover f from |V g f (X)| where X ⊆ R 2 .Clearly, this will only be possible up to a global phase.In order to decide whether two functions agree up to global phase one can study their tensor product f ω as defined in (4).

Proposition 3.1. Let f, h : R → C be two complex-valued maps and denote by F, H
(1) If p ∈ R such that F (p, 0) ̸ = 0 then for every ω ∈ R we have Then there exists a τ ∈ T with the following property: Whenever t ∈ [p 1 , p J ] with t = p j + ω where p j is the largest value not exceeding t then In particular, f can be reconstructed on the interval [p 1 , p J ] from the values F (X) where X is the union of the diagonal line segments The first statement then follows directly from the equations The second statement follows directly from the first and the third follows from the second.It remains to show Part 4. Let L j and ν j be defined as above.We have by definition Therefore, |f (p 1 )| which further shows that f can be reconstructed on [p 1 , p J ] up a global phase from the union of the diagonal line segments D j .□ Part 4 of the preceding proposition can be interpreted as follows: if the tensor product is known on the diagonal line segments D j , 1 ≤ j ≤ J − 1, then f can be reconstructed on the interval [p j , p j+1 ] up to a global phase and this reconstruction is given by L j .If on each such interval the local reconstruction is multiplied by the phase τ ν 1 ν 2 • • • ν j−1 then the resulting function defined on the union of all the intervals [p j , p j+1 ] agrees up to the phase factor τ with f .We call the multiplication by ν 1 ν 2 • • • ν j−1 a phase synchronization.The phase synchronization is possible whenever f (p j ) ̸ = 0.This motivates Definition 3.2.Let p 1 , . . ., p J ∈ R and γ > 0. We say that f : R → C satisfies condition (P) if Suppose now that f is an element of the Gaussian shift-invariant space V ∞ β (φ) where β > 0 is an arbitrary step-size.An important property of V ∞ β (φ) is that the shift-invariant structure is invariant under the tensor product operation.Proposition 3.3.Let φ = φ σ be a Gaussian window function with variance σ 2 and let V ∞ β (φ) be the corresponding shift-invariant space with generator φ and step-size The tensor product of f is given by Using the product formula shows that f ω takes the form The statement follows at once from the bound □ We aim to reconstruct f from its tensor-product f ω using the abstract reconstruction result given in Proposition 3.1 (4).Since f ω is generated by φ ω one could recover f ω from the inner products ⟨f ω , T β 2 n φ ω ⟩, n ∈ Z, using a biorthogonal expansion as discussed in Section 2.2.We already verified the Riesz basis property of (T β 2 n φ ω ) n in Corollary 2. 4. It remains open how to gain access to the inner products ⟨f ω , T β 2 n φ ω ⟩ ∈ C.This information is encoded in the spectrogram of f .Proposition 3.4.For every f ∈ L ∞ (R) and every g ∈ L 2 (R) one has Proof.Under the above assumptions we have f T x g ∈ L 2 (R) and therefore The convolution F 2 p x * p x can be written as where we used the change of variables k → k − x 2 in the first equality.It follows that which yields the assertion.□ Denoting the dual-generator of φ ω by φ ω , the following result is an immediate consequence of Proposition 3.4.
is the biorthogonal expansion of f ω w.r.t. the dual generator φ ω then for every n ∈ Z.In particular, every f ω is determined uniquely by |Gf 20) is a direct consequence of Proposition 3.3 and Proposition 3.4 above.In particular, this implies that f ω is determined uniquely by |Gf ( β 2 Z × R)|.Since ω was arbitrary, Proposition 3.1 shows that f is determined up to a global phase factor by |Gf f is real-valued, then the Gabor phase retrieval problem can be interpreted as a reconstruction of f from its modulus (up to a variance-change in the generator).For if the frequency variable in the Gabor transform is set to zero then ( 21) ) and a result due to Gröchenig shows that h is determined by |h(X)| whenever X ⊆ R has lower Beurling density This result is sharp and does not hold for complex-valued maps [19, Section 2].Corollary 3.5 shows that in the complex case, (c n ) n ⊂ C, uniqueness can be guaranteed if we set ε = 0 and extend the grid β 2 Z×{0} to parallel lines β 2 Z × R in the time-frequency plane.In order to transform the uniqueness part of Corollary 3.5 into a uniqueness result from measurements lying on a lattice, one may suppose further properties on the step-size β as in [22,Theorem 3.6].We conclude the present section with an explicit reconstruction formula from spectrogram measurements.
Then there exists a unimodular constant τ ∈ T such that for every ω ∈ R and Proof.Combining Corollary 3.5 with Proposition 3.1 and fact that the dual generator φ ω is real-valued yields the statement.□ Theorem 3.6 implies that in a Gaussian shift-invariant setting it suffices to know the spectrogram on parallel vertical lines in order to reconstruct f from its tensor product.In our further analysis it will be useful to work with an explicit formula for the dual generator of the tensor product φ ω .Such a formula can be readily established.The product formula .

3.2.
Stability.This section is devoted to the stability analysis of the biorthogonal expansion of the tensor product f ω as well as the explicit inversion formula presented in Theorem 3.6.Let the mixed norm ∥ • ∥ α,p be defined as in equation (6).Using this norm, stability estimates for tensor products of functions in where Proof.By Hölder's inequality we have Combining the previous two inequalities yields the first part of the statement.To obtain the second assertion, we observe that the 1-periodization of ).According to Theorem 2.1, the value A(σ, β, ω) is then given as the minimum of the map It was shown by Janssen [27, p. 178] that for every c ∈ (0, 1) the theta function ϑ 3 (•, c) attains its minima at t ∈ π 2 Z, thereby proving the statement. □ By virtue of the RKHS structure it is evident that an analogue version of Corollary 3.7 can be derived where the L 2 -norm on the left-hand side is replaced by the (φ ω ) and (φ ω ) is a RKHS.Therefore, the point evaluation functionals (φ ω ) and there exists Combining the fact that the generator is a Gaussian with an explicit form of K x (see, for instance, [4, Theorem 4.1]), we observe that is bounded on R which in turn yields the bound for some constant c > 0.Moreover, for a fixed ω ∈ R, Corollary 3.7 shows that the tensor product is stably determined by spectrogram measurements and the stability constant it follows that the stability constant for reconstructing the modulus f 0 = |f | 2 is the best among the constants {A(σ, β, ω) : ω ∈ R} and the stability deteriorates rapidly in ω.Classical examples for instabilities are of the form and they show that the stability constant grows like e n 2 [3].Thus, one cannot expect to stably reconstruct arbitrary non-positive functions without imposing further assumptions.We surpass instabilities arising from signals of the form (24) [1,23,24].Such situations can be excluded by the requirement that f satisfies condition (P): and γ > 0 with |f (p j )| ≥ γ for all j then the distance of the areas where |f | is concentrated is determined by the value r = max 1≤j≤J−1 p j+1 − p j .In this setting we expect that the quotient distance min τ ∈T ∥f − τ g∥ can be controlled by the norm of the difference of the corresponding spectrograms up to factor which depends on e r 2 and 1 γ .This is indeed the case as we shall elaborate in the following.First, we define a constant C(σ, β) depending only on the variance of the generator as well as the step-size of the underlying Gaussian shift-invariant space by (25) C(σ, β) := sup where Λ is the map which characterizes the dual generator, see equation (22).We start with local stability estimates around points p j ∈ R, 1 ≤ j ≤ J, J ∈ N where f (p j ) ̸ = 0 for all j.
Proof.In view of Corollary 3.5, both maps |f | 2 and |g| 2 can be expanded in terms of the dual generator of φ 0 = φ 2 .Hence, we have for every t ∈ R the estimate g(p j ) ̸ = 0 the explicit reconstruction formula from tensor products shows that The choice of τ j implies that Setting A := |f ω (p j + ω) − g ω (p j + ω)| and employing the reconstruction formula from Theorem 3.6 and the explicit form of the dual generator gives The assertion follows from that.□ For j ∈ N, let τ j , I j and c j be defined analogue to Lemma 3.8.The following elementary Lemma will be useful in the proof of the main stability result of the present section.
Suppose that c j , τ j , τ m are defined as in Lemma 3.8.If we define then we have the implication Using the definition of τ j we can re-write the left-hand side of the previous inequality as Since ) we obtain with Lemma 3.8 .
and for a fixed j ∈ {1, . . ., J} we have By the assumption on the distance between the p j 's and the definition of the I j 's we have p j ∈ I j−1 for every j ∈ {2, . . ., J}.Therefore, Lemma 3.9 shows that where B is defined as in (26).Using standard estimates it follows that, which implies in combination with equation ( 28) that 2C(σ, β) then clearly for every τ ∈ T we have the estimate The statement follows from case 1 and case 2. □ Additional assumptions on σ and β yield explicit upper bounds on the constant C = C(σ, β).To achieve such bounds we first observe that according to Section 2.2 the map F −1 Λ satisfies an exponential decay.Thus, there are constants In an exemplary way we can derive a result of the following form which will result in explicit upper bounds on C(σ, β).
Remark 3.12 (On the dependence on ∥f ∥ L ∞ (I) and ∥g∥ L ∞ (I) in the stability estimte ( 27)).Shift-invariant spaces , generated by a Gaussian with standard deviation u > 0 satisfy the norm equivalence

Moreover, recent sampling results for shift-invariant spaces with totally positive generator of Gaussian type show that
with a universal constant depending only on σ and β.With this choice of a sampling set Ω, the right-hand side of the stability estimate (27) can be made dependent only on γ, σ, β, |Gf 4. Algorithmic approximation in V ∞ β (φ) In the final section of the present article we shall combine our previous results to derive an algorithmic approximation of functions f ∈ V ∞ β (φ) from finitely many spectrogram samples |Gf (X)|, X ⊂ R 2 , |X| < ∞.Assume that the samples X are given on a finite grid of the form (30) X = β 2 {−N, . . ., N } × h{−H, . . ., H} where N, H ∈ N and h, β > 0. We call 1 h the sampling density of X in frequency direction.If η ∈ R (2N +1)×(2H+1) denotes a noise matrix then the given sampling set S takes the form To quantify the contribution of the noise η we use the ℓ ∞ -operator norm of η, which is the maximum absolute row sum of η.Our only assumption on the noise η is that we have control over the norm ∥η∥ ∞ ≤ δ for some noise-level δ > 0. Otherwise, the error is arbitrary and can be adversarial.
local reconstructive functions L j by and phases ν 0 , . . ., ν J−1 ∈ T by Note that c j , L j and ν j depend only on the sampling set S. Assuming that L j and ν j are well-defined, i.e. c j > 0 and L j (p j+1 − p j ) ̸ = 0, gives rise to the definition of a numerical approximation routine.
Definition 4.1.Suppose that f : R → C is a measurable and bounded map which satisfies condition (P) with γ > 0 and Let the grid X be defined as in (30) and suppose that L j and ν j are well-defined.The numerical approximation routine R : In order to derive approximation results by means of the numerical approximation routine R, we define for j = 1, . . ., J − 1 the error term This error term allows us to state conditions under which the functions L j and the phases ν j are well-defined.More importantly, it yields an upper bound on the quotient distance of f and R.
Step 1: Upper bounding min τ ∈T ∥f −τ R∥ L ∞ (I) .Assume for the timebeing that the phases ν j and the constants c j are well-defined.Let t ∈ I such that t ∈ (p j , p j+1 ] with t = p j + ω and 2 ≤ j ≤ J − 1.Then The term A can be bounded by where the third inequality follows by induction.Each summand in the last term of the sum above satisfies Combining the inequalties ( 32), ( 33) and (34) gives Observing that Step 2. Well-definedness of c j .If the map S j is defined by .
By assumption we have In particular, the constants c j are well-defined.
Step 3: Upper bounding E j (ω).We have shown that the constants c j are well-defined.Recalling that if f (p j ) ̸ = 0 then the identity 1 holds for every ω ∈ R. We can continue with the following estimate: Using the bound obtained in equation (36) together with the fact that Step 4: Well-definedness of the phases ν j .We use the upper bound on E j to show that the phases ν j are well-defined.First, observe that it follows it follows from the upper bound obtained in (37) that |f (p j+1 )− In particular, this shows that (39) Thus, the phases ν j are well-defined for all j.
If Q is given as above then from inequality (38) we obtain :=Ξ Inequality ( 35) and the lower bound 39 implies that min In particular, the assumption on ε and the fact that min{γ, γ 2 , γ 4 , γ 5 } = min{γ, γ 5 } and max{1, ∥f Application to Gaussian shift-invariant spaces.Recall that the inner products ⟨f ω , T x φ ω ⟩ were given in terms of the Fourier integral To allow an approximation of functions in f ∈ V ∞ β (φ) from |Gf (X)| where X is a set of finitely many sampling points, we discretize the above integral using a suitable quadrature formula.If σ is the standard deviation of the Gaussian φ then we define σ ′ := 1 2πσ .The choice of a quadrature formula is based on the following factorization formula of the spectrogram.
where S x (t) = ℓ b ℓ (x)e πiβℓt is a trigonometric series whose coefficients satisfy the Gaussian bound In particular, for every x ∈ R the map z → |Gf (x, z)| 2 extends from R to an entire function on C. The extension is given by the map C ∋ z → πσ 2 φ σ ′ (z)S x (z).
Proof.See the Appendix 5.4.□ Recall that for h > 0 and H ∈ N∪{∞} the trapezoidal rule approximation of a map W : R → C is defined by ( 42) Applying this quadrature rule to the Fourier integral (40) gives ( 43) The quantity 2M in the numerator is optimal.
Truncating the series which corresponds to the biorthogonal expansion of f ω and using the approximation (43) of the inner products given in (40), results in the map ( 44) which is precisely the function √ c j L j (with p = p j ) appearing in the error term E j defined in equation (31).The pointwise distance between the map ω → f ω (p + ω) and S can be upper bounded as follows.
⌉ + m for some m ∈ N and some r, s > 0. Let h > 0 and H ∈ N. Then for every p ∈ [−s, s] and every ω ∈ [0, r] we have the error bound where K = K(σ, β) and ν = ν(σ, β) denote the decay constants of the dual generator as defined in (29).
Proof.See the Appendix 5.5.□ Lemma 4.5 becomes meaningful if we interpret it in a qualitative way: recall that the constants K and ν depend only on σ and β.Suppose that f satisfies condition (P) and that the variance σ 2 and the step-size β are fixed, i.e. both the window function and the underlying signal space is fixed.Under these assumption, it follows from inequality (45) and the definition of the error terms E j that (46) where a and b are constants depending only on σ and β.The constant r is given by Let the universal constant in inequality (46) be denoted by D = D(σ, β, r).
Then the following approximation result holds.
. Suppose that the sampling density and the measurements are given on the grid

. , H} such that the grid size characterizing parameters N, H satisfy
and Then the following holds: if the noise level does not exceed where R denotes the numerical approximation routine.
Proof.Denote the four summands on the right-hand side of inequality (46) by Ξ 1 , . . ., Ξ 4 .It follows that if Elementary manipulations show that the condition DΞ i ≤ ε 4(J−1) for i = 1, . . ., 4 is satisfied by the above assumptions on h, H, N and ∥η∥ ∞ .The statement is therefore a consequence of Theorem 4.2.□ Based on the previous statement we can study the growth of the number of samples needed to achieve the bound (48) for an increasing length of the reconstruction interval I = [−s, s].Consider the noiseless case, η = 0. Suppose that condition (P) holds on the entire real-line, i.e. there exists a partition so that the maximal distance between two consecutive points p j and p j+1 remains upper bounded by r, i.e. sup j∈Z (p j+1 − p j ) ≤ r.By going to a sub-partition we can always assume that for an interval [−s, s] there exists p j , p k ∈ P with p j ≤ −s, p k ≥ s and k − j ≤ 4s r + O(1).Therefore, the assumptions on the bounds on h, H and N as stated in Theorem 4.6 can be re-written as where D ′ is a constant depending only on σ, β and r.A direct consequence is Corollary 4.7.Under the assumptions above, the following holds: in order to achieve a reconstruction of f ∈ V ∞ β (φ) with defining sequence c ∈ ℓ ∞ (Z) on an interval I = [−s, s] of length 2s up an error bound of the form (48) using the numerical approximation routine R, the number N (s) of spectrogram samples satisfies the upper growth estimate Assumptions on σ and β turn the qualitative results above into quantitative results by means of explicit upper bounds on a, b and D. This can be done in a similar fashion as at the end of Section 3 and exploiting Lemma 3.11.
Proof.This follows by applying the bounds on K and ν derived in Lemma 3.11 to each term on the right-hand side of inequality (45).□ 4.3.Algorithm.In a natural way, the foregoing results lead to a reconstruction algorithm whenever f ∈ V ∞ β (φ) satisfies condition (P).This is an assumption on the signal itself which is unknown if only spectrogram samples are accessible.However, we have shown that the map f 0 = |f | 2 can be recovered in a globally stably way by means of a biorthogonal expansion.Let the grid X and the sampling set S be given as above and set uniformly on compact intervals.This motivates the following 4-step approximation procedure.
) for all j ∈ {1, . . ., J} (3) Define for j ∈ {1, . . ., J} local reconstructive functions L j by and phases ν j by Suppose that the samples are given on the grid X with Since p was arbitrary, this estimate holds for every p ∈ [−s, s].Let p j ∈ I be one of the detected points.Then F (p j ) ≥ γ = 3 2 γ 2 and therefore Since j ∈ {1, . . ., J} was arbitrary and J ≥ 2, we conclude from the inequality above that the points p 1 < • • • < p J satisfy condition (P). (

Select all points t
and do this as long as there are no points to remove anymore.Call the resulting points u i 1 < • • • < u i K .Let i k be the first index such that u i k+1 − u i k ≤ r and J ∈ N the largest value such that u i k+ℓ − u i k+ℓ−1 ≤ r for all ℓ = 1, . . ., J. Finally, set p ℓ = u i k+ℓ−1 , ℓ = 1, . . ., J. Remark 4.11.Suppose for simplicity that β, γ, r ≈ 1 and that the noise is Gaussian of the form N (0, σ 2 ).Empirical simulations highlight that if the samples are taken on a grid which satisfies the conditions of Theorem 4.9 then Algorithm 1 reconstructs functions in V ∞ β (φ) with almost no visible error provided that the noise level satisfies σ ≤ 0.01.If the noise level exceeds 0.01 then the algorithm usually terminates earlier since no point p ∈ R with |f (p)| ≥ γ could be found.
. Compactly supported signals.Algorithm 1 was designed for functions in shift-invariant spaces with Gaussian generator.One could inquire about the justification of applying Algorithm 1 to different signal classes such as compactly supported maps.To give a meaningful answer to this question we define for a, h > 0 the map This is essentially the local reconstructive function defined in Section 4.1 with (noiseless) samples given on the infinite lattice X = β 2 Z × hZ.Proof.General frame theory implies that the orthogonal projection P ω (g) of a map g ∈ L 2 (R) onto V 2 β 2 (φ ω ) is given by To visualize the approximation, we multiplied R with a global phase (the phase of f at p 1 , where p 1 arises from step 2 of Algorithm 1).
Hence, the map t → |Gf (x, t)| 2 is band-limited with bandwidth at most 4a.The assumption on h implies that h ≤ 1 4a and Shannon's sampling theorem [26,Theorem 6.13] Setting x = β 2 n and applying the Fourier transform in the second argument of the spectrogram results in and this proves the statement.□ This result shows that if we run Algorithm 1 with a compactly supported function rather than a function in V ∞ β (φ) then in each step the algorithm approximates the projection of f ω onto V 2 β 2 (φ ω ).Heuristically, compactly supported functions that can be well-approximated by a linear combination of Gaussians and which satisfy condition (P) can be reconstructed with this methodology.In Figure 7 we visualize the performance of Algorithm 1 assuming that f is a continuous function with compact support (but is not differentiable).
Band-limited signals.Suppose now that f has compact support in the frequency domain, i.e. f is a band-limited signal belonging to the Paley-Wiener space a (R) from noisy samples as soon as the standard deviation σ of the Gaussian window φ lies below a constant depending on the bandwidth a of f .This follows from two observations.Firstly, f is generated by a sinc function (Shannon's sampling theorem) and a sinc can be well-approximated by a linear combination of equally spaced Gaussians with standard deviation σ ≈ 1 a , see Figure 8 for a visualization.Secondly, if f has bandwidth a then f ω has bandwidth 2a for every ω ∈ R. Hence, P W 2 a (R) exhibits a similar invariance of the tensor product operation as in the Gaussian shift-invariant setting (Proposition 3.3).Since φ ω is (up to a constant) a Gaussian with standard deviation σ √ 2 for every ω ∈ R, we observe that under the assumption σ ≈ 1 a every f ω can be well-approximated by projecting f ω onto the Gaussian shift-invariant space V 2 β 2 (φ ω ).Clearly, general frame theory shows that the projection of (φ ω ) can be derived as soon as the inner products ⟨f ω , T x φ ω ⟩ are available.Since (Proposition 3.4), the inner products can be recovered up to a small error using a suitable quadrature rule.In case of a band-limited signal f ∈ P W 2 a (R) we further observe that .The spectrogram samples are taken on the grid X = 0.04{50, . . ., 50} × β 2 {−120, . . ., 120}.To each spectrogram sample Gaussian noise N (0, σ 2 ) is added with σ = 0.0001.To visualize the approximation, we multiplied R with a global phase (the phase of f at p 1 , where p 1 arises from step 2 of Algorithm 1).
where we used the convolution theorem for the Fourier transform.Thus, the map t → |Gf (x, t)| satisfies a Gaussian decay as in Theorem 4.3.Figure 9 depicts the output of Algorithm 1 for an input function f ∈ P W 2 a (R) in a Gaussian noise regime.
The foregoing discussion highlights that a combination of approximation properties of Gaussian shift-invariant spaces, numerical integration theory and the abstract theory presented in Section 4.1 point towards new approximation results from phaseless samples for function spaces beyond V ∞ β (φ).We leave this research direction open for future work.The inner products ⟨ θ, M −kβ φ⟩ can be written as where If φ(t) ̸ = 0 then Ψ β (t) ̸ = 0 and therefore we have 1 D φ = φ which implies that FSθ = 1 β φ. □ 5.2.Proof of Proposition 2.5.
be the Fourier coefficients of the reciprocal theta function as given in equation ( 16) and ( 17).Then .
The modulus of the coefficients |a n | is upper bounded by It follows that where the right-hand side is even in t.In addition, the series on the right of the previous inequality satisfies .
If t ̸ = 0 then the term A(t) can be upper bounded via where in the rightmost inequality we used the assumption that β 4 ≤ σ ≤  Step 2: Lower bounding ξ(c).For every 0 < c < 1 we have Proof.
For t, y, x, ω ∈ R let W be defined by W (t + iy) := |Gf (x, t + iy)| 2 e −2πiω(t+iy) .From this we see that inside the strip U the estimate holds and the theorem of Trefethen and Weidemann, Theorem 4.4, implies that Step 2: Estimating the cut-off error.For the integrand W as defined above and an H ∈ N we estimate the cut-off error |I ∞ h (W ) − I H h (W )|.For every ω ∈ R we have Using a classical Gaussian tail bound of the form ∞ q e −pt 2 dt ≤ π 4p e −pq 2 , p, q > 0 with p = 2π 2 σ 2 and q = Hh results in and the triangle inequality implies that

Theorem 2 . 3 .
Suppose that (T βn ϕ) n is a Riesz basis for V 2 β (ϕ) with frame operator S. Let D = {t ∈ R : Ψ β (t) ̸ = 0} and let 1 D be the indicator function of D. If the map θ is defined as 3. □Note that the zeros of ϑ 3 are bounded away from the real axis.Hence, the reciprocal theta function 1 ϑ 3 defines an analytic and periodic function in an open strip containing the real axis.It follows that the Fourier coefficients a n = a n (c) of 1 ϑ 3 decay at least exponentially.If (15) 1 ϑ 3 (πβt, c) = n∈Z a n e 2πiβnt
in the following way: A closer look on the map f := f n,+ shows that on the interval I = [−βn, βn] it holds that |f n,+ (βn)| = |f n,+ (−βn)| ≈ 1 whereas |f | ≈ 0 on a large part of the interior of I. Therefore, |f | (and also |Gf |) is concentrated on two almost disjoint components which is a classical indicator for instabilities

Figure 4 .
Figure 4.The contour plot depicts the spectrogram of a function f ∈ V ∞ β (φ) and the white dots visualize the location of the samples S. The values 2 β and 1 h are the densities in time (x) resp.frequency (ω) direction.

4. 1 .
The numerical approximation routine R. We start by presenting our reconstruction approach in an abstract setting, leaving the ansatz open for extensions to other signal classes (see Section 4.4).Assume for the timebeing that the function f has only the property of being measurable and bounded.Let p 1 < • • • < p J and γ > 0 such that f satisfies condition (P).For simplicity we assume that −p 1 = p J =: s and set I = [−s, s].The aim is to recover f on the interval I. Generalizations to other intervals are straightforward and can be achieved by a simple translation.Motivated by Proposition 3.1 and Theorem 3.6 we define constants c j by

Theorem 4 . 4 (
hk)| 2 e −2πiωhk , which essentially represents the discrete Fourier transform.We apply a classical result on the exponentially convergent trapezoidal rule [34, Theorem 5.1].Trefethen-Weideman).Let a > 0 and suppose that W is an analytic map in the strip U a = {z ∈ C : |Im(z)| < a}.Suppose further that W (x) → 0 uniformly as |x| → ∞ in the strip.If M (W ) := sup t+iy∈Ua ∞ −∞ |W (t + iy)| dt < ∞ then, for any h > 0, the trapezoidal rule approximation I ∞ h (W ) of R W exists and satisfies

( 1 )
2 s εr and D, a, b are defined as above.If the noise level satisfies ∥η∥ ∞ ≤ εr 16hDs then the following holds.Let p 1 < • • • < p J ∈ I be the points detected in Step 2 of Algorithm 1 with input parameters r = r, γ = 3 2 γ 2 , s = s.Then p 1 < • • • < p J satisfy condition (P) with constants r and γ (2) If p 1 < • • • < p J ∈ I are the points detected in Step 2 of Algorithm 1 with input parameters r = r, γ = 3 2 γ 2 , s = s and if R is the output function of Algorithm 1 then

2 ∞ 1 e 1 σh − 1
s, s] are the points detected in Step 2 of Algorithm 1 then the condition p j+2 − p j ≥ r implies that 4s r + 1 ≥ J and therefore J −1 ≤ 4s r .Let p ∈ [−s, s], let S be defined as in equation (44) and let F be the function defined in Step 2 of Algorithm 1. Then S(0) = F (p). Combining the inequalities (45) and (46) with the assumptions on h, H, N and arguing in an analogue fashion as in the proof of Theorem 4.6 yields ||f (p)| 2 − F (p)| = |f 0 (p) − S(0)| ≤ D ∥c∥ 2 ∞ e −am + ∥c∥

β 2 ≤ 1 .
By making use of the elementary bound e −bt 2 ≤ e b 4 e −b|t| (t, b ∈ R) and the condition on σ and β we deduce the chain of inequalities