*2.5. Dispersion Curve Estimation with Matrix Pencil Method*

The wavenumber–frequency relations can be calculated based on a set of wave propagation signals collected in a number of points spaced evenly along a measurement line. One of the popular methods is the use of the two-dimensional fast Fourier transform (2D-FFT [48,54,55]). However, the curves obtained are in the form of 2D maps which can be problematic when comparing different results. To solve this problem, a simple and robust algorithm called the Matrix Pencil (MP) technique [56–58] can be used. The method is based on the calculation of the one-dimensional fast Fourier transform (FFT). Its advantage is a good elimination of noise effects on the obtained results. Let us have the set of time-domain signals measured in the *m* points distributed along *z* axis with the interval Δ*z*. The Fourier coefficients *X*(*z*,ω) for each signal *x*(*z*,*t*) can be calculated with respect to the formula:

$$X(z,\omega) \, := \int\_{-\infty}^{\infty} x(z,t)e^{-i\omega t}dt\tag{3}$$

where ω is a circular frequency. Then, the wavenumbers are estimated using the forward/backward averaging technique. The sequence *xn* of *m* complex values at each fixed frequency ω<sup>0</sup> is assumed to be the sum of *p* complex exponentials representing wavenumbers, as follows:

$$\mathbf{x}\_{\mathbb{H}} = \mathbf{x}\_{\mathbb{H}}(\omega\_0) = X(z\_{\mathbb{H}}, \omega\_0) \ = \sum\_{j=1}^{p} a\_j e^{-i\mathbf{k}\_j z} \qquad \text{or} \quad = 1, 2, \dots, m. \tag{4}$$

In the above formula, *p* is the model order, so-called pencil parameter that should be chosen based on the unknown number of searched signal modes *q* and satisfy the limitation:

$$q \le p \le m - q.\tag{5}$$

The algorithm requires the construction of a parent Hankel matrix **H** (*m* − *p* by *p* + 1) from the sequence *xn* as presented below:

$$\mathbf{H} \cdot = \begin{bmatrix} \mathbf{x}\_1 & \mathbf{x}\_2 & \cdots & \mathbf{x}\_{p+1} \\ \mathbf{x}\_2 & \mathbf{x}\_3 & \cdots & \mathbf{x}\_{p+2} \\ \vdots & \vdots & \ddots & \vdots \\ \mathbf{x}\_{m-p} & \mathbf{x}\_{m-p+1} & \cdots & \mathbf{x}\_m \end{bmatrix} \tag{6}$$

Two submatrices of **H** need to be created, the **H**<sup>0</sup> formed by deleting the first column of **H** and **H**<sup>1</sup> formed by deleting the last column of **H**. The sets of forward λ*<sup>f</sup>* and backward λ*<sup>b</sup>* exponential estimates can be calculated by solving the two eigenvalue problems, respectively:

$$\left(\mathbf{H}\_0^+\mathbf{H}\_1-\lambda\_f\mathbf{I}\right)\mathbf{\mathcal{e}} = 0, \quad \left(\mathbf{H}\_1^+\mathbf{H}\_0-\lambda\_b\mathbf{I}\right)\mathbf{\mathcal{e}} = 0\tag{7}$$

where (.)<sup>+</sup> indicates the Moore-Penrose generalized inverse and **I** is the identity matrix. The sets of forward *kf* and backward *kb* complex wavenumbers are computed as follows:

$$k\_f = \frac{\ln \lambda\_f}{i \Lambda z}, \quad k\_b = -\frac{\ln \lambda\_b}{i \Lambda z}. \tag{8}$$

The final wavenumber *k* values can be calculated by the averaging technique presented in [56]. Firstly, the sets *kf* and *kb* are sorted and matched in pairs in the order of sorting. If the relative difference between two values of a certain pair is within a specified tolerance *l*, the pair is averaged arithmetically, if not, the pair is discarded.
