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

${\displaystyle {\hat {f}}(t)=\sum _{i=1}^{M}A_{i}e^{\sigma _{i}t}\cos(2\pi f_{i}t+\phi _{i})}$

to the observed $f(t)$. After some manipulation utilizing Euler's formula, the following result is obtained. This allows more direct computation of terms.

${\displaystyle {\begin{aligned}{\hat {f}}(t)&=\sum _{i=1}^{M}A_{i}e^{\sigma _{i}t}\cos(2\pi f_{i}t+\phi _{i})\\&=\sum _{i=1}^{M}{\frac {1}{2}}A_{i}e^{\pm j\phi _{i}}e^{\lambda _{i}t}\end{aligned}}}$

where:

• $\lambda _{i}=\sigma _{i}\pm j\omega _{i}$ 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$ 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 as

${\displaystyle F_{n}={\hat {f}}(\Delta _{t}n)=\sum _{m=1}^{M}\mathrm {B} _{m}e^{\lambda _{m}\Delta _{t}n},\quad n=0,\dots ,N-1.}$

If ${\hat {f}}(t)$ happens to consist of damped sinusoids, then there will be pairs of complex exponentials such that

${\begin{aligned}\mathrm{B} _{a}&={\frac {1}{2}}A_{i}e^{{\phi _{i}j}},\\\mathrm{B} _{b}&={\frac {1}{2}}A_{i}e^{{-\phi _{i}j}},\\\lambda _{a}&=\sigma _{i}+j\omega _{i},\\\lambda _{b}&=\sigma _{i}-j\omega _{i},\end{aligned}}$

where

${\begin{aligned}\mathrm{B} _{a}e^{{\lambda _{a}t}}+\mathrm{B} _{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}).\end{aligned}}$

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

${\displaystyle {\hat {f}}(\Delta _{t}n)=-\sum _{m=1}^{M}{\hat {f}}[\Delta _{t}(n-m)]P_{m},\quad n=M,\dots ,N-1.}$

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

${\displaystyle z^{M}+P_{1}z^{M-1}+\dots +P_{M}=\prod _{m=1}^{M}\left(z-e^{\lambda _{m}}\right).}$

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

1) Construct and solve the matrix equation for the $P_m$ values:

${\displaystyle {\begin{bmatrix}F_{M}\\\vdots \\F_{N-1}\end{bmatrix}}={\begin{bmatrix}F_{M-1}&\dots &F_{0}\\\vdots &\ddots &\vdots \\F_{N-2}&\dots &F_{N-M-1}\end{bmatrix}}{\begin{bmatrix}P_{1}\\\vdots \\P_{M}\end{bmatrix}}.}$

Note that if $N\neq 2M$, 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

${\displaystyle z^{M}+P_{1}z^{M-1}+\dots +P_{M},}$

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 that may be used to solve for the $\mathrm{B} _{m}$ values:

${\begin{bmatrix}F_{{k_{1}}}\\\vdots \\F_{{k_{M}}}\end{bmatrix}}={\begin{bmatrix}(e^{{\lambda _{1}}})^{{k_{1}}}&\dots &(e^{{\lambda _{M}}})^{{k_{1}}}\\\vdots &\ddots &\vdots \\(e^{{\lambda _{1}}})^{{k_{M}}}&\dots &(e^{{\lambda _{M}}})^{{k_{M}}}\end{bmatrix}}{\begin{bmatrix}\mathrm{B} _{1}\\\vdots \\\mathrm{B} _{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 $\lambda_m$ will yield ambiguities, since only $e^{{\lambda _{m}}}$ was solved for, and $e^{{\lambda _{m}}}=e^{{\lambda _{m}\,+\,q2\pi j}}$ for an integer $q$. This leads to the same Nyquist sampling criteria that discrete Fourier transforms are subject to:

$\left|Im(\lambda _{m})\right|=\left|\omega _{m}\right|<{\frac {1}{2\Delta _{t}}}.$

Example

Notes

1. Hauer, J.F.; Demeure, C.J.; Scharf, L.L. (1990). "Initial results in Prony analysis of power system response signals". IEEE Transactions on Power Systems. 5: 80–89. doi:10.1109/59.49090.

References

• Carriere, R.; Moses, R.L. (1992). "High resolution radar target modeling using a modified Prony estimator". IEEE Transactions on Antennas and Propagation. 40: 13–18. doi:10.1109/8.123348.