Prony's method explained

Prony analysis (Prony's method) was developed by Gaspard Riche de Prony in 1795. However, practical use of the method awaited the digital computer.^[1] Similar to the Fourier transform, Prony's method extracts valuable information from a uniformly sampled signal and builds a series of damped complex exponentials or damped sinusoids. This allows the estimation of frequency, amplitude, phase and damping components of a signal.

The method

Let

f(t) be a signal consisting of N evenly spaced samples. Prony's method fits a function\hat{f}(t) = \sum ^{N} _{i =1} A _{i} e ^{\sigma _{i} t} \cos(\omega _{i} t + \phi _{i)}

to the observed

f(t) . After some manipulation utilizing Euler's formula, the following result is obtained, which allows more direct computation of terms:\begin{align} \hat{f}(t) & = \sum ^{N} _{i =1} A _{i} e ^{\sigma _{i} t} \cos(\omega _{i} t + \phi _{i)} \\ & = \sum ^{N} _{i =1} \tfrac{ 1 }{ 2 } A _{i} \left( e ^{j\phi _{i}} e ^{λ ^{+} _{i} t} + e ^{- j\phi _{i}} e ^{λ ^{-} _{i} t} \right), \end{align}

where

λ ^{\pm} _{i} = \sigma _{i} \pm j \omega _{i} are the eigenvalues of the system,\sigma _{i} = - \omega _{0,i} \xi _{i} are the damping components,\omega _{i} = \omega _{0,i} \sqrt{ 1- \xi ^{2 }} _{i} are the angular-frequency components,\phi _{i} are the phase components,A _{i} are the amplitude components of the series,j is the imaginary unit (j ^{2} = -1 ).

Representations

Prony's method is essentially a decomposition of a signal with

M complex exponentials via the following process:

Regularly sample

\hat{f}(t) so that the n -th of N samples may be written asF _{n} = \hat{f}(\Delta _{t} n) = \sum ^{M} _{m =1} \Beta _{m} e ^{λ _{m} \Delta _{t} n} , n =0 , ... ,N -1 .

If

\hat{f}(t) happens to consist of damped sinusoids, then there will be pairs of complex exponentials such that\begin{align} \Beta _{a} & = \tfrac{ 1 }{ 2 } A _{i} e ^{\phi _{i} j} , \\ \Beta _{b} & = \tfrac{ 1 }{ 2 } A _{i} e ^{- \phi _{i} j} , \\ λ _{a} & = \sigma _{i} + j \omega _{i,} \\ λ _{b} & = \sigma _{i} - j \omega _{i, \end{align}} where\begin{align} \Beta _{a} e ^{λ _{a} t} + \Beta _{b} e ^{λ _{b} t} & = \tfrac{ 1 }{ 2 } A _{i} e ^{\phi _{i} j} e ^{(\sigma _{i} + j \omega _{i)} t} + \tfrac{ 1 }{ 2 }A _{i} e ^{- \phi _{i} j} e ^{(\sigma _{i} - j \omega _{i)} t} \\ & = A _{i} e ^{\sigma _{i} t} \cos(\omega _{i} t + \phi _{i). \end{align}}

Because the summation of complex exponentials is the homogeneous solution to a linear difference equation, the following difference equation will exist:

\hat{f}(\Delta _{t} n) = \sum ^{M} _{m =1} \hat{f}[\Delta _{t} (n - m)] P _{m,} n = M, ... ,N -1 .

The key to Prony's Method is that the coefficients in the difference equation are related to the following polynomial:

z ^{M} - P _{1} z ^{M-1} - ... - P _{M} = \prod ^{M} _{m =1} \left(z - e ^{λ _{m}} \right).

These facts lead to the following three steps within Prony's method:

1) Construct and solve the matrix equation for the

P _{m} values:\begin{bmatrix} F _{M} \\ \vdots \\ F _{N-1} \end{bmatrix} = \begin{bmatrix} F _{M-1} & ... & F _{0} \\ \vdots & \ddots & \vdots \\ F _{N-2} & ... & F _{N-M-1} \end{bmatrix} \begin{bmatrix} P _{1} \\ \vdots \\ P _{M \end{bmatrix}.}

Note that if

N \ne 2 M , a generalized matrix inverse may be needed to find the values P _{m} .

2) After finding the

P _{m} values, find the roots (numerically if necessary) of the polynomialz ^{M} - P _{1} z ^{M-1} - ... - P _{M.}

The

m -th root of this polynomial will be equal to e ^{λ _{m}} .

3) With the

e ^{λ _{m}} values, the F _{n} values are part of a system of linear equations that may be used to solve for the \Beta _{m} values:\begin{bmatrix} F _{k _{1}} \\ \vdots \\ F _{k _{M}} \end{bmatrix} = \begin{bmatrix} (e ^{λ _{1}} ) ^{k _{1}} & ... & (e ^{λ _{M}} ) ^{k _{1}} \\ \vdots & \ddots & \vdots \\ (e ^{λ _{1}} ) ^{k _{M}} & ... & (e ^{λ _{M}} ) ^{k _{M}} \end{bmatrix} \begin{bmatrix} \Beta _{1} \\ \vdots \\ \Beta _{M \end{bmatrix},} where M unique values k _{i} are used. It is possible to use a generalized matrix inverse if more than M samples are used.

Note that solving for

λ _{m} will yield ambiguities, since only e ^{λ _{m}} was solved for, and e ^{λ _{m}} = e ^{λ _{m} + q 2 \pi j} for an integer q . This leads to the same Nyquist sampling criteria that discrete Fourier transforms are subject to\left|\operatorname{Im}( λ _{m)\right|} = \left|\omega _{m\right|} < (\pi)/(\Delta _{t}) .

See also

• Generalized pencil-of-function method
• Computation of Prony decomposition using Autoregression analysis
• Application of Prony decomposition in Time-frequency analysis

References

• 13–18. 10.1109/8.123348. High resolution radar target modeling using a modified Prony estimator. IEEE Transactions on Antennas and Propagation. 40. 1992. Carriere. R.. Moses. R. L..
• 35–39. Interpretation of the Proni method for solving long-range problems. Radioelectronics and Communications Systems. 41. 1998. Slyusar. V. I.. 1.

Notes and References

1. 10.1109/59.49090 . 80–89 . Initial results in Prony analysis of power system response signals . IEEE Transactions on Power Systems . 5 . 1990 . Hauer . J. F. . Demeure . C. J. . Scharf . L. L. . 1 . 1990ITPSy...5...80H . 10217/753 . free.