/
Home |

Estimation of sinusoid and exponential signals using eigen-analysis of covariance matrices.

- Introduction
- Least Squares
- Prony's Method
- Modified Prony Methods
- Real Exponential Fitting
- Frequency Estimation
- Articles
- Software
- Links

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

(1) ∑ | p+1 | ξ_{k} D^{k-1}μ(t) = 0
for all t |

k=1 |

where D is the differential operator and the ξ_{k} are
unknown coefficients. A typical solution to (1) is a sum of exponential functions

(2) μ(t) = ∑ | p | α_{j} exp(β_{j}t) |

j=1 |

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

(3) μ(t) = ∑ | p/2 | α_{j} sin(β_{j}t
+ φ_{j}) |

j=1 |

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,

(4) μ(t) = ∑ | q | α_{j} exp(β_{j}t)
sin(ω_{j}t) |

j=1 |

which combines transient and periodic behaviour. An example would be a pendulum which swings which an exponentially decreasing amplitude.

In practice we don't observe μ(t) exactly, because there is always experimental error.
Instead we observe y_{i} = μ(t_{i}) + ε_{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,

∑
[y_{i} - μ(t_{i})]^{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.

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 Baron Gaspard Riche de Prony. The original reference is

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, volume1, cahier 22, 24-76.

A contemporary treatment of modern Prony methods can be found in

Marple, S. L. (1987).

Digital Spectral Analysis with Applications. Prentice-Hall.

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, y_{i} = μ(t_{i}) + ε_{i}
, are made at equi-spaced times t_{i}, 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 μ(t_{i}) satisfy an exact
difference equation when the t_{i} 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).

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.

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.

- Osborne, M. R. (1975). Some special nonlinear least squares problems.
*SIAM Journal of Numerical Analysis*,**12**, 571-592. - Kay, S. M., and Marple, S. L. (1981). Spectrum analysis - a modern perspective.
*Proc. IEEE*,**69**, 1380-1419. - Smyth, G. K. (1985).
*Coupled and Nested Iterations in Nonlinear Estimation*. PhD Thesis, Australian National University, Canberra. - 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. - 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. - 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. (PDF) - 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. - Kundu, D. (1993). Estimating the parameters of undamped exponential signals.
*Technometrics*,**35**, 215-218. - 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. - Kundu, Debasis. (1994). Estimating the parameters of complex-valued exponential signals.
*Comput. Statist. Data Anal.***18**, no. 5, 525-534. - Kannan, N., and Kundu, D. (1994). On modified EVLP and ML methods for estimating
superimposed exponential signals.
*Signal Processing*,**39**, 223-233. - 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. (PDF) - Smyth, G. K., and Hawkins, D. M. (1997). Robust frequency estimation using elemental
sets.
*Computing Science and Statistics*,**28**, 659-662. - Kundu, Debasis, and Kannan, Nandini. (1997). Constrained maximum likelihood estimators
for superimposed exponential signals.
*Comm. Statist. Simulation Comput.***26**, no. 2, 733-764. - 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. - Smyth, G. K., and Hawkins, D. M. (2000). Robust frequency estimation using elemental
sets.
*Journal of Computational and Graphics Statistics***9**, 196-214. (Abstract - Software - PDF) - Smyth, G. K. (2000). Employing constraints for improved frequency estimation by
eigenanalysis methods.
*Technometrics***42**, 277-289. (Abstract - PDF)

- StatBox. 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.
- pronyfreq. An S-Plus function for frequency estimation which implements the ORA algorithm of Kahn et al (1992).
- robfreq. An S-Plus function for frequency estimation which implements the multistage algorithm of Smyth and Hawkins (1999).

Home - About Us -
Contact Us Copyright © Gordon Smyth |