Non-stationary Prony method

Equation 7 can be written as:

$$\sum_{m=1}^{M}\hat{a}_m x[n-m] = x[n].$$ (9)

If the $\hat{a}_m$ in equation 9 are time dependent, then we have:

$$\sum_{m=1}^{M}\hat{a}_m[n] x[n-m] \approx x[n],$$ (10)

which is an under-determined linear system. There are many methods for solving under-determined linear system, such as Tikhonov method (Tikhonov, 1963). In this paper, we apply shaping regularization (Fomel, 2009,2007) to regularize the under-determined linear system, and obtain (for details see Appendix):

$$\mathbf{\hat{a}} = \mathbf{F}^{-1}\mathbf{\eta},$$ (11)

where $\mathbf{\hat{a}}$ is a vector composed of $\hat{a}_m [n]$, the elements of vector $\mathbf{\eta}$ are $\eta_i[n] = \mathbf{S}[x_i^* [n] x[n]]$, where $x_i[n] = x[n-i]$, $x_i^*[n]$ stands for the complex conjugate of $x_i[n]$ and $\mathbf{S}$ is the shaping operator. The elements of matrix $\mathbf{F}$ are:

$${F}_{ij}[n] = \sigma^2 \delta_{ij} + \mathbf{S}[x_i^*[n]x_j[n] - \sigma^2 \delta_{ij}],$$ (12)

where $\sigma$ is the regularization parameter. Solving equation 11, we obtain the coefficients vector $\hat{a}_m [n]$ and form a polynomial below:

$$\mathbf{P}(z) = z^M + \hat{a}_1[n] z^{M-1} + \cdots + \hat{a}_M[n].$$ (13)

For the roots computation $\hat{z}_m[n], m=1,2,\cdot,M$ of the above polynomial, we use the method proposed by Toh and Trefethen (1994). The instantaneous frequency of each different component is derived from the following equation:

$$f_m[n] = \Re \left[ \arg \left( \frac{\hat{z}_m[n]}{2\pi \Delta t} \right) \right].$$ (14)

From the instantaneous frequency, we compute the local phase according to the following equation:

$$\Phi_m[n] = 2\pi\sum_{k=0}^{n}f_m[k]\Delta t.$$ (15)

Solving the following equation using regularized non-stationary regression method (Fomel, 2013):

$$x[n] = \sum_{m=1}^{M} \hat{A}_m[n]e^{j \Phi_m[n]} = \sum_{m=1}^{M} c_m[n].$$ (16)

Finally the narrow-band intrinsic mode functions $c_m[n]$ are computed based on equation 16