Exact reconstruction of sparse non-harmonic signals from their Fourier coefficients

In this paper, we derive a new reconstruction method for real non-harmonic Fourier sums, i.e., real signals which can be represented as sparse exponential sums of the form f(t)=∑j=1Kγjcos(2πajt+bj)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$f(t) = \sum _{j=1}^{K} \gamma _{j} \, \cos (2\pi a_{j} t + b_{j})$$\end{document}, where the frequency parameters aj∈R\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$a_{j} \in {\mathbb {R}}$$\end{document} (or aj∈iR\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$a_{j} \in {\mathrm i} {\mathbb {R}}$$\end{document}) are pairwise different. Our method is based on the recently proposed numerically stable iterative rational approximation algorithm in Nakatsukasa et al. (SIAM J Sci Comput 40(3):A1494–A1522, 2018). For signal reconstruction we use a set of classical Fourier coefficients of f with regard to a fixed interval (0, P) with P>0\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$P>0$$\end{document}. Even though all terms of f may be non-P-periodic, our reconstruction method requires at most 2K+2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$2K+2$$\end{document} Fourier coefficients cn(f)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$c_{n}(f)$$\end{document} to recover all parameters of f. We show that in the case of exact data, the proposed iterative algorithm terminates after at most K+1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$K+1$$\end{document} steps. The algorithm can also detect the number K of terms of f, if K is a priori unknown and L≥2K+2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$L \ge 2K+2$$\end{document} Fourier coefficients are available. Therefore our method provides a new alternative to the known numerical approaches for the recovery of exponential sums that are based on Prony’s method.


