Prony Method

Prony method extracts damped complex exponential functions (or sinusoids) from a given signal by solving a set of linear equations (Lobos et al., 2003; Mitrofanov and Priimenko, 2015; Peter and Plonka, 2013; Prony, 1795). The Prony method allows for estimation of frequencies, amplitudes and phases of a signal (For details see Appendix). Assume we want to solve the problem:

$\displaystyle x[n] = \sum_{k=1}^{M}A_k e^{(\alpha_k + j\omega_k)(n-1)\Delta t + j\phi_k},$ (2)

if let $h_k = A_k e^{j\phi_k}, z_k = e^{(\alpha_k + j\omega_k)\Delta t}$, we derive the concise form

$\displaystyle x[n] = \sum_{k=1}^{M}h_k z_k^{n-1}.$ (3)

The above M equations can be written in a matrix form:

$\displaystyle \left[ \begin{array}{cccc}
z_1^0 & z_2^0 & \cdots & z_M^0\\
z_1^...] =
\left[ \begin{array}{c} x[1] x[2] \vdots  x[M]
\end{array} \right].$ (4)

The above $z_k, k=1,2, \cdots , M$ of equation 4 can be computed by solving a polynomial of the form:

$\displaystyle \mathbf{P}(z) = \prod_{k=1}^{M}(z-z_k),$ (5)

equation 5 can also be written in a form:

$\displaystyle \mathbf{P}(z) = a_0 z^M + a_1 z^{M-1} + \cdots + a_{M-1}z + a_M.$ (6)

The coefficients $a_k$ of the polynomial can be computed by solving the following equation:

$\displaystyle \sum_{m=0}^{M}a_mx[n-m] = 0.$ (7)

We use the method proposed by Toh and Trefethen (1994) to compute the roots $z_k$ of equation 6. If the roots are solved, the $h_k$ can be computed using equation 3. Finally, the components are computed based on the equation below(For details see Appendix):

$\displaystyle c_k[n] = h_kz_k^{n-1} = h_kz_k^{n-1}, k=1,2,\cdots,M.$ (8)