Appendix A: Prony method

Prony method can extract damped complex exponential signals from a given data by solving a set of linear equations (Prony, 1795; Peter and Plonka, 2013; Mitrofanov and Priimenko, 2015; Lobos et al., 2003). Assume the N complex data samples $x[1],x[2],\cdots x[N]$, we approximate the data by $M$ exponential functions:

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

where $A_k$ is the amplitude, $\Delta t$ is the sampling period, $\alpha_k$ is the damping factor, $\omega_k$ is the angular frequency, $\phi_k$ is the initial phase. If we let $h_k = A_k e^{j\phi_k}, z_k = e^{(\alpha_k + j\omega_k)\Delta t}$, we then derive the concise form below:

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

The approximation problem above can be solved based on the error minimization:

$$\mathbf{min} \sum_{n=1}^{N} \left\vert\epsilon[n]\right\vert^2 = ... ...{min}\sum_{n=1}^{N}\left \vert x[n] - \sum_{k=1}^{M}h_k z_k^{n-1}\right\vert^2.$$ (19)

This turns to be a nonlinear problem. It can be solved using Prony method that utilizes linear equation solutions. If there are as many data samples as parameters of the approximation problem, the above M equations 18 can be expressed:

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

 20 can be written in a matrix form as below:

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

Prony proposed to define the polynomial that has the above $z_k, k=1,2, \cdots , M$ as its roots (Prony, 1795):

$$\mathbf{P}(z) = \prod_{k=1}^{M}(z-z_k).$$ (22)

Equation  22 can be rewritten in the form below:

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

Shifting the index on equation 20 from $n$ to $n-m$ , and multiplying by parameter $a[m]$, then we derive:

$$\sum_{m=0}^{M}a_mx[n-m] = \sum_{k=1}^{M}h_kz_k^{n-M-1} \sum_{m=0}^M a[m]z_{k}^{M-m}.$$ (24)

Notice $z_k, k=1,2, \cdots , M$ are roots of equation 23, then equation 24 be written as:

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

Solving equation 25 for the polynomial coefficients. In subsequent steps we compute the frequencies, damping factors and the phases according to Algorithm 1. After all the parameters are computed, we then compute the components of the input signal. For details see Algorithm 1 as follows: