Exact Reconstruction of Extended Exponential Sums using Rational Approximation of their Fourier Coefficients

In this paper we derive a new recovery procedure for the reconstruction of extended exponential sums of the form $y(t) = \sum_{j=1}^{M} \left( \sum_{m=0}^{n_j} \, \gamma_{j,m} \, t^{m} \right) {\mathrm e}^{2\pi \lambda_j t}$, where the frequency parameters $\lambda_{j} \in {\mathbb C}$ are pairwise distinct. For the reconstruction we employ a finite set of classical Fourier coefficients of $y$ with regard to a finite interval $[0,P] \subset {\mathbb R}$ with $P>0$. Our method requires at most $2N+2$ Fourier coefficients $c_{k}(y)$ to recover all parameters of $y$, where $N:=\sum_{j=1}^{M} (1+n_{j})$ denotes the order of $y$. The recovery is based on the observation that for $\lambda_{j} \not\in \frac{{\mathrm i}}{P} {\mathbb Z}$ the terms of $y$ possess Fourier coefficients with rational structure. We employ a recently proposed stable iterative rational approximation algorithm in [12]. If a sufficiently large set of $L$ Fourier coefficients of $y$ is available (i.e., $L>2N+2$), then our recovery method automatically detects the number $M$ of terms of $y$, the multiplicities $n_{j}$ for $j=1, \ldots , M$, as well as all parameters $\lambda_{j}$, $j=1, \ldots , M$ and $ \gamma_{j,m}$ $j=1, \ldots , M$, $m=0, \ldots , n_{j}$, determining $y$. Therefore our method provides a new stable alternative to the known numerical approaches for the recovery of exponential sums that are based on Prony's method.