Introduction
Classical Fourier analysis methods provide for any real square integrable signal f (t) a Fourier series representation on a given interval (0, P), P > 0, of the form with Fourier coefficients c n ( f ) = 1 P P 0 f (t) e −2π int/P dt, n ∈ Z, and α n = 2 |c n ( f )|, β n = atan2(Im c n ( f ), Re c n ( f )) for n ≥ 0, where atan2 denotes the modified inverse tangent, see e.g. [23,Chapter 1]. If f is P-periodic and differentiable, then its Fourier series (1.1) converges uniformly to f . However, if f is smooth but not P-periodic, then the P-periodization of f "forced" by the Fourier series representation in (1.1) usually leads to a discontinuity at the interval boundaries t = 0 and t = P, respectively, and thus to a slow decay of the Fourier coefficients. In applications, it frequently happens that a signal is only given on an interval of length P, where it appears to be non-periodic, even if f may be periodic with a certain period P 1 which is not of the form P/n for some positive integer n. Considering for example the signal f (t) = cos(2π √ 2t) + cos(2π √ 3t), (1.2) which contains only two different frequency parameters, we observe that this signal is non-periodic with regard to any interval (0, P). For P = 1, the corresponding Fourier series is given by Thus, the question occurs, how to reconstruct a non-harmonic Fourier sum, i.e., how to compute the much more informative representation (1.2) directly from suitable measurements of f . Contents of this paper The goal of this paper is to reconstruct non-harmonic Fourier sums f of the form from a finite number of its classical Fourier coefficients c n ( f ) corresponding to a Fourier series of f on (0, P). Here, we assume that K ∈ N, γ j ∈ (0, ∞), (a j , b j ) ∈ (0, ∞)×[0, 2π), and that the frequency parameters a j are pairwise distinct. As we will show in the sequel, the restrictions made for K , γ j , a j , and b j will ensure uniqueness of the presentation (1.3). Note that f in (1.3) admits real (nonnegative) frequency parameters a j and therefore essentially generalizes usual trigonometric polynomials. The example in (1.2) with frequencies √ 2 and √ 3 is covered by our model (1.3) taking K = 2, γ 1 = γ 2 = 1, a 1 = √ 2, a 2 = √ 3, and b 1 = b 2 = 0. Observe that the function f in (1.3) is only P-periodic for some P > 0 if all parameters a j can be written in the form a j = c q j with a positive constant c ∈ R and rational numbers q j ≥ 0, i.e., only in this case, there exists a real number P > 0 such that f (t + P) = f (t) for all t ∈ R.
In Sect. 2 we show that our model (1.3) is well-defined, i.e., that all parameters K , a j , b j , γ j , j = 1, . . . , K , are (with the given restrictions) uniquely determined for a non-harmonic Fourier sum f . If all terms f j of f are non-P-periodic, i.e., if all frequency parameters a j > 0 in (1.3) satisfy that a j / ∈ 1 P N, then we can show that the modified Fourier coefficientsc n ( f ) := Re c n ( f ) + i n Im c n ( f ), n > 0, have a special structure of the form r (n 2 ), where r (z) is a rational function of type (K − 1, K ). Conversely, r (z) already provides all information to find the parameters determining f in (1.3).
Section 3 is devoted to the new reconstruction method. Using a modification of the recently proposed AAA algorithm for iterative rational approximation, see [20], we compute r (z) from a set of (at least) 2K + 1 classical Fourier coefficients of f . Then a partial fraction decomposition together with a non-linear bijective transform provides the wanted parameters in (1.3). Numerical stability of the rational approximation algorithm is ensured using a barycentric representation of the numerator and the denominator polynomial of r (z). Compared to other rational interpolation algorithms, a further important advantage of the employed modified AAA algorithm is that we do not need a priori knowledge on the number K of terms in (1.3) but can determine K in the iteration process, supposed that L ≥ 2K + 1 Fourier coefficients are available.
We show in Sect. 4, that a signal f with K non-P-periodic terms f j as in (1.3) can theoretically be determined from 2K Fourier coefficients c n ( f ) with n ∈ ⊂ N if K is known beforehand. If L ≥ 2K + 1 Fourier coefficients are available, then our method based on the AAA algorithm always provides the wanted rational function r (z) after K iteration steps (and using 2K +1 modified Fourier coefficients). Thereby, also the number K of terms in (1.3) is determined.
In Sect. 5, our new reconstruction method is generalized to the case that f in (1.3) also contains P-periodic terms f j with frequencies a j ∈ 1 P N. It turns out that there is no a priori information needed about possibly occurring P-periodic terms of f . In this case, we first compute the rational function r (z) that determines the non-Pperiodic part of f , where we again employ the modified AAA algorithm from Sect. 3. Afterwards, the P-periodic terms of f can be found in a post-processing step, if all c n j ( f ) with n j = Pa j are contained in the given set of Fourier coefficients. In particular, we show that f can be always completely recovered from L ≥ 2K + 2 Fourier coefficients, i.e., our method determines the number K of terms as well as all parameters γ j , a j and b j , j = 1, . . . , K , and automatically recognizes frequencies a j ∈ 1 P N that correspond to P-periodic terms. In Sect. 6, we generalize the model in (1.3). Beside supposing (a j , b j ) ∈ (0, ∞) × [0, 2π), we can also admit parameters (a j , b j ) ∈ i(0, ∞)×iR. By cos(ix) = cosh(x), this leads to terms of the form γ j cosh((−i)(2πa j + b j ) in (1.3). The considered generalization still admits the same rational structure of Fourier coefficients and can therefore be treated similarly as (1.3) with the proposed reconstruction method.
Finally we provide some numerical experiments. The Matlab implementation of our reconstruction algorithm is provided at the Software section of our homepage http://na.math.uni-goettingen.de. Related literature Our model (1.3) can be viewed as a sparse expansion into exponentials with 2K terms via Euler's identity, i.e., The problem of parameter identification of exponential sums from suitable input values appears in many applications as e.g. sparse array of arrival estimation [18], image super-resolution [32], nondestructive testing [5], system theory for parametrized model reduction [15], and phase retrieval [1]. Exponential sums have been extensively studied within the last years, based on Prony's method and its relatives, see e.g. [2,4,7,11,21,[24][25][26][27]29,31,34]. To reconstruct an exponential sum via Prony's method, one usually employs equidistant function values f (t 0 + h ), = 0, . . . , L, and the number of given values should be at least 2M, where M is the number of terms in the exponential sum. In our case, the number of exponential terms is M = 2K , but the symmetry properties can be exploited such that the samples f (h ), = 0, . . . , 2K − 1, are theoretically sufficient for the recovery of f in the noise-free case, see e.g. [24]. However, Prony's method involves Hankel or Toeplitz matrices with possibly high condition numbers, and therefore requires a very careful numerical treatment. Our new method for reconstruction of signals of the form (1.3) in this paper is based on rational approximation and can be seen as a good alternative to the Prony reconstruction approach. Another way to look at the model (1.3) is to view it as a special case of a signal decomposition into so-called intrinsic mode functions (IMFs) in adaptive data analysis, see [14]. Empirical mode decomposition (EMD) is based on a model that decomposes the signal f into K IMFs, with nonnegative envelope functions γ j (t) and so-called instantaneous phase functions φ j (t), see e.g. [14]. As already pointed out in [8], despite certain restrictions, as e.g. that γ j (t) and φ j (t) are smooth with γ j (t) ≥ 0 and φ j (t) ≥ 0 for t ∈ R, a representation of the form (1.4) is far from being unique. For example, the function f (t) in (1.2) has the form (1.4) with K = 2, constant functions γ 1 (t), γ 2 (t), and with φ 1 (t) = 2π √ 2t and φ 2 (t) = 2π √ 3t. However, f (t) in (1.2) can also be written as a single IMF in − The non-uniqueness of the model (1.4) often prevents a simple interpretation of the obtained decomposition. Compared to (1.4), the main advantages of the non-harmonic Fourier sum (1.3) are that the representation of f is unique, and, that the model (1.3) has a direct physical interpretation, similarly to classical Fourier sums.
There are also other approaches to represent signals by adaptive generalized Fourier sums using the so-called Takenaka-Malmquist basis, an adaptive orthonormal basis, see [22,28]. While the greedy algorithm in [28] only slightly improves the signal approximation compared to classical Fourier sums, it has been shown in [22], that strong decays of adaptive Fourier expansions can be achieved, if the sequence of classical Fourier coefficients of a signal can be well approximated using a short exponential sum. Our approach in the current paper is somehow vice versa, the (modified) Fourier coefficients of f are represented by rational functions in order to reconstruct the special exponential sum f .
While we focus on signal reconstruction in the current paper, there remains the question of (almost) optimal signal approximation by non-harmonic Fourier sums, which we will study in the future. Obviously, each square integrable signal in (0, P) can be arbitrarily well approximated be a non-harmonic Fourier sum (1.3) for K → ∞, since it is a direct generalization of classical Fourier sums. However, approximation rates for signals in certain function classes are not completely known so far. Research on non-harmonic Fourier series particularly focussed on functional analytic questions, see e.g. [33]. In particular, it has been shown that {e 2π ia j · } j∈Z forms a Riesz basis in L 2 ([0, 1]) for a given increasing sequence {a j } j∈Z if and only if |a j − j| < 1/4 for j ∈ Z, see [16], while completeness of this function system is ensured for |a j | ≤ | j| + 1/4, where a j can even be complex, see [19]. Note that for finite non-harmonic Fourier sums as in (1.3) we do not need any further assumption on the distribution of (pairwise distinct) frequencies to ensure the uniqueness of the presentation.

Unique representation of non-harmonic signals
We will show that the model in (1.3) is well-defined, since all occurring parameters K , γ j , a j and b j , j = 1, . . . , K , are uniquely determined for a function f given on an interval with positive length. More precisely, we can show the following: . Then h(t) = 0 for all t ∈ T , and h(t) has the structure By assumption, the number L of distinct frequency parameters x j in the representation (2.1) of h satisfies L ≥ max{K , M}, and we can rewrite h(t) in the form wherex ∈ {x j : j = 1, . . . , K + M} are now pairwise distinct. Ifx occurs only once in the set {x j : j = 1, . . . , K + M}, sayx = x j , then α = μ j cos(y j ) β = μ j sin(y j ).
Remark 2.2 1. Theorem 2.1 also shows that a function f of the form (1.3) with the given restrictions on the parameters γ j , a j , b j cannot vanish on any interval T ⊂ R with positive length. 2. Since cos(2πa j t) and thus also the functions f , g and h in the proof of Theorem 2.1 are analytic functions, the assumption h(t) = 0 for t ∈ T already implies h(t) = 0 for t ∈ C.
3. The linear independence of the system {e 2π ix t , e −2π ix t : = 1, . . . , L} in the proof of Theorem 2.1 also follows from the fact that an exponential sum of the form h(t) in (2.4) can appear as a general solution of a linear difference equation of order 2L with constant coefficients, see e.g. [3]. 4. Observe that the function model can simply be extended by adding a constant component f 0 (t) = ±γ 0 = γ 0 cos(2πa 0 t + b 0 ) with γ 0 > 0, a 0 = 0 and either b 0 = 0 for f 0 (t) > 0 or b 0 = π for f 0 (t) < 0. This extended model also satisfies the assertion of Theorem 2.1. If we admit beside (a j , b j ) ∈ (0, ∞) × [0, 2π) also (a j , b j ) = (0, 0) and (a j , b j ) = (0, π), the proof of Theorem 2.1 can be suitably modified.
Proof Since φ(t) = γ cos(2πat + b) is a differentiable function, its restriction onto the interval (0, P) is also differentiable. Thus the Fourier expansion of φ converges pointwise for all t ∈ (0, P), see [23,Chapter 1]. For the real part of c n (φ) we obtain with cos x cos y = 1 2 (cos(x + y) + cos(x − y)) Assuming that a = n P it follows with sin x − sin y = 2 sin( x−y 2 ) cos( x+y 2 ) For a = n P > 0, the function φ is P-periodic, and we simply find This is also achieved from (2.6) by taking the limit with the rule of L' Hospital, lim a→ n P Re c n (φ) = lim a→ n P Pγ a π(a 2 P 2 − n 2 ) sin(aπ P) cos(aπ The formula (2.7) for the imaginary part of c n (φ) can be derived analogously.

Representation of Fourier coefficients by rational functions
We consider now where the a j are assumed to be pairwise distinct. As shown in Theorem 2.3, we have for n ∈ Z and a j /

9)
A j := − Pγ j a j π sin(a j π P) cos(a j π P + b j ), (2.10) Note that C j is real and positive for real values a j . We show that A j , B j , C j , j = 1, . . . , K , completely determine all Fourier coefficients c n ( f ) and thus f .
Proof We can assume that C j > 0. Then (2.9) implies that a j = √ C j P > 0. Further, taking the weighted sum A 2 j + a 2 j P 2 B 2 j using (2.10) and (2.11) we obtain Since γ j > 0, we can determine γ j uniquely. Inserting the found representations for a j and γ j into (2.10) and (2.11), we conclude for as well as and thus sign(A j −B j C j cot( C j π)) = sign(sin(b j )). If B j = 0 then sin(a j π P+ b j ) = 0 and thus b j ∈ {−π C j mod π, −π C j mod π + π }.
Remark 2.6 1. Since we had assumed that a j / ∈ 1 P N and in particular a j = 0, it follows that C j = 0. As seen from Theorem 2.3, we always have C j > 0 for the considered model. In Sect. 6, we will generalize the model to treat also the case C j < 0 which leads to generalized expansions involving also cosine hyperbolic terms.
2. Please note that the function acot used in our Matlab implementation is defined differently, it maps to [−π/2, π/2).
In particular, the real part Re c n ( f ) and the imaginary part Im c n ( f ) can for all n ∈ Z be written as Re is a monic polynomial of degree K , and where is a (complex) algebraic polynomial of degree (at most) K − 1. In other words, for i.e., the modified Fourier coefficientc n ( f ) can be represented by a rational function of type (K − 1, K ), evaluated at z = n 2 , where q K (z) in (2.13) and p K −1 (z) in (2.14) are coprime.

