Mathematical basis

A transient electromagnetic signal can be represented as:^[3]

${\displaystyle y(t)=x(t)+n(t)\approx \sum _{i=1}^{M}R_{i}e^{s_{i}t}+n(t);0\leq t\leq T,}$

where

${\displaystyle y(t)}$ is the observed time-domain signal,

${\displaystyle n(t)}$ is the signal noise,

${\displaystyle x(t)}$ is the actual signal,

${\displaystyle R_{i}}$ are the residues (${\displaystyle R_{i}}$),

${\displaystyle s_{i}}$ are the poles of the system, defined as ${\displaystyle s_{i}=-\alpha _{i}+j\omega _{i}}$,

${\displaystyle z_{i}=e^{(-\alpha _{i}+j\omega _{i})T_{s}}}$ by the identities of Z-transform,

${\displaystyle \alpha _{i}}$ are the damping factors and

${\displaystyle \omega _{i}}$ are the angular frequencies.

The same sequence, sampled by a period of ${\displaystyle T_{s}}$, can be written as the following:

${\displaystyle y[kT_{s}]=x[kT_{s}]+n[kT_{s}]\approx \sum _{i=1}^{M}R_{i}z_{i}^{k}+n[kT_{s}];k=0,...,N-1;i=1,2,...,M}$,

Generalized pencil-of-function estimates the optimal ${\displaystyle M}$ and ${\displaystyle z_{i}}$'s.^[4]

Noise-free analysis

For the noiseless case, two ${\displaystyle (N-L)\times L}$ matrices, ${\displaystyle Y_{1}}$ and ${\displaystyle Y_{2}}$, are produced:^[3]

${\displaystyle [Y_{1}]={\begin{bmatrix}x(0)&x(1)&\cdots &x(L-1)\\x(1)&x(2)&\cdots &x(L)\\\vdots &\vdots &\ddots &\vdots \\x(N-L-1)&x(N-L)&\cdots &x(N-2)\end{bmatrix}}_{(N-L)\times L};}$ ${\displaystyle [Y_{2}]={\begin{bmatrix}x(1)&x(2)&\cdots &x(L)\\x(2)&x(3)&\cdots &x(L+1)\\\vdots &\vdots &\ddots &\vdots \\x(N-L)&x(N-L+1)&\cdots &x(N-1)\end{bmatrix}}_{(N-L)\times L}}$

where ${\displaystyle L}$ is defined as the pencil parameter. ${\displaystyle Y_{1}}$ and ${\displaystyle Y_{2}}$ can be decomposed into the following matrices:^[3]

${\displaystyle [Y_{1}]=[Z_{1}][B][Z_{2}]}$

${\displaystyle [Y_{2}]=[Z_{1}][B][Z_{0}][Z_{2}]}$

where

${\displaystyle [Z_{1}]={\begin{bmatrix}1&1&\cdots &1\\z_{1}&z_{2}&\cdots &z_{M}\\\vdots &\vdots &\ddots &\vdots \\z_{1}^{(N-L-1)}&z_{2}^{(N-L-1)}&\cdots &z_{M}^{(N-L-1)}\end{bmatrix}}_{(N-L)\times M};}$ ${\displaystyle [Z_{2}]={\begin{bmatrix}1&z_{1}&\cdots &z_{1}^{L-1}\\1&z_{2}&\cdots &z_{2}^{L-1}\\\vdots &\vdots &\ddots &\vdots \\1&z_{M}&\cdots &z_{M}^{L-1}\end{bmatrix}}_{M\times L}}$

${\textstyle [Z_{0}]}$ and ${\textstyle [B]}$ are ${\textstyle M\times M}$ diagonal matrices with sequentially-placed ${\textstyle z_{i}}$ and ${\textstyle R_{i}}$ values, respectively.^[3]

If ${\textstyle M\leq L\leq N-M}$, the generalized eigenvalues of the matrix pencil

${\displaystyle [Y_{2}]-\lambda [Y_{1}]=[Z_{1}][B]([Z_{0}]-\lambda [I])[Z_{2}]}$

yield the poles of the system, which are ${\displaystyle \lambda =z_{i}}$. Then, the generalized eigenvectors ${\displaystyle p_{i}}$ can be obtained by the following identities:^[3]