Introduction
Recently, we have proposed a new reconstruction method to recover real functions of the form from a limited number of classical Fourier coefficients of y from its Fourier expansion on a given fixed interval [0, P ], see [17].
This paper continues and strongly generalizes our research started in [17]. We present a new reconstruction method to recover complex extended exponential sums, i.e., sums being of polynomial exponential form, γ j,m t m e λ j t , γ j,m ∈ C, γ j,n j = 0, λ j ∈ C.
First we introduce the main objects studied in the paper. For N ∈ N, we consider the set of proper exponential sums We call n = n(y) the order of the exponential sum y(t) and M = M (y) its length. Obviously, Y 0 N is a subset of Y N and we have n(y) ≥ M (y) in the general case and n(y) = M (y) for y ∈ Y 0 N . In fact, Y 0 N is a dense subset of Y N , and Y N is closed with respect to the maximum norm in C [a, b] for any compact interval [a, b]. Moreover, Y N is an existence set for the space C[a, b] of continuous functions given on a compact interval [a, b], see, [23] or [5,Chapter VI]. Approximation with extended exponential sums has a long history. We refer to [10], [26] that are dedicated to the question of approximation of classes of smooth functions L p [a, b], 1 ≤ p < ∞, and C[a, b] by extended exponential sums. Further, it is well-known for a long time that there is a close connection between approximation with exponential sums and rational approximation, see [5,24].
The extended exponential sums y(t) in Y N of order n ≤ N are solutions of homogeneous linear differential equations of order n with constant coefficients a 0 , . . . , a n−1 ∈ C of the form y (n) + a n−1 y (n−1) + . . . + a 1 y + a 0 y = 0. (1. 3) The corresponding characteristic polynomial is given as p(λ) = λ n + a n−1 λ n−1 + . . . + a 1 λ + a 0 = M j=1 (λ − λ j ) n j +1 , i.e., p(λ) has M distinct roots λ j with multiplicity n j + 1, j = 1, . . . , M . In particular, the functions t m e λ j t , m = 0, . . . , n j , j = 1, . . . , M , are linearly independent and form a basis of the space of the solutions of the differential equation (1.3). Similarly, it can be shown that Y N is the solution space for homogeneous linear difference equations of order n ≤ N , see e.g. [4], and therefore Y N is closely related to the characterization of Hankel operators of finite rank, [9,13].
Extended exponential sums appear in many applications in system identification and sparse approximation, see e.g. [1,14]. For a comprehensive study of extended exponential sums from an algebraic point of view, we refer to [11].
Our goal in this paper is to reconstruct the exponential sums in Y 0 N and Y N . This problem has been extensively studied using Prony's method and its generalizations, see for example [6], [7], [16], [21], [19], [20]. However, the known studies mainly focussed on the recovery of proper exponential sums. The extended model is less often treated, see [3,2,15,22,24,25]. Since Prony's method involves computations with Hankel or Toeplitz matrices with possibly high condition numbers, it requires a very careful numerical treatment. This is particularly true if the distance between two distinct frequency parameters λ j 1 and λ j 2 is very small. The extended exponential sums appear if such frequencies collide.
Our new method for reconstruction of signals in form of (extended) exponential sums is based on rational approximation and can be seen as a good alternative to the Prony reconstruction approach. As input information for our new algorithm we use a finite set of their Fourier coefficients. More precisely, we consider the Fourier expansion of y ∈ Y N on [0, P ] for some P with 0 < P < ∞ of the form y(t) = k∈Z c k (y) e 2πikt/P with Fourier coefficients c k (y) := 1 P P 0 y(t) e −2πikt/P dt. We will derive algorithms to reconstruct y from a sufficiently large set of given Fourier coefficients c k (y). Although Y 0 N is a subset of Y N we will separate the reconstruction of proper exponential sums as a special case due to the fact that in many applications only the set Y 0 N is considered. If we have −iλ j ∈ 1 P Z for all frequencies λ j , i.e., if y(t) in (1.1) or (1.2) does not possess any P -periodic terms, then we will employ the special property that the classical Fourier coefficients c k (y) of y ∈ Y 0 N and y ∈ Y N of full order N can usually be represented by a rational function r N of type (N − 1, N ), i.e., c k (y) = r N (k) for k ∈ Z. Our reconstruction approach is then based on the recovery of this rational function r N using a modification of the AAA algorithm that has been recently proposed in [12]. Numerical stability of this algorithm is ensured by barycentric representation of the constructed approximant. The second important property of the AAA algorithm is that it works iteratively thereby enlarging the polynomial degree of the numerator and denominator of the rational approximant step by step. The modified AAA algorithm will terminate if the needed order N of the exponential sum is reached. This gives us the opportunity to reconstruct the order N of the exponential sum, supposed that a sufficiently large set of L ≥ 2N + 1 Fourier coefficients is given.
If y ∈ Y N is an extended exponential sum, then the corresponding rational function r N has M multiple poles with multiplicities n j + 1 such that M j=1 (1 + n j ) = N . Having reconstructed the rational function r N (t), we will show, how all wanted parameters that determine y ∈ Y N can be uniquely computed from r N . More exactly, beside N , M , and n j , j = 1, . . . , M , we can determine all frequencies λ j and all polynomial coefficients γ j,m j = 1, . . . , M , m = 0, . . . , n j , from r N .
Further we will extend the algorithm for the case when P -periodic components also appear in the expansions (1.1) and (1.2). In this case we first reconstruct the rational function that determines the non-P -periodic part of y(t) using modified AAA-algorithm and recover all parameters of the non-P -periodic part of y(t) as before. In a second step, the P -periodic part of y(t) is determined. Note that this is only possible if the given set of Fourier coefficients contains all c k j (y) that are related to the occurring frequencies that have to be recovered, i.e., we need knowledge about c k j (y) if λ j = ik j P is a frequency parameter that occurs in y(t). For proper exponential sums with periodic terms it can be simply shown that only a finite number of Fourier coefficients of y(t) loses its special rational structure. The recovery of P -periodic parts for (truly) extended exponential sums requires a special treatment, since in this case the polynomial coefficients corresponding to P -periodic terms still influence all Fourier coefficients of y.
Outline. In Section 2 we prove that the signals y ∈ Y N (and consequently y ∈ Y 0 N ) are uniquely determined by given parameters M , n j , λ j , γ j,m for j, . . . , M , m = 0, . . . , n j . In Section 3 we shortly describe the main ideas of the needed modified AAA-algorithm for rational approximation of functions. In Section 4 we consider the recovery of proper exponential sums. We separate two possible cases. In Section 4.1 we study the recovery of proper exponential sums that contain only non-P -periodic terms, and Section 4.2 is dedicated to the case when y ∈ Y 0 N contains also P -periodic components. In Section 5 we consider the recovery of extended exponential sums. In Section 5.1 we show that the Fourier coefficients of y ∈ Y N can be represented by a rational function, and we show, how the wanted parameters can be reconstructed from this rational representation. Again we separate the consideration of the following two cases: In Section 5.2 we study the reconstruction of extended exponential sums containing only non-P -periodic components. In Section 5.3, we assume that extended exponential sums contain also P -periodic terms. The proposed reconstruction algorithms in Sections 4 and 5 are illustrated by numerical examples. The corresponding software can be found on our homepage http://na.math.unigoettingen.de We conclude the paper with some final remarks in Section 6.
Notation. As usual N, Z, R and C are reserved for natural, integer, real and complex numbers and let N 0 := N ∪ {0}.

