Parameter estimation and signal reconstruction

In frequency analysis an often appearing problem is the reconstruction of a signal from given samples. Since the samples are usually noised, pure interpolating approaches are not recommended and appropriate approximation methods are more suitable as they can be interpreted as a kind of denoising. Two approaches are widely used. One uses the reflection coefficients of a finite sequence of Szegő polynomials and the other one the zeros of the so called Prony polynomial. We show that both approaches are closely related. As a kind of inverse problem, it’s not surprising that they have in common that both methods depend very sensitive on sampling errors. We use known properties of the signal to estimate the positions of the zeros of the corresponding Szegő- or Prony-like polynomials and construct adaptive algorithms to calculate these ones. Hereby, we get the corresponding parameters in the exponential parts of the signal, too. Then, the coefficients of the signal (as a linear combination of such exponential functions) can be obtained from a system of linear equations by minimizing the residuals with respect to a suitable norm as a kind of denoising.

A problem in signal theory is to get the signal h or an approximation of it from given samples. From the sampling theorem it is known that 2m samples are necessary for the reconstruction, i.e. to get the j , j , j = 1, … , m.
The classical approach of Prony [35] uses 2m samples h k to calculate a monic polynomial m in its monomial representation whose zeros are the desired z j , i.e. m (z) = ∏ m j=1 (z − z j ) . By construction, m has only simple zeros. If we know the z j we get the corresponding j via j = Log z j (complex logarithm).
If all j are known then the coefficients j , j = 1, … , m , can be obtained from (1) as the solution of a system of linear equations, using the samples as interpolation conditions, see Sect. 9.
Usually, the measured values h k contain some errors which may also distort the Prony polynomial. From the numerical analysis it is well known that even small perturbations of the monomial coefficients may cause essential changes for the zeros of m . This can already be the case for small values of m. Thus, Prony's original method is not recommended for general signal reconstruction.
Nevertheless, methods which are based on this approach are widely used by engineers, in particular for special applications, e.g. in the medical, radiophysical and electromechanical context, see [9-11, 15, 22, 23, 28] for some few examples.
Another classical way constructs the Prony polynomial not in its monomial representation, see e.g. [1,2,13,[17][18][19][20][21]. Using given samples, some authors construct moments of a weight function. By a Levinson-like algorithm, an appropriate finite sequence of Szegő polynomials (s ) m =0 resp. its reflection coefficients can be obtained. Here, s m coincides up to normalization with the Prony polynomial m if all zeros of m are in the interior or all of them are on the boundary of the unit disc.
The reflection coefficients are the entries of a Hessenberg matrix whose eigenvalues are the zeros of s m . But this method has some drawbacks, especially if some zeros are in the interior of the unit disc and some of them are unimodular ones. Then, the convergence of the zeros of s m to the unimodular ones of m is not guaranteed [29,30].
Besides the least squares approaches, one can use the l 1 -approximation. For example, in [36, p. 255] it was noted, that this approximation appears to be markedly superior among other l p -approximations, p > 1 , in presence of outliers or wild points. In the next section we sketch an algorithm, based on the method given in [42], to calculate an approximate Prony polynomial by using the l 1 -norm for the overdetermined case (i.e., if we have at least 2m samples).
In this article we focus on the overdetermined case with N ≥ 2m samples: Theoretically the signal can be reconstructed in the absence of noised samples and numerical uncertainties. But in practice, due to noising and occuring numerical errors it can only be approximated.
Nevertheless, if we have too few samples (i.e. an underdetermined system), we notice that an l 1 -solution is in most cases the sparsest solution [7], i.e. the obtained approximation of the Prony polynomial has in most cases the minimal number of non zero coefficients in its monomial representation, see e.g. [4,5,24]. Thus, an l 1 -solution for the underdetermined analogues of (4) resp. of (21) yields usually sparse vectors and hence a mimimal value of m and therefore a sparse approximation for h.
Usually, m is not known a priori and will be estimated. If m is initially overestimated by ≥ m , by the l 1 -norm approach we get a linear programming problem. Depending on the kind of the underlying noise, the number of basis variables in the final tableau gives for the real case an 'optimal' value of m which is the numerical rank of the matrix of the overdetermined system of linear equations [42]. (For the complex case see there, too.) For 'moderate' noise this coincides with the true value of m. If we have too few samples, the underdetermined system yields an approximation of the calculated signal in the above mentioned sparse sense.
We see that there exist various methods to calculate a Prony polynomial m or an approximation of it. A possible approach is sketched in the next section.
However, the main objective of this article is not the calculation of an (approximate) Prony polynomial but the determination of its zeros and thus the parameters ( j ) m j=1 of the signal. In principle, each zero finder for polynomials can be used. But using properties of this polynomial, which can be derived from known characteristics of the signal, maybe a more appropriate approach.
Here, we use knowledges on possible positions of the zeros of m to get more insights in the structure of this polynomial (Sect. 3). For example, we use connections between the Szegő-and the Prony-polynomials in Sects. 4 and 5 to find a factorization of m as the product of two polynomials m = s m− , whose zeros are in the interior resp. on the boundary of the unit disc.
More precisely, all the zeros of s are simple and lie in [[e − , 1)) . In Sect. 5 we construct a sequence of Szegő polynomials (s ) =0 which we use in Sect. 6 to determine their reflection coefficients. With them we built up a Hessenberg matrix whose eigenvalues are the zeros of s .
In Sect. 7 we derive an algorithm to calculate the zeros of m− . We use the property that all of its zeros are on the unit circle for the design of this algorithm.
From the calculated zeros z j of m we obtain j as the complex logarithms j = Log z j . After calculating the corresponding coefficients j in Sect. 9, the signal h, given as in (1), is reconstructed respectively approximated from the given samples.
Numerically, it is difficult to distinguish between zeros on the unit circle and 'near' to it, especially for a noised signal. Thus, this case deserves a closer attention in Sect. 8.
In total, we get an algorithm to reconstruct or approximate a given signal from possible noised samples. Some tests and numerical remarks are given in Sect. 10.

Determining an approximate prony polynomial
Prony [35] has shown that the coefficients of the monomial representation of m , can be obtained from 2m samples h ∶= h( ), = 0, … , 2m − 1 , where h is the exponential sum in (1), by solving the system of linear equations Unfortunately, the zeros of m depend very sensitive on these coefficients: For a polynomial in monomial representation it is known that already slight perturbations of some coefficients may vary the position of its zeros substantially (a famous example is Wilkinson's polynomial), even for small m. Here, the coefficients depend on the samples, which are often distorted in the real world, e.g. by sampling errors, and the value of m is usually not known a priori. Thus, it is recommended to use known properties of the signal for its reconstruction/approximation. Algorithms for calculating zeros of a polynomial or the eigenvalues of its companion matrix with satisfactory stability behaviours are for example given in [3].
In [42] we modified Prony's method and used that p 0 ≠ 0 (instead of p m = 1 ). This is possible because we have assumed that all zeros lie in [[e − , 1]] , and thus 0 is not a zero of m . We may renormalise m so that p 0 becomes an arbitrary non zero value, e.g. p 0 ∶= 1 , which can be done by defining p ∶= With the residual vector we developed an algorithm which minimizes ‖r(p)‖ † in [42]. The norm ‖ ⋅ ‖ † is defined by ‖z‖ † ∶= ‖ Re z‖ 1 + ‖ Im z‖ 1 , z ∈ ℂ N− , which coincides with the l 1 -norm ‖ ⋅ ‖ 1 for real or pure complex vectors. Without additional work, we get a reasonable value of m, which takes possible noised values of the samples into account, even if H N− , has singular values with small moduli. We can interpret m (real case) resp. 2m (complex case) as the number of basis variables in the final tableau of the algorithm described in [42]. Of course, the numerical rank of H N− , depends on some tolerance parameters in this algorithm.
It can be shown that such a solution p † is a near best approximation for which the l 1 -error within a factor √ 2 of its minimal value ‖r(p ⋆ )‖ 1 . The solution (p 0 , … , p m ) T of (3) can be used as a solution of (4) by p = as an approximation to the 'true' Prony polynomial that particularly takes erroneous samples into account. Especially its zeros can be used as approximations for zeros of the Prony polynomial which can be obtained by suitable root finding algorithms.
In [42] the main focus is the construction of the Prony-like polynomial m and the calculation of the coefficients j , j = 1, … , m . Here, we set the key aspect in determining the zeros of m , using special properties of the signal. In the subsequent sections we don't distinguish between the Prony(-like) polynomial and its numerical calculated approximation and focus on the goal to calculate the zeros of m which maybe not monic.
Although Szegő polynomials are widely used in frequency analysis in classical algorithms, see e.g. [17][18][19][20][21], it is also known from these papers that their construction is sensitive not only due to erroneous samples but as well due to the number of samples. Furthermore, the convergences of the algorithms are not always guaranteed [21,29,30]. These articles suggest that the moments used there depend sensitive on the samples, too.
We remark that most of the just before mentioned articles make the assumption that in their model all j have real parts zero, i. e. that |z j | = 1 for all j.
Apart from a degeneracy case, Szegő polynomials have all its zeros in the interior of the unit disc, so it seems that their usage would be inappropriate in that context.
We consider the more general case 0 < e − ≤ |z j | ≤ 1 for all j. Starting from m , in Sect. 4 we construct Szegő polynomials (s ) =0 . In Sect. 6 we use their reflection coefficients to find the zeros in the interior of the unit disc and a polynomial m− , whose zeros are all unimodular, m = s m− . Here, the number of zeros of m resp. s in [[e − , 1)) , will be calculated by our algorithm, too.
By numerical reasons we cannot distinguish exactly between zeros |z j | < 1 and |z j | = 1 . In fact, we distinguish between zeros z j of a Szegő polynomial s with �z j � < 1 − √ , 0 < ≪ 1 , and zeros of a polynomial m− whose zeros are in . The latter ones can be obtained by another approach which calculates zeros with a fixed modulus, where we vary the moduli by a bisection strategy. Both polynomials, s and m− , will be calculated in Sect. 5. 1]] . (Note, that we have normalized p 0 = 1 there.) But even for erroneous polynomials the calculated zeros in the interior of the unit disc are usually close to the exact ones in our tests. Thus, we don't add this restriction. Only zeros on [ [1,1]] or close to it are sometimes more affected in our tests, see example 1. These ones can be improved e.g. by iterative methods like the Newton method. Of course, this requires that erroneous samples don't affect the calculation of m respectively the positions of its zeros too much.
Besides the case |p m | = |p 0 | , i.e. that all zeros are unimodular, which we examine in Sect. 7, we have to consider the case |p m | − |p 0 | > 0 , i.e. that at least one zero has modulus < 1 . Then we also have By numerical reasons we choose an ∈ ℝ , 0 < ≪ 1 , and assume that �p m � > √ . Although the leading coefficient can be chosen arbitrarily ( ≠ 0 ) it shouldn't be chosen to have a too small modulus.
Furthermore, we distinguish in the subsequent sections between the cases �z j � (there exist some zeros in the interior of the unit disc; but it may be that there are unimodular zeros ), besides possible unimodular zeros, all (at least one) zeros ). In the next sections we generate a sequence of polynomials ( m− ) =0 , where m− has only unimodular zeros and is a common divisor of all m− , = 0, … , . Especially, its zeros are the unimodular ones of m . The zeros of m− will be calculated by the algorithm, derived in Sect. 7.
By construction, s j ∶= j+m− ∕ m− , j = 0, … , , have all their zeros in the interior of the unit disc. We will see that it is a sequence of Szegő polynomials, whose coefficients are the entries of a Hessenberg matrix. Its eigenvalues are the zeros of s resp. of m in [[e − , 1)) . Usually, the zeros of s have moduli less that 1 − √ for sufficient small . Nevertheless, there maybe zeros of m in [[1 − √ , 1)). By a modification of the algorithm for calculating the unimodular zeros, the number of such zeros can be determined and calculated by a bisection strategy. This will be done in Sect. 8.
It should be noted that

