Improving the Accuracy of the Pencil of Function Method Increasing Its Matrix Polynomial Degree

0.1. Background

The accurate representation of a system’s transient signal as a sum of damped complex harmonic functions presents a significant numerical challenge that is essential for solving various engineering and physics problems. A critical aspect of this challenge is estimating the natural complex resonance frequencies (NCRF) of the linear system under investigation. Numerous techniques have been proposed to address this issue, including the Fast Fourier Transform, Prony’s method, Total Least Squares Prony’s method, and more recently, the Generalized Pencil of Functions (GPOF) method, also known as the Matrix Pencil method. The GPOF method, introduced by Hua and Sarkar in 1989 [1], was enhanced in 1995 to improve its performance in noisy environments [2]. Due to its computational efficiency and robustness against noise [3], the GPOF method has been widely adopted for various practical applications, including power system harmonics monitoring, biomedical monitoring, remote vital sign surveying, magnetic resonance imaging, super-resolution microscopy, noninvasive blood glucose monitoring, gravitational wave detection, radar, identification of concealed firearms, military aircraft recognition, sonar, seismic data processing, meta-material design, and Coriolis mass flow meter signal processing [4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21]. In the geosciences and remote sensing fields, the GPOF method has proven useful for tasks such as localization and identification of buried objects, soil characterization, and topographic studies through radar techniques [22,23,24].

In this context, some applications may demand greater precision than the conventional GPOF algorithm can offer in order to extract valuable information from transient signals. This is particularly relevant for engineering challenges associated with radar and spectroscopy system design, as well as image super-resolution processing in medical imaging, microscopy, and astronomy [25,26,27,28].

In scenarios where applications require the discrimination of closely spaced complex natural frequencies, it may be also necessary to employ a more precise algorithm. An example of this is found in super-resolution microscopy, where very similar complex natural frequencies need to be estimated [11].

Conventionally, the GPOF method utilizes a linear matrix pencil, which involves a polynomial matrix of degree one [1,2]. In this work, we develop a GPOF method based on a polynomial matrix pencil of degree greater than one. We demonstrate that this new nonlinear GPOF approach can enhance computational accuracy for applications requiring improved precision, provided that the transient response of the system is measured and sampled with low noise levels.

0.2. Conventional Generalized Pencil of Function (GPOF) Method Revisited

Let $y\left(t\right)$ a signal that represents the transient response of a linear system, which can be expressed as a sum of damped complex harmonics plus noise $n\left(t\right)$:

$$y\left(t\right)=\sum _{n=1}^N{r}_{n}exp\left(s_{n}t\right)+n\left(t\right).$$

(1)

Then, the data set $y\left[0\right],y\left[1\right],\dots ,y\left[m\right],\dots ,y[M-1]$ is generated by means of an uniform sampling of $y\left(t\right)$, as follows:

$$y\left[m\right]=y\left(m{t}_s\right)=\sum _{n=1}^N{r}_{n}exp\left(m\left(s_n{t}_s\right)\right)=\sum _{n=1}^N{r}_n{z}_n^m,$$

(2)

where $t_s$ is the sampling interval, N is the number of complex natural resonance frequencies (CNRFs) or poles, $s_n$ is the $n^{\mathrm{th}}$ system NCRF, $r_n$ is the $n^{\mathrm{th}}$ complex residue, and $z_n=exp\left(s_n{t}_s\right)$ is the $n^{\mathrm{th}}$ discrete complex pole. Here, M is the total number of samples, and m is an integer such that $1\le m<M$. It is important to note that ${\sigma }_n=-\Re \left(s_n\right)$ represents the damping factor of the $n^{\mathrm{th}}$ harmonic, while ${\omega }_n=\Im \left(s_n\right)$ corresponds to the natural oscillation frequency of resonance for that harmonic.

Next, we highlight some fundamental aspects of the linear algebra associated with the conventional GPOF method. A matrix pencil of degree $\Omega$ is defined as a polynomial function of the form:

$$\mathcal{L}\left(\lambda \right)=\sum _{k=0}^{\Omega }{\lambda }^k{\mathbf{A}}_k={\mathbf{A}}_0+\lambda {\mathbf{A}}_1+{\lambda }^2{\mathbf{A}}_2+\dots +{\lambda }^{\Omega }{\mathbf{A}}_{\Omega },$$

(3)

where the polynomial coefficients ${\mathbf{A}}_0,{\mathbf{A}}_1,\dots ,{\mathbf{A}}_{\Omega }$ are real rectangular matrices and $\Omega$ is a positive integer. Here, $\lambda$ is a complex scalar [29]. The conventional linear GPOF method is derived from a matrix pencil function of degree $\Omega =1$:

$$\mathcal{L}\left(\lambda \right)={\mathbf{A}}_0+\lambda {\mathbf{A}}_1,$$

(4)

where ${\mathbf{A}}_0$ and ${\mathbf{A}}_1$ are Hankel matrices containing samples from the transient response $y\left(t\right)$. The objective of this approach is to solve the generalized eigenvalue problem:

$$\mathcal{L}\left(\lambda \right)\mathbf{x}=({\mathbf{A}}_0+\lambda {\mathbf{A}}_1)\mathbf{x}=0,$$

(5)

where each computed value of $\lambda$ corresponds to each one of the N possible complex eigenvalue-eigenvector pairs represented by $z_n$ , for $n=1,2,\dots ,N$ [1].

1. Proposed Generalized Pencil of Function Method with an Increased Matrix Polynomial Degree

1.1. Definitions

In this work, we propose a matrix pencil polynomial of the form