Uniqueness of the Representation of Extended Exponential Sums
In order to be able to reconstruct functions y(t) from Y N uniquely, we need to ensure that all parameters of can be uniquely determined. For this purpose we assume the following restrictions on the parameters: (1) M ∈ N, M < ∞; (2) n j ∈ N 0 , n j < ∞; (3) γ j,m ∈ C for j = 1, . . . , M and m = 0, . . . , n j , and γ j,n j = 0 for j = 1, . . . , M ; (4) λ j ∈ C are pairwise distinct for j = 1, . . . , M .
Conditions (3) and (4) do not restrict the set Y N with N ≥ M j=1 (1 + n j ). If γ j,n j = 0, then the corresponding term in (2.1) can be removed, and if λ j 1 = λ j 2 for two frequency parameters in (2.1), the corresponding terms can be combined into one term. Note that the third condition implies that the polynomial p n j (t) has exactly the degree n j . In the special case y ∈ Y 0 N , the restriction (2) simplifies to n j = 0 for j = 1, . . . , M , and restriction (3) reads γ j = γ j,0 = 0 for j = 1, . . . , M . We show now that with the restrictions above, all parameters of y(t) are uniquely determined. Theorem 2.1. Let two extended exponential sums be of the form and assume that the parameters for y 1 (t) and y 2 (t) satisfy the restrictions (1) − (4) given above. If y 1 (t) = y 2 (t) pointwise on a finite interval with positive length, then M 1 = M 2 =: M and (after suitable permutation of summands) n j = n j =: n j , λ j = λ j for j = 1, . . . , M , and γ j,m = γ j,m for j = 1, . . . , M , m = 0, . . . , n j .
Proof. We consider the function Let y 1 (t) = y 2 (t) for t ∈ I, where I ⊂ R is some interval with finite positive length. Then y ≡ 0 on I.
Let M denote the number of pairwise distinct frequency parameters λ j in the representation (2.2) of the function y(t). According to the restriction (4) on the frequency parameters we have M ≥ max{M 1 , M 2 }. Then we can rewrite Each parameterλ can occur once or twice in the set Λ. If it occurs once, for examplẽ λ = λ j , then (2.4) implies that the corresponding polynomial coefficient in p n (t) = n m=0γ ,m t m completely vanishes, i.e.,γ ,m = 0 for m = 0, . . . , n . But this contradicts the assumption (3). Therefore, eachλ occurs twice in the set Λ and we can conclude that M 1 = M 2 = M . Further, for each = 1, . . . , M , we find j 1 , j 2 ∈ {1, . . . , M } such that λ = λ j 1 = λ j 2 . Then the polynomial coefficient corresponding to eλ t in (2.3) satisfies Since by restriction (3), γ j 1 ,n j 1 = 0 and γ j 2 ,n j 2 = 0, we conclude that n j 1 = n j 2 =: n j and γ j 1 ,m = γ j 2 ,m for m = 0, . . . , n j .

The modified AAA Algorithm for Rational Approximation
Our reconstruction algorithms for exponential sums in Sections 4 and 5 are essentially based on the following slight modification of the AAA algorithm in [12]. For the convenience of the reader we shortly summarize this algorithm with the needed slight modification for rational functions of type (N − 1, N ). For a similar modification for special point sets we refer to [17].
For a given sufficiently large finite set of pairwise distinct points Γ ⊂ R and a corresponding set of values {f (z) : z ∈ Γ}, the (modified) AAA algorithm iteratively computes a rational function r N of type (N − 1, N ) such that r N (z) = f (z) for z ∈ S, where S ⊂ Γ is a subset of N + 1 given points, and such that the error |r N (z) − f (z)| is small for the remaining points z ∈ Γ \ S.
At the iteration step J ≥ 1, we proceed as follows. Assume that we have the set S J+1 = {z 1 , . . . , z J+1 } ⊂ Γ, where we want to interpolate f (z j ), and let Γ J+1 := Γ \ S J+1 be the point set where we will approximate. We introduce the corresponding vectors The rational function r J (z) is constructed in barycentric form r J (z) : where w j ∈ C, j = 1, . . . , J + 1, are weights. This representation already implies that the interpolation conditions r J (z j ) = f (z j ) are satisfied for w j = 0, j ∈ {1, . . . , J + 1}. The vector of weights w := (w 1 , . . . , w J+1 ) T is now chosen such that r(z) approximates the remaining data and additionally satisfies the side conditions The first condition is a normalization condition. The second condition in (3.2) ensures that r J (z) is of a type (J − 1, J). To compute w, we consider the restricted least-squares problem At the first iteration step J = 1, w ∈ C 2 is already completely fixed by the two side conditions. For J > 1, we define the matrices and rewrite the term in (3.3) as Now, the minimization problem in (3.3) takes the form To find the solution vector w = w J+1 of (3.4) approximately, we compute the right (normalized) singular vectors v 1 and v 2 of the matrix A J+1 corresponding to the two smallest singular values σ 1 ≤ σ 2 of A J+1 and take such that w J+1 2 2 = 1 and w T J+1 f S J+1 = 0. Having determined the weight vector w J+1 , the rational function r J is completely fixed by (3.1). Finally we consider the errors |r J (z) − f (z)| for all z ∈ Γ J+1 , where we do not interpolate. The algorithm terminates if max z∈Γ J+1 |r J (z) − f (z)| < for a predetermined bound or if J reaches a predetermined maximal degree. Otherwise, we find the next point for interpolation as • Compute the singular vectors v 1 and v 2 corresponding to two smallest singular values of A j+1 ; compute w : end (for) Output: N = j the order of the rational function r N S = (z j ) N +1 j=1 ∈ C N +1 is the vector of points with the interpolation property f S = (f (z j )) N +1 j=1 ∈ C N +1 is the vector of the corresponding interpolation function values w = (w j ) N +1 j=1 ∈ C N +1 is the weight vector. (3.5) which are determined by the output parameters of this algorithm.

Remark 3.2. Observe that the interpolation condition
is a "non-achievable" value for this rational interpolation problem. We will use this property of the modified AAA algorithm to detect and to reconstruct also P -periodic terms in exponential sums.
In order to rewrite r N (z) in (3.5) in the form of a partial fraction decomposition, we need to determine g 1 , . . . , g N and ρ 1 , . . . , ρ N from the output of Algorithm 3.1.
The zeros of denominatorq N (z) are the poles ρ j of r N (z) and can be computed by solving an (N + 2) × (N + 2) generalized eigenvalue problem (see [12] or [17]), Two eigenvalues of this generalized eigenvalue problem are infinite and the other N eigenvalues are the wanted zeros ρ j ofq N (z) (see [17] for more detailed explanation). We apply the following Algorithm 3.3 to the output of Algorithm 3.1.