Construction of a recurrence for the prony polynomial
In this section we consider the case From this we conclude Since 2 < < √ , should be chosen, so that 2 is larger that the machine-epsilon and 1∕ 2 smaller than the largest finite floating point number.
For a polynomial m ∈ Π m (space of (complex) algebraic polynomials with highest degree m) we define its reciprocal polynomial ⋆ m by If m is given in its monomial representation (2), this can be written as Since m (0) = p 0 ≠ 0 ≠ p m we see that ⋆ m has degree m, too. We take the approach and want to find the coefficients p (1) of a polynomial m−1 , , such that (7) is satisfied. This approach is based on the construction of the recurrence of Szegő polynomials. We start with m , generate a sequence m−1 , m−2 , … , m− of polynomials for a suitable ≤ m (this will be specified later) to get a factorization m = m− s , where all zeros of s are in [[e − , 1)) and all zeros of m− are unimodular.
From (7) we immediately see Combining this identity with (7) we get On the other hand, using the monomial representations of m and ⋆ m , we obtain The right side can be simplified if we choose the normalization p (1) m−1 ∈ ℝ�{0} for the leading coefficient of (1) m−1 . This has no effects for the zeros of the polynomials. Usually, we may assume p (1) m−1 = 1 , i.e. that m−1 is monic. For a real valued leading coefficient of m−1 we get ⋆ m (z) = p 0 z m + ⋯ + p m .  (7) holds. Besides its normalization m−1 is uniquely determined.
We get a sequence of polynomials m− , since all zeros of m− have modulus 1 for ≤ m − 1.
In the inner loop of Algorithm 1 the product is obviously independent of and can be calculated recursively by We obtain Algorithm 2 as a more elaborated reformulation of Algorithm 1. It is not necessary that m is given in its monomial representation. (14) can be rewritten as from which we get in the same manner by (13) In summary, we obtain Applying these recurrences, starting with m and a m = m (0) , the polynomials can be generated 'backwards'. Again, if is unknown, these recurrences stop if the denominator 2 − |a | 2 is close to zero with = m − , i.e. with m− which has only unimodular zeros for < m. If m is given in monomial representation, the construction of the monomial coefficients of in Algorithm 2 corresponds with the Schur-Cohn algorithm/ Schur-Cohn test [14][Chapter 6].
As in (10)   Lemma 1 Each s j , defined in (18), has all its zeros in the interior of the unit disc and 0 is no zero of them. Furthermore, s j and s j−1 , j = 1, … , , have no common zeros.
Proof By (18) it can be seen that the leading coefficient of Since renormalising a polynomial doesn't change its zeros, we consider the monic polynomials  (17). ◻ Obviously, the zeros of m in [[e − , 1)) are these ones of s (which have moduli ≥ e − by the given assumptions), and its unimodular zeros are these ones of m− .
The elements (s j ) j=0 satisfy the recurrence (19), whose reflection coefficients have moduli less than one. By [12, Theorem 4.1 and the sentence before] it follows that these polynomials are orthogonal with respect to a sequence of complex numbers, respectively with respect to a linear functional. More precisely, the examinations before this cited Theorem 4.1 show that these polynomials are orthogonal on the (complex) unit circle [ [1,1]] with respect to a bounded nondecreasing measure with at least + 1 points of increase on Thus, the designation 'Szegő-polynomials' for the elements of both sequences (s j ) j=0 and (s j ) j=0 is justified, because all s j and s j differ only by a normalisation factor. Especially, both polynomials have the same zeros for j = 1, … , , which do not affect the orthogonality properties of both sequences of polynomials.
But this article does not focus on the orthogonality of these polynomials, rather we need the recurrence and the reflection coefficients.