Modified AAA algorithm for sparse signal representation
We want to exploit the special structure of the Fourier coefficients of functions f which are built by function atoms of the form f j (t) = γ j cos(2πa j t + b j ) in order to study the following problem: How can we reconstruct a function f of the form (1.3) from a given set of its Fourier coefficients in a numerically stable and efficient way? We assume here that only the structure of f is known, i.e., we need to recover the number K of terms in (1.3) as well as the parameters γ j , a j , b j for j = 1, . . . , K .
For the reconstruction process we need to keep in mind that the rational representation of Fourier coefficients in (2.12), or in (2.15) respectively, is only valid for the non-P-periodic terms f j of f , i.e., for a j / ∈ 1 P N 0 . If f contains components f j (t) = γ j cos(2πa j t + b j ) with a j P = n j ∈ N 0 , then, as shown in Theorem 2.5, these components will provide only one non-zero Fourier coefficient c n j ( f ) (with nonnegative index), which destroys the rational function structure (2.15) at z = n 2 j . Therefore, our approach consists of two parts. In the first step, we will reconstruct the non-P-periodic part of f , and in a second step, we will determine possible Pperiodic terms of f that can be obtained from the set of Fourier coefficients.
To reconstruct the non-P-periodic part of f in (1.3) we will extensively use the structure of the Fourier coefficients c n ( f ) found in (2.15) and employ a modification of the recently proposed AAA algorithm in [20]. Differently from other rational interpolation algorithms, the modified AAA algorithm provides essentially higher numerical stability and enables us to determine also the order of the rational approximant which is at the same time the number K of terms in (1.3). The AAA algorithm can be seen as a method for rational approximation, where certain values of a given function are interpolated, while other given values are approximated using a least squares approach. With this algorithm, we will determine the rational function in (2.15) by interpolation or approximation of the given modified Fourier coefficientsc n ( f ). The algorithm works iteratively, where at each iteration step the degree of the polynomials determining the rational function grows by 1, and a next interpolation value is chosen at the point, where the error of the rational approximation found so far is maximal. The algorithm terminates if either the error at all given points considered for approximation is less than a given bound or if a certain fixed degree of the rational function is reached. If the rational function is found, we can extract the parameters A j , B j , C j and then finally obtain the wanted representation with parameters γ j , a j , b j of the non-P-periodic part of f from Theorem 2.5.
Numerical stability of the AAA algorithm is ensured by a barycentric representation, as already considered in [13,30] and exploited also in [12], to compute rational minimax approximations. In [20], the AAA-algorithm is presented for rational functions r (z) = p(z)/q(z), where the polynomials p and q have the same degree. Therefore, we need to modify the approach for our purpose, similarly as proposed in [12]. A side effect of the linearization procedure within the algorithm is that unattainable interpolation points lead to vanishing weight components, see [30]. This behavior of the algorithm enables us to determine possible P-periodic terms of f in a postprocessing step, where we need to inspect all Fourier coefficients that cannot be well approximated by the found rational function.
If f does not contain P-periodic terms and the given Fourier coefficients of f are exact, then we will be able to determine f uniquely from 2K + 1 Fourier coefficients. This will be shown in Sect. 4. Otherwise, we will need 2K + 2 Fourier coefficients to recover f , where all P-periodic terms are simply determined in a postprocessing step, see Sect. 5.