Algorithm 3.3 (Reconstruction of parameters g j and ρ j of partial fraction representation).
Input: S ∈ Z N +1 , f S ∈ C N +1 , w ∈ C N +1 the output vectors of Algorithm 3.1.
• Build the matrices in (3.7) and solve this eigenvalue problem to find the vector

Recovery of Proper Exponential Sums
In this section we study the recovery of proper exponential sums y ∈ Y 0 N of the form where λ j are assumed to be pairwise distinct. Note that we use frequencies 2πλ j instead of λ j just for convenience. We want to recover this exponential sum from a small set of Fourier coefficients of y obtained from the Fourier series expansion on a finite interval [0, P ] ⊂ R of length P > 0. First we study the special structure of the Fourier coefficients of y(t).
the Fourier coefficients c k (φ) for k ∈ Z are given by Proof. Since the function φ(t) is differentiable on R, the Fourier coefficients of φ and the Fourier series expansion are well-defined (see, for example, [18, Chapter 1]). A simple computation yields for k = −iλP

Functions Containing only Non-P -periodic Terms
Let us assume first that −iλ j ∈ 1 P Z for j = 1, . . . , N in (4.1). We introduce Then, the Fourier coefficients of y(t) in (4.1) can by Lemma 4.1 be written as In other words, the sequence of Fourier coefficients c k (y) is already determined by a rational function , and moreover, there is a bijection between the parameters sets γ j , λ j , j = 1, . . . , N , determining y(t) in (4.1) and .
We obtain Then y is uniquely determined by 2N of these Fourier coefficients and Algorithm 3.1 (with Γ = (k) k∈Γ and f := (c k (y)) k∈Γ with c k (y) in (4.3)) terminates after N steps taking N + 1 interpolation points and provides a rational function degree at most N − 1 and q N (z) of degree exactly N , is already completely determined by 2N (independent) interpolations conditions r N (k) = c k (y), if we can assume that the rational interpolation problem is solvable at all. But solvability can be assumed since we know that the coefficients c k (y) possess the structure given in (4.4). Linear independence of the conditions follows also from (4.4), since the coefficients c k (y) cannot be presented by a rational function of smaller type than (N − 1, N ). Algorithm 3.1 chooses at the N th step a set S N +1 ⊂ Γ of N +1 indices for interpolation such that f (k) := c k (y) for k ∈ S N +1 (i.e. S = (k) k∈S N +1 and f S = (c k (y)) k∈S N +1 in Algorithm 3.1), and builds the matrix Using the known structure of c k (y) in (4.3) we find the factorization The matrix A N +1 has exactly the rank N , since all three matrix factors have full rank N . Thus, there is a right (normalized) singular vector v 1 of A N +1 to the singular value σ 1 = 0, i.e., A N +1 v 1 = 0, and the factorization above implies that also are linearly independent, it follows that all components of v 1 are nonzero. We conclude that the choice of a weight vector w = v 1 ensures that the rational function r N (z) in barycentric form (3.5) constructed by Algorithm 3.1 satisfies all interpolation conditions r N (k) = c k (y) for k ∈ S N +1 and moreover r N (n) = c n (y) for n ∈ Γ \ S N +1 by A N +1w = 0. Since this rational function r N (z) satisfies 2N +1 interpolation conditions, it follows that it coincides with the rational function r N (z) in (4.4). But r N (z) in (4.4) is of type (N − 1, N ), it follows thatw also satisfies the second side condition of minimization problem (3.4). Thus, Algorithm 3.1 will provide the weight vector w =w = v 1 at the N th iteration step.
The proof of Theorem 4.2 implies that in the considered case the kernel vector w = v 1 of the Loewner matrix A N +1 already satisfies the two side conditions. Therefore, the modification of the original AAA-algorithm that ensures that the resulting rational approximant is of type (N − 1, N ) is not needed in this case. The reconstruction of y(t) in (4.1) can now be summarized as follows.

