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}^{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.

$\begin{align} \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{align}$

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

$F_n = \hat{f}(\Delta_t n) = \sum_{m=1}^{M} \Beta_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{align} \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, \end{align}$

where

$\begin{align} \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). \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=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:

$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:

$\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 \ne 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

$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 $\Beta_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} \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 $\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 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. et al. (1990). "Initial Results in Prony Analysis of Power System Response Signals". IEEE Transactions on Power Systems, 5, 1, 80-89.

References

• 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.