Rational interpolation using barycentric representation
Let us assume now that we are given a set of Fourier coefficients c n ( f ), n ∈ ⊂ N of the function f of the form (1.3) with L := # ≥ 2K +1. We assume first that all terms of f are non-P-periodic, such that we obtain the rational structure ofc n =c n ( f ) as given in (2.15). We want to find a rational function r K are satisfied. Assuming that the given modified Fourier coefficientsc n of f in (1.3) are exact, we will show in Sect. 4 that r K (z) will be the wanted rational function in (2.15) that determines f . As in [13,17,20], we will use the barycentric representation of r K (z) = where w j , j = 1, . . . , K + 1, are nonzero weights, and where n j ∈ for j = 1, . . . , K + 1. Here, n 2 j , j = 1, . . . , K + 1, cannot occur as poles of r K (z), since the poles C j in (2.13) satisfy C j = a 2 j P 2 / ∈ N by assumption. It can be simply observed thatp K (z)/q K (z) is indeed a rational function of type (K , K ). In order to ensure that r K (z) is of the wanted type (K −1, K ), we require the additional condition Let S K +1 := {n 1 , . . . , n K +1 } be the index set, where for nonzero weights w k the interpolation conditions are already satisfied. To determine r K (z) using (3.1) we still need to fix the normalized weight vector w = (w 1 , . . . , w K +1 ) T . According to [20], this is done by solving a least squares problem in order to minimize the error where beside w 2 = 1, in our case we will incorporate the side condition K +1 j=1 w jcn j = 0 to ensure that r K (z) is of type (K − 1, K ). The algorithm is described in the next two sections and closely follows the approach in [20] with the modification that we want to get a rational function of type (K − 1, K ) instead of type (K , K ).

Initialization of the modified AAA algorithm
We start by initializing the modified AAA-algorithm as follows. First, we choose the two given modified Fourier coefficientsc n 1 ,c n 2 with largest modulus (where n 1 = n 2 and n 1 , n 2 ∈ ) for interpolation and compute a rational function r 1 (z) of type (0, 1) that interpolatesc n 1 at z = n 2 1 andc n 2 at z = n 2 2 . We determine the rational function , such that r 1 (n 2 1 ) =c n 1 and r 1 (n 2 2 ) =c n 2 holds. A barycentric form of r 1 (z) as in (3.1) is given by with the (complex) weights satisfying |w 1 | 2 + |w 2 | 2 = 1 and w 1cn 1 + w 2cn 2 = 0. Obviously,p 1 (z) andq 1 (z) are themselves rational functions of type (at most) (1, 2). The condition w 1cn 1 +w 2cn 2 = 0 ensures that the polynomial is only constant (and not linear).
In order to decide, which interpolation point should be taken at the next iteration step, we consider the error |r 1 (n 2 ) −c n | for all n ∈ \ {n 1 , n 2 }. Following the lines of [20], we use the notation S 2 := {n 1 , n 2 } ⊂ and 2 := \ S 2 . Let the Cauchy matrix C 2 be given by C 2 := 1 n 2 −n 2 j n∈ 2 ,n j ∈S 2 with 2 columns and L − 2 rows.
. As in the original AAA algorithm, the remaining freedom to choose w is now used in order to approximate the (modified) Fourier coefficientsc n by r J −1 (n 2 ) for n ∈ J applying a (linearized) least squares approach. Observing that r J −1 (z)q J −1 (z) =p J −1 (z), we consider the minimization problem Similarly to [20], we define the matrices Then we can write such that the minimization problem in (3.5) takes the form To solve the minimization problem (3.7) approximatively, we compute the right (normalized) singular vectors v 1 and v 2 of A J corresponding to the two smallest singular values σ 1 ≤ σ 2 of A J and take a linear combination w J = μ 1 v 1 + μ 2 v 2 such that w J 2 = 1 and w T Jc S J = 0. These conditions are satisfied for Having determined the weight vector w J , the rational function r J −1 is completely fixed by (3.3). Now, we consider the errors |r J −1 (n 2 ) −c n | for all n ∈ J , where we do not interpolate. The algorithm terminates if max n∈ J |r J −1 (n 2 ) −c n | < for a predetermined bound or if J − 1 reaches a predetermined maximal degree. Otherwise, we find the next index for interpolation as The values r J −1 (n 2 ) =p J −1 (n 2 ) q J −1 (n 2 ) can be simply computed by the vectors as suggested in [20], where . * indicates pointwise multiplication. We summarize the modified AAA algorithm to compute the vectors w K +1 and the index vector S K +1 := (n 1 , . . . , n K +1 ) T determining the rational function r K (z) = 3) that interpolates the given modified Fourier coefficientsc n if n is a component of S K +1 , and approximatesc n otherwise. 1. Find the two componentsc n 1 ,c n 2 ofc with largest absolute values.
Update S,c S , ,c by adding n 1 , n 2 as components of S and deleting these components in , addingc n 1 ,c n 2 inc S and deleting them inc.
) v 2 and normalize w = 1 w 2 w. 5. Compute p = C m (w. * c S ) , q = C m w, and r = p./q ∈ C L−m . 6. If r −c ∞ < tol then stop. end(for) where S is the vector of indices,c S the corresponding coefficient vector of interpolation values, and w the weight vector to determine the rational function r K via (3.3).
As we will show in Sect. 4, it will be sufficient to employ L = 2K + 1 Fourier coefficients if the function f has no P-periodic components f j . In Sect. 5 we will prove that the algorithm can be also applied if f contains P-periodic terms f j with frequencies a j ∈ 1 P N. In this case we need L = 2K + 2 Fourier coefficients, and the rational output function will only determine the Fourier coefficients of the nonperiodic part of f . In both cases, for L ≥ 2K + 1 the Algorithm 3.2 will terminate after K or K + 1 iteration steps. Algorithm 3.2 involves SVDs of the matrices A m of size (L − m) × m, m = 3, . . . , K + 1 (or m = 3, . . . , K + 2), to compute their singular vectors v 1 and v 2 yielding an overall complexity of O(L K 3 ) flops, see also [20]. This is modest since K is usually small.