Functions Containing also P -periodic Terms
Now, we assume that the exponential sum (4.1) contains beside non-P -periodic terms y j (t) = γ j e 2πλ j t with −iλ j ∈ 1 P Z also P -periodic terms with −iλ j ∈ 1 P Z. As seen in Lemma 4.1, each P -periodic term provides only one non-zero Fourier coefficient, i.e., c k (γ j e 2πλ j t ) = γ j δ k,−iλ j P . Therefore, we assume that the index set Γ of given Fourier coefficients c k (y) contains all integers {−iλ j P : j = 1, . . . , N } ∩Z. If c k (y) with k = −iλ j P is not provided, then the term y j (t) = γ j e 2πλ j t cannot be identified from the given data. Now the function y in (4.1) can be written as y = y (1) + y (2) , where (4.6) and N 1 < N . The part y (1) (t) is non-P -periodic with the Fourier coefficients Then The reconstruction of y(t) is now based on the observation that all but c k (y), k ∈ Σ, still have the structure of a rational function r N 1 , i.e., c k (y) = r N 1 (k) for k ∈ Z \ Σ, and y (1) (t) can be reconstructed by Algorithm 3.1 while the P -periodic part y (2) (t) can be determined in a post-processing step.
Fourier coefficients of the Fourier expansion of y on the finite interval [0, P ] ⊂ R with P > 0, and assume that the (unknown) index set Σ of non-zero Fourier coefficients of y (2) (t) is a subset of Γ. Then y (1) and y (2) can be uniquely recovered from this set of Fourier coefficients. Algorithm 3.1 (with Γ = (k) k∈Γ and f := (c k (y)) k∈Γ ) terminates after at most N + 1 steps and provides a rational function Proof. The proof employs similar ideas as the proof of Theorem 5.1 in [17] despite the different context. We therefore only sketch the main ideas of the proof.
At the (N + 1)-th iteration step, Algorithm 3.1 has chosen a set S ⊂ Γ of N + 2 indices used for interpolation. Let Γ S := Γ \ S. Then we obtain the matrix We show that A N +2 has at most rank N . Let s be the number of indices of Σ being contained in S. We consider the partial matrix A 11 of A N +2 obtained by deleting the rows and columns corresponding to indices in Σ. Then, A 11 has at least N + 2 − s ≥ N 1 + 2 columns and, similarly as in the proof of Theorem 4.8, it follows from a factorization argument that A 11 has rank N 1 . Therefore, the submatrix of A N +2 built with the columns of A N +2 that correspond to the N + 2 − s columns of A 11 has at most rank N 1 + (N − N 1 − s) = N − s, since there are at most N − N 1 − s rows of A N +2 which are deleted in A 11 . We conclude that the full matrix A N +2 has at most rank N .
Thus, Algorithm 3.1 finds two vectors v 1 and v 2 with A N +2 v 1 = A N +2 v 2 = 0 and thus, there is a weight vector w with A N +2 w = 0 satisfying the side conditions (c k ) T k∈S w = 0 and w 2 = 1. Observe that c k (y) = c k (y (1) ) for all k ∈ Γ \ Σ. Since any N 1 columns out of the columns of A 11 are linearly independent, w contains at least N 1 + 1 nonzero components corresponding to columns of A 11 and therefore, the rational function r N 1 obtained by Algorithm 3.1 indeed interpolates all Fourier coefficients of y (1) . Therefore, Algorithm 3.1 determines y (1) . In a post processing step we find all nonzero Fourier coefficients of y (2) by inspecting c k (y (2) ) = c k (y) − r N 1 (k), k ∈ Γ, and can determine y (2) . Remark 4.5. Comparing Theorems 4.2 and 4.4 we see that the reconstruction procedure may require N + 1 instead of N steps if the exponential sum also contains P -periodic components, where N is the order of the corresponding exponential sum. In practice, Algorithm 3.1 often terminates after N steps, even if P -periodic terms appear. The Fourier coefficients that are deteriorated by the P-periodic part can however also produce a Froissart doublet. Actually, if Algorithm 3.1 stops, then all indices of Fourier coefficients that corresponds to the P -periodic components, have been taken into the set S, i.e., Σ ⊂ S. Otherwise, these Fourier coefficients would cause an error in the approximation step, since they do not have the wanted rational structure.
To reconstruct the function y(t) in (4.1), we can again just apply Algorithm 4.3 to reconstruct all parameters of the non-periodic part y (1) (t). The set Σ can be found by determining all integer poles C j that are found in the second step of Algorithm 4.3. The coefficients γ j for j = N 1 + 1, ..., N can be now reconstructed via (4.7) Example 4.6. We consider the following proper exponential sum, see Figure 1,  We use P = 6, and employ the 59 Fourier coefficients c k (y 1 ), k = −29, . . . , 29, to recover y 1 (t). The function y 1 (t) = y 1 (t) contains one 6-periodic term and 5 non-6periodic terms, y Algorithm 3.1 iteratively employs the 7 Fourier coefficients c 8 (y 1 ), c −12 (y 1 ), c 9 (y 1 ), c −13 (y 1 ), c −16 (y 1 ), c 20 (y 1 ) and c 0 (y 1 ) (in this order) for interpolation before it stops after 6 iteration steps with error 3.6 · 10 −16 . Here, the Fourier coefficient c −12 (y 1 ), which contains information about y (2) 1 , is already taken. The obtained rational function r(z) is already completely determined by the remaining 6 Fourier coefficients c k (y 1 ) = c k (y 1 ) for k = −12. Indeed we observe that the second component of w ∈ C 7 vanishes, indicating that c −12 (y) is a "non-achievable" point for this rational interpolation. After skipping this vanishing term in w and in S correspondingly, we obtain a rational function r 5 (z) in the barycentric form of type (4,5) that is determined by , and the vector of Fourier coefficients with indices corresponding to the index vector S. To reconstruct the non-periodic part y 1 (t) of y(t), we apply Algorithm 4.3 as described in the previous subsection. Finally, we reconstruct the periodic part y whereλ andγ denote the reconstructed parameter vectors.

Recovery of Real Proper Exponential Sums with Real Frequencies
In this section, we consider the recovery of real proper exponential sums from a small number of its Fourier coefficients obtained for the Fourier series expansion of y on a given finite interval [0, P ]. In this case, we can derive a special algorithm in real arithmetic. We start with studying the structure of the Fourier coefficients of y.
Proof. The proof follows the same lines as the proof of Lemma 4.1.
Now, the function y in (4.9) can be rewritten as y(t) = N j=1 y j (t) with y j (t) := γ j e 2πα j t , and, according to Lemma 4.7, its Fourier coefficients satisfy the representation with real parameters C j := α 2 j P 2 , A j := γ j P α j π e α j πP sinh(πα j P ) B j := γ j π e α j πP sinh(πα j P ), for all j = 1, . . . , N . Conversely, the coefficients A j and B j , j = 1, . . . , N , uniquely determine the parameters α j and γ j of the exponential sum (4.9). We have and . (4.12) Obviously, we also have α 2 j = C j /P 2 , and C j can hence be written as We consider now the following modification of Fourier coefficients c k (y) in (4.10), which can be seen as the sample values of a rational function The reconstruction algorithm is now based on this observation. We obtain The proof of Theorem 4.8 can be derived analogously as for Theorem 4.2. In particular, Algorithm 3.1 provides a reconstruction algorithm for the rational function r N (z) that determines the modified Fourier coefficients of y. For reconstruction of the parameters of y(t) in (4.9), we thus need again to proceed with the following steps. Algorithm 4.9 (Reconstruction of the parameters α j and γ j from (4.9)). Input: Γ = (k 2 ) k∈Γ ⊂ Z, f := (c k (y)) k∈Γ withc k (y) in (4.13)), 1) Apply Algorithm 3.1 to compute a rational function r N (z) of type (N − 1, N ) from a set of L ≥ 2N + 1 Fourier coefficients with nonnegative index. Use the input data Γ = (k 2 ) k∈Γ and f := (c k (y)) k∈Γ withc k (y) in (4.13). Algorithm 3.1 then provides the vector of used (squared) indices S = (z j ) N +1 j=1 with z j = k 2 j , k j ∈ Γ ⊂ N, the vector of used modified Fourier coefficients f S = (c k j (y)) N +1 j=1 and the weight vector w = (w j ) N +1 j=1 to determine r N (z) of the form (3.5).

2) Rewrite r N in the form
i.e., extract A j , B j , C j from S, f S , and w. This is done by employing Algorithm 3.3 with the output g j = A j + iB j and ρ j = C j .

Recovery of Extended Exponential Sums
In this section we study the recovery of extended exponential sums y ∈ Y N ,

Representation of Fourier Coefficients via Rational Functions
First we study the structure of Fourier coefficients of functions y(t) in (5.1). We start with the following result regarding the expansion of one component in (5.1) into a Fourier series on [0, P ] with some given P > 0.
If −iλ ∈ 1 P Z, i.e., if there exists an n ∈ Z with −iλP = n, then Thus, for k ∈ Z, Let now −iλP ∈ Z. Then there is an n ∈ Z such that −iλP = n, and we obtain Thus we also obtain with n j ∈ N 0 , λ j ∈ C, and γ j,m ∈ C with γ j,n j = 0. Then the Fourier coefficients of y j (t) with respect to the Fourier series expansion on [0, P ] are of the following form: where for j = 1, . . . , M , and = 0, . . . , n j ,

5)
where for j = 1, . . . , M , and = 0, . . . , n j − 1, Proof. Using the formula (5.2) in Theorem 5.1, it follows for −iλ j P ∈ Z that For the case −iλ j P = k j ∈ Z the proof is similar to the one above where we use (5.3) instead of (5.2).

