Accuracy of Algebraic Fourier Reconstruction for Shifts of Several Signals
PPreprint submitted toSAMPLING THEORY IN SIGNAL AND IMAGE PROCESSING DVI file produced onMay 20, 2019
Accuracy of Algebraic Fourier Reconstruction forShifts of Several Signals
Dmitry Batenkov † (cid:93) , Niv Sarig † , Yosef Yomdin † { dima.batenkov,niv.sarig,yosef.yomdin } @weizmann.ac.ilThis research was supported by the Adams Fellowship Program of the Israeli Academyof Sciences and Humanities, ISF Grant No. 779/13, and by the Minerva foundation. † Department of Mathematics, Weizmann Institute of Science, Rehovot 76100, Israel (cid:93)
Department of Mathematics, Ben-Gurion University of the Negev, Beer-Sheva 84105,Israel
Abstract
We consider the problem of “algebraic reconstruction” of linear combi-nations of shifts of several known signals f , . . . , f k from the Fourier sam-ples. Following [5], for each j = 1 , . . . , k we choose sampling set S j to be asubset of the common set of zeroes of the Fourier transforms F ( f (cid:96) ) , (cid:96) (cid:54) = j ,on which F ( f j ) (cid:54) = 0. It was shown in [5] that in this way the reconstructionsystem is “decoupled” into k separate systems, each including only one ofthe signals f j . The resulting systems are of a “generalized Prony” form.However, the sampling sets as above may be non-uniform/not “denseenough” to allow for a unique reconstruction of the shifts and amplitudes.In the present paper we study uniqueness and robustness of non-uniformFourier sampling of signals as above, investigating sampling of exponentialpolynomials with purely imaginary exponents. As the main tool we applya well-known result in Harmonic Analysis: the Tur´an-Nazarov inequality([18]), and its generalization to discrete sets, obtained in [12]. We illustrateour general approach with examples, and provide some simulation results. Key words and phrases : non-uniform sampling, Tur´an-Nazarov inequality,exponential fitting, Prony systems — 94A20, 65T40 a r X i v : . [ m a t h . C A ] A p r D. BATENKOV AND N. SARIG AND Y. YOMDIN
In this paper we investigate robustness of Fourier reconstruction of signals ofthe following a priori known form: F ( x ) = k (cid:88) j =1 q j (cid:88) q =1 a jq f j ( x − x jq ) , (1.1)with a jq ∈ C , x jq ∈ R . We assume that the signals f , . . . , f k : R → R are known(in particular, their Fourier transforms F ( f j ) are known), while a jq , x jq are theunknown signal parameters, which we want to find from Fourier samples of F .Practical importance of signals of the form (1.1) is well recognized in the lit-erature. For instance, they appear in digital processing of neuronal signals,bioimaging, image processing and ultrawideband communications [11, 13]. Theyare of relevance also in inverse moment problems, an important subject in math-ematical physics [14, 21].We follow a general line of the “Algebraic Sampling” approach (see [6, 23, 11] andreferences therein), i.e. we reconstruct the values of the unknown parameters,solving a system of non-linear equations, imposed by the measurements. Theequations in this system appear as we equate the “symbolic” expressions of theFourier samples, obtained from (1.1), to their actual measured values.Our specific strategy, as suggested in [5, 24], is as follows: we choose a samplingset S j ⊂ R , j = 1 , . . . , k, in a special way, in order to “decouple” the reconstruc-tion system, and to reduce it to k separate systems, each including only one ofthe signals f j . To achieve this goal we take S j to be a subset of the sets W j ofcommon zeroes of the Fourier transforms F ( f (cid:96) ) , (cid:96) (cid:54) = j . It was shown in [5] thatthe decoupled systems turn out to be exactly the same as those which appearin the fitting of exponential polynomials on sets S j (systems (2.2) in Section 2below).In this paper we restrict ourselves to one-dimensional case. A presentationof the Fourier Decoupling method in several variables, as well as some initialuniqueness results, can be found in [24, 5]. On the other hand, we explicitlyassume here that k ≥
2. So the usual methods which allow one to solve thisproblem “in closed form” in the case of shifts of a single function (see [11, 4, 24])are not directly applicable. Still, as it was shown in [5], in many cases anexplicit reconstruction from a relatively small collection of Fourier samples of F is possible. Let us also stress that the decoupling method of [24, 5] in dimensionone can “generically” be applied only to the shifts of at most two different signals. ccuracy of Algebraic Fourier Reconstruction for Shifts of Several Signals S j are the intersections ofat least two different discrete sets (the sets of zeroes of the Fourier transforms F ( f (cid:96) ) , (cid:96) (cid:54) = j ), so “generically” S j are empty. However, in many important“non-generic” situations of k > k .If the points s j(cid:96) ∈ S j , (cid:96) = 1 , , . . . , form an arithmetic progression, the recon-struction systems (2.2) are very closely related to the standard Prony system(see, for instance, [7] and discussion therein). However, the sampling sets S j ,being subsets of the sets W j of zeroes of the Fourier transforms F ( f (cid:96) ) , (cid:96) (cid:54) = j ,are completely defined by the original signals f (cid:96) , and cannot be altered in orderto make sampling more stable. These sets usually are non-uniform, thereforethe standard methods for robust solution of Prony systems cannot be applied.Even if S j forms an arithmetic progression, it may turn out to be “insufficientlydense” to allow a robust reconstruction of the shifts and amplitudes (see anexample in Section 3 below). Because of these reasons, we restrict ourselves toonly one solution method for system (2.2) - that of the least squares fitting,mainly because of its relative insensitivity to the specific geometry of the sam-pling set. Accordingly, we do not consider in this paper other approaches, whichcan be more efficient in certain specific circumstances. Let us only mention thatnon-uniform sampling is an active area of research, see e.g. [1, 17] and referencestherein.The main goal of the present paper is to study uniqueness and robustness of theFourier decoupling method. We define a “metric span” ω ( S ) of sampling sets S ,which is a simple geometric quantity, taking into account both the geometry of S , as well as the maximal shifts allowed in the signal F (which are the maximalfrequencies of the exponential polynomials appearing in the Fourier transformof F ). Our main results - Theorem 2.1 and Corollary 2.1 below - provide, interms of the metric span ω a “density-like” geometric condition on the commonsets W j of zeroes of the Fourier transforms F ( f (cid:96) ) , (cid:96) (cid:54) = j , which, in the case of nonoise, guarantees uniqueness of the least square reconstruction via the decoupledsystems. In the noisy case Theorem 2.1 provides an upper bound for the maximalerror of the least square reconstruction. The proof of these results relies on astability estimate for non-uniform sampling of exponential polynomials, whosederivation constitutes the bulk of the paper (Section 4). The principal resultthere, which might be of independent interest, is Theorem 4.4, whose proof is, inturn, based on the classical Tur´an-Nazarov inequality [18], and its generalizationto discrete sets, obtained recently in [12].We hope that our results will provide useful criteria for applicability of FourierDecoupling method to specific models of the form (1.1) in specific applications, D. BATENKOV AND N. SARIG AND Y. YOMDIN and a guiding principle for designing relevant sampling strategies. While the sta-bility bounds in Theorem 2.1 increase exponentially in the number of shifts q j ,this seems to be consistent with the general principle that reconstruction meth-ods based on sparsity are poorly conditioned with respect to model complexity,see [10].The paper is organized as follows: in Section 2 the method of Fourier decouplingof [24, 5] is presented in some detail, next we define the metric span ω and giveour main results. In Section 3 one specific example is considered in detail,illustrating, in particular, the importance of the frequency bound in the generalresults of Section 2. In Section 4 we study uniqueness and robustness of non-uniform sampling of exponential polynomials. Finally, in Section 5 some resultsof numerical simulations are presented. We consider signals of the form (1.1): F ( x ) = k (cid:88) j =1 q j (cid:88) q =1 a jq f j ( x − x jq ) , a jq ∈ C , x jq ∈ R . Here f j are known, while a jq , x jq are the unknown signal parameters, which wewant to find from Fourier samples F ( F )( s ) of F at certain sample points s ∈ R .Let F ( f j ) be the (known) Fourier transforms of f j .For F of the form (1.1) and for any s ∈ R we have for the sample of the Fouriertransform F ( F ) at s F ( F )( s ) = k (cid:88) j =1 q j (cid:88) q =1 a jq e − πisx jq F ( f j )( s ) . (2.1)In the case k = 1 we could divide the equation (2.1) by F ( f )( s ) and obtaindirectly a Prony-like equation. However, for k ≥ f , . . . , f k using the freedom in the choice of the sample set S . Let Z (cid:96) = (cid:8) x ∈ R , F ( f (cid:96) )( x ) = 0 (cid:9) denote the set of zeroes of the Fourier transform F ( f (cid:96) ). For each j = 1 , . . . , k ccuracy of Algebraic Fourier Reconstruction for Shifts of Several Signals S j to be a subset of the set W j = W j ( f , . . . , f k ) = (cid:18)(cid:92) (cid:96) (cid:54) = j Z (cid:96) (cid:19) \ Z j of common zeroes of the Fourier transforms F ( f (cid:96) ) , (cid:96) (cid:54) = j , but not of F ( f j ). Forsuch S j all the summands in (2.1) vanish, besides those with the index j . Hencewe obtain: Proposition 2.1. ([5]) Let for each j = 1 , . . . , k the sampling set S j satisfy S j = { s j , . . . , s jm j } ⊂ W j . Then for each j the corresponding system of equations (2.1) on the sample set S j takes the form q j (cid:88) q =1 a jq e − πix jq s j(cid:96) = c j(cid:96) , (cid:96) = 1 , , . . . , s j(cid:96) ∈ S j , (2.2) where c j(cid:96) = c j(cid:96) ( F ) = F ( F )( s j(cid:96) ) / F ( f j )( s j(cid:96) ) . These decoupled systems are exactly the same as the fitting systems for expo-nential polynomials H j ( s ) = (cid:80) q j q =1 a jq e − πix jq s on sets S j . So various methodsof exponential fitting can be applied (see, for example, [15, 20, 21, 25] and refer-ences therein.) As it was mentioned above, the main problem is that the samplesets W j may be non-uniform, and/or not sufficiently dense to provide a robustfitting. Indeed, the zeroes sets Z (cid:96) of the Fourier transforms F ( f (cid:96) ) may be anyclosed subsets G (cid:96) of R : it is enough to take F (cid:96) to be smooth rapidly decreasingfunctions on R with zeroes exactly on G (cid:96) , and to define the signals f (cid:96) as theinverse Fourier transforms of F (cid:96) . In particular, as a typical situation, Z (cid:96) maybe arbitrary finite sets or discrete sequences of real points.The main results of this paper provide a simple “density” condition on the sets W j ( f , . . . , f k ) as above, which guarantees a robust least square reconstructionof the signal F as in (1.1). We need some definitions:Let S be a bounded subset of R , and let I = [0 , R ( S )] be the minimal intervalcontaining S . Let λ ∈ R + be fixed. We put M = M ( N, λ, R ( S )) = N − (cid:98) λR ( S ) π (cid:99) , where for a real A , (cid:98) A (cid:99) denotes the integer part of A . Definition 2.1.
For N ∈ N , λ ∈ R + , the ( N, λ )-metric span of S is defined as ω N,λ ( S ) = max { , sup (cid:15)> (cid:15) [ M ( (cid:15), S ) − M ( N, λ, R ( S ))] } , D. BATENKOV AND N. SARIG AND Y. YOMDIN where M ( (cid:15), S ) is the (cid:15) -covering number of S , i.e. the minimal number of (cid:15) -intervals covering S ∩ I . Definition 2.2.
For each j = 1 , . . . , k the maximal frequency η j of the j -thequation in the decoupled system (2.1) is defined by η j = max q =1 ,...,q j π | x jq | . The minimal gap σ j of the j -th equation in (2.1) is defined by σ j = min ≤ p For each j = 1 , . . . , k the minimal divisor κ j = κ j ( S j ) of the j -th equation in the decoupled system (2.1) on S j is defined by κ j = min s ∈ S j |F ( f j )( s ) | . The sample gap ρ j of the j -th equation in (2.1) on I j is defined by ρ j = 3 R j σ j πq j ( q j + 1) f or η j R j ≤ πq j , and ρ j = 2 σ j η j q j ( q j + 1) otherwise. Theorem 2.1. Assume that for each j = 1 , . . . , k we have ω j := ω q j ,η j ( S j , R j ) > . Then the parameters a jq , x jq , q = 1 , . . . , q j , j = 1 , . . . , k of the signal F as in(1.1) can be uniquely reconstructed via the least square solution of the equations(2.2) on the sample sets S j , assuming that the measured samples of F ( F )( s ) atall the sample points are exact.In the case of noisy measurements, with the maximal error of the sample F ( F )( s j(cid:96) ) for s j(cid:96) ∈ S j being at most δ j (sufficiently small), we have the following boundsfor the reconstruction errors ∆ a jq , ∆ x jq of a jq , x jq : ∆ a jq ≤ κ j · (cid:18) R j ρ j ω j (cid:19) q j · δ j , (2.3)∆ x jq ≤ | a jq | κ j · (cid:18) R j ρ j ω j (cid:19) q j · δ j . (2.4) ccuracy of Algebraic Fourier Reconstruction for Shifts of Several Signals Proof: This theorem follows directly from Theorem 4.4 below, which estimatesthe accuracy of the least square sampling of exponential polynomials with purelyimaginary exponents on a given sampling set S . The only adaptation we haveto make is that the right hand sides c j(cid:96) of the equations (2.1) are given by c j(cid:96) = c j(cid:96) ( F ) = F ( F )( s j(cid:96) ) / F ( f j )( s j(cid:96) ), and hence the Fourier sampling error ismagnified by F ( f j )( s j(cid:96) ) . Consequently, the minimal divisor κ j = κ j ( S j ) of the j -th equation in (2.2) on S j appears in the denominator of (2.3) and (2.4). (cid:3) As a corollary we show that if the Fourier zeroes sets W j are “sufficiently dense”then the decoupling approach provides a robust reconstruction of the signal F .The notion of density we introduce below is a very restricted one. Much moreaccurate definition, involving not only the asymptotic behavior of W j ∩ [0 , R ]as R tends to infinity, but also its finite geometry, can be given. We planto present these results separately. Notice also a direct connection with theclassical Sampling Theory, in particular, with Beurling theorems of [9, 16]. Seealso [17, 19] and references therein. Definition 2.4. Let S be a discrete subset of R + . The “central density” D ( S )of S is defined as D ( S ) = lim sup R →∞ | S ∩ [0 ,R ] | R . Corollary 2.1. If for j = 1 , . . . , k we have D ( W j ) > η j π then the decouplingprocedure on appropriate sampling sets S j ⊂ W j provides a robust reconstructionof the signal F . Proof: If D ( W j ) > η j π then for arbitrarily big R j we have | W j ∩ [0 , R j ] | >M (2 q j , η j , R j ) , so taking sufficiently small (cid:15) > ω q j ,η j ( S j , R j ) is strictly positive. Application of Theorem 2.1 completes theproof. (cid:3) Some examples of Fourier decoupling have been presented in [24, 5]. In thepresent paper we consider one of these examples in more detail, stressing thequestion of robust solvability of the resulting decoupled systems. As everywherein this paper, we restrict ourselves to the case of one-dimensional signals. Someinitial examples in dimension two can be found in [24, 5].Let f be the characteristic function of the interval [ − , , while we take f ( x ) = D. BATENKOV AND N. SARIG AND Y. YOMDIN δ ( x − 1) + δ ( x + 1) . So we consider signals of the form F ( x ) = N (cid:88) q =1 [ a q f ( x − x q ) + a q f ( x − x q )] . (3.1)We allow here the same number N of shifts for each of the two signals f and f . Easy computations show that F ( f )( s ) = (cid:114) π sin ss and F ( f )( s ) = (cid:114) π cos s. So the zeros of the Fourier transform of f are the points πn, n ∈ Z \ { } andthose of f are the points ( + n ) π, n ∈ Z . These sets do not intersect, so we have W = { πn } , and W = { ( + n ) π } , and we can take as S , S any appropriatesubsets of these sets. Notice that the central density of W , W , according toDefinition 2.4, is π . By Corollary 2.1, if both the maximal frequencies η , η arestrictly smaller than 1 then the decoupling procedure on appropriate samplingsets S j ⊂ W j , j = 1 , , provides a robust reconstruction of the signal F . Let usshow that this condition on the frequencies η , η is sharp.The decoupled systems, given by (2.2) above, take the form N (cid:88) q =1 a q e − + n ) π ix q = N (cid:88) q =1 α q e iφ q n = c n , n ∈ Z , (3.2) N (cid:88) q =1 a q e − nπ ix q = N (cid:88) q =1 β q e iψ q n = c n , n ∈ Z , (3.3)where α q = a q e − iπ x q , φ q = − π x q , β q = a q , ψ q = − π x q ,c n = F ( F )(( + n ) π ) / F ( f )(( + n ) π ), c n = F ( F )( nπ ) / F ( f )( nπ ).Now we put in the equations (3.2), (3.3) α = i , α = − i , α q = 0 for q = 3 , . . . , N , φ = π, φ = − π , ccuracy of Algebraic Fourier Reconstruction for Shifts of Several Signals β = i , β = − i , β q = 0 for q = 3 , . . . , N , ψ = π, ψ = − π .This corresponds to the following shifts and amplitudes in the signal F : x = − π , x = π , a = − i a = i , and thus η = 1. x = − π , x = π , a = − a = , η = 1.For this specific signal F the exponential polynomials (3.2), (3.3) are both equalto sin( πs ), so they both vanish identically at all the sampling points s = n ∈ Z .Thus, allowing η = η = 1 we cannot reconstruct uniquely our signals from thesamples on the sets W , W , no matter how many sampling points we take.On the other hand, put η = η = η < 1. Let us take as S m (respectively, S m )the set of points of the form ( + n ) π (respectively, nπ ) for n = 0 , . . . , m . Wehave R ( S m ) = ( + m ) π, R ( S m ) = mπ . The number of the sample points in eachcase is m + 1 . So in computing ω N,η,mπ ( S m ) we have M = 4 N − (cid:98) ηmππ (cid:99) ∼ N + ηm − 1. So we have ω N,λ,mπ ( S m ) ∼ sup (cid:15)> (cid:15) [ M ( (cid:15), S m ) − N + 1 − ηm ] . Substituting here (cid:15) < π tending to π , we get M ( (cid:15), S m ) = m + 1, so ω ( S m ) ∼ [(1 − η ) m − N + 2] π. Essentially the same expression we get for ω ( S m ). So themetric spans of S m and S m are positive for m > N − − η . Applying Theorem 2.1we conclude that for such m the least square sampling on the sets S m , S m , iswell posed, and get explicit estimates for its accuracy. It would be very desirableto check the sharpness of this conclusion. Our numerical simulations, presentedbelow, provide an initial step in this direction. The decoupling method of [24, 5], presented in Section 2 above, reduces theFourier reconstruction problem for signals of the form (1.1) to a system of de-coupled equations (2.2), which are, for each j = 1 , . . . , k, the sampling equationsfor exponential polynomials of the form H j ( s ) = (cid:80) q j q =1 a jq e − πix jq s on samplingsets S j . So from now on we deal with sampling of exponential polynomials,not returning any more to the original problem of the Fourier reconstruction oflinear combinations of shifts of several signals.0 D. BATENKOV AND N. SARIG AND Y. YOMDIN We study robustness of sampling of exponential polynomials on the real line.Let H ( s ) = N (cid:88) j =1 a j e λ j s , a j , λ j ∈ C , s ∈ R , (4.1)be an exponential polynomial of degree N . We consider the following problem. Given a sampling set S ⊂ R , can an exponential polynomial H of degree N bereconstructed (i.e. its coefficients a j , λ j be recovered) from its known values on S ? If so, how robust can this reconstruction procedure be with respect to noisein the data? Let us now elaborate some assumptions we keep below.1. In this paper we deal with the case of only purely imaginary exponents λ j = ıφ j , φ j ∈ R . This assumption, which is satisfied in the case of Fourierreconstruction of the linear combinations of shifts of several signals, i.e. forthe fitting problem (2.2) above, strongly simplifies the presentation. Weplan to describe the general case of arbitrary complex exponents separately.2. We restrict ourselves to the least square reconstruction method, and donot consider other possible reconstruction schemes.3. In the noisy setting, we investigate the case of sufficiently small noise levels(for a more detailed explanation of this assumption see Theorem 4.3 below,and [3, 7, 26, 8]).As it was shown above, in order to ensure well-posedness of the (even noiseless)reconstruction problem, a certain “density” of the sampling set S with respectto the frequency set { φ , . . . , φ N } must be assumed. Accordingly, we assume anexplicit upper bound λ on the frequencies φ j and incorporate this bound intothe definition of the metric span of the sampling sets (compare Definition 2.1above). We shall also assume a lower bound on the minimal distance betweenthe frequencies: | φ j − φ i | ≥ ∆ . Without this assumption we cannot bound theaccuracy of the reconstruction of the amplitudes a j : indeed, as the exponents λ j of the exponential polynomial H ( s ) as in (4.1) collide, while the amplitudes a j tend to infinity in a pattern of divided finite differences, H ( s ) remains bounded ccuracy of Algebraic Fourier Reconstruction for Shifts of Several Signals λ > , ∆ > { φ , . . . , φ N } ≤ λ, min i 3. Now we choose inside the interval I a certain arithmetic progression ofpoints ¯ S = { s , s , . . . , (2 N − s } . The reconstruction problem on ¯ S isreduced to the standard Prony system. The right hand side of this Pronysystem, i.e. the values of ˜ H on ¯ S , would deviate from those of H not morethan allowed by the estimate of the previous step. Then, the deviationsof the reconstructed parameters ˜ a j , ˜ φ j from the original ones can finallybe estimated by the Lipschitz constant of the inverse Prony mapping, aspresented in [7] (see Theorem 4.3 below). An appropriate choice of s ispossible if we assume (as we do) that the exponents φ j do not collide.The rest of this section is organized as follows. The discrete Tur´an-Nazarovinequality is presented in Subsection 4.3. The stability estimates for the inverseProny mapping are reproduced in Subsection 4.4. The formulation of the finalestimate and its proof using the above steps are presented in Subsection 4.5. Let I = [0 , R ( S )] be the minimal interval containing S . Let N ∈ N and λ ∈ R + be fixed. We recall that the metric span ω N,λ ( S ) was defined asmax { , sup (cid:15)> (cid:15) [ M ( (cid:15), S ) − M ( N, λ, R ( S )] } , where M ( N, λ, R ) = N − (cid:98) λRπ (cid:99) .In this section we shall prove the following special case of the main result of [12]: Theorem 4.1. Let H ( s ) = (cid:80) Nj =1 a j e λ j s be an exponential polynomial, where a j ∈ C , λ j = iφ j , φ j ∈ R , λ = max j =1 ,...,N | φ j | . Let S, I be as above, with ω N,λ ( S ) > . Then we have sup I | H ( s ) | ≤ (cid:18) R ( S ) ω N,λ ( S ) (cid:19) N − · sup S | H ( s ) | . (4.3)We follow the lines of proof of Theorem 1.3 of [12]. We shall use the followingtwo results from [18]. Lemma 4.1 (Lemma 1.3 in [18], Langer’s lemma) . Let p ( z ) = (cid:80) nk =1 c k e iλ k z bean exponential polynomial ( < λ < · · · < λ n = λ ) not vanishing identically.Then the number of complex zeros of p ( z ) in an open vertical strip x < (cid:60) z Let ρ := sup S | H ( s ) | and consider the sublevel set V ρ := { t ∈ I : | H ( s ) | < ρ } . Define p ( z ) = H ( z ) − ρ . It is an exponential polynomial of degree at most N with purely imaginary exponents whose absolute values are bounded from aboveby 2 λ . By Lemma 4.1 above the number of solutions to the equation p ( t ) = 0in I (which is equivalent to | H ( s ) | = ρ ) is at most N − (cid:22) λR ( S ) π (cid:23) = M ( N, λ, R ( S )) . Therefore, the set V ρ consists of at most M ( N, λ, R ( S )) subintervals ∆ i . Nowfix (cid:15) > (cid:15) -covering number M ( (cid:15), V ρ ). In order to cover each ofthe ∆ i ’s, we need at most µ (∆ i ) (cid:15) + 1 (cid:15) -intervals. Overall, we get M ( (cid:15), V ρ ) (cid:54) M ( N, λ, R ( S )) + µ ( V ρ ) (cid:15) . (4.5)By definition of ρ , we obviously have that S ⊆ V ρ , therefore M ( (cid:15), S ) (cid:54) M ( (cid:15), V ρ ).Substituting this into (4.5), multiplying by (cid:15) and taking supremum w.r.t. (cid:15) weobtain µ ( V ρ ) (cid:62) ω N,λ ( S ). Now we just apply Theorem 4.2 with p = H and E = V ρ . (cid:3) Notice that the result of Theorem 4.1 does not depend at all on the minimaldistance between the exponents of H , which is crucial in the rest of our estimates,and does not imply any bound on the amplitudes a j of H . Let x , . . . , x N be pairwise distinct complex numbers, and let a , . . . , a N benonzero complex numbers. In [7] we introduced the “Prony map”, P : C N → C N , defined by P ( x , . . . , x N , a , . . . , a N ) = ( m , . . . , m N − ) , m k = N (cid:88) j =1 a j x kj . D. BATENKOV AND N. SARIG AND Y. YOMDIN This mapping can be considered as the sampling operator for the exponentialpolynomial H ( s ) = (cid:80) Nj =1 a j x sj on the integer points s ∈ { , , . . . , N − } . In[7] we provided local perturbation estimates for P , as follows. Theorem 4.3. Let x , . . . , x N be pairwise distinct complex numbers, and let a , . . . , a N be nonzero complex numbers. Let x = ( m , . . . , m N − ) be the imageof the point ( x , . . . , x N , a , . . . , a N ) under the Prony map P . Let δ > besufficiently small, so that the inverse map P − is defined in the δ -neighborhood U of x . Let ˜ x be some point in this neighborhood: ˜ x = ( m + δ , . . . , m N − + δ N − ) , | δ i | < δ. Then the image of ˜ x under P − satisfies | a j − ˜ a j | (cid:54) C ( x , . . . , x N ) δ, | x j − ˜ x j | (cid:54) C ( x , . . . , x N ) | a j | − δ, (4.6) where C ( x , . . . , x N ) depends only on the configuration of the nodes x , . . . , x N . In fact, as we show in [3], in the case that x , . . . , x N belong to the unit circle,the constant C can be bounded from above by C (cid:54) · (cid:18) (cid:19) N , (4.7)where Λ = min i Let H ( s ) = (cid:80) Nj =1 a j e ıφ j s be an a-priori unknown exponentialpolynomial satisfying max | φ j | ≤ λ, min | φ i − φ j | ≥ ∆ for some fixed λ, ∆ . Letthere be given the noisy samples h ( s ) of H ( s ) on a finite set S = { s , . . . , s n } ⊂ R , with the noise bounded by δ, i.e. | h ( s (cid:96) ) − H ( s (cid:96) ) | ≤ δ, (cid:96) = 1 , . . . , n. Assumethat ω ( S ) := ω N,λ ( S ) > . Then, for sufficiently small δ , the amplitudes ˜ a j and ccuracy of Algebraic Fourier Reconstruction for Shifts of Several Signals the frequencies ˜ φ j of the least square fitting exponential polynomial ˜ H ( s ) satisfy,for j = 1 , . . . , N, the following inequalities: | a j − ˜ a j | ≤ √ n · (cid:18) R ( S ) ρω ( S ) (cid:19) N · δ, (4.8) | φ j − ˜ φ j | ≤ √ n | a j | − · (cid:18) R ( S ) ρω ( S ) (cid:19) N · δ, (4.9) where ρ = R ( S )∆2 πN ( N +1) for λR ( S ) ≤ πN , and ρ = λN ( N +1) otherwise. In par-ticular, in the case of zero noise, the least square reconstruction of ˜ H ( s ) on S ,under the constraints as above, is unique, up to a transposition of the indices. Proof: First of all, let us establish the following easy bound: Lemma 4.2. For s ∈ S we have | ˜ H ( s ) − H ( s ) | ≤ √ nδ . Proof: Indeed, the quadratic deviation σ ( H, h ) of H from h on S does notexceed nδ , where n , as above, denotes the number of elements in S . Since˜ H ( s ) is the exponential polynomial of the least square deviation from h , wehave σ ( ˜ H, h ) ≤ nδ , which directly implies σ ( ˜ H, H ) ≤ nδ and hence | ˜ H ( s ) − H ( s ) | ≤ √ nδ, for each s ∈ S . (cid:3) Now we get directly the following bound: Corollary 4.1. For H, S and I = [0 , R ( S )] the minimal interval containing S we have sup I | ˜ H ( s ) − H ( s ) | ≤ (cid:18) R ( S ) ω N,λ ( S ) (cid:19) N − · √ nδ. (4.10) Proof: We notice that by Lemma 4.2 we have sup S | ˜ H ( s ) − H ( s ) | ≤ √ nδ .Substituting into Theorem 4.1 (which is applied to the exponential polynomial H ( s ) − ˜ H ( s ) of degree 2 N with purely imaginary exponents, bounded in absolutevalue by λ ), we get the required bound. (cid:3) The bound of Corollary 4.1 does not imply by itself any bound on the amplitudes a j . They may tend to infinity, as the exponents collide, following the patternof divided finite differences (see [26, 8]). So the continuation of the proof in-corporates the a priori known lower bound ∆ on the differences between theexponents of H . We get estimates of the reconstruction accuracy of a j and λ j D. BATENKOV AND N. SARIG AND Y. YOMDIN via solving an appropriate auxiliary Prony system, and applying Theorem 4.3above.Let I = [0 , R ( S )] be as above. Fix certain s ∈ (0 , R N ] and consider the points s , s , . . . , (2 N ) s ∈ I . We denote ν k = H ( ks ) , k = 0 , , . . . , the values of H at the points ks . We get H ( ks ) = N (cid:88) j =1 a j e (cid:96) j ks = N (cid:88) j =1 a j x kj = ν k , k ∈ Z , (4.11)where x j = e λ j s = e iφ j s . So for each choice of s ∈ (0 , R N ] we obtain a Pronysystem N (cid:88) j =1 a j x kj = ν k , k = 0 , . . . , N − , (4.12)which is satisfied by a j and x j = e iφ j s , j = 1 , . . . , N . It is well known thatif x i (cid:54) = x j for i (cid:54) = j, then the solution a j , x j = e iφ j s , j = 1 , . . . , N of (4.12)is unique, up to a permutation of the indices. Moreover, the robustness of thesolutions of (4.12) with respect to the perturbations of the right-hand side, isdetermined by the mutual distances | x i − x j | , i (cid:54) = j (see [7, 6] and Subsection 4.4).So our next goal is to choose s ∈ (0 , R N ] in such a way that Λ = min i (cid:54) = j | x i − x j | be sufficiently large. To achieve this goal we have to find s such that all theangles ∆ i,j s are separated from the integer multiples 2 πm, m ∈ Z , where∆ i,j = | φ j − φ j | . An easy example shows that there are “bad” choices of s : assume that thefrequencies φ j in H are of the form φ j = s · πm j , with s ∈ R , m j ∈ Z , m i (cid:54) = m j for i (cid:54) = j . Then for s = s we have x = x = . . . = x N . The next lemma showsthat most choices of s are good, assuming that ∆ = min i Let ¯ R, q , . . . q r ∈ R + be given, with q = min q (cid:96) , Q = max q (cid:96) There exists s ∈ (0 , ¯ R ] such that all the angles q (cid:96) s , (cid:96) = 1 . . . , r, are separatedfrom the integer multiples πm, m ∈ Z by at least ¯ h , defined as ¯ h = ¯ Rq r ≤ π for Q ¯ R ≤ π , and as ¯ h = πq rQ ≤ π for Q ¯ R > π . Proof: By assumptions we have q (cid:96) ≤ Q . Hence for each (cid:96) = 1 , . . . , r the interval q l · (0 , ¯ R ] contains at most Q ¯ R π + 1 integer multiples 2 πm . For h > U ( h )denote the h -neighborhood of these points. Denote by µ the standard Lebesguemeasure on R . We have µ ( U ( h )) ≤ h [ Q ¯ Rπ + 2]. Now let V (cid:96) ( h ) denote the set ofthose s ∈ (0 , ¯ R ] for which q (cid:96) s ∈ U ( h ). We conclude that µ ( V (cid:96) ( h )( h )) ≤ hq (cid:96) [ Q ¯ Rπ + ccuracy of Algebraic Fourier Reconstruction for Shifts of Several Signals ≤ hq [ Q ¯ Rπ + 2]. Finally, denoting V ( h ) the set of the points s ∈ (0 , ¯ R ] for which q (cid:96) s ∈ U ( h ) for at least one index (cid:96) = 1 , . . . , r , we get µ ( V ( h )) ≤ rhq [ Q ¯ Rπ + 2]. Iffor some h we have µ ( V ( h )) < | (0 , ¯ R ] | = ¯ R , then there exists s ∈ (0 , ¯ R ] suchthat all the angles q l s , (cid:96) = 1 , . . . , r, are separated from the integer multiples2 πm, m ∈ Z at least by h .Now we consider two cases: Q ¯ R ≤ π and Q ¯ R > π . In the first case µ ( V ( h )) ≤ rhq , and the inequality µ ( V ( h )) ≤ ¯ R is valid with h = ¯ h = ¯ Rq r ≤ π . Inthe second case µ ( V ( h )) ≤ rhq Q ¯ Rπ , and the inequality µ ( V ( h )) ≤ ¯ R is valid,starting with h = ¯ h = πq rQ ≤ π . This completes the proof of the lemma. (cid:3) In our case of the angles ∆ i,j and s ∈ (0 , R N ] we have, respectively, ¯ R = R N , r = N ( N +1)2 , Q = 2 λ, q = ∆. Applying Lemma 4.3 we obtain the following result: Corollary 4.2. There exists s ∈ (0 , R N ] such that all the angles ∆ i,j · s , πN . Accordingly, the minimal distance min | x i − x j | , i (cid:54) = j, betweenthe points x j = e iφ j s , j = 1 , . . . , N in (4.12) is at least ρ = π ¯ h, which is ρ = R ∆2 πN ( N +1) for λR ≤ πN , and ρ = λN ( N +1) for λR > πN . Proof: The result on the separation of the angles follows directly from Lemma4.3. The result for the distances follows from the fact that always ¯ h ≤ π . (cid:3) Now we can complete the proof of Theorem 4.4. We fix s whose existence isguaranteed by Lemma 4.3, and form the Prony system, which is satisfied by theparameters of H : N (cid:88) j =1 a j x kj = ν k , k = 0 , . . . , N − , x j = e λ j s , (4.13)with ν k = H ( ks ) the values of H at the points ks , k = 0 , . . . , N − 1. Noticethat these values are not exactly known. However, by Corollary 4.1 we knowthat sup J | ˜ H ( s ) − H ( s ) | ≤ (cid:18) Rω N,λ ( S ) (cid:19) N − · √ nδ, (4.14)where ˜ H ( s ) = (cid:80) Nj =1 ˜ a j e ˜ λ j s is the polynomial of the least square approximationon S . In particular, denoting ˜ ν k = ˜ H ( ks ) , k = 0 , . . . , N − , the values of ˜ H D. BATENKOV AND N. SARIG AND Y. YOMDIN at the points ks , we get | ˜ ν k − ν k | ≤ (cid:16) Rω N,λ ( S ) (cid:17) N − · √ nδ. Now the parameters ˜ a j , ˜ λ j of ˜ H satisfy the Prony system N (cid:88) j =1 ˜ a j ˜ x kj = ˜ ν k , k = 0 , . . . , N − , ˜ x j = e ˜ λ j s . (4.15)Finally we apply Theorem 4.3 to Prony system (4.13) and its perturbation (4.15),taking into account the expression (4.7) for the constant C in Theorem 4.3.Noticing that the distances between the nodes x j of the unperturbed Pronysystem (4.13) are bounded from below by ρ via Corollary 4.2, we arrive at (4.8)and (4.9). Uniqueness of reconstruction for δ = 0 follows directly from (4.8) and(4.9). This completes the proof of Theorem 4.4. (cid:3) ω N,λ ( S ) : some examples The metric span ω N,λ ( S ) can be explicitly computed in many important cases.In particular, we have the following simple result: Proposition 4.1. Let N, λ be fixed. Assume that a subset S ⊂ R with R ( S ) = R contains M ( N, λ, R ) + 1 points, and let η be the minimal distance between theneighboring points in S . Then ω N,λ ( S ) = η. Proof: For (cid:15) ≥ η we have M ( (cid:15), S ) − M ( N, λ, R ) ≤ 0. For (cid:15) < η this differenceis 1. Hence the supremum in Definition 2.1 is achieved as (cid:15) tends to ρ from theleft. (cid:3) Corollary 4.3. Let N, λ be fixed. Assume that a subset S ⊂ R contains M ( N, λ, R ) + 1 points. Then ω N,λ ( S ) ≤ RM ( N,λ,R ) , and this value is achievedonly for S consisting of M ( N, λ, R ) + 1 points at the distance RM ( N,λ,R ) one fromanother. Proof: For the equidistant configuration the minimal distance η between theneighboring points in S is RM ( N,λ,R ) . Otherwise η is strictly smaller. (cid:3) .Now let us consider equidistant configurations with a larger number of samplingpoints. ccuracy of Algebraic Fourier Reconstruction for Shifts of Several Signals Proposition 4.2. Let N, λ be fixed. Assume that a subset S ⊂ R contains m + 1 ≥ M ( N, λ, R ) + 1 points at the distance Rm from one another. Then ω N,λ ( S ) = (cid:18) Rm (cid:19) [ m + 1 − M ( N, λ, R )] . (4.16) Proof: For the equidistant configuration S the minimal distance η betweenthe neighboring points in S is Rm . On the other hand, for each (cid:15) < η we have M ( (cid:15), S ) = m + 1 , while for kη ≤ (cid:15) ≤ ( k + 1) η, k = 1 , , ..., we have M ( (cid:15), S ) = m +1 k . An easy computation then shows that the supremum of (cid:15) [ M ( (cid:15), S ) − M ( N, λ, R )]is achieved for (cid:15) tending to η from the left, and it is equal to ( Rm )[ m + 1 − M ( N, λ, R )] . (cid:3) Remark 4.1. As substituted into the expression of Theorem 4.4, the resultsabove imply the corresponding bound for the accuracy of the least square re-construction on S . In particular, the expression (4.16) above seems to pro-vide a non-trivial recommendation for the choice of the number of equidistantsample points inside a given interval I . Indeed, for m = M ( N, λ, R ) we get ω ( S ) = RM ( N,λ,R ) . But for m = 2 M ( N, λ, R ) we get ω ( S ) = R M ( N, λ, R ) (cid:2) M ( N, λ, R ) + 1 (cid:3) , which is approximately R for large M ( N, λ, R ) - improvement by M times.For m tending to infinity ω ( S ) tends to R , so we do not achieve any essentialimprovement any more. Thus the recommendation may be to take m of order KM d with K between, say 2 and 5. Remark 4.2. Combining Theorem 4.4 and Proposition 4.1 we can also predictthe rate of the degeneration of the reconstruction problem on S as two points of S collide. By the same method we can analyse also the cases of more complicatedcollisions between the sampling points. In this section we present results of initial numerical experiments. Our goal inthese very preliminary simulations has been to numerically investigate the qual-itative dependence of the reconstruction error on the geometry of the samplingset S . Our results below are indeed qualitatively consistent with the bounds ofTheorem 4.4.0 D. BATENKOV AND N. SARIG AND Y. YOMDIN e rr o r Figure 1: In this experiment, we changed the mutual distance d between thesubsequent points of S , while keeping the two endpoints fixed. ε = 10 − , λ =1 , R = 60 , N = 2. The size of S is n = 35. The error is plotted versus thevalue of d in red. For comparison, the value d is plotted in blue. 20 30 40 50 60 700 101e-71e-61e-51e-41e-31e-2 n e rr o r Test Figure 2: In this experiment, we increased n , the number of points in S , keepingthe range (i.e. the value of R ) fixed. ε = 10 − , λ = 0 . , R = 10 , N = 2. Here M (2 N, λ, R ) = 15. The error is plotted versus the value of n . ccuracy of Algebraic Fourier Reconstruction for Shifts of Several Signals H ( s ), and modified the samplingset S according to the description of each experiment below. The sampling values { H ( s i ) , s i ∈ S } have been perturbed by the (random) amount ε ∼ − .Subsequently, the least-squares approximation to H ( s ) has been obtained bythe standard sequential quadratic programming algorithm (implemented by thefunction sqp in GNU Octave environment). The initial values for the algorithmhave been taken to be equal to the true values perturbed by the (random)amount ε , specified in each experiment below. We have plotted the recoveryerror for one of the frequencies (specifically, | ∆ φ | ).In the first experiment we changed the distance d between s , . . . , s n − , whilekeeping the endpoints s , s n (and thereby the value of R ) fixed. The number ofpoints was chosen to be exactly n = M (2 N, λ, R ) + 1. According to Proposition4.1, in this case we have ω ( S ) = d . As can be seen in Figure 1, the error isroughly proportional to ω ( S ) .In the second experiment, we have kept the endpoints of the set S fixed (0 and R ), while increasing the number n of (equispaced) points in S . According toFigure 2, a significant improvement in accuracy appears when the number ofsamples passes M (2 N, λ, R ) which is 15 in this case. References [1] B. Adcock, M. Gataric, and A.C. Hansen On stable reconstructions fromunivariate nonuniform Fourier measurements. Preprint. Arxiv: 1310.7820 [2] D. Batenkov. Complete Algebraic Reconstruction of Piecewise-SmoothFunctions from Fourier Data. To appear in Mathematics of Computation. [3] D. Batenkov. Decimated Generalized Prony systems. Preprint.arXiv:1308.0753. [4] D. Batenkov, N. Sarig, and Y. Yomdin. An “algebraic” reconstruction ofpiecewise-smooth functions from integral measurements. Functional Differ-ential Equations , 19(1-2):9–26, 2012.[5] D. Batenkov, N. Sarig, and Y. Yomdin. Decoupling of Reconstruction Sys-tems for Shifts of Several Signals. Proc. of Sampling Theory and Applica-tions (SAMPTA) , 2013.[6] D. Batenkov and Y. Yomdin. Algebraic reconstruction of piecewise-smoothfunctions from Fourier data. Proc. of Sampling Theory and Applications(SAMPTA) , 2011.2 D. BATENKOV AND N. SARIG AND Y. YOMDIN [7] D. Batenkov and Y. Yomdin. On the accuracy of solving confluent Pronysystems. SIAM J. Appl. Math. , 73(1):134–154, 2013.[8] D. Batenkov and Y. Yomdin. Geometry and Singularities of the PronyMapping. To appear in Proceedings of 12th International Workshop onReal and Complex Singularities , 2013.[9] A. Beurling. Balayage of Fourier-Stiltjes Transforms. The collected Worksof Arne Beurling, Vol.2, Harmonic Analysis . Birkhauser, Boston, 1989.[10] D.L. Donoho. Superresolution via sparsity constraints. SIAM Journal onMathematical Analysis, IEEE Transactions on Signal Processing, Vol. 55, Nr. 5, Part 1, pp. 1741-1757, 2007.[12] O. Friedland and Y. Yomdin. An observation on Tur´an-Nazarov inequality. Studia Mathamatica , 218(1), pp. 27–39, 2013. DOI 10.4064/sm218-1-2 [13] K. Gedalyahu, R. Tur, and Y.C. Eldar. Multichannel sampling of pulsestreams at the rate of innovation. IEEE Transactions on Signal Processing ,59(4):1491–1504, 2011.[14] B. Gustafsson, C. He, P. Milanfar and M. Putinar. Reconstructing planardomains from their moments. Inverse Problems , 16(4):1053–1040, 2000.[15] Liviu Gr Ixaru and Guido Vanden Berghe. Exponential Fitting . Springer,May 2004.[16] H. Landau. Necessary density conditions for sampling and interpolation ofcertain entire functions. Acta Mathematica , 117(1):37–52, 1967.[17] F. Marvasti. Nonuniform sampling: theory and practice . Springer, 2001.[18] F.L. Nazarov. Local estimates of exponential polynomials and their appli-cations to inequalities of uncertainty principle type. St Petersburg Mathe-matical Journal , 5(4):663–718, 1994.[19] A. Olevski, A. Ulanovski. Near critical density irregular sampling inBernstein spaces. Mathematisches Forschungsinstitut Oberwolfach gGmbH Exponential data fitting and its appli-cations . Bentham Science Publishers, January 2010. ccuracy of Algebraic Fourier Reconstruction for Shifts of Several Signals SIAM Journal on Scientific Computing ,33(4):1920, 2011.[22] B.D. Rao and K.S. Arun. Model based processing of signals: A state spaceapproach. Proceedings of the IEEE , 80(2):283–309, 1992.[23] N. Sarig and Y. Yomdin. Signal Acquisition from Measurements via Non-Linear Models. Mathematical Reports of the Academy of Science of theRoyal Society of Canada , 29(4):97–114, 2008.[24] N. Sarig. Algebraic reconstruction of ”shift-generated” signals from integralmeasurements . PhD thesis, Weizmann Institute of Science, 2010.[25] P. Stoica and R.L. Moses. Spectral analysis of signals . Pearson/PrenticeHall, 2005.[26] Y. Yomdin. Singularities in algebraic data acquisition.