$$\mathcal{L}\left(\lambda \right)={\mathbf{A}}_0+{\lambda }^{\Omega }{\mathbf{A}}_{\Omega },$$

(6)

where $\Omega$ is a positive integer. This leads to the solution of a generalized nonlinear eigenvalue problem given by

$$\mathcal{L}\left(\lambda \right)\mathbf{x}=({\mathbf{A}}_0+{\lambda }^{\Omega }{\mathbf{A}}_{\Omega })\mathbf{x}=0,$$

(7)

which is aimed at determining the complex natural frequencies of the linear system under study from its transient response $y\left(t\right)$. To achieve this, we follow a procedure analogous to the conventional linear matrix pencil method.

First, we organize the data samples appropriately. A new Hankel matrix, denoted as $\mathbf{Y}$, is defined with dimensions $(M-L-\Omega +1)\times (L+\Omega )$:

$$\mathbf{Y}={\left(\begin{array}{cccc}y\left[0\right]& y\left[1\right]& \cdots & y[L+\Omega -1]\\ y\left[1\right]& y\left[2\right]& \cdots & y[L+\Omega ]\\ ⋮& ⋮& \ddots & ⋮\\ y[M-L-\Omega -1]& y[M-L-\Omega ]& \cdots & y[M-2]\\ y[M-L-\Omega ]& y[M-L-\Omega +1]& \cdots & y[M-1]\end{array}\right)}_{(M-L-\Omega +1)\times (L+\Omega )}.$$

(8)

This matrix $\mathbf{Y}$ includes the entire dataset $y\left[0\right],y\left[1\right],\dots ,y[M-1]$, where L is a positive integer known as the pencil parameter, satisfying the condition $M>(L+\Omega )$. Additionally, the $j^{th}$ column of $\mathbf{Y}$ is represented by the vector ${\mathbf{y}}_j$, defined as

$${\mathbf{y}}_j={(y[j-1],y\left[j\right],\dots ,y[M-L-\Omega +j-2],y[M-L-\Omega +j-1])}^H,$$

(9)

for $j=1,2,\dots ,(M-L-\Omega +1)$, where the superscript H denotes the conjugate transpose. Thus, the Hankel matrix $\mathbf{Y}$ can be expressed as

$$\mathbf{Y}={({\mathbf{y}}_0,{\mathbf{y}}_1,\dots ,{\mathbf{y}}_{L+\Omega -2},{\mathbf{y}}_{L+\Omega -1})}_{(M-L-\Omega +1)\times (L+\Omega )}.$$

(10)

To implement the proposed method, we need to construct a generic matrix ${\mathbf{Y}}_k$ by extracting $\Omega$ columns from $\mathbf{Y}$. The proposed structure for ${\mathbf{Y}}_k$ is given by

$${\mathbf{Y}}_k={({\mathbf{y}}_{k-1},{\mathbf{y}}_k,{\mathbf{y}}_{k+1},\dots ,{\mathbf{y}}_{k+L-2})}_{(M-L-\Omega +1)\times L},$$

(11)

for $k=1,2,\dots ,(\Omega +1)$. We also define a set of $\Omega +1$ Hankel matrices derived from Eq. (11):

$${\mathbf{Y}}_1,{\mathbf{Y}}_2,{\mathbf{Y}}_3,\dots ,{\mathbf{Y}}_k,\dots ,{\mathbf{Y}}_{\Omega +1}.$$

(12)

Some specific elements of this set can be represented as follows:

$${\mathbf{Y}}_1={({\mathbf{y}}_0,{\mathbf{y}}_1,\dots ,{\mathbf{y}}_{L-2},{\mathbf{y}}_{L-1})}_{(M-L-\Omega +1)\times L},$$

(13)

$${\mathbf{Y}}_2={({\mathbf{y}}_1,{\mathbf{y}}_2,\dots ,{\mathbf{y}}_{L-1},{\mathbf{y}}_L)}_{(M-L-\Omega +1)\times L},$$

(14)

$${\mathbf{Y}}_3={({\mathbf{y}}_2,{\mathbf{y}}_3,\dots ,{\mathbf{y}}_L,{\mathbf{y}}_{L+1})}_{(M-L-\Omega +1)\times L},$$

(15)

$${\mathbf{Y}}_{\Omega }={({\mathbf{y}}_{\Omega -1},{\mathbf{y}}_{\Omega },\dots ,{\mathbf{y}}_{L+\Omega -3},{\mathbf{y}}_{L+\Omega -2})}_{(M-L-\Omega +1)\times L},$$

(16)

$${\mathbf{Y}}_{\Omega +1}={({\mathbf{y}}_{\Omega },{\mathbf{y}}_{\Omega +1},\dots ,{\mathbf{y}}_{L+\Omega -2},{\mathbf{y}}_{L+\Omega -1})}_{(M-L-\Omega +1)\times L}.$$

(17)

It is worth highlighting that ${\mathbf{Y}}_1$ comprises the first L columns of the matrix $\mathbf{Y}$ and that ${\mathbf{Y}}_{\Omega +1}$ comprises the last L columns of the matrix $\mathbf{Y}$.

1.2. Properties of the Hankel Matrices Set

In this section, we leverage certain intrinsic properties of uniformly sampled transient exponential signals and Hankel matrices in order to propose some expressions that are useful to develop the proposed method. As in [1], it is convenient to write that

$${\mathbf{Y}}_1={\mathbf{Z}}_A\mathbf{R}{\mathbf{Z}}_B,$$

(18)

$${\mathbf{Y}}_2={\mathbf{Z}}_A\mathbf{R}{\mathbf{Z}}_0{\mathbf{Z}}_B,$$