Partial fraction representation of the rational approximant
Assume that we have found the rational approximant r K (z) after K iteration steps. In other words, we have now given the vector of interpolation indices (n 1 , . . . , n K +1 ) T and the weight vector As we will show in Sect. 4, r K (z) determined by Algorithm 3.2 coincides with the desired rational function in (2.15) for exact data, and, in particular, w j = 0 for j = 1, . . . , K + 1. To extract the wanted parameters A j , B j , and C j in (2.12) from {n 1 , . . . , n K +1 } and w K +1 , we need to rephrase r K (z) in the form The parameters C j , j = 1, . . . , K , are the zeros of the rational functionq K (z) in (3.9), since z = n 2 j cannot occur as poles of r K (z) due to the assumption C which is by (2.9) equivalent with to a j / ∈ 1 P N 0 . To compute the zeros ofq K (z) we again draw from the results in [20] or [17] and consider the generalized eigenvalue problem Observe that lim z→±∞q (z) = 0 causes two infinite eigenvalues that we are not interested in. The other K eigenvalues are the wanted zeros ofq K (z). This can be simply seen by taking the eigenvectors v λ corresponding the eigenvalues λ of the form Having found the K zeros C j of this eigenvalue problem, we obtain from the interpolation conditions r K (n 2 ) =c n for = 1, . . . , K + 1 with (3.10) the linear equation system in order to determine A j + iB j , j = 1, . . . , K . We summarize the reconstruction of the parameter vectors (A j ) K j=1 , (B j ) K j=1 , and (C j ) K j=1 from the output of Algorithm 3.2 in Algorithm 3.3.

Exact reconstruction using Algorithm 3.2
Assume that the Fourier coefficients c n ( f ) of a function f in (1.3) with K components f j are given for n ∈ ⊂ N with # = L ≥ 2K + 1, where we suppose that a j / ∈ 1 P Z for all components f j of f . We show that all parameters determining f can be uniquely reconstructed from 2K Fourier coefficients c n ( f ) if K is known beforehand.
Moreover, if L ≥ 2K + 1 Fourier coefficients c n ( f ) are given, then Algorithm 3.2 terminates after K steps (i.e., taking K + 1 interpolation points). Thereby, also K is determined if it is not a priori known.

