Prony's method

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 sinusoids. This allows for 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\_\{i=1\}^\{N\}\; A\_i\; e^\{sigma\_i\; t\}\; cos(2pi\; f\_i\; t\; +\; phi\_i)$ to the observed $f(t)$. After some manipulation utilizing Euler's formula, the following result is obtained. This allows for more direct computation of terms.:$hat\{f\}(t)\; =\; sum\_\{i=1\}^\{N\}\; A\_i\; e^\{sigma\_i\; t\}\; cos(2pi\; f\_i\; t\; +\; phi\_i)=\; sum\_\{i=1\}^\{N\}\; frac\{1\}\{2\}\; A\_i\; e^\{phi\_i\; j\}e^\{lambda\_i\; t\}$ where::$lambda\_i\; =\; (-sigma\_i\; pm\; j\; omega\_i\; )t$ are the eigenvalues of the system, $sigma\_i$ are the damping components, $phi\_i$ are the phase components, $f\_i$ are the frequency components, $A\_i$ are the amplitude components of the series, and $j=sqrt\{-1\}$.

Example

References

[1] Hauer, J.F. et al (1990). "Initial Results in Prony Analysis of Power System Response Signals". "IEEE Transactions on Power Systems", 5, 1, 80-89.

How to

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 as follows::$F\_n=hat\{f\}(Delta\_t\; n)\; =\; sum\_\{m=1\}^\{M\}\; Beta\_m\; e^\{Lambda\_m\; t\}$

If $hat\{f\}(t)$ happens to be consist of dampened sinusoids then there will be pairs of complex exponentials such that :$Beta\_a\; =\; frac\{1\}\{2\}\; A\_i\; e^\{phi\_i\; j\}$:$Beta\_b\; =\; frac\{1\}\{2\}\; A\_i\; e^\{-phi\_i\; j\}$:$Lambda\_a\; =\; sigma\_i\; +\; j\; omega\_i$ :$Lambda\_b\; =\; sigma\_i\; -\; j\; omega\_i$where:$Beta\_a\; e^\{Lambda\_a\; t\}\; +\; Beta\_b\; e^\{Lambda\_b\; t\}\; =\; frac\{1\}\{2\}\; A\_i\; e^\{phi\_i\; j\}\; e^\{(sigma\_i\; +\; j\; omega\_i)\; t\}\; +\; frac\{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)$

??Because the sumation of complex exponentials is the homogeneous solution to a linear differential equation the following difference equation will exist??::$hat\{f\}(Delta\_t\; n)\; =\; -sum\_\{m=1\}^\{M\}\; hat\{f\}(Delta\_t\; (n-m))\; P\_m$The key to Prony's Method is that the coefficients in the difference equation are related to the following polynomial:$sum\_\{m=1\}^\{M+1\}\; P\_m\; x^\{m-1\}\; =\; prod\_\{m=1\}^\{M\}\; (\; x\; -\; e^\{Lambda\_m\}\; )$

These facts lead to the following three steps to Prony's Method

1) Construct and solve the matrix equation for the $P\_m$ values:

Note that if $N$ ≠ $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 polynomial :$sum\_\{m=1\}^\{M+1\}\; P\_m\; x^\{m-1\}$

The $m^\{th\}$ root of this polynomial will be equal to $e^\{Lambda\_m\}$.

3) With the $e^\{Lambda\_m\}$ values the $F\_n$ values are part of a system of linear equations which may be used to solve for the $Beta\_m$ values:

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 $Lambda\_m$ will yield ambiguities since only $e^\{Lambda\_m\}$ was solved for, and $e^\{Lambda\_m\}=e^\{Lambda\_m+q\; 2\; pi\; j\}$ for and integer $q$. This leads to the same nyquist sampling criteria that discrete fourier transforms are subject to: :

reference:Rob Carriere and Randolph L. Moses, “High Resolution Radar Target Modeling Using a Modified Prony Estimator,” IEEE Trans. Antennas Propogat., vol.40, pp. 13-18, January 1992.