Recovery of Extended Exponential Sums with Non-P -periodic Terms
We consider first the simpler case where y(t) possesses only non-P-periodic terms.
where C j and A j, , j = 1, . . . , M , = 0, . . . , n j , are given as in (5.4). Then we have r N (k) = c k (y) for k ∈ Z. Moreover, the parameters λ j and γ j,m , j = 1, . . . , M , m = 0, . . . , n j , of y(t) are given by (5.8) and recursively by Proof. Using Corollary 5.2, we obtain for all k ∈ Z. Thus, the Fourier coefficients of y(t) in (5.1) can be represented by r N (z) in (5.7) such that r N (k) = c k (y) for k ∈ Z. Note that the polynomials q N (z) and p N −1 (z) determining r N (z) are coprime. In particular, y(t) is completely determined by r N (z). From (5.4) we obtain λ j = iC j P for j = 1, . . . , M , and, taking the vectors A j := (A j,0 , . . . , A j,n j ) T and γ j := (γ j,0 , . . . , γ j,n j ) T for j = 1, . . . , M , we conclude from (5.4) (5.11) The assumption C j = −iP λ j ∈ Z yields invertibility of the upper triangular matrix in the formula above, and moreover, we can compute γ j recursively by backward substitution to get (5.9).
Therefore, it suffices to determine the parameters of r N (z) to reconstruct y(t). We obtain the following generalization of Theorem 4.2.
Theorem 5.4. Let y be of the form (5.1) with N = M j=1 (1 + n j ) ∈ N, pairwise distinct λ j ∈ C \ i P Z and γ j,m ∈ C for j = 1, . . . , M , m = 0, . . . , n j , with γ j,n j = 0. Let {c k (y) : k ∈ Γ} with Γ ⊂ Z be a set of L ≥ 2N +1 Fourier coefficients of the Fourier expansion of y on the finite interval [0, P ] ⊂ R with P > 0. Then y is uniquely determined by 2N of these Fourier coefficients and Algorithm 3.1 (with Γ = (k) k∈Γ and f := (c k (y)) k∈Γ ) terminates after N steps taking N + 1 interpolation points and provides a rational function r N (z) that satisfies c k (y) = r N (k) for all k ∈ Z.
Proof. As shown in Lemma 5.3, the Fourier coefficients of y(t) in (5.1) admit a representation of the form (5.10), i.e., there exists a rational function r N (z) of type (N − 1, N ) such that we have r N (k) = c k (y). This rational function is already determined by 2N interpolation conditions. To show that Algorithm 3.1 provides this rational function (in barycentric form) after N iteration steps, we again have to inspect the generalized Cauchy , where S ⊂ Γ denotes the index set for interpolation. Using similar arguments as in the proof of Theorem 4.2, we can show that A N +1 possesses exactly rank N and that the vector w ∈ C N +1 satisfying A N +1 w = 0 is the weight vector determining the rational function r N (z) in the barycentric form (3.5).
To reconstruct y(t) in (5.1) we now can proceed as follows.
Step 1. First, we apply Algorithm 3.1 with P > 0, the period for the computation of Fourier coefficients, Γ ∈ Z L , the vector of indices of given Fourier coefficients, and f = c = (c k (y)) k∈Γ ∈ C L , the vector of given Fourier coefficients. After N iteration steps, we obtain r N (z) in barycentric form (3.5) that is determined by the vector of interpolation indices S = (k 1 , . . . , k N +1 ) T , the vector of corresponding Fourier coefficients f S = (c k (y)) N +1 =1 , and the weight vector w = (w j ) N +1 j=1 .
Step 2. The rational function r N (z) can be rewritten as a partial fraction decomposition To recover the parameters C j and A j, , we need to apply a modified version of Algorithm 3.3. As before, the parameters C j are the poles of r N (z), i.e., the zeros ofq N (z). This time, the zeros C j may come with multiplicity n j + 1 ≥ 1. Afterwards, the parameters A j, are computed using the interpolation conditions c k (y) = r N (k ), = 1, . . . , N + 1.
We summarize the modified Algorithm 5.5.

Algorithm 5.5 (Reconstruction of a partial fraction representation).
Input: S ∈ Z N +1 , c S ∈ C N +1 , w ∈ C N +1 the output vectors of Algorithm 3.1.