Theorem 4.1 Let f be of the form (1.3)
with K ∈ N, γ j ∈ (0, ∞), and (a j , b j ) ∈ (0, ∞) ×[0, 2π) j = 1, . . . , K . Further, let a j be pairwise different and a j / ∈ 1 P N for a given P > 0. Assume that we have a given set of classical Fourier coefficients c n ( f ), n ∈ ⊂ N (with regard to period P) with L = # ≥ 2K + 1. Then f is uniquely determined by 2K of these Fourier coefficients. If L ≥ 2K + 1, then Algorithm 3.2 terminates after K steps, takes K + 1 interpolation points and determines the rational function r K (z) = K j=1 A j +iB j z−C j in (2.15) that interpolates allc n ( f ), n ∈ N, exactly. Moreover, K as well as γ j , a j , and b j , j = 1, . . . , K , are uniquely determined by r K . Proof 1. From (2.15) it follows that there exists a rational function r K (z) = p K −1 (z) q K (n 2 ) for all n ∈ N with q K (z) in (2.13) and p K −1 (z) in (2.14). In particular, p K −1 (z) and q k (z) are coprime.
First, we show that r K (z) is uniquely determined by 2K coefficientsc n j ( f ), n j ∈ , j = 1, . . . , 2K . By (2.14), p K −1 (z) has at most degree K − 1. Further, q K (z) has exactly degree K by (2.13). We use the notation p K −1 (z) = K −1 r =0 p r z r and q K (z) = z K + K −1 r =0 q r z r . Then the interpolation conditions and with (p T , q T ) = ( p 0 , . . . , p K −1 , q 0 , . . . q K −1 , 1) T ∈ R 2K +1 . The kernel of W has at least dimension 1, and by construction, the vector (p T , q T ) generating the rational function r K (z) = p K −1 (z)/q K (z) satisfies (4.1). We show that the kernel of W has exactly dimension 1. Suppose to the contrary that there exists another vector in the kernel of W being linearly independent of (p T , q T ). Then we also find a kernel vector, whose last component vanishes, i.e., of the form (p T ,q T ) = (p 0 , . . . ,p K −1 ,q 0 , . . .q K −2 ,q K −1 , 0) T . Thus, there exist polynomials p K −1 andq K −1 of at most degree K − 1 satisfying Using the known structure ofc n j ( f ) = p K −1 (n 2 j )/q K (n 2 j ) in (2.15), we obtain Since the degree of the involved polynomial products is at most 2K − 1, it follows that is a monic polynomial of degree K and the polynomials p K −1 (z), q K (z) are coprime, and we conclude thatq K −1 (z) possesses all K linear factors of q K (z). This leads to a contradiction, sinceq K −1 (z) has degree at most K − 1. Therefore, there exists only one normalized solution vector of the form (4.1), which is already uniquely defined by 2K modified Fourier coefficients of f , and this solution vector (p T , q T ) T determines the rational polynomial r K (z) = p K −1 (z) q K (z) that satisfies all 2K interpolation conditions. 2. We show now that Algorithm 3.2 leads to this unique solution r K (z) after K steps. Assume that contains L ≥ 2K + 1 indices. At the K -th iteration step we have chosen a set S K +1 of K + 1 pairwise different indices n ∈ for interpolation and start with the ansatz such that the interpolation conditions r (n 2 ) =c n are already satisfied for = 1, . . . , K + 1, if w = 0. Let K +1 := \ S K +1 . Now, Algorithm 3.2 determines the weight vector w = (w ) K +1 =1 as a linear combination of the two right singular vectors v 1 and v 2 of corresponding to the two smallest singular values σ 1 ≤ σ 2 with side conditions w 2 = 1 and K +1 =1 w c n = 0. From (2.12) and (2.15) it follows that where the two Cauchy matrices have full rank K and where the entries −(A j + iB j ) of the diagonal matrix do not vanish for all j = 1, . . . , K . Thus, A K +1 has exactly rank K and therefore a kernel of dimension 1. Let v 1 be the normalized right singular vector of A K +1 to σ 1 = 0, i.e., A K +1 v 1 = 0. The factorization of A K +1 also implies We observe that, if v 1 = 0 had one or more vanishing components, then K columns of the Cauchy matrix 1 n 2 −C j j=1,...,K ,n ∈S K +1 would be linearly dependent. But this is not possible, since the interpolation points n ∈ S K +1 are pairwise distinct and therefore any K columns of this Cauchy matrix are linearly independent. Thus, all components of v 1 are nonzero.
If we determine the rational polynomial r (z) in (4.2) with w := v 1 , then it follows r (n 2 ) =c n for n ∈ S K +1 by construction, since all weight components are nonzero. Moreover, the condition i.e., it follows that r (n 2 ) =c n for all n ∈ K +1 . Thus, the first part of the proof implies that the obtained rational function r (z) in (4.2) coincides with r K (z) = p K −1 (z) q K (z) defined by (2.13)-(2.15), and therefore has to be of the wanted type (K − 1, K ) since it is already uniquely defined by the interpolation conditions. We conclude that w = v 1 already satisfies the side condition K +1 =1 w c n = 0 and is therefore the weight vector computed at the K -th iteration step of Algorithm 3.2.

How to proceed if the function contains P-periodic terms
Let us assume that the function there are also periodic components with a j ∈ 1 P N. Now, we will study the question, how Algorithm 3.2 proposed in Sect. 3 behaves in this case and how we can reconstruct f (t).
We assume that the index set of given Fourier coefficients of f contains the integer n j , if a j = is non-P-periodic such that the modified Fourier coefficientsc n (φ 1 ) = Re c n (φ 1 ) + i n Im c n (φ 1 ) are of the form similarly as in (2.15). The P-periodic part of f (t), possesses only K 2 := K − K 1 nonzero Fourier coefficients with non-negative index. Let := {P a K 1 +1 , . . . , P a K } ⊂ N denote the corresponding index set. Then K 2 = # , andc Sincec n ( f ) =c n (φ 1 ) +c n (φ 2 ) for n ∈ N, it follows that only K 2 modified Fourier coefficients of f are not of the form as in (2.15) while allc n ( f ) with n / ∈ satisfỹ c n ( f ) =c n (φ 1 ) and can be reconstructed by a rational function of type (K 1 − 1, K 1 ). Let us now examine, how to reconstruct f (t) in this setting.

Theorem 5.1 Assume that f (t)
is of the form f (t) = φ 1 (t) + φ 2 (t) with φ 1 (t) and φ 2 (t) in (5.1) and (5.3), where φ 2 possesses the K 2 nonzero Fourier coefficients c n (φ 2 ), n ∈ ⊂ N. Assume that we have given a set c n ( f ), n ∈ ⊂ N of L ≥ 2K + 2 Fourier coefficients of f , where the unknown set is contained in . Then φ 1 and φ 2 , i.e., K 1 , K 2 as well as all parameters γ j , a j , b j determining φ 1 and φ 2 can be completely recovered from this set of Fourier coefficients. In particular, Algorithm 3.2 terminates after at most K + 1 steps and provides a rational function r K 1 (z) of type Proof Let as beforec n ( f ) := Re c n ( f ) + i n Im c n ( f ). Assume that we have taken a set S K +2 ⊂ of K + 2 given indices as interpolation points at the (K + 1)-th iteration step of Algorithm 3.2. We will show that the matrix A K +2 = c n −c k n 2 −k 2 n∈ K +2 ,k∈S K +2 ∈ C L−K −2×K +2 has rank K and possesses a kernel vector w ∈ C K +2 that satisfies the side condition w Tc S K +2 = 0. This in turn will imply that the weight vector w found in Algorithm 3.2 determines the rational function r K 1 (z) which interpolatesc n (φ 1 ) for all n ∈ N.