Zeros in the interior of the unit disc
In this section we determine the zeros of s . We abbreviate j ∶= and (18)  Here, we have used Since is similar to , the eigenvalues of both matrices are the same and s is the characteristic polynomial of and hence of , too, c.f. [2,Theorem 3.2], [14]. From (10) we see, as expected, that the entries of the non degenerated Hessenberg matrix are independent of the representation of the s j , because there appear only terms

Zeros on the unit circle
In this section we consider (z) = ∑ =0 p ( ) z on [ [1,1]]. Although its zeros can be determined with standard zero finders, see e.g. [26,27], we sketch an algorithm which uses the knowledge that has only unimodular zeros.
We follow ideas given in [25] and write z = Thus, the common zeros t j of the trigonometric polynomials u and v in (− , ] yield the zeros z j of m in [ [1,1]] and j = i t j , j = 1, … , . By the given assumptions there are simple ones. In the following we consider the case that has only real coefficients p ( ) , = 0, … , (the case that it has only imaginary coefficients can be treated in the same manner). For the general case one may consider the real polynomial | | 2 . From its 2 zeros z j , z j only these ones with Re z j ≥ 0 needed to be calculated. The other ones can be obtained by conjugation. Then, verify which of them are the zeros of .
For p ( ) ∈ ℝ for all = 0, … , , the functions u , v can be written as Furthermore, we define the even trigonometric polynomial ṽ −1 by j k j k Since ([ [1,1]]) is symmetric to the real axis, it is sufficient to consider (e i t ) for t ∈ [0, ] instead for − < t ≤ . As the transformation x ∶= cos t, t ∈ [0, ], x ∈ [−1, 1] , is bijective we have where T (resp. U ) denotes the th Chebyshev polynomial of the first (resp. second) kind. Since a possible zero of in ±1 can be detected and divided out by the Horner algorithm, we assume in the following that (±1) ≠ 0 and especially that is even because is a real polynomial with only non real and simple zeros. In this case, a zero of on [ [1,1]] can be obtained from the common zeros of r and r −1 in (−1, 1) , i.e. as the zeros of the greatest common divisor ( gcd ) of r and r −1 .
As has zeros on |z| = 1 and by the symmetry of its zeros with respect to the real axis, there are 2 zeros for t ∈ (0, ) . Thus, r ∕2 ∶= gcd(r , r −1 ) ∈ Π ∕2 can be calculated by the Euclidean algorithm for Chebyshev expansions which we have introduced in [25]. This algorithm yields r ∕2 as an U -expansion. It has all its ∕2 (simple) zeros in (−1, 1) , which can be obtained successively e.g. by Newton's method with deflation. The deflation can be done by the Euclidean algorithm for Chebyshev expansions [25], too, or by the Clenshaw algorithm in the version which we have described in [41, p. 605 f.], where q (1) n−1 from there is the deflated polynomial and y from there denotes the zero. Alternatively, the algorithm in [6] calculates the roots as eigenvalues of a matrix, whose coefficients depend on the recurrence coefficients of the U and on the coefficients of r ∕2 in terms of the U .

