CINXE.COM

<html> <head> <meta http-equiv="content-type" content="text/html;charset=utf-8"> <link rel="stylesheet" type="text/css" href="../res.css"> <meta name="description" content="Notes on Prony's and related methods for sinusoid and exponential signal estimation."> <meta name="keywords" content="Prony's method, exponential fitting, frequency estimation, differential equations, iterated eigenvalue methods, elemental sets, robust estimation, nonlinear least squares, high breakdown, high efficiency."> <title>Prony Estimation</title> </head> <body>  <table border="0" cellpadding="0" cellspacing="0"> <tr> <td><img border="0" src="../images/statsci24.png" width="148" height="27"></td> <td>  / <a href="../index.html">Home</a></td> </tr> </table>  <h1>Prony Estimation</h1> Estimation of sinusoid and exponential signals using eigen-analysis of covariance matrices. <hr> <h2>Table of Contents</h2> <ol> <li>Introduction</li> <li>Least Squares</li> <li>Prony's Method</li> <li>Modified Prony Methods</li> <li>Real Exponential Fitting</li> <li>Frequency Estimation</li> <li>Articles</li> <li>Software</li> <li>Links</li> </ol> <h2>Introduction</h2> Many functions, used to model physical systems in engineering or physics, arise as the solution to a homogeneous differential equation. Let μ(t) be the concentration or intensity of the physical process at time t. It is often the case that there are physical reasons for supposing that μ(t) satisfies a differential equation like <table border="0" cellspacing="0" cellpadding="0"> <tr> <td rowspan="3">(1)         ∑</td> <td>p+1</td> <td rowspan="3">ξk Dk-1μ(t) = 0           for all t</td> </tr> <tr> <td></td> <td></td> <td></td> </tr> <tr> <td>k=1</td> <td></td> <td></td> </tr> </table> where D is the differential operator and the ξk are unknown coefficients. A typical solution to (1) is a sum of exponential functions <table border="0" cellspacing="0" cellpadding="0"> <tr> <td rowspan="3">(2)           μ(t) = ∑</td> <td>p</td> <td rowspan="3">αj exp(βjt)</td> </tr> <tr> <td></td> <td></td> <td></td> </tr> <tr> <td>j=1</td> <td valign="top"></td> <td></td> </tr> </table> where the βj are rate constants (usually negative) and the αj are amplitudes. Such a solution is appropriate for a transient system which dies to zero as time goes on, such as the concentration of an intermediary compound in a chemical reaction.  Another typical solution to (1) is a sum of sinusoids <table border="0" cellspacing="0" cellpadding="0"> <tr> <td rowspan="3">(3)           μ(t) = ∑</td> <td>p/2</td> <td rowspan="3">αj sin(βjt + φj)</td> </tr> <tr> <td></td> <td></td> <td></td> </tr> <tr> <td>j=1</td> <td valign="top"></td> <td></td> </tr> </table> Such a solution is persistent and periodic, and is appropriate for electrical or physical systems driven by persistent regular oscillating forces. Signals in speech recognition are often of this type. A third common solution is the damped sinusoid, <table border="0" cellspacing="0" cellpadding="0"> <tr> <td rowspan="3">(4)           μ(t) = ∑</td> <td>q</td> <td rowspan="3">αj exp(βjt) sin(ωjt)</td> </tr> <tr> <td></td> <td></td> <td></td> </tr> <tr> <td>j=1</td> <td valign="top"></td> <td></td> </tr> </table> which combines transient and periodic behaviour.  An example would be a pendulum which swings which an exponentially decreasing amplitude. <h2>Least Squares</h2> In practice we don't observe μ(t) exactly, because there is always experimental error. Instead we observe yi = μ(ti) + εi, where the εi are random observation errors. If the errors are approximately Gaussian, then we will try to estimate or identify the unknown parameters in the exponential or sinusoid signal by minimizing the sum of squared errors,            ∑ [yi - μ(ti)]2 This is a nonlinear least squares problem in the unknown parameters αj and βj. The basic idea of Prony estimation is to represent the signal μ(t) not in terms of the αj and βj but in terms of the coefficients of the differential equation (1). The coefficients are identified by computing the eigenvectors of a suitably calculated covariance matrix. <h2>Prony's Method</h2> Prony's method is a technique for extracting the sinusoid or exponential signals by solving a set of linear equations for the coefficients of the recurrence equation that the signals satisfy. It is closely related to Pisarenko's method, which finds the smallest eigenvalue of an estimated covariance matrix. The classical method of Count de Prony models a sequence of 2p observations made at equally spaced times by a linear combination of p exponential functions. Prony's ingenious method converted the problem to a system of linear equations. See the MacTutor biography of <a href="http://www-groups.dcs.st-and.ac.uk/history/Mathematicians/De_Prony.html">Baron Gaspard Riche de Prony</a>. The original reference is <blockquote> de Prony, Baron Gaspard Riche (1795). Essai éxperimental et analytique: sur les lois de la dilatabilité de fluides élastique et sur celles de la force expansive de la vapeur de l'alkool, à différentes températures. Journal de l'école Polytechnique, volume 1, cahier 22, 24-76. </blockquote> A contemporary treatment of modern Prony methods can be found in <blockquote> Marple, S. L. (1987). Digital Spectral Analysis with Applications. <a href="http://www.prenhall.com">Prentice-Hall</a>. </blockquote> <h2>Modified Prony Methods</h2> Unfortunately, Prony's method is well known to perform poorly when the signal is imbedded in noise; Kahn et al (1992) show that it is actually inconsistent. The Pisarenko form of the method is consistent but inefficient for estimating sinusoid signals and inconsistent for estimating damped sinusoids or exponential signals. A modified Prony algorithm that is equivalent to maximum likelihood estimation for Gaussian noise was originated by Osborne (1975). It was generalized in Smyth (1985) and Osborne and Smyth (1991) to estimate any function which satisfies a difference equation with coefficients linear and homogeneous in the parameters. Osborne and Smyth (1991) considered in detail the special case of rational function fitting, and proved that the algorithm is asymptotically stable in that case. The modified Prony algorithm for exponential fitting will estimate, for fixed p, any function μ(t) that solves a constant coefficient differential equation (1). Perturbed observations, yi  =   μ(ti) + εi , are made at equi-spaced times ti, i  =  1, ..., n, where the εi are independent with mean zero and variance σ2. The solutions to (1) include complex exponentials, damped and undamped sinusoids and real exponentials, depending on the roots of the polynomial with the ξk as coefficients. The modified Prony algorithm has the great practical advantage that it will estimate any of these functions according to which best fits the available observations. The modified Prony algorithm uses the fact that the μ(ti) satisfy an exact difference equation when the ti are equally spaced. The algorithm directly estimates the coefficients, γk say, of this difference equation. The residual sum of squares after estimating the αj can be written in terms of the γk. The derivative with respect to γ = (γ1, ..., γp+1)T can then be written as 2B(γ)γ, where B is a symmetric matrix function of γ. The modified Prony algorithm finds the eigenvector of B(γ)γ = λγ corresponding to λ = 0 by the fixed point iteration in which γk+1 is the eigenvector of B(γk) with eigenvalue nearest zero. The eigenvalue λ is the Lagrange multiplier for the scale of γ in the homogeneous difference equation. Inverse iteration proves very suitable for the actual computations. Although the algorithm estimates all functions in the same way, the practical considerations and asymptotic arguments differ depending on whether the signals are periodic or transient, real or complex. The algorithm has been applied to real exponential signals by Osborne and Smyth (1991, 1995), to real sinusoidal signals by Mackisack et al (1991) and Kahn et al (1992) and to exponentials with imaginary exponents in complex noise by Kundu (1993). <h2>Real Exponential Fitting</h2> When the polynomial with coefficients ξk has real distinct roots, the solution to (1) is a sum of real exponential functions (2). If rate coefficients are negative, then this represents a function which is transient, i.e., decays to zero as time goes on. This represents the outcome of a chemical experiment or other process in which the reaction completes in a finite amount of time. Real exponential fitting is one of the most important, difficult and frequently occurring problems of applied data analysis. Applications include radioactive decay, compartment models and atmospheric transfer functions. Estimation of the αj and βj is well known to be numerically difficult. General purpose algorithms often have great difficulty in converging to a minimum of the sums of squares. This can be caused by difficulty in choosing initial values, ill-conditioning when two or more βj are close, and other less important difficulties associated with the fact that the ordering of the βj is arbitrary. The modified Prony algorithm solves the problem of ordering the βj and is relatively insensitive to starting values. It also solves the ill-conditioning problem as far as convergence of the algorithm is concerned, but may return a pair of damped sinusoids in place of two exponentials which are coalescing. In some applications the restriction to positive coefficients αj is natural. A convex cone characterization is then possible, and special algorithms have been proposed. We prefer to treat the general problem with freely varying coefficients since this is appropriate for most compartment models. A common attempt to reduce the difficulty of the general problem has been to treat it as a separable regression, i.e., to estimate the coefficients by linear least squares conditional on the rate constants βj. Another approach has been suggested by Ross (Nonlinear Regression, 1990, Section 3.1.4), who suggests that the coefficients of the differential equation (1) comprise a more "stable" parametrization of the problem than do the parameters of (2). Both of these strategies are part of the modified Prony algorithm. <h2>Frequency Estimation</h2> When the polynomial with coefficients ξk has distinct complex roots on the unit circle, the solution to (1) is a sum of sinusoidal signals (3) with p/2 distinct frequencies. This represents a function which is periodic and which can be observed on an ongoing basis. The frequency estimation problem is mathematically similar to the real exponential problem, but has its own difficulties and features. The least squares frequency estimators are known to be consistent at rate O(n-3/2) while the coefficient estimators are consistent at rate O(n-1/2). The sum of squares objective function is not convex as a function of the frequencies. Rather it has minima which are approximately 2π/n apart. General purpose least squares algorithms have no trouble converging when used to fit sinusoidal models to data, but they typically converge to a local minima rather than to the global minimum of the sum of squares. To achieve efficient estimation of the frequencies it is necessary to resolve the frequencies to the correct global minimum. One algorithm which achieves this is the ORA Prony algorithm from Kahn et al (1992). The Prony methods can be constrained to give periodic solutions when the solution is know apriori to be a sum of sinusoids. See Kannan and Kundu (1994), Kundu and Kannan (1997) and Smyth (1998). Another method which can be used to handle the non-convexity of the sum of squares is to use elemental sets. An elemental set is the minimum number of observations which are sufficient to identify the parameters in the model, 3p/2 for the frequency estimation problem. The elemental set method consists of using Prony's method to find identify the frequencies and coefficients for a large number of randomly selected elemental sets. This methodology not only allows the global minimum of the sum of squares to be resolved, but also provides a way to handle outliers in the data. Smyth and Hawkins (1997, 1999) show how to obtain frequency estimators with both high efficiency and high breakdown properties. <h2>Articles</h2> <ul> <li>Osborne, M. R. (1975). Some special nonlinear least squares problems. SIAM Journal of Numerical Analysis, 12, 571-592.</li> <li>Kay, S. M., and Marple, S. L. (1981). Spectrum analysis - a modern perspective. Proc. IEEE, 69, 1380-1419. </li> <li>Smyth, G. K. (1985). Coupled and Nested Iterations in Nonlinear Estimation. PhD Thesis, Australian National University, Canberra. </li> <li>Bresler, Y., and Macovski, A. (1986). Exact maximum likelihood parameter estimation of superposed exponential signals in noise. IEEE Trans. Acoust., Speech, Signal Processing, 34, 1081-1089. </li> <li>Osborne, M.R. and Smyth, G.K. (1986). An algorithm for exponential fitting revisited. Contribution to Essays in Time Series and Allied Processes: Papers in honour of E. J. Hannan. J. Gani and M. B. Priestley (Editors). Applied Probability Trust, Sheffield, 419-430.</li> <li>Osborne, M.R. and Smyth, G.K. (1991). A modified Prony algorithm for fitting functions defined by difference equations. SIAM Journal of Scientific and Statistical Computing, 12, 362-382. (<a href="../smyth/pubs/pronyr.pdf">PDF</a>)</li> <li>Kahn, M., Mackisack, M. S., Osborne, M. R., and Smyth, G. K. (1992). On the consistency of Prony's method and related algorithms. Journal of Computational and Graphical Statistics, 1, 329-349.</li> <li>Kundu, D. (1993). Estimating the parameters of undamped exponential signals. Technometrics, 35, 215-218.</li> <li>Mackisack, M.S., Osborne, M.R., and Smyth, G.K. (1994). A modified Prony algorithm for estimating sinusoidal frequencies. Journal of Statistical Computation and Simulation, 49, 111-124.</li> <li>Kundu, Debasis. (1994). Estimating the parameters of complex-valued exponential signals. Comput. Statist. Data Anal. 18, no. 5, 525-534.</li> <li>Kannan, N., and Kundu, D. (1994). On modified EVLP and ML methods for estimating superimposed exponential signals. Signal Processing, 39, 223-233.</li> <li>Osborne, M. R., and Smyth, G. K. (1995). A modified Prony algorithm for fitting sums of exponential functions. SIAM Journal of Scientific Computing, 16, 119-138. (<a href="../smyth/pubs/pronyexp.pdf">PDF</a>) </li> <li>Smyth, G. K., and Hawkins, D. M. (1997). Robust frequency estimation using elemental sets. Computing Science and Statistics, 28, 659-662. </li> <li>Kundu, Debasis, and Kannan, Nandini. (1997). Constrained maximum likelihood estimators for superimposed exponential signals. Comm. Statist. Simulation Comput. 26, no. 2, 733-764.</li> <li>Mitra, Amit, and Kundu, Debasis. (1997). Consistent method for estimating sinusoidal frequencies: a non-iterative approach. J. Statist. Comput. Simulation 58, no. 2, 171-190.</li> <li>Smyth, G. K., and Hawkins, D. M. (2000). Robust frequency estimation using elemental sets. Journal of Computational and Graphics Statistics 9, 196-214. (<a href="../smyth/pubs/element.html">Abstract</a> - <a href="../s/robfreq.html">Software</a> - <a href="../smyth/pubs/element.pdf">PDF</a>) </li> <li>Smyth, G. K. (2000). Employing constraints for improved frequency estimation by eigenanalysis methods. Technometrics 42, 277-289. (<a href="../smyth/pubs/constrai.html">Abstract</a> - <a href="../smyth/pubs/constrai.pdf">PDF</a>)</li> </ul> <h2>Software</h2> <ul> <li><a href="../matlab/statbox.html">StatBox</a>. A statistics toolbox for Matlab. Includes a modified Prony (least squares) algorithm for estimating solutions of homogeneous differential equations. Gordon Smyth, Walter and Eliza Hall Institute of Medical Research.</li> <li><a href="../s/pronyfre.html">pronyfreq</a>. An S-Plus function for frequency estimation which implements the ORA algorithm of Kahn et al (1992).</li> <li><a href="../s/robfreq.html">robfreq</a>. An S-Plus function for frequency estimation which implements the multistage algorithm of Smyth and Hawkins (1999).</li> </ul>   <hr> <div align="center"><center> <table border="0" cellpadding="0" cellspacing="0"> <tr> <td height="30"></td> </tr> <tr> <td align="center"><form align="center" action="http://www.statsci.org/cgibin/searchwg.pl" method="get" name="SearchForm"> <input type="text" size="32" name="Terms"> <input type="submit" value="Search"> <a href="../search.html">Help</a> </form> </td> </tr> <tr> <td align="center" valign="bottom" height="50"> <a href="../index.html">Home</a> - <a href="../about.html">About Us</a> - <a href="../contact.html">Contact Us</a> Copyright © <a href="../smyth/index.html">Gordon Smyth</a></td> </tr> </table> </center></div>  </body> </html>