Let
∪ = with ∩ = ∅ be the partition of such that ⊂ S K +2 and ⊂ K +2 , and let K 2 := # and K 2 := # denote the numbers of elements of and , such that K 1 + K 2 + K 2 = K . Then the K + 2 − K 2 indices in S K +2 \ correspond to modified Fourier coefficients with the rational structurẽ q K 1 (n 2 ) = r K 1 (n 2 ) as in (5.2), and the same is true for the Assume that the rows and columns of the matrix A K +2 are ordered such that the first K +2 − K 2 columns of A K +2 correspond to S K +2 \ , while the last K 2 columns correspond to the index set . Similarly, we suppose that the rows of A K +2 are ordered such that the first L − K − 2 − K 2 rows correspond to the indices K +2 \ , while the remaining K 2 rows correspond to indices in . In other words, we obtain Since A 11 is only composed of the modified Fourier coefficients of the non-periodic function φ 1 , it follows similarly as in the proof of Theorem 4.1 that A 11 possesses rank K 1 . More exactly, we have with (5.2) the matrix factorization where the two Cauchy matrices and the diagonal matrix have full rank K 1 . Therefore, A 11 possesses a kernel of dimension K + 2 − K 1 − K 2 = K 2 + 2. Since A 21 contains only K 2 ≤ K 2 rows, it follows that A 11 A 21 has at most rank K 1 + K 2 and therefore possesses a kernel of dimension at least 2. Thus, A K +2 has at most rank 2. We will prove that A K +2 has exactly rank K by showing that rank (A 11 , A 12 ) = K 1 + K 2 and similarly, that rank A 11 A 21 = K 1 + K 2 . As in the proof of Theorem 4.1, we can always find a linear combination of K 1 columns of A 11 to represent the columns for all k ∈ . Indeed for each of these columns we have can be written as linear combination of the columns in andc k (φ 2 ) = 0 for k ∈ . Therefore, it suffices to show that the concatenation of the first matrix factor of A 11 and 1 n 2 −k 2 n∈ K +2 \ ,k ∈ , i.e., 1 n 2 − C j n∈ K +2 \ , j=1,...,K 1 , 1 n 2 − k 2 n∈ K +2 \ ,k ∈ has full rank K 1 + K 2 . This is obviously true since this Cauchy matrix has L − K − 2 − K 2 ≥ K − K 2 = K 1 + K 2 rows and the values C j , j = 1, . . . , K 1 and k 2 , k ∈ are pairwise distinct and also distinct from n 2 with n ∈ K +2 \ . Similarly we can show that rank A T 11 , A T 21 can be simplified to and has rank K 1 + K 2 . 4. Thus rank A K +2 = K , i.e., the dimension of the kernel of A K +2 is 2. Therefore, Algorithm 3.2 always finds a vector w in the kernel of A K +2 which satisfies also the side condition w Tc S K +2 = 0. Moreover, it follows from the previous observations that any vector w in the kernel of A K +2 , is of the form w = (w T , 0 T ) ∈ C K +2 , wherẽ w ∈ C K +2−K 2 is in the kernel of A 11 A 21 and particularly in the kernel of A 11 . Thus, it follows from Theorem 4.1 thatw has at least K 1 + 1 nonzero components and provides the rational function r K 1 (z) that interpolates all modified Fourier coefficients of φ 1 . However, differently from the proof of Theorem 4.1,w may possess more than K 1 + 1 nonzero components, and the computation of r K 1 (z) may involve the removal of Froissart doublets. 5. Having determined r K 1 (z) to interpolate all modified Fourier coefficients of φ 1 , we can find φ 2 of the form (5.3) by capturing all modified Fourier coefficientsc n ( f ) withc n ( f ) =c n (φ 1 ), n ∈ . For all n ∈ , we computec n (φ 2 ) =c n ( f ) − r K 1 (n 2 ). Then φ 2 can be reconstructed from (5.3) and (5.4), where is found as the set of indices n ∈ withc n (φ 2 ) = 0, and K 2 = # .

Remark 5.2 We can also reconstructf
where γ 0 = 0 is a constant. Then γ 0 can be seen as a periodic component of the function, and we have c n (f ) = c n ( f ) for all n = 0. Thus, if beside the set of Fourier coefficients c n ( f ), n ∈ , in Theorem 5.1 also c 0 (f ) is known, then the constant part γ 0 can be reconstructed, too.

Generalization of the model
The model (1.3) for signals considered in the previous sections can be generalized.
In particular, the Fourier coefficients c n (φ) do not vanish for all n ∈ N.
Thus, the Fourier coefficients of the signal in model (6.1) have still the same structure as found for the model (1.3) A j = γ j |a j |P π sinh(π |a j |P) cosh((−i(πa j P + b j |)) = − γ j a j P π sin(πa j P) cos(πa j P + b j ), Hence, the obtained parameters A j , B j , C j have exactly the same form as in (2.9)-(2.11). Moreover, we can reconstruct the parameters γ j , a j , b j of f in (6.1) from A j , B j , C j via a generalization of Theorem 2.5.

Corollary 6.3 Let f (t) be given as in
Then, there is a bijection between the parameters γ j , a j , b j , j = 1, . . . , K , determining f (t) and the parameters A j , B j , C j in (2.9)-(2.11), for j = 1, . . . , K . For C j > 0 we obtain a j , γ j , b j via Theorem 2.5. For C j < 0, we find Proof For C j < 0 it follows that a j = i √ |C j | P . Further, we obtain The parameter γ j > 0 is thus uniquely defined, since we always have A 2 j +C j B 2 j > 0. Finally, inserting the found parameters γ j and a j into (2.10) and (2.11), we obtain where arccosh is the inverse of cosh and maps onto [0, ∞). Note that for C j < 0, we necessarily have A j > 0.
Therefore, Algorithm 3.2 can also be applied to a set of Fourier coefficients of the generalized model (6.1) to obtain a rational function that approximates the Fourier coefficients of the non-periodic part of f . Then, we apply Algorithm 3.3 as before to find the partial fraction decomposition of the rational function as described in Sect. 3.4, and can reconstruct the wanted parameters for the nonperiodic part of f in (6.1) using Corollary 6.3. Finally, if f contains a periodic part φ 2 as studied in Sect. 5, i.e., if there are parameters a j ∈ 1 P N, then φ 2 can be reconstructed via Theorem 5.1.

Remark 6.4
The model (6.1) for real non-periodic functions f is the most general model, such that Fourier coefficients of f can be written as in (2.12). In particular, complex values for C j cannot occur in (2.12) for real functions, since we always have c −n ( f ) = c n ( f ) for n ∈ N.

Numerical experiments
In this section we present two numerical experiments, which show that the considered reconstruction scheme provides very good reconstruction results even for small frequency gaps, if P is chosen suitably. We will compare our results for parameter reconstruction with the ESPRIT method/matrix pencil method (a stabilized variant of Prony's method) as it is for example summarized in [25,Section 3], see also [10,27]. We remark that a direct comparison of our method with Prony-like methods is difficult, since these methods require different input values. Our model (1.3) can be rewritten as an exponential sum of the form where For the matrix pencil method we will employ equidistant function samples f (hk), k = 1, . . . , L, L ≥ 4K . To choose a suitable step size h, we need to assume that a good estimate of the upper bound of occurring frequencies is known beforehand, since we need to take h −1 > max{2 a j : j = 1, . . . , K } to avoid aliasing. In this way we ensure that 2πλ j h is in [−π, π) such that λ j can be uniquely determined from e 2π iλ j h .
Since all parameters γ j , b j , a j , j = 1, . . . , K , can already be reconstructed from the shorter exponential sumf we will also apply the matrix pencil method to the equidistant function samplesf (hk), k = 1, . . . , L, L ≥ 2K . Further, for the application of Prony's method, we will take 2K (or K ) as a known parameter in the matrix pencil algorithm, while our new method will detect the parameter K itself.
In a second example we consider the function (7.4) i.e., f (t) is given by the parameter vectors One Froissart doublet occurs in r (z). This is due to the fact that the Fourier coefficient c 4 corresponding to the periodic part of f has been chosen for interpolation by Algorithm 3.2 only in the last iteration step. According to the proof of Theorem 5.1, we therefore need 7 iteration steps to generate a kernel of A 8 of dimension 2. Application of Algorithm 3.3 then leads to 6 finite eigenvalues C 1 , . . . , C 6 of (3.11), while the equation system at the second step of Algorithm 3.3 yields a vector (A j + iB j ) 6 j=1 with one vanishing component. This component and the corresponding component C j (which is in N up to numerical error) are removed to obtain the rational function of type (4, 5) determining the non-periodic part φ 1 of f . We reconstruct the parameter vectorsã,b andγ according to Theorem 2.5 and Theorem 5.1 with (5.4) and errors a −ã ∞ = 9.6 · 10 −13 , b −b ∞ = 2.7 · 10 −12 , γ −γ ∞ = 3.4 · 10 −12 .  For comparison, the matrix pencil method applied to function values f (kh) with h = 0.05 (since we need h −1 > 2 √ 89 to avoid aliasing) and k = 1, . . . , L leads here also to very good recovery results since in this case, we have no small frequency gaps, and Prony's method does not need to care for frequency parameters being close to 1 P Z. We obtain for f in (7.4) and the corresponding functionf the results in Tables 3  and 4. One needs to keep in mind here that the almost optimal step size h = 0.05 has been taken here, which needs a priori knowledge about the size of the maximal frequency. Already h = 0.04 and L = 80 gives for the recovery of f only the errors a −ã ∞ = 2.90 · 10 −11 instead of 1.57 · 10 −14 , b −b ∞ = 3.12 · 10 −10 instead of 1.97 · 10 −13 , and γ −γ ∞ = 4.85 · 10 −10 instead of 5.25 · 10 −13 .

Conclusions
In this paper we have proposed a new method to recover functions f of the form (1.3) from their Fourier coefficients. This method exploits the special structure of the Fourier coefficients of f on some interval [0, P). More precisely, if all frequencies a k satisfy a k / ∈ 1 P N, then c k ( f ) can be represented as function values of a special rational function of type (K − 1, K ). In turn, we can apply an algorithm for rational approximation to recover all Fourier coefficients of f and to reconstruct all parameters a j , b j and γ j determining f . We need however to pay attention if frequency parameters of the form k P with k ∈ N occur which lead to P-periodic terms in f . The method in our paper can be generalized to the recovery of extended exponential sums of the form with γ j,m , λ j ∈ C and γ j,n j = 0, see [9].
Our method is substantially different from Prony-like methods for the recovery of exponential sums of the form (1.3). In particular, the usual Prony method employs equidistant function values while our method requires the knowledge of Fourier coefficients. The performance of Prony-based algorithms like ESPRIT and matrix pencil method strongly depends on the chosen step size for the input function values. This step size needs to be taken very small if f contains large frequencies. In such cases, only more sophisticated numerical methods that exploit aliasing effects can provide satisfying recovery results, [7]. For our method, no direct knowledge is needed about the highest occurring frequency. Besides the different approach to the recovery of exponential sums that we have taken here, there might be close connections to Prony's method that still need to be explored.
The numerical stability of our algorithm stems from the numerical stability of the AAA algorithm which is the main ingredient of our method. However, we need to study the numerical stability of this algorithm (in comparison to Prony-like methods) more exactly in the future. In particular, there are two main issues. First, what happens if the gap between neighboring frequencies gets small? Our numerical experiments in [9] show, that the AAA-algorithm still works well while the Prony-like methods break down in this case. Second, how does the algorithm behave for frequencies close to the set 1 P N?
The recovery of parameters of exponential sums from function values (or Fourier coefficients) is an ill-posed problem. Even for small changes in the input data the recovery methods may lead to parameter vectors that are far away from the parameter vectors determining the original signal. This is however not a failure of the recovery method, since the new parameter vectors often lead to exponential sums that approximate the given noisy values better than the exponential sum with the original parameters, see e.g. [34] for examples.
Our future research will focus on the problem of approximation of signals by (short) exponential sums. Signals as in (1.3) can be seen as truly adaptive counterparts of Fourier sums. A general answer to the question, how well smooth functions can be approximated by exponential sums, is open. We are still in need of reliable algorithms to recover digital signals using a small number of structural parameters.