Let ${\displaystyle f(t)}$  be a signal consisting of ${\displaystyle N}$  evenly spaced samples. Prony's method fits a function

${\displaystyle {\hat {f}}(t)=\sum _{i=1}^{N}A_{i}e^{\sigma _{i}t}\cos(\omega _{i}t+\phi _{i})}$ 

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

${\displaystyle {\begin{aligned}{\hat {f}}(t)&=\sum _{i=1}^{N}A_{i}e^{\sigma _{i}t}\cos(\omega _{i}t+\phi _{i})\\&=\sum _{i=1}^{N}{\frac {1}{2}}A_{i}\left(e^{j\phi _{i}}e^{\lambda _{i}^{+}t}+e^{-j\phi _{i}}e^{\lambda _{i}^{-}t}\right),\end{aligned}}}$ 

where

${\displaystyle \lambda _{i}^{\pm }=\sigma _{i}\pm j\omega _{i}}$  are the eigenvalues of the system,

${\displaystyle \sigma _{i}=-\omega _{0,i}\xi _{i}}$  are the damping components,

${\displaystyle \omega _{i}=\omega _{0,i}{\sqrt {1-\xi _{i}^{2}}}}$  are the angular-frequency components,

${\displaystyle \phi _{i}}$  are the phase components,

${\displaystyle A_{i}}$  are the amplitude components of the series,

${\displaystyle j}$  is the imaginary unit (${\displaystyle j^{2}=-1}$ ).

Prony's method is essentially a decomposition of a signal with ${\displaystyle M}$  complex exponentials via the following process:

Regularly sample ${\displaystyle {\hat {f}}(t)}$  so that the ${\displaystyle n}$ -th of ${\displaystyle 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 ${\displaystyle {\hat {f}}(t)}$  happens to consist of damped sinusoids, then there will be pairs of complex exponentials such that

${\displaystyle {\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

${\displaystyle {\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 within Prony's method:

1) Construct and solve the matrix equation for the ${\displaystyle 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 ${\displaystyle N\neq 2M}$ , a generalized matrix inverse may be needed to find the values ${\displaystyle P_{m}}$ .

2) After finding the ${\displaystyle P_{m}}$  values, find the roots (numerically if necessary) of the polynomial

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

The ${\displaystyle m}$ -th root of this polynomial will be equal to ${\displaystyle e^{\lambda _{m}}}$ .

3) With the ${\displaystyle e^{\lambda _{m}}}$  values, the ${\displaystyle F_{n}}$  values are part of a system of linear equations that may be used to solve for the ${\displaystyle \mathrm {B} _{m}}$  values:

${\displaystyle {\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 ${\displaystyle M}$  unique values ${\displaystyle k_{i}}$  are used. It is possible to use a generalized matrix inverse if more than ${\displaystyle M}$  samples are used.

Note that solving for ${\displaystyle \lambda _{m}}$  will yield ambiguities, since only ${\displaystyle e^{\lambda _{m}}}$  was solved for, and ${\displaystyle e^{\lambda _{m}}=e^{\lambda _{m}\,+\,q2\pi j}}$  for an integer ${\displaystyle q}$ . This leads to the same Nyquist sampling criteria that discrete Fourier transforms are subject to

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

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