Zeros near the unit circle
Until now we have distinguished between finding zeros z j with modulus less than one or equal one. By numerical reasons, we have chosen a small value , 0 < ≪ 1 , ) and especially that the zeros of s resp. the eigenvalues of are there, too. Consequently, Algorithm 2 calculates a sequence of polynomials whose all zeros are there.
Furthermore, for m− we have calculated its unimodular zeros. If we would have exact samples and no numerical inaccuracies, we could choose = 0 and would obtain all zeros of m on the unit circle (i.e. that all m − zeros of m− are there).    21) has an unique solution. The overdetermination of this system can be used to regularize the effects of numerical errors during the calculations. This can be done by minimizing ‖r( )‖ , where r( ) ∶= h − E N,m and ‖ ⋅ ‖ denotes a suitable norm, usually the l 1 , l 2 ones, or the norm ‖ ⋅ ‖ † from Sect. 2. In a last step, all j e j ⋅ of h with | j | < for a suitable small positive will be removed, and we get m ∶= #{j ∶ | j | ≥ } , so that a possible overestimated value of m can be corrected.

Numerical remarks
Example 1 To test the algorithm for various values of m and we have randomly generated j , j = 1, … , m, and thus the zeros of a polynomial m which we expand in its monomial representation, using Mathematica. Its zeros in the interior of the unit disc are computed by the algorithm in Sect. 6 and compared with the exact zeros by the l 2 -norm of the residual vectors. The results are given in Table 1, where the residual-norms are mean values of the considered examples.
Even in cases with larger norms of the residuals, most zeros are accurately determined and only some few 'outliers' appeared, see Fig. 1 for an example (randomly generated). Thus, it may be appropriate, especially for large m, to use an iterative (e.g. Newton's) method to improve the accuracy of the numerically calculated zeros.  with a Hessenberg matrix 20 6 ∈ ℂ 14,14 , whose superdiagonal elements have moduli > √ . The first element of has an absolute value greater than one. The element on superdiagonal position (1, 2) has a modulus close to zero, indicating that is too small.
Obviously, the element on position (1, 1) of 20 5 , i.e. 6 5 ∈ ℂ , is an approximate eigenvalue, but is outside of the unit disc and close to a zero 0.998768 − 0.0496163i of m on the unit circle, see Fig. 2. If would be correct, the eigenvalues of 20 5 would be the zeros of a Szegő polynomial s 15 which are, due to Lemma 1, all in the interior of the unit disc. This indicates that is chosen too large resp. too small.
If we choose = 10 −7 we obtain = 20 6 , whose eigenvalues are the zeros of s 14 . They are all in the interior of the unit disc. Since 20 5 is nearly degenerated, the eigenvalues of 20 6 are closely to the corresponding ones of 20 5 in the interior of the unit disc but without its incorrect one outside there. The l 2 -norm of the vector of the differences of the calculated zeros of s 14