(19)

$${\mathbf{Z}}_0=exp\left(t_s\mathbf{S}\right),$$

(20)

where the matrices ${\mathbf{Z}}_A$, ${\mathbf{Z}}_B$, $\mathbf{R}$, $\mathbf{S}$, and ${\mathbf{Z}}_0$ are defined as follows:

$${\mathbf{Z}}_A={\left(\begin{array}{ccccc}1& 1& \cdots & 1& 1\\ z_1& z_2& \cdots & z_{N-1}& z_N\\ z_1^2& z_2^2& \cdots & z_{N-1}^2& z_N^2\\ ⋮& ⋮& \ddots & ⋮& ⋮\\ z_1^{M-L-\Omega }& z_2^{M-L-\Omega }& \cdots & z_{N-1}^{M-L-\Omega }& z_N^{M-L-\Omega }\end{array}\right)}_{(M-L-\Omega +1)\times N},$$

(21)

$${\mathbf{Z}}_B={\left(\begin{array}{ccccc}1& z_1& z_1^2& \cdots & z_1^{L-1}\\ 1& z_2& z_2^2& \cdots & z_2^{L-1}\\ ⋮& ⋮& ⋮& ⋮& ⋮\\ 1& z_{N-1}& z_{N-1}^2& \cdots & z_{N-1}^{L-1}\\ 1& z_N& z_N^2& \cdots & z_N^{L-1}\end{array}\right)}_{N\times L},$$

(22)

$$\mathbf{R}=\mathrm{diag}[r_1,r_2,\dots ,r_{N-1},r_N],$$

(23)

$$\mathbf{S}=\mathrm{diag}[s_1,s_2,\dots ,s_{N-1},s_N],$$

(24)

$${\mathbf{Z}}_0=\mathrm{diag}[z_1,z_2,\dots ,z_{N-1},z_N]=exp\left(t_s\mathbf{S}\right).$$

(25)

Next, we define ${\mathbf{Z}}_B^{†}$ as the Moore-Penrose pseudo-inverse of ${\mathbf{Z}}_B$. A relationship between ${\mathbf{Y}}_2$ and ${\mathbf{Y}}_1$ can then be established from Eqs. (18) and (19):

$${\mathbf{Y}}_2={\mathbf{Y}}_1\left({\mathbf{Z}}_B^{†}{\mathbf{Z}}_0{\mathbf{Z}}_B\right).$$

(26)

Note that ${\mathbf{Y}}_1$ in Eq. (13) consists of the first L columns of the matrix $\mathbf{Y}$. Moreover, due to the properties of Hankel matrices, ${\mathbf{Y}}_2$ can be constructed by eliminating the first column of ${\mathbf{Y}}_1$ (denoted as ${\mathbf{y}}_0$) and appending the column ${\mathbf{y}}_{L+1}$ as the last column. This implies that, relative to $\mathbf{Y}$, ${\mathbf{Y}}_2$ is shifted by one column compared to ${\mathbf{Y}}_1$. Consequently, each element of the matrix ${\mathbf{Y}}_2$ is shifted by one sampling interval compared to its corresponding element in ${\mathbf{Y}}_1$.

Similarly, it can be inferred from Eqs. (14) and (15) that ${\mathbf{Y}}_3$ is also shifted by one column with respect to ${\mathbf{Y}}_2$. Therefore, by induction, the relationship between ${\mathbf{Y}}_3$ and ${\mathbf{Y}}_2$ mirrors that between ${\mathbf{Y}}_2$ and ${\mathbf{Y}}_1$. Hence, we can express:

$${\mathbf{Y}}_3={\mathbf{Y}}_2\left({\mathbf{Z}}_B^{†}{\mathbf{Z}}_0{\mathbf{Z}}_B\right).$$

(27)

Then, a generalized relationship can be derived between ${\mathbf{Y}}_{k+1}$ and ${\mathbf{Y}}_k$:

$${\mathbf{Y}}_{k+1}={\mathbf{Y}}_k\left({\mathbf{Z}}_B^{†}{\mathbf{Z}}_0{\mathbf{Z}}_B\right).$$

(28)

Substituting the expression for ${\mathbf{Y}}_2$ from Eq. (19) into Eq. (27), it results:

$${\mathbf{Y}}_3={\mathbf{Z}}_A\mathbf{R}{\mathbf{Z}}_0^2{\mathbf{Z}}_B.$$

(29)

Finally, by examining Eqs. (18), (19), and (29), we can propose a generalized expression for ${\mathbf{Y}}_{k+1}$:

$${\mathbf{Y}}_{k+1}={\mathbf{Z}}_A\mathbf{R}{\mathbf{Z}}_0^k{\mathbf{Z}}_B.$$

(30)

1.3. Defining the Eigenvalue Problem

Since the diagonal matrix ${\mathbf{Z}}_0$ contains the values of $z_n$ to be determined, it is important to note that the product ${\mathbf{Y}}_1^{†}{\mathbf{Y}}_{k+1}$, derived from Eq. (30) and the inversion of Eq. (18), can be expressed as

$${\mathbf{Y}}_1^{†}{\mathbf{Y}}_{k+1}={\mathbf{Z}}_B^{†}{\mathbf{Z}}_0^k{\mathbf{Z}}_B,$$

(31)

where ${\mathbf{Z}}_0^k$ is a diagonal matrix containing the eigenvalues of ${\mathbf{Y}}_1^{†}{\mathbf{Y}}_{k+1}$. Consequently, Eq. (31) can also be reformulated as

$${\mathbf{Y}}_{k+1}-{\mathbf{Y}}_1{\mathbf{Z}}_B^{†}{\mathbf{Z}}_0^k{\mathbf{Z}}_B=0.$$

(32)

This leads us to define a new generalized eigenvalue problem as follows:

$$({\mathbf{Y}}_{k+1}-z_n^k{\mathbf{Y}}_1){\mathbf{x}}_n=0,$$

(33)

where each $z_n^k$ represents the eigenvalues of the linear system to be computed, and ${\mathbf{x}}_n$ denotes their associated eigenvectors. Notably, if $k=\Omega$, the expression in parentheses in Eq. (33) assumes the form of the matrix pencil of degree $\Omega$ defined in Eq. (6).

Furthermore, if we left-multiply Eq. (33) by ${\mathbf{Y}}_1^{†}$ and set $k=\Omega$, denoting $z_n^{\Omega }=\xi$, we can express it as

$$({\mathbf{Y}}_1^{†}{\mathbf{Y}}_{\Omega +1}-{\xi }_n\mathbf{I}){\mathbf{x}}_n=0,$$

(34)

where $\mathbf{I}$ is an $L\times L$ identity matrix.

1.4. Computing the Values of the Complex Natural Frequencies

In this work, the eigenvalue problem defined in Eq. (34) is not solved directly. Instead, to mitigate the adverse effects of noise present in $y\left(t\right)$, we decompose the matrix $\mathbf{Y}$ using its singular value decomposition (SVD) as described in [2]:

$$\mathbf{Y}=\mathbf{U}\Sigma {\mathbf{V}}^H,$$

(35)

where H denotes the conjugate transpose operator, and $\Sigma$ is a diagonal matrix whose diagonal elements, denoted as ${\varsigma }_n$, are the singular values of $\mathbf{Y}$ for $n=1,2,\dots ,N$. These singular values are organized in $\Sigma$ in descending order, such that ${\varsigma }_1>{\varsigma }_2>\dots >{\varsigma }_N$.

Next, we compute the ratios ${\displaystyle \frac{{\varsigma }_n}{{\varsigma }_1}}$ for $n=1,2,\dots ,N$, and define a truncation criterion as follows:

$${\displaystyle \frac{{\varsigma }_T}{{\varsigma }_1}}\approx {10}^{-p},$$

(36)

where T is an integer satisfying $N\ge T\ge 1$, and p represents the number of significant decimal digits in the data.

We then define truncated versions of the matrices $\mathbf{U}$, $\Sigma$, and $\mathbf{V}$ as ${\mathbf{U}}_T$, ${\Sigma }_T$, and ${\mathbf{V}}_T$, respectively. The matrices ${\mathbf{U}}_T$ and ${\mathbf{V}}_T$ are formed by taking the first T columns of $\mathbf{U}$ and $\mathbf{V}$, while ${\Sigma }_T$ consists of the first T rows and columns of $\Sigma$ [2].

To filter out noise, we obtain a less noisy but approximate matrix $\mathbf{Y}$ as follows:

$$\mathbf{Y}\approx {\mathbf{U}}_T{\Sigma }_T{\mathbf{V}}_T^H.$$

(37)

Following a similar methodology to that outlined in [2], which enhances the accuracy of the GPOF method in noisy environments, both matrices ${\mathbf{Y}}_1$ and ${\mathbf{Y}}_{\Omega +1}$ must retain L columns.

Thus, we can express the singular value decompositions of ${\mathbf{Y}}_1$ and ${\mathbf{Y}}_{\Omega +1}$ as follows:

$${\mathbf{Y}}_1\approx {\mathbf{U}}_T{\Sigma }_T{\mathbf{V}}_1^H,$$

(38)

$${\mathbf{Y}}_{\Omega +1}\approx {\mathbf{U}}_T{\Sigma }_T{\mathbf{V}}_{\Omega +1}^H,$$

(39)

where ${\mathbf{V}}_1$ is constructed by retaining the first L rows of ${\mathbf{V}}_T$, and ${\mathbf{V}}_{\Omega +1}$ is formed from the last L rows of ${\mathbf{V}}_T$.

From Eqs. (34), (38), and (39), we can then define and solve the following eigenvalue problem computationally:

$$\left[{\left({\mathbf{V}}_1^H\right)}^{†}{\mathbf{V}}_{\Omega +1}^H-z_n^{\Omega }\mathbf{I}\right]{\mathbf{x}}_n=0.$$

(40)

Letting ${\xi }_n=z_n^{\Omega }$ represent the eigenvalues of ${\left({\mathbf{V}}_1^H\right)}^{†}{\mathbf{V}}_{\Omega +1}^H$ for $n=1,2,\dots ,N$, we note that there are multiple possible solutions for any $z_n$. Since $z_n=exp\left(s_n{t}_s\right)$ with $s_n=-{\sigma }_n+ı{\omega }_n$, and ${\xi }_n=\left|{\xi }_n\right|exp[ı(arg\left({\xi }_n\right)+2\pi k^{\prime })]$ for any integer $k^{\prime }$, the multiplicity of possible values for $z_n$ and the imaginary part of $s_n$ can be established by rewriting the equality ${\xi }_n=z_n^{\Omega }$ as follows:

$$|{\xi }_n|exp(ı(arg\left({\xi }_n\right)+2\pi k^{\prime }))=exp(-{\sigma }_n{t}_s\Omega )exp(ı{\omega }_n{t}_s\Omega ).$$

(41)

Next, we can solve Eq. (41) for the damping factor ${\sigma }_n$ and the angular resonance frequency ${\omega }_n$:

$${\sigma }_n=-\Re \left\{s_n\right\}=-{\displaystyle \frac{ln|{\xi }_n|}{\Omega t_s}},$$

(42)

$${\omega }_n\left(k^{\prime }\right)=\Im \left\{s_n\left(k^{\prime }\right)\right\}={\displaystyle \frac{arg\left({\xi }_n\right)}{\Omega t_s}}+{\displaystyle \frac{2\pi k^{\prime }}{\Omega t_s}}.$$

(43)

While ${\sigma }_n$, as given by Eq. (42), is unique, ${\omega }_n$, as given by Eq. (43), is not, warranting further investigation. According to the Nyquist sampling criterion, $|{\omega }_n\left(k^{\prime }\right)|<{\displaystyle \frac{\pi }{t_s}}$, constraining all possible values of $k^{\prime }$ to the interval:

$${\displaystyle \frac{\Omega -arg\left({\xi }_n\right)/\pi }{2}}>k^{\prime }>{\displaystyle \frac{-\Omega -arg\left({\xi }_n\right)/\pi }{2}}.$$

(44)

Thus, $z_n\left(k^{\prime }\right)$, $s_n\left(k^{\prime }\right)$, and ${\omega }_n\left(k^{\prime }\right)$ are functions of the integer variable $k^{\prime }$. This solution multiplicity is intrinsic to any matrix pencil polynomial of degree $\Omega >1$. Designating that $s_n\left(k_0^{\prime }\right)$ is the true value among the set of possibilities of $s_n\left(k^{\prime }\right)$, it is useful to define the scalar $\chi \left(k^{\prime }\right)$ as the absolute value of the determinant of the matrix pencil expressed in Eq. (34), evaluated for $\Omega =1$ (the conventional linear case):

$$\chi \left(k^{\prime }\right)=\left|det\left({\mathbf{Y}}_1^{†}{\mathbf{Y}}_2-z_n\left(k^{\prime }\right)\mathbf{I}\right)\right|=\left|det\left({\mathbf{Y}}_1^{†}{\mathbf{Y}}_2-exp\left(s_n\left(k^{\prime }\right)t_s\right)\mathbf{I}\right)\right|.$$

(45)

According to linear algebra theory, in the absence of noise, a true complex natural frequency $s_n\left(k_0^{\prime }\right)$ and its related eigenvalue $z_n\left(k_0^{\prime }\right)$ yield $\chi \left(k^{\prime }\right){|}_{k^{\prime }=k_0^{\prime }}=0$.

In practical cases, however, noise is always present in the signal $y\left(t\right)$. Thus, the correct solution for $z_n\left(k_0^{\prime }\right)$ and $s_n\left(k_0^{\prime }\right)$ is found at $k^{\prime }=k_0^{\prime }$, where $k^{\prime }$ belongs to the set of integer values defined by the restriction (44). The value $k^{\prime }=k_0^{\prime }$ must minimize the function $\chi \left(k^{\prime }\right)$ (see Figure 1). Consequently, the appropriate values of ${\omega }_n\left(k_0^{\prime }\right)$, $s_n\left(k_0^{\prime }\right)$, and $z_n\left(k_0^{\prime }\right)$, which characterize the studied signal $y\left(t\right)$, can be estimated. This process must be repeated to determine each one of the N values of $z_n$, ${\omega }_n$, and $s_n$, since $n=1,2,\dots ,N$.

If necessary, the residues of the signal can also be estimated as the diagonal elements of the matrix $\mathbf{R}$ obtained from Eq. (18):

$$\mathbf{R}=\left({\mathbf{Z}}_B^{†}{\mathbf{Y}}_1{\mathbf{Z}}_A^{†}\right).$$

(46)

In this way, the analytic approach to the system transient signal $y\left(t\right)$ defined in Eq. (1) is finally established.

1.5. Relation Between the Proposed Method and the Conventional GPOF Method

The well-known conventional GPOF method is a specific case of the previously proposed method, obtained by setting the pencil degree to one ($\Omega =1$). To illustrate this, we can evaluate the matrices $\mathbf{Y}$, ${\mathbf{Y}}_1$, ${\mathbf{Y}}_{\Omega +1}$, ${\mathbf{Z}}_A$, and ${\mathbf{Z}}_B$ at $\Omega =1$. The matrix $\mathbf{Y}$ given by Eq. (8) then takes the form:

$$\mathbf{Y}{|}_{\Omega =1}={\left(\begin{array}{cccc}y\left[0\right]& y\left[1\right]& \dots & y\left[L\right]\\ y\left[1\right]& y\left[2\right]& \dots & y[L+1]\\ ⋮& ⋮& & ⋮\\ y[M-L]& y[M-L-1]& \dots & y[M-2]\\ y[M-L-1]& y[M-L]& \dots & y[M-1]\end{array}\right)}_{(M-L)\times (L+1)},$$

(47)

which corresponds to the Hankel matrix $\mathbf{Y}$ proposed in [1]. Similarly, the matrix set defined in Eq. (12) for $\Omega =1$ consists of only two matrices, ${\mathbf{Y}}_1$ and ${\mathbf{Y}}_2$, which are the same as those defined for the conventional GPOF method in [1]. Furthermore, for $\Omega =1$, the auxiliary matrices ${\mathbf{Z}}_A$ and ${\mathbf{Z}}_B$ presented in Eqs. (21) and (22) become equal to the matrices ${\mathbf{Z}}_1$ and ${\mathbf{Z}}_2$, respectively, both of which are also defined for the conventional GPOF method in [1]:

$${\mathbf{Z}}_{\mathbf{A}}{|}_{\Omega =1}={\mathbf{Z}}_{\mathbf{1}}={\left(\begin{array}{ccccc}1& 1& \dots & 1& 1\\ z_1& z_2& \dots & z_{N-1}& z_N\\ z_1^2& z_2^2& \dots & z_{N-1}^2& z_N^2\\ ⋮& ⋮& & ⋮\\ z_1^{M-L-1}& z_2^{M-L-1}& \dots & z_{N-1}^{M-L-1}& z_N^{M-L-1}\end{array}\right)}_{(M-L)\times N},$$

(48)

$${\mathbf{Z}}_{\mathbf{B}}{|}_{\Omega =1}={\mathbf{Z}}_{\mathbf{2}}={\left(\begin{array}{ccccc}1& z_1& z_1^2& \dots & z_1^{L-1}\\ 1& z_2& z_2^2& \dots & z_2^{L-1}\\ ⋮& ⋮& ⋮& & ⋮\\ 1& z_{N-1}& z_{N-1}^2& \dots & z_{N-1}^{L-1}\\ 1& z_N& z_N^2& \dots & z_N^{L-1}\end{array}\right)}_{N\times L}.$$

(49)

Therefore, for $\Omega =1$, the eigenvalue problem defined by Eq. (34) is simplified to the conventional form:

$$({\mathbf{Y}}_1^{†}{\mathbf{Y}}_2-z_n\mathbf{I}){\mathbf{x}}_n=0,$$

(50)

as described in [1].

1.6. Analytical Estimation of the Variation in Complex Natural Frequencies

We now present a first-order approximation of the numerical error incurred when computing the complex natural frequency $s_n$ using the proposed method (Eq. (6)) compared to the conventional method (Eq. (5)). First, given that $z_n^{\Omega }=exp(\Omega t_s{s}_n)={\xi }_n$, we can establish:

$$\delta s_n(\Omega )={\displaystyle \frac{1}{t_s\Omega }}{\displaystyle \frac{\delta {\xi }_n}{{\xi }_n}},$$

(51)

where $\delta$ is the first-order differential operator, and $\delta s_n(\Omega )$ represents the error associated with the proposed algorithm in computing $s_n$ from ${\xi }_n$.

Notably, for the particular case of $\Omega =1$, we have:

$$\delta s_n\left(1\right)={\displaystyle \frac{1}{t_s}}{\displaystyle \frac{\delta z_n}{z_n}},$$

(52)

where $\delta s_n\left(1\right)$ indicates the perturbation in $s_n$ when it is computed using the conventional GPOF method, resulting from the propagation of error in ${\xi }_n$ [1].

A relationship between both errors can be established as follows. First, observe that the conventional algorithm with SVD decomposition estimates the values of $z_n$ numerically as the eigenvalues of the matrix ${\left({\mathbf{V}}_1^{†}\right)}^H{\mathbf{V}}_2^H$ (Eq. (40) for $\Omega =1$). Second, the proposed algorithm estimates $z_n^{\Omega }$ from the eigenvalues of the matrix ${\left({\mathbf{V}}_1^{†}\right)}^H{\mathbf{V}}_{\Omega +1}^H$ (Eq. (40) for $\Omega >1$) in a similar manner. Thus, the precision of the estimate of $z_n$ appears to be of the same order as that of $z_n^{\Omega }$. Consequently, the absolute errors for both estimates can be approximated as:

$$|\delta {\xi }_n|\approx |\delta z_n|.$$

This leads to the ratio of errors $\delta s_n(\Omega )$ and $\delta s_n\left(1\right)$ being approximated as:

$${\displaystyle \frac{|\delta s_n(\Omega )|}{|\delta s_n\left(1\right)|}}\approx {\displaystyle \frac{1}{\Omega }}{\displaystyle \frac{|z_n^{\Omega }|}{|z_n|}}.$$

(53)

Under the condition that $|z_n|\approx 1$, which can typically be achieved by selecting a sufficiently small sampling period such that $t_s\to 0$ implies $|z_n|\to 1$, we can define an Improvement Accuracy Ratio (IAR) as follows:

$$\mathrm{IAR}(\Omega )\approx {\displaystyle \frac{|\delta s_n(\Omega )|}{|\delta s_n\left(1\right)|}}\approx {\displaystyle \frac{1}{\Omega }}.$$

(54)

From Eq. (54), we observe that the IAR decreases as the polynomial matrix degree $\Omega$ increases, suggesting that the precision of the proposed method may surpass that of the conventional method, since $\Omega >1$.

2. Results

To investigate the performance of the proposed method, we designed and conducted four experiments. The first three experiments utilized single precision numeric format, while the fourth employed a double precision one.

For this study, we used an artificial test transient signal $y\left(t\right)$ defined as follows:

$$\begin{array}{c}y\left(t\right)=r_{1}exp(-{\sigma }_{1}t)cos\left({\omega }_{1}t\right)+r_{2}exp(-{\sigma }_{2}t)cos\left({\omega }_{2}t\right)\hfill \\ \hfill \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}+r_{3}exp(-{\sigma }_{3}t)cos\left({\omega }_{3}t\right)+r_{4}exp(-{\sigma }_{4}t)cos\left({\omega }_{4}t\right)+n\left(t\right),\end{array}$$

(55)

where $n\left(t\right)$ represents additive Gaussian white noise. The complex natural frequencies to be estimated are defined as $s_1=-{\sigma }_1+j{\omega }_1$, $s_2=-{\sigma }_2+j{\omega }_2$, $s_3=-{\sigma }_3+j{\omega }_3$, $s_4=-{\sigma }_4+j{\omega }_4$, $s_5=-{\sigma }_1-j{\omega }_1$, $s_6=-{\sigma }_2-j{\omega }_2$, $s_7=-{\sigma }_3-j{\omega }_3$, and $s_8=-{\sigma }_4-j{\omega }_4$.

The details of these experiments are presented in Sections 3.1, 3.2, 3.3, and 3.4 below.

2.1. A First Comparison of Accuracy from Proposed and Conventional Methods

In this first case study, the parameters of the signal $y\left(t\right)$ were defined as follows: $r_1=10$, ${\sigma }_1=1.010$, ${\omega }_1=2\pi ·1.251$; $r_2=7$, ${\sigma }_2=1.510$, ${\omega }_2=2\pi ·2.561$; $r_3=3$, ${\sigma }_3=2.010$, ${\omega }_3=2\pi ·3.901$; and $r_4=1$, ${\sigma }_4=3.010$, ${\omega }_4=2\pi ·6.112$. The number of elements N on the diagonal of the matrix ${\mathbf{Z}}_0$, the pencil parameter L, the pencil order $\Omega$, and the truncation factor T were set to $N=20$, $L=42$, $\Omega =6$, and $T=8$, respectively.

The signal $y\left(t\right)$ was sampled from $t_0=0$ to $t_f=t_{s}M$ with a sampling frequency of $f_s=40$, corresponding to a sampling time of $t_s={\displaystyle \frac{1}{f_s}}={\displaystyle \frac{1}{40}}$ and $M=100$ samples.

Additionally, $y\left(t\right)$ was contaminated with noise to achieve a signal-to-noise ratio (SNR) of 150 dB. The errors $\Delta s_1,\Delta s_2,\dots ,\Delta s_8$ were defined as follows:

$$\Delta s_n={\widehat{s}}_n-s_n=-({\widehat{\sigma }}_n-{\sigma }_n)+j({\widehat{\omega }}_n-{\omega }_n)=-\Delta {\sigma }_n+j\Delta {\omega }_n,$$

(56)

where ${\widehat{s}}_n$, ${\widehat{\sigma }}_n$, and ${\widehat{\omega }}_n$ represent the estimated values of $s_n$, ${\sigma }_n$, and ${\omega }_n$, respectively.

Finally, the errors were computed for both algorithms. Figure 2 illustrates the resulting errors from the proposed method (for $\Omega =6$) and the conventional GPOF method, plotted in the complex plane for an SNR of 150 dB.

2.2. Average Accuracy from the Proposed and Conventional Methods in Presence of Different Levels of Noise

In this case study, the parameters of the signal $y\left(t\right)$ were defined as follows: $r_1=10$, ${\sigma }_1=1.010$, ${\omega }_1=2\pi ·1.251$; $r_2=7$, ${\sigma }_2=1.510$, ${\omega }_2=2\pi ·2.561$; $r_3=3$, ${\sigma }_3=2.010$, ${\omega }_3=2\pi ·3.901$; and $r_4=1$, ${\sigma }_4=3.010$, ${\omega }_4=2\pi ·6.112$. The number of elements N on the diagonal of the matrix ${\mathbf{Z}}_0$, the pencil parameter L, the pencil order $\Omega$, and the truncation factor T were set to $N=20$, $L=42$, $\Omega =20$, and $T=8$, respectively.

The signal $y\left(t\right)$ was sampled from $t_0=0$ to $t_f=t_{s}M$ with a sampling frequency of $f_s=57$, resulting in a sampling time of $t_s={\displaystyle \frac{1}{f_s}}={\displaystyle \frac{1}{57}}$ and $M=100$ samples.

Both methods were then applied to $K=15$ realizations of $y\left(t\right)$ in the presence of noise, over 24 different values of SNR, ranging from 30 dB to 170 dB. From the estimated values ${\widehat{s}}_1$ (i.e., ${\widehat{\sigma }}_1$ and ${\widehat{\omega }}_1$), the root mean square errors $\mathrm{MSE}\left({\widehat{\sigma }}_1\right)$ and $\mathrm{MSE}\left({\widehat{\omega }}_1\right)$ were calculated for all SNRs as follows:

$$\mathrm{MSE}\left({\widehat{\sigma }}_1\right)={\displaystyle \frac{{\sum }_{k=1}^K{({\widehat{\sigma }}_{1,k}-{\sigma }_1)}^2}{K}},$$

(57)

$$\mathrm{MSE}\left({\widehat{\omega }}_1\right)={\displaystyle \frac{{\sum }_{k=1}^K{({\widehat{\omega }}_{1,k}-{\omega }_1)}^2}{K}},$$

(58)

where ${\widehat{\sigma }}_{1,k}$ and ${\widehat{\omega }}_{1,k}$ denote the estimated values of ${\sigma }_1$ and ${\omega }_1$ for any algorithm in the $k^{th}$ realization.

The obtained results are presented in Figure 3 and Figure 4. In these figures, the inverses of the normalized mean square errors ${\left({\displaystyle \frac{\mathrm{MSE}\left({\widehat{\sigma }}_1\right)}{{\sigma }_1^2}}\right)}^{-1}$ and ${\left({\displaystyle \frac{\mathrm{MSE}\left({\widehat{\omega }}_1\right)}{{\omega }_1^2}}\right)}^{-1}$ are plotted as functions of the signal-to-noise ratio.

2.3. Accuracy from the Proposed Method as a Function of Different Pencil Degrees

In this case study, the parameters of the signal $y\left(t\right)$ were set as follows: $r_1=10$, ${\sigma }_1=0.5$, ${\omega }_1=2\pi ·1.251$, and $r_2=r_3=r_4=0$. The number of elements N on the diagonal of the matrix ${\mathbf{Z}}_0$, the pencil parameter L, and the truncation factor T were set to $N=10$, $L=42$, and $T=2$, respectively. The signal $y\left(t\right)$ was sampled from $t_0=0$ to $t_f=t_{s}M$ with a sampling frequency of $f_s=25$, corresponding to a sampling time of $t_s={\displaystyle \frac{1}{f_s}}={\displaystyle \frac{1}{25}}$ and $M=100$ samples.

Both methods were then applied to $K=500$ realizations of $y\left(t\right)$ in the presence of noise with a signal-to-noise ratio (SNR) of 125 dB, for 20 distinct values of $\Omega$, in the range comprised from $\Omega =1$ to $\Omega =20$.

An average error in the estimates of ${\widehat{s}}_1$ from both methods was defined as:

$$\overline{|\Delta s_1(\Omega )|}={\displaystyle \frac{1}{K}}\sum _{k=1}^K|{\widehat{s}}_{1,k}(\Omega )-s_1|,$$

(59)

where ${\widehat{s}}_{1,k}(\Omega )$ denotes the estimated value for the $k^{th}$ realization. Subsequently, this mean error was computed using Eq. (59). The mean error values were later used to compute the estimated improvement accuracy ratio (IAR) as follows:

$$\widehat{IAR}={\displaystyle \frac{\overline{|\Delta s_1(\Omega )|}}{\overline{|\Delta s_1\left(1\right)|}}}.$$

(60)

Figure 5 illustrates the inverse of the numerically estimated value of $\widehat{IAR}$. For comparison, it is shown the inverse of the analytical IAR value, which is given by ${\displaystyle \frac{1}{\Omega }}$ according to Eq. (54), as well as the inverse of the IAR value estimated from Eq. (60).

2.4. A Comparison Between the Accuracy of the Proposed Method and the Conventional One When the Values of Complex Natural Frequencies Are Close to Each Other

Finally, we examine the performance of both the proposed and conventional algorithms applied to a signal $y\left(t\right)$ where the complex natural frequencies $s_n$ are very close to each other.

In this case study, the signal $y\left(t\right)$ consists of three damped harmonic functions characterized by the following parameters: $r_1=10$, $r_2=9$, $r_3=8$, $r_4=0$, ${\sigma }_1=2.002$, ${\sigma }_2=2.004$, ${\sigma }_3=2.003$, ${\omega }_1=1.0320$, ${\omega }_2=1.0325$, and ${\omega }_3=1.0332$. No Gaussian noise was added to $y\left(t\right)$ in this instance, meaning only numerical noise is present. The number of elements N on the diagonal of the matrix ${\mathbf{Z}}_0$, the pencil parameter L, the pencil degree $\Omega$, and the truncation factor T were set to $N=20$, $L=\left[{\displaystyle \frac{5M}{12}}\right]=42$, $\Omega =11$, and $T=6$, respectively. The signal $y\left(t\right)$ was sampled from $t_0=0$ to $t_f=t_{s}M$ with a sampling frequency of $f_s=1$, corresponding to a sampling time of $t_s=1$ and $M=100$ samples.

The estimations of ${\widehat{s}}_1$, ${\widehat{s}}_2$, and ${\widehat{s}}_3$ from both methods for one realization are plotted in Figure 6.

3. Discussion

Based on the results presented above, the proposed method accurately calculates the complex natural frequencies of the system $s_n$ from a given transient response $y\left(t\right)$. Furthermore, it demonstrates greater accuracy than the conventional linear method, particularly in the presence of low noise levels. Figure 2, Figure 3 and Figure 4 illustrate that the proposed method yields lower estimation errors for $s_n$ compared to the conventional method when the signal-to-noise ratio (SNR) exceeds 100 dB, especially when using single-precision numeric format. However, for SNRs below 100 dB, the precision of the proposed method becomes comparable to that of the conventional method (see Figure 3 and Figure 4).

The improvement in precision of the proposed method under low noise conditions becomes increasingly evident as $\Omega$ rises, as shown in Figure 5. The improvement factor $\mathrm{IAR}(\Omega )$ indicates that the estimation error for ${\widehat{s}}_1$ decreases if $\Omega$ increases. A good agreement between numerically and analytically estimated IAR values is observed in Figure 5.

Another notable advantage of the proposed method is its performance when analyzing transient signals with closely spaced complex natural frequencies. When employing double precision and considering only numerical noise, the estimated values obtained using the proposed method are significantly more precise than those derived from the conventional matrix pencil method, as illustrated in Figure 6.

4. Conclusion

In this work, we propose a nonlinear matrix pencil method with a pencil degree $\Omega$, such as $\Omega \ge 1$, in order to estimate the natural complex frequencies of resonance in a given linear system. The proposed method has been thoroughly investigated and compared with the conventional GPOF method.

An analytical approximation for the estimation error of these complex frequencies was developed and validated through simulations. The results demonstrate that the proposed method not only accurately estimates the complex natural frequencies of the system but also achieves higher precision under low noise condition compared to the GPOF method. Specifically, when using single precision, the improvement in accuracy is observed for signal-to-noise ratios (SNR) greater than 100 dB, and this enhancement becomes more pronounced as the nonlinear pencil degree increases.

Moreover, the findings indicate that the proposed method is particularly effective in estimating complex natural frequencies when they are closely spaced under low noise condition.