${\displaystyle [Y_{1}]^{+}[Y_{1}]p_{i}=p_{i};}$    ${\displaystyle i=1,...,M}$

${\displaystyle [Y_{1}]^{+}[Y_{2}]p_{i}=z_{i}p_{i};}$    ${\displaystyle i=1,...,M}$

where the ${\displaystyle ^{+}}$ denotes the Moore–Penrose inverse, also known as the pseudo-inverse. Singular value decomposition can be employed to compute the pseudo-inverse.

Noise filtering

If noise is present in the system, ${\textstyle [Y_{1}]}$ and ${\textstyle [Y_{2}]}$ are combined in a general data matrix, ${\textstyle [Y]}$:^[3]

${\displaystyle [Y]={\begin{bmatrix}y(0)&y(1)&\cdots &y(L)\\y(1)&y(2)&\cdots &y(L+1)\\\vdots &\vdots &\ddots &\vdots \\y(N-L-1)&y(N-L)&\cdots &y(N-1)\end{bmatrix}}_{(N-L)\times (L+1)}}$

where ${\displaystyle y}$ is the noisy data. For efficient filtering, L is chosen between ${\textstyle {\frac {N}{3}}}$ and ${\textstyle {\frac {N}{2}}}$. A singular value decomposition on ${\textstyle [Y]}$ yields:

${\displaystyle [Y]=[U][\Sigma ][V]^{H}}$

In this decomposition, ${\textstyle [U]}$ and ${\textstyle [V]}$ are unitary matrices with respective eigenvectors ${\textstyle [Y][Y]^{H}}$ and ${\textstyle [Y]^{H}[Y]}$ and ${\textstyle [\Sigma ]}$ is a diagonal matrix with singular values of ${\textstyle [Y]}$. Superscript ${\textstyle H}$ denotes the conjugate transpose.^[3][4]

Then the parameter ${\textstyle M}$ is chosen for filtering. Singular values after ${\textstyle M}$, which are below the filtering threshold, are set to zero; for an arbitrary singular value ${\textstyle \sigma _{c}}$, the threshold is denoted by the following formula:^[1]

${\displaystyle {\frac {\sigma _{c}}{\sigma _{max}}}=10^{-p}}$,

${\textstyle \sigma _{max}}$ and p are the maximum singular value and significant decimal digits, respectively. For a data with significant digits accurate up to p, singular values below ${\textstyle 10^{-p}}$ are considered noise.^[4]

${\textstyle [V_{1}']}$ and ${\textstyle [V_{2}']}$ are obtained through removing the last and first row and column of the filtered matrix ${\textstyle [V']}$, respectively; ${\textstyle M}$ columns of ${\textstyle [\Sigma ]}$ represent ${\textstyle [\Sigma ']}$. Filtered ${\textstyle [Y_{1}]}$ and ${\textstyle [Y_{2}]}$ matrices are obtained as:^[4]

${\displaystyle [Y_{1}]=[U][\Sigma '][V_{1}']^{H}}$

${\displaystyle [Y_{2}]=[U][\Sigma '][V_{2}']^{H}}$

Prefiltering can be used to combat noise and enhance signal-to-noise ratio (SNR).^[1] Band-pass matrix pencil (BPMP) method is a modification of the GPOF method via FIR or IIR band-pass filters.^[1][5]

GPOF can handle up to 25 dB SNR. For GPOF, as well as for BPMP, variace of the estimates approximately reaches Cramér–Rao bound.^[3][5][4]

Calculation of residues

Residues of the complex poles are obtained through the least squares problem:^[1]

${\displaystyle {\begin{bmatrix}y(0)\\y(1)\\\vdots \\y(N-1)\end{bmatrix}}={\begin{bmatrix}1&1&\cdots &1\\z_{1}&z_{2}&\cdots &z_{M}\\\vdots &\vdots &\ddots &\vdots \\z_{1}^{N-1}&z_{2}^{N-1}&\cdots &z_{M}^{N-1}\end{bmatrix}}{\begin{bmatrix}R_{1}\\R_{2}\\\vdots \\R_{M}\end{bmatrix}}}$