Recovery of Extended Exponential Sums Containing also P -periodic Terms
Now let us assume that the function (5.1) contains also P -periodic components. In this case it can be represented as y = y (1) + y (2) , where (1 + n j ) = N − N 1 denote the length and the order of y (2) (t). We define the set Σ := {−iλ j P : j = M 1 + 1, . . . , M } corresponding to the M 2 periodic terms in y (2) .
2. We consider now the representation for the Fourier coefficients of the P -periodic part y (2) (t). We denote C j := −iλ j P = k j ∈ Z for j = M 1 + 1, . . . , M , then we have Σ = {k M 1 +1 , . . . , k M }. According to (5.5), with A * j, as in (5.6). Thus all Fourier coefficients c k (y (2) ), k ∈ Z \ Σ, still have a rational structure. We consider the rational function Note that the polynomials p * N 2 −M 2 −1 and q * N 2 −M 2 do not have common zeros and r * . Taking into account (5.18) and (5.19) we conclude that 20) and thus also c k (y) = r † N −M 2 (k) for k ∈ Z \ Σ.
Finally, we study the recovery of extended exponential sums y in (5.1) that can be written as y = y (1) + y (2) , where y (1) is the non-P -periodic part defined by (5.13) and y (2) is the P -periodic part defined by (5.14). If y (2) is not just a proper exponential sum then this recovery problem essentially differs from the recovery of proper exponential sums in Section 4.2, since y (2) also possesses infinitely many nonzero Fourier coefficients if we have some n j > 0 for j ∈ {M 1 + 1, . . . , M }, see (5.18). But by Lemma 5.7, all but M 2 Fourier coefficients still have the structure of a rational function. Using this information, we can now reconstruct the function y = y (1) + y (2) as follows.
Theorem 5.8. Let y in (5.1) be of the form y = y (1) + y (2) as in (5.13) and (5.14), with γ j,m ∈ C, γ j,n j = 0, and where λ j ∈ C are pairwise distinct. Further, let M 1 < M , and let M 2 := M − M 1 be the number of P -periodic components in y (2) . Denote by Σ := {−iλ j P : j = M 1 + 1, . . . , M } ⊂ Z the index set corresponding to the frequencies of the P -periodic part y (2) . Let {c k (y) : k ∈ Γ}, with Γ ⊂ Z be a set of L ≥ 2N + 2 Fourier coefficients of the Fourier expansion of y on the finite interval [0, P ] ⊂ R with P > 0. Assume that Σ ⊂ Γ. Then y can be uniquely recovered from this set of Fourier coefficients. Algorithm 3.1 (with Γ = (k) k∈Γ and f := (c k (y)) k∈Γ ) terminates after at most N + 1 steps and provides a rational function Theorem 5.8 can be proved along the lines of Theorem 4.4. To recover y, we can now proceed as follows.
Step 1. We apply Algorithm 3.1 to Γ = (k) k∈Γ and f := (c k (y)) k∈Γ and find after at most N + 1 iteration steps the (sub)vectors S = (z j ) N +1 j=1 = (k ) N +1 =1 , c S = (c k (y)) N +1 =0 ∈ C N +1 and w = (w j ) N +1 j=1 ∈ C N +1 . If S contains integer indices k j ∈ Z of the form k j = −iλ j P , j ∈ {M 1 + 1, . . . , M }, then the corresponding components of w vanish, since these interpolation points do not possess a rational function structure. Therefore, we simply remove the zero components of w and the corresponding components of S and c S to obtain the parameter vectors that determine the rational function r † N −M 2 (z) of type (N − M 2 − 1, N − M 2 ) in barycentric form (3.5).
Step 2. We apply now a modification of Algorithms 5.5 as follows. From the definition of r † N −M 2 (z) in Lemma 5.7 it follows that we find N − M 2 poles (counting also multiplicities): each C j ∈ Z is a pole with multiplicity n j + 1 and belongs to rational function r N 1 (z) that corresponds to the non-P -periodic part y (1) . Each C j ∈ Z is a pole with multiplicity n j and belongs to the rational function r * N 2 −M 2 (z) which corresponds to the P -periodic part y (2) . Taking into account this information we find M pairwise distinct C j that are the poles of the r † N −M 2 (z) with multiplicities n j + 1 for non-P -periodic components and n j for P -periodic components. The set Σ can be easily determined by Σ = Step 4. We compute the coefficients γ j,m , j = M 1 + 1, . . . , M , m = 0, . . . , n j , by solving the linear system (5.22) by using values A * j, , j = M 1 + 1, . . . , M , = 0, . . . , n j − 1, and the vectorc Σ = (c k (y)) k∈Σ , withc k (y) in (5.21), which can be determined using When the exponential sum contains a proper P -periodic part, the corresponding component −iλ j P does not appear as a pole because in this case n j = 0. Therefore, we apply a technique similar to the one that we used in Section 4.2 in order to detect the corresponding frequency λ j , namely we compere the Fourier coefficients and the values of the rational function constructed by Algorithm 3.1. The coefficient γ j can be found then via (4.7).
4 of y 4 . After omitting this first term in w and in the corresponding index in the vector S, we obtain the rational function r 9 (z) of type (8,9)  We compute the vector C of poles, the number M of pairwise distinct poles, the vector n of multiplicities of poles and the parameters A j, for = 0, . . . , n j , j = 1, . . . , M 1 , A * j, for = 0, . . . , n j − 1, j = M 1 + 1, . . . , M , as it is explained in Steps 2 and 3 above. We obtain the poles of r 9 (z) We assume that two poles C j 1 and C j 2 , j 1 = j 2 are equal if |C j 1 − C j 2 | < 0.001 holds, and find M = 3 different poles: C 1 = −5.839999999992412 + 0.800000000001007i with multiplicity n 1 = 3, C 2 = −25.436980952935279 − 0.399999999999682i with multiplicity n 2 = 2, and C 3 = 12.000000000000027 − 0.000000000000030ii with multiplicity n 3 = 4. The poles C 1 and C 2 corresponds to the non-8-periodic part of the exponential sum y 4 , therefore they appear in C with multiplicities n 1 + 1 and n 2 + 1 respectively. The pole C 3 corresponds to the 8-periodic part of y 4 therefore it comes with multiplicity n 3 . To determine the values for the poles C 1 , C 2 and C 3 we have taken the average values. Finally, we reconstruct λ j , j = 1, 2, 3, and γ j, , j = 1, 2, = 0, . . . , n j , as described in Step 3 and The recovery errors are λ − λ ∞ = 9.56 · 10 −13 , γ − γ ∞ = 3.11 · 10 −10 .

Conclusions
In Sections 4 and 5 we have considered the recovery of proper and extended exponential sums. We have shown that sums of the form (1 + n j ) = N , γ j,m ∈ C and pairwise distinct λ j ∈ C can always be recovered from at most 2N + 2 Fourier coefficients of a Fourier expansion on a finite interval. Numerical stability of this procedure essentially depends on the numerical stability of the underlying AAA algorithm for rational approximation of these Fourier coefficients. Observe that the considered model also covers sums of the form y 0 (t) = with δ m ∈ C, where the polynomial term occurs for λ = 0. Obviously, λ 0 = 0 ∈ i P Z for any P > 0, and we need to apply the procedure described in Section 5.3 for its recovery. Our models also cover real signals of the form with γ j ∈ R \ {0}, α j ∈ R and b j ∈ [0, 2π) considered in [17]. The algorithms in [17] for recovery of y(t) in (6.1) are also based on rational approximation of the Fourier coefficients, but have essentially used the additional information that all parameters are real, and are therefore different from the algorithms for the complex case considered here. Moreover, because of a different representation of the rational function in form of partial fraction decomposition the algorithm in [17] requires to use modified coefficientsc k (f ) = Re c k (f )+ i k Im c k (f ) instead of c k (f ). Our approach can now be also applied to the recovery of where γ j (t) are polynomial of finite degree.