# Multi-output subspace identification of complex Bloch wavenumbers in 1D periodic structures

Alvaro Gavilán-Rojas, Qinghua Zhang, Olivier Robin, Christophe Droz

► To cite this version:

Alvaro Gavilán-Rojas, Qinghua Zhang, Olivier Robin, Christophe Droz. Multi-output subspace identification of complex Bloch wavenumbers in 1D periodic structures. Mechanical Systems and Signal Processing, 2026, 242, pp.113577. 10.1016/j.ymssp.2025.113577. hal-05365706

HAL Id: hal-05365706

https://hal.science/hal-05365706v1

Submitted on 14 Nov 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Multi-output subspace identification of complex Bloch wavenumbers in 1D periodic structures

Alvaro Gavilán Rojas<sup>a,b</sup>, Qinghua Zhang<sup>a</sup>, Olivier Robin<sup>b</sup>, Christophe Droz<sup>a</sup>

<sup>a</sup> Univ. Gustave Eiffel, Inria, COSYS-SII, 14S, Rennes, France

<sup>b</sup> Centre de Recherche Acoustique-Signal-Humain, Université de Sherbrooke, Sherbrooke, Canada

## Abstract

The experimental characterization of complex dispersion curves is challenging in phononic crystals, composites, periodic, architected or metamaterials. Recent studies have highlighted the importance of subspace identification methods in determining wave propagation properties through complex wavenumbers and, consequently, in characterizing a complex structure experimentally. Still, such methods have not yet been adapted for 1D periodic structures with periodic sampling limitations. This work introduces a Subspace-based complex Bloch WAveNumber identification method (SWAN) which can take advantage of full-field vibration measurements (i.e., multiple data points per unit cell) to statistically reduce the negative impact of having a limited number of unit cells. The SWAN method is based on a state-space representation of the wave finite element method. A symplectic state-space model is formulated and mathematically proved to represent the original system. Eventually, the proposed method enhances complex wavenumber estimates when a small number of unit cells is available. In addition, a general-purpose, adaptive spectral mask is introduced to reject physically irrelevant identification results, enabling straightforward denoising of the identified dispersion curves. The proposed approach is validated through numerical and experimental applications.

**Keywords:** Complex wavenumber, Subspace identification, Full-field measurements, WFEM, Phononic crystals, Symplectic state-space model

## Notation

<table><thead><tr><th>Symbol</th><th>Description</th></tr></thead><tbody><tr><td>κ</td><td>complex wavenumber</td></tr><tr><td>λ</td><td>propagation constant</td></tr><tr><td>Δ</td><td>unit cell length</td></tr><tr><td>ω</td><td>angular frequency</td></tr><tr><td>N</td><td>number of unit cells of the waveguide</td></tr><tr><td>N<sub>o</sub></td><td>number of unit cells considered for identification</td></tr><tr><td>A</td><td>state-transition matrix</td></tr><tr><td>n<sub>w</sub></td><td>model order – number of waves to identify</td></tr><tr><td>x<sub>k</sub></td><td>state vector</td></tr><tr><td>C</td><td>output matrix</td></tr><tr><td>y<sub>k</sub></td><td>output vector</td></tr><tr><td>l</td><td>number of outputs (sensors)</td></tr><tr><td>H</td><td>block-Hankel measurements matrix of size (pl) × s</td></tr><tr><td>p</td><td>spatial horizon</td></tr><tr><td>s</td><td>number of columns of H</td></tr><tr><td>O</td><td>extended observability matrix</td></tr><tr><td>Y</td><td>output vectors matrix of size l × N<sub>o</sub></td></tr></tbody></table>
# 1. Introduction

Periodic structures exhibit bandgap regions, i.e., frequency bands with strongly attenuated elastic waves. This phenomenon has been utilized mainly for Noise, Vibration and Harshness (NVH) mitigation – e.g., through phononic crystals and local resonance [1, 2, 3, 4]. Beyond mere acoustic comfort, the reduction of NVH is a global health concern [5].

Dispersion curves provide the essential information to characterize bandgap regions and the functional properties of periodic structures. In the case of 1D waveguides – which are addressed in this work –, they represent the evolution of wavenumbers $κ \in C$ through the frequency domain. Bragg bandgap regions can be identified at the conceptual design stage, and the degree of effectiveness in attenuating waves and vibrations can be evaluated through the imaginary part of the wavenumber and vibration response computations, respectively. Next, in the experimental validation stage, vibration measurements on a real prototype can confirm vibration attenuation in bandgap regions [6, 7, 8], and identification methods can also be applied to vibration measurements to estimate the complex dispersion curves [9, 10]. This can also verify the match between the mechanical model from the design stage and the real prototype from the experimental validation stage. Moreover, dispersion curve characterization is often seen as an intermediate step towards structural parameter identification [11, 12], which is also relevant to identify effective properties (e.g., Young's modulus, structural loss factor) of complex periodic media or their underlying partial differential equations [13].

Several methods in the literature address the problem of identifying wavenumbers. A very first one to mention is the spatial Fourier transform [14]: by analyzing time-harmonic vibration fields at a fixed frequency, underlying space-harmonic information (spatial wavelengths $Λ \in R$) is retrieved in the form of real wavenumbers $κ = 2π/Λ$. Despite the ease of implementation, two major drawbacks remain: the imaginary part of the wavenumber is not identified – and therefore, wave attenuation cannot be evaluated – and the spatial resolution of the measurements limits the wavenumber resolution. However, the spatial Fourier transform is still advantageous for frequency-wavenumber domain analysis of guided wavefields, not only to characterize guided waves but also for damage characterization [15, 16]. Prony's method can be used to identify complex wavenumbers. Like the Fourier transform, Prony's method searches for damped sinusoids in a signal. This method has been applied to identify complex wavenumbers in homogeneous structures [17, 18]. Recently, Ribeiro et al. [10] proposed a method for identifying complex Bloch wavenumbers in periodic structures using Prony's method and extended it to multiple measurements per unit cell: this is the condition referred to as *multi-output* in this work, which is linked to the spatial richness of data in periodic waveguides. Their work elucidates the need for spatially periodic measurements with the same periodicity of the unit cell. Another relevant method is the *INverse CONvolution Method* (INCOME) proposed by Boukadia et al. [19], where complex wavenumber identification in periodic structures can be achieved by taking into account a symplectic property of the waves propagating in the structure, which improves the identification results by finding roots of a polynomial that represent both right- and left-going waves. It is worth noting that INCOME can accommodate cases with multiple displacements per unit cell (see (53) in [19]) and has been applied to create artificial unit cells in homogeneous structures to enhance estimation results. Some other methods that can be mentioned are the *Inhomogeneous Wave Correlation* (IWC), developed by Berthaut and co-authors [20, 21, 9], which is essentially for mono-modal identification: a single physical mode is identified – e.g., bending mode. It has been adapted for multi-modal identification with a subtraction procedure in [22]: the *multi-modal* term is used in this work to refer to modal richness of identification results, when several physical modes are identified – e.g., bending, torsion or high-order modes. The IWC was applied for bandgap characterization in 1D phononic crystals in [23], and an extended version that takes into account the source location was proposed in [24]. The McDaniel method [25] was also proposed for multi-mode wavenumber identification. However, these two methods rely on nonlinear optimization and assume that spatial harmonic functions describe the vibration field – instead of spatial Bloch-harmonic functions required for periodic media. Another wavenumber identification method for periodic structures, utilizing the *Wave Superposition Method*, was developed by Junyi et al. [26, 27]. This approach enables bandgap identification with a relatively small number of unit cells and sensors, and relies on single-output identification. Time-domain subspace identification was recently proposed in [28]; however, the method has not yet been extended to multi-output identification, nor have multi-modal results been reported. Another recent contribution to note is the family of algebraic wavenumber identification methods: Li et al. proposed an *Algebraic Wavenumber Identification* (AWI) [29], which can be applied to periodic structures. In a recent version (*Generalized algebraic K-space identification* (GAKSI) [30]), they extended the approach to multi-output measurements. Despite being a promising method, the authors of these papers have reported only mono-modal identification cases [29, 30]. After identifying the inconveniences for successfully implementing the method for multi-modal identification, this method is not considered for comparison in the present paper.

P. Margerit and co-authors made remarkable contributions using subspace methods for the wavenumber identification problem [31, 32, 33]. Based on the subspace method *Estimation of Signal Parameters via Rotational*
<table><caption>Table 1: List of the reviewed wavenumber identification methods and their capabilities to perform multi-output – which in the current study means to have several measurements inside the unit cell – and multi-modal identification. The orange ✓(n.r.) sign means that the methodology is theoretically capable of handling multi-modal identification, but no numerical or experimental validations were reported, to the authors’ knowledge.</caption><thead><tr><th>Wavenumber identification method</th><th>Complex</th><th>Multi-output</th><th>Multi-modal</th></tr></thead><tbody><tr><td>Spatial Fourier [14]</td><td>✗</td><td>✗</td><td>✓</td></tr><tr><td>Time-domain subspace identification [28]</td><td>✓</td><td>✗</td><td>✗</td></tr><tr><td>Wave superposition method [26]</td><td>✓</td><td>✗</td><td>✓</td></tr><tr><td>ESPRIT / HRWA [31]</td><td>✓</td><td>✗</td><td>✓</td></tr><tr><td>Inhomogeneous Wave Correlation [20]</td><td>✓</td><td>✗</td><td>✓</td></tr><tr><td>Mc Daniel [25]</td><td>✓</td><td>✗</td><td>✓</td></tr><tr><td>Algebraic Wavenumber Identification [30]</td><td>✓</td><td>✓</td><td>✓(n.r.)</td></tr><tr><td>Prony's method [10]</td><td>✓</td><td>✓</td><td>✓(n.r.)</td></tr><tr><td>Inverse Convolution Method [19]</td><td>✓</td><td>✓</td><td>✓</td></tr></tbody></table>

*Invariance Techniques* (ESPRIT) [34], initially developed in the signal processing community for direction of arrival estimation, multi-modal identification is successfully achieved in 1D [31] and 2D [32] structures. Margerit’s work appears to be inspired by the application of ESPRIT to modal analysis, as proposed by Ege et al. [35]. These two applications – wave propagation and structural vibrations – are closely related within the identification problem. In structural vibrations and modal analysis, the *modal parameters* are (i) natural frequencies, (ii) damping ratios, and (iii) normal modes. The equivalent to these modal parameters for elastic wave propagation are (i) real part of the wavenumber, (ii) imaginary part of the wavenumber, and (iii) wave modes (or Bloch modes, for periodic waveguides). However, from the three modal parameters, only the frequencies and damping ratios are usually identified with ESPRIT, and even the reconstruction of normal modes relies on the individual estimation of amplitudes and phases – a signal reconstruction – at each measurement point [35]. This is because the underlying assumption in the ESPRIT method is an invariance property of the sensor array geometry, and this invariance is considered just once (both in 1D [31] and 2D [32]). In simpler terms, if several sets of measurements exhibit the same invariance property – e.g., several measurement points inside a unit cell that are periodically repeated – there is no other choice than to identify properties from each set individually. This means that, for 1D periodic waveguides, the ESPRIT method [31] follows a single-output, line-wise measurement formalism (leading to Hankel data matrices), as most of the methods described in the previous paragraph. A summary of all the described methods for wavenumber identification is provided in Table 1. To conclude this review, the only methods capable of addressing the problem of multi-output and multi-modal identification are INCOME, Prony’s method, and GAKSI. Consequently, INCOME and Prony’s method will be compared to the methodology proposed in this paper.

To enhance the estimation of signal parameters and overcome the limitation described above, multiple invariances (e.g., for uniform linear sensor arrays and uniformly sampled time series) can be collectively included in a single identification procedure, and this has been proposed by Swindlehurst et al. as a *Multiple Invariance* ESPRIT (MI-ESPRIT) [36]. This method fully exploits the invariance structure shared by the abovementioned measurements. To the authors’ knowledge, Swindlehurst [36] is also the first one to point out the links between subspace methods originated from signal processing (e.g., Pisarenko [37], MUSIC [38], ESPRIT [39], Matrix Pencil [40]) and those originated in the system identification community [41]. In the latter, exploitation of the entire invariance structure of *state-space models* is considered through an *extended observability matrix*, and the identification problem is easily solved by least-squares minimization, as in the original ESPRIT version. Further details on these links are found in [42]. Using subspace methods that rely on state-space models intrinsically considers multiple measurements per unit cell (leading to block-Hankel data matrices). In that sense, overcoming the limitation of line-wise measurement formalism will be relevant to improving the wavenumber identification problem addressed in this work. Although the reasoning behind this change from single-output to multi-output measurement leads to a change from Hankel to block-Hankel data matrices – as it will be shown in the paper –, this trivial consequence implies a subspace identification framework change. As a final remark, it is worth mentioning that the separation between signal and noise subspaces often relies on a *model order parameter*. Both Ege and Margerit propose using a model order selection criterion based on an ESTimation of ERror (ESTER) criterion [43], which is rooted in the statistics of the collected data. In the system identification community, stabilization diagrams are more commonly used, where identification results are examined against some model order evolution [44]. However, most of the reviewed methods for wavenumber identification fix this model order parameter – usually denoted by $n_w$, the number of waves to estimate – supposing its *a priori*
knowledge, which is in practice the case for most wavenumber identification applications.

In this paper, a data-driven subspace method is proposed for the identification of complex wavenumbers in periodic structures. This method takes advantage of multiple measurement points within the periodic unit cells, which are not necessarily arranged on a regular grid, as is often the case with complex unit cell geometries. Multiple measurements can be obtained from a full-field measurement technique (e.g., laser Doppler vibrometry [45], optical deflectometry [46]) to compensate relatively few unit cells in most experiments [27, 10] to characterize periodic waveguides. After deriving the state-space model of the 1D wave propagation problem in a periodic waveguide, two different approaches for subspace-based identification are proposed: (i) a classical identification scheme within subspace methods for system identification, and (ii) a new symplectic identification scheme, that takes into account the simultaneous existence of both right- and left-going waves. Additionally, after the multi-modal identification of the $n_w$ waves at each frequency – the $n_w$ eigenvalues conforming the spectrum of an estimated matrix, which are the propagation constants $\lambda = \exp(i\kappa\Delta) \in \mathbb{C}$, where $\Delta \in \mathbb{R}$ is the unit cell length –, a frequency-dependent spectral mask acting on the propagation constants is proposed to retain in the identified dispersion curves only the physically meaningful part of the spectrum – e.g., eigenvalues located near the unit circle or Brillouin edges. The two subspace methods and the spectral mask presented in this work have potential applications beyond 1D periodic waveguides, since 1D homogeneous waveguides can be comprised of artificial unit cells, which are built just for the sake of estimation accuracy improvement. This is demonstrated in the present work.

The remainder of this paper is organized as follows: Section 2 outlines the theoretical framework of the proposed subspace-based complex wavenumber identification method. Section 3 provides the frequency-dependent spectral mask, a practical strategy for a posteriori wavenumber selection based on the identified propagation constants. Section 4 presents a numerical validation of the proposed identification method. Section 5 provides two experimental applications that demonstrate the effectiveness of the methods under various scenarios: single- and multi-output measurements, and artificial unit cells. Conclusions are presented in Section 6.

## 2. Subspace-based wavenumber identification

### 2.1. 1D-periodic waveguide system model in state-space form

To present the approach of this paper for wavenumber identification from sensor measurements, we express in this section the wave propagation equations and their connection to measurements in a state-space form. In practice, for wavenumber identification, a detailed knowledge of these equations is not required – only their existence is assumed. The model of a periodic waveguide is based on the wave finite element method (WFEM) [47, 48], which expresses the dynamic motion equations of the periodic waveguide in terms of its smallest repetitive element, the $\Delta$-periodic unit cell (UC). After finite element discretization, the time-harmonic motion equation at some angular frequency $\omega \in \mathbb{R}$ for a periodic unit cell is given by:

$$ \mathbf{D}\mathbf{u}_k = \mathbf{f}_k \quad \forall k \in [1, N] \quad \text{with} \quad \mathbf{D} = \mathbf{D}(\omega) = -\omega^2\mathbf{M} + i\omega\mathbf{C} + \mathbf{K}, \quad (1) $$

where $N$ is the number of unit cells in the waveguide, $\mathbf{D}(\omega) \in \mathbb{C}^{n_d \times n_d}$ represents the frequency-dependent dynamic stiffness matrix, formed with the classical finite element mass, damping and stiffness matrices $\mathbf{M}$, $\mathbf{C}$ and $\mathbf{K}$, all being real symmetric and of size $n_d$, which is equal to the number of degrees of freedom of the unit cell. The nodal displacements and forces are $\mathbf{u}_k \in \mathbb{C}^{n_d}$ and $\mathbf{f}_k \in \mathbb{C}^{n_d}$, respectively. The time dependence of these displacements and force vectors have been simplified in (1), which are $\mathbf{u}_k(t) = \mathbf{u}_k \exp(i\omega t)$ and $\mathbf{f}_k(t) = \mathbf{f}_k \exp(i\omega t)$.

Let the set of indices $\{1, 2, \dots, n_d\}$ be partitioned into three subsets, namely, $L \in \mathbb{N}^n$ for the left boundary of the unit cell, $R \in \mathbb{N}^n$ for the right boundary and $I \in \mathbb{N}^m$ for the inner degrees of freedom, where $n$ is the number of interface degrees of freedom and $m$ is the number of inner degrees of freedom. The corresponding unit cell vectors and matrices are partitioned in the same manner. This partition is depicted in Figure 1(a).

<figure><img src="image_2.png" /><figcaption>Figure 1: (a) Unit cell of size $\Delta$ and its partitioned degrees of freedom: $L$ (left), $I$ (inner), and $R$ (right). (b) Periodic waveguide formed by an assembly of $N$ unit cells in the waveguide direction.</figcaption></figure>
To form a periodic waveguide, as depicted in Figure 1(b), $N > 2$ unit cells are assembled in one direction, so that the discretization grids of the k-th and (k+1)-th unit cells are coherent at their interface. The displacement continuity can be written as $\mathbf{u}_k^R = \mathbf{u}_{k+1}^L$ and the dynamic force equilibrium as $\mathbf{f}_k^R + \mathbf{f}_{k+1}^L = \mathbf{0}$. The equation of harmonic motion (1) is then developed with partitioned degrees of freedom as:

$$ \begin{bmatrix} \mathbf{D}_{LL} & \mathbf{D}_{LI} & \mathbf{D}_{LR} \\ \mathbf{D}_{IL} & \mathbf{D}_{II} & \mathbf{D}_{IR} \\ \mathbf{D}_{RL} & \mathbf{D}_{RI} & \mathbf{D}_{RR} \end{bmatrix} \begin{Bmatrix} \mathbf{u}_k^L \\ \mathbf{u}_k^I \\ \mathbf{u}_{k+1}^L \end{Bmatrix} = \begin{Bmatrix} \mathbf{f}_k^L \\ \mathbf{0} \\ -\mathbf{f}_{k+1}^L \end{Bmatrix} \quad \forall k \in [1, N]. \quad (2) $$

This equation involving both inner and boundary nodes of each unit cell is then reduced by *dynamic condensation* [48], based on the fact that no external forces act in the inner region of the unit cell ($\mathbf{f}_k^I = \mathbf{0}$), yielding an expression of $\mathbf{u}_k^I$ in terms of $\mathbf{u}_k^L$ and $\mathbf{u}_{k+1}^L$:

$$ \mathbf{u}_k^I = -\mathbf{D}_{II}^{-1} (\mathbf{D}_{IL}\mathbf{u}_k^L + \mathbf{D}_{IR}\mathbf{u}_{k+1}^L) \quad \forall k \in [1, N]. \quad (3) $$

The inner nodes are then eliminated by substituting (3) into (2), yielding the *condensed equation of motion*:

$$ \begin{Bmatrix} \mathbf{f}_k^L \\ -\mathbf{f}_{k+1}^L \end{Bmatrix} = \begin{bmatrix} \tilde{\mathbf{D}}_{LL} & \tilde{\mathbf{D}}_{LR} \\ \tilde{\mathbf{D}}_{RL} & \tilde{\mathbf{D}}_{RR} \end{bmatrix} \begin{Bmatrix} \mathbf{u}_k^L \\ \mathbf{u}_{k+1}^L \end{Bmatrix} \quad \text{where} \quad \tilde{\mathbf{D}}_{mn} = \mathbf{D}_{mn} - \mathbf{D}_{mI}\mathbf{D}_{II}^{-1}\mathbf{D}_{In} \quad \forall m, n \in \{L, R\}. \quad (4) $$

Finally, the reduced equation (4) is rearranged as a *state-space model* [49], which is essential for the subspace identification approach presented in the next subsections. It is composed of a state equation and an output equation, as follows:

$$ \begin{cases} \mathbf{x}_{k+1} = \mathbf{A}\mathbf{x}_k \\ \mathbf{y}_k = \mathbf{C}\mathbf{x}_k + \mathbf{n}_k \end{cases} \quad \text{with} \quad \mathbf{x}_k = \begin{Bmatrix} \mathbf{u}_k^L \\ \mathbf{f}_k^L \end{Bmatrix} \quad \forall k \in [1, N]. \quad (5) $$

where the state vector $\mathbf{x}_k \in \mathbb{C}^{2n}$ consists of the nodal displacements $\mathbf{u}_k^L$ and forces $\mathbf{f}_k^L$, the *output vector* $\mathbf{y}_k \in \mathbb{C}^l$ contains the $l$ $\Delta$-periodic measurements inside the k-th unit cell and on its left boundary, and the measurement noise vector $\mathbf{n}_k \in \mathbb{C}^l$ is assumed to be a zero-mean stationary white noise sequence. The involved measurements are either harmonic displacements or harmonic velocities. While displacement measurements are directly included in the measurement vector $\mathbf{y}_k$, velocities are converted to displacements through the relationship $\mathbf{u}_k = \dot{\mathbf{u}}_k/(i\omega)$ before being included. In the following, $\mathbf{y}_k$ will refer to displacement measurements, including converted velocity measurements. The transfer matrix $\mathbf{A} \in \mathbb{C}^{2n \times 2n}$ [48] resulting from the rearrangement of (4) is

$$ \mathbf{A} = \begin{bmatrix} -\tilde{\mathbf{D}}_{LR}^{-1}\tilde{\mathbf{D}}_{LL} & \tilde{\mathbf{D}}_{LR}^{-1} \\ (\tilde{\mathbf{D}}_{RR}\tilde{\mathbf{D}}_{LR}^{-1}\tilde{\mathbf{D}}_{LL} - \tilde{\mathbf{D}}_{RL}) & -\tilde{\mathbf{D}}_{RR}\tilde{\mathbf{D}}_{LR}^{-1} \end{bmatrix}, \quad (6) $$

and the matrix $\mathbf{C} \in \mathbb{C}^{l \times 2n}$ relating the output vector to the state vector will be discussed later, together with the configuration of the sensor measurements contained in $\mathbf{y}_k$.

The eigenvalues of $\mathbf{A}$ are the propagation constants $\lambda \in \mathbb{C}$, which are related to the wavenumbers of interest in the considered context. These are the eigenvalues characterizing the system that are aimed to identify through subspace methods. Note that $\mathbf{A}$ has a symplectic structure. This fundamental property was previously used by the authors to improve the accuracy and computational efficiency for the finite element resolution of $\lambda$ in the WFEM [50]. This property will be used again in Section 2.3 for subspace identification.

The output vector $\mathbf{y}_k$, in the simplest case, contains harmonic displacement measurements located on the left boundary of the k-th unit cell, corresponding to part of the components of $\mathbf{u}_k^L$. In this simple case, since $\mathbf{u}_k^L$ is a sub-vector of the state vector $\mathbf{x}_k$ defined in (5), the matrix $\mathbf{C}$ relating $\mathbf{y}_k$ to $\mathbf{x}_k$ contains only 0's and 1's. This simple case is too restrictive. In practice, it is essential to make measurements *inside unit cells* to enhance the information contained in sensor data. Then the output vector $\mathbf{y}_k$ is no longer a sub-vector of the condensed state vector $\mathbf{x}_k$ defined in (5). This connection between $\mathbf{y}_k$ and $\mathbf{x}_k$ is the key step for using multiple-point measurements per unit cell in wavenumber identification. It is detailed in Appendix A and summarized below.

In Appendix A, three different strategies to include inner measurements in the output vector $\mathbf{y}_k$ are presented. Firstly, two *augmented state* solutions are described in Appendix A.1 and Appendix A.2: to relate $\mathbf{y}_k$ to $\mathbf{x}_k$, an obvious way is to complete the state vector $\mathbf{x}_k$ to include the inner nodal displacements, so that $\mathbf{y}_k$ becomes a sub-vector of $\mathbf{x}_k$. These two strategies share the same drawback of a large state vector and the associated large transfer matrix $\mathbf{A}$. A third and better solution, avoiding augmented states and large transfer matrices, is proposed below and detailed in Appendix A.3.
It is in fact possible to keep the condensed state vector $\mathbf{x}_k$ as defined in (5), yet to relate $\mathbf{y}_k$ containing both inner and boundary measurements to $\mathbf{x}_k$, so that

$$ \mathbf{y}_k = \mathbf{C}\mathbf{x}_k + \mathbf{n}_k. \qquad (7) $$

This is the same output equation as in (5), but now with a matrix $\mathbf{C}$ filled with complex numbers, instead of the 0's and 1's in the previous solutions. The key idea is to express the inner nodal displacements $\mathbf{u}_k^I$ as a linear function of $\mathbf{x}_k$. In (3) $\mathbf{u}_k^I$ is already expressed in terms of $\mathbf{u}_k^L$ and $\mathbf{u}_{k+1}^L$. The former is part of $\mathbf{x}_k$ as defined in (5), but not the latter. However, $\mathbf{u}_{k+1}^L$ is part of $\mathbf{x}_{k+1}$, which is equal to $\mathbf{A}\mathbf{x}_k$ according to the state equation in (5).

As a final note about $\mathbf{y}_k$ and $\mathbf{C}$, in each unit cell, it is not necessary to choose regularly spaced measurement points -i.e., periodic measurements inside each unit cells are not required, but the measurement positions must be the same in different unit cells, so that the matrix $\mathbf{C}$ relating $\mathbf{y}_k$ to $\mathbf{x}_k$ is independent of $k$. In this sense, the measurements contained in $\mathbf{y}_k$ are $\Delta$-periodic, with $\Delta$ being the unit cell length of the waveguide. This measurement scheme is depicted in Figure 2.

<figure><img src="image_3.png" /><figcaption>Figure 2: Periodic measurements along the waveguide: each set of measurement points $\mathbf{y}_k$ is at a distance $\Delta$ (the unit cell size) from its neighboring sets of measurement points, i.e., $\mathbf{y}_{k-1}$ and $\mathbf{y}_{k+1}$. Measurement points corresponding to different output vectors are represented with different markers.</figcaption></figure>

## 2.2. Subspace method for system identification

Subspace identification methods [49] can be applied to state-space models, such as (5) derived in the previous subsection, to extract the eigenvalues of the transfer matrix $\mathbf{A}$ from sensor measurements contained in $\mathbf{y}_k$. More specifically, the problem formulated in this study is well addressed by stochastic subspace identification algorithms [41], and considering the relatively small number of unit cells for the experimental characterization of waveguides addressed in this work, some considerations of the ESPRIT [34] are introduced. The short introduction hereinafter follows the work presented in Ref. [51].

Consider a waveguide consisting of $N$ unit cells. To limit boundary effects of the finite length waveguide, the first and the last unit cells are excluded when sensor data are collected. Then the effective number of unit cells is<fnref n="1" /> $N_o = N - 2$. Accordingly, these unit cells are numbered with $k = 1, 2, \dots, N_o$. Each unit cell involves $l$ measurement points, collected in $\mathbf{y}_k \in \mathbb{C}^l$ for the $k$-th cell.

The whole spatial measurement data of $N_o$ cells is arranged in a $l \times N_o$ matrix as

$$ \vec{\mathbf{Y}} = \begin{bmatrix} \mathbf{y}_1 & \mathbf{y}_2 & \dots & \mathbf{y}_{N_o} \end{bmatrix}, \qquad (8) $$

where $\vec{\cdot}$ indicates that the numbering of the $N_o$ unit cells is made in a direction chosen by default – e.g., from left to right in Figure 1(b).

The first step of subspace identification is to build a block-Hankel matrix $\mathcal{H} \in \mathbb{C}^{(pl) \times s}$ from the measurements $\vec{\mathbf{Y}}$ as follows:

$$ \mathcal{H} = \begin{bmatrix} \mathbf{y}_1 & \mathbf{y}_2 & \dots & \mathbf{y}_s \\ \mathbf{y}_2 & \mathbf{y}_3 & \dots & \mathbf{y}_{s+1} \\ \vdots & \vdots & \ddots & \vdots \\ \mathbf{y}_p & \mathbf{y}_{p+1} & \dots & \mathbf{y}_{N_o} \end{bmatrix}, \qquad (9) $$

where the integer $p$ is the spatial horizon, and $s = N_o - p + 1$. By repeatedly applying the state and output equations in (5), it is straightforward to show (see, e.g., [49] for the details) that

$$ \mathcal{H} = \mathcal{O}\mathcal{X} + \mathcal{N}, \qquad (10) $$

<footnote n="1">¹More generally, $N_o = N - 2d$ is considered to discard the $d$ extreme unit cells at both ends. This approach permits to eventually handle real-world structures including complex boundary conditions.</footnote>
with

$$ \mathcal{O} = \begin{bmatrix} \mathbf{C} \\ \mathbf{CA} \\ \vdots \\ \mathbf{CA}^{p-1} \end{bmatrix}, \quad \mathcal{X} = \begin{bmatrix} \mathbf{x}_1 & \mathbf{x}_2 & \dots & \mathbf{x}_s \end{bmatrix} \quad \text{and} \quad \mathcal{N} = \begin{bmatrix} \mathbf{n}_1 & \mathbf{n}_2 & \dots & \mathbf{n}_s \\ \mathbf{n}_2 & \mathbf{n}_3 & \dots & \mathbf{n}_{s+1} \\ \vdots & \vdots & \ddots & \vdots \\ \mathbf{n}_p & \mathbf{n}_{p+1} & \dots & \mathbf{n}_{N_o} \end{bmatrix}. \quad (11) $$

The matrix $\mathcal{O}$ can be estimated through the singular value decomposition (SVD) of $\mathcal{H}$ as<fnref n="2" />:

$$ \mathcal{H} = \begin{bmatrix} \mathbf{U}_1 & \mathbf{U}_2 \end{bmatrix} \begin{bmatrix} \mathbf{S}_1 & \mathbf{0} \\ \mathbf{0} & \mathbf{S}_2 \end{bmatrix} \begin{bmatrix} \mathbf{V}_1^H \\ \mathbf{V}_2^H \end{bmatrix} \approx \mathbf{U}_1 \mathbf{S}_1 \mathbf{V}_1^H \quad \text{and} \quad \hat{\mathcal{O}} = \mathbf{U}_1, \quad (12) $$

where $\mathbf{S}_1$ contains the $n_w$ singular values mainly due to $\mathcal{O}\mathcal{X}$, $\mathbf{S}_2$ contains the remaining small singular values due to the perturbation by the noise term $\mathcal{N}$, the approximation $\mathbf{U}_1\mathbf{S}_1\mathbf{V}_1^H$ neglects the small singular values, and thus $\hat{\mathcal{O}} = \mathbf{U}_1$ is an estimate of $\mathcal{O}$. This SVD truncation is done up to the model order $n_w$ – where $n_w$ stands for number of waves – which is tied physically to the guided waves propagating in the waveguide, it is usually known beforehand and rarely greater than 4 (in the case of bi-modal identification [31]).

When subspace identification is applied to periodic waveguides, due to the limited number of unit cells, it may be recommended to use $p = 2$, as in ESPRIT [36], to favor the parameter $s$ and therefore the column snapshots taken for the SVD procedure and the $\hat{\mathcal{O}}$ estimate. The impact of the number of unit cells is discussed in Appendix B.

Let $\hat{\mathcal{O}}_\uparrow$ and $\hat{\mathcal{O}}_\downarrow$ denote the parts of $\hat{\mathcal{O}}$ without the last block-row and without the first block-row, respectively. An estimate of $\mathbf{A}$ is then obtained by solving for $\hat{\mathbf{A}}$ the relationship $\hat{\mathcal{O}}_\uparrow\hat{\mathbf{A}} \approx \hat{\mathcal{O}}_\downarrow$ in the least-squares sense:

$$ \hat{\mathbf{A}} = (\hat{\mathcal{O}}_{\uparrow}^{\mathrm{H}}\hat{\mathcal{O}}_{\uparrow})^{-1}\hat{\mathcal{O}}_{\uparrow}^{\mathrm{H}}\hat{\mathcal{O}}_{\downarrow}. \quad (13) $$

In this approach, the estimate of $\mathcal{O}$ is not unique, nor is the estimate of $\mathbf{A}$. The estimates $\hat{\mathcal{O}}$ and $\hat{\mathbf{A}}$ given in (12) and (13) correspond to one of the possible solutions. It is important to note that this least-squares estimation could benefit from a larger spatial horizon $p > 2$, but this condition could also be detrimental for the number of snapshots $s$ for the SVD, as mentioned earlier. Nevertheless, all the possible solutions of $\mathbf{A}$ have the same set of eigenvalues [49]. These $n_w$ eigenvalues obtained from the eigendecomposition of $\hat{\mathbf{A}}$, represented in the vector $\hat{\lambda}$, are the estimated propagation constants, which are related to the estimated wavenumbers $\hat{\kappa}$ through $\hat{\lambda} = \exp(i\hat{\kappa}\Delta)$.

## 2.3. Symplectic system subspace identification

The general subspace method recalled in the previous subsection from the system identification literature does not take into account that the transfer matrix $\mathbf{A}$ is symplectic: its spectrum comes in reciprocal pairs $(\lambda, 1/\lambda)$ and its inverse $\mathbf{A}^{-1}$ is symplectic too. The physical meaning of this is that every guided wave mode – e.g., torsional, flexural – is propagated in two directions: from right to left, and inversely. In this particular case, the same identification procedure applied to the measurements arranged in the reverse direction (from the right side to the left side of the waveguide) yields, in principle, the same spectrum, but the numerical results are different in practice. A simple way to consider this fact would be to execute the identification procedure twice, each processing the same dataset arranged in one of the two opposite directions, and then to average their results. Though this averaging strategy improves the results, it is delicate to match the elements in the results of the two identification procedures in order to average them. Alternatively, to avoid the drawback of the averaging procedure and further improve the accuracy, a new identification procedure is proposed here by performing once subspace identification on the concatenated data arranged in the two directions, as detailed below.

The reversed measurement matrix $\overleftarrow{\mathbf{Y}}$ is built by collecting the spatial data in the opposite direction as

$$ \overleftarrow{\mathbf{Y}} = [\mathbf{y}_{N_o} \ \mathbf{y}_{N_o-1} \ \dots \ \mathbf{y}_1], \quad (14) $$

and is associated to the reversed state equation

$$ \mathbf{x}_k = \mathbf{A}^{-1} \mathbf{x}_{k+1}, \quad (15) $$

<footnote n="2">²Within the family of subspace identification methods, it is common to perform projections of the block-Hankel data matrices to asymptotically suppress the noise term $\mathcal{N}$, leading to covariance-driven and data-driven algorithms [41]. In the present derivation, no projection is done and the block-Hankel data matrix is used directly to estimate $\mathcal{O}$. This is in accordance with the ESPRIT formulation that was proposed either with autocovariances or data matrix [34]. The small number of unit cells compared to the large number of measurements per unit cell justify this choice.</footnote>
which, with $\mathbf{z}_k = \mathbf{x}_{N_o-k+1}$, can be written as a standard state equation

$$ \mathbf{z}_{k+1} = \mathbf{A}^{-1}\mathbf{z}_k. \qquad (16) $$

From (8) and (14), a concatenated measurement matrix $\mathbf{Y}$ is built as follows:

$$ \mathbf{Y} = \begin{bmatrix} \vec{\mathbf{Y}} \\ \left\langle \overleftarrow{\mathbf{Y}} \right\rangle \end{bmatrix} = \begin{bmatrix} \mathbf{y}_1 & \mathbf{y}_2 & \cdots & \mathbf{y}_{N_o} \\ \mathbf{y}_{N_o} & \mathbf{y}_{N_o-1} & \cdots & \mathbf{y}_1 \end{bmatrix}. \qquad (17) $$

Based on this concatenated measurement matrix, it would be quite natural to consider the identification problem of the augmented state-space system combining the state equations of (5) and (16):

$$ \left\{ \begin{aligned} \mathbf{x}_{k+1} &= \mathbf{A}\mathbf{x}_k \\ \mathbf{y}_k &= \mathbf{C}\mathbf{x}_k + \mathbf{n}_k \end{aligned} \right. \qquad (18) $$

with, for $\forall k \in \{1, 2, \dots, N_o - 1\}$,

$$ \mathbf{A} = \begin{bmatrix} \mathbf{A} & \mathbf{0} \\ \mathbf{0} & \mathbf{A}^{-1} \end{bmatrix}, \quad \mathbf{C} = \begin{bmatrix} \mathbf{C} & \mathbf{0} \\ \mathbf{0} & \mathbf{C} \end{bmatrix}, \quad \mathbf{x}_k = \left\{ \begin{array}{l} \mathbf{x}_k \\ \mathbf{x}_{N_o-k+1} \end{array} \right\}, \quad \mathbf{y}_k = \left\{ \begin{array}{l} \mathbf{y}_k \\ \mathbf{y}_{N_o-k+1} \end{array} \right\}, \quad \mathbf{n}_k = \left\{ \begin{array}{l} \mathbf{n}_k \\ \mathbf{n}_{N_o-k+1} \end{array} \right\}. \qquad (19) $$

This augmented system approach would lead to an estimate of the $2n \times 2n$ matrix $\mathbf{A}$ containing information about both $\mathbf{A}$ and $\mathbf{A}^{-1}$. Such a result would not be fundamentally different from the two identification procedures separately applied to $\vec{\mathbf{Y}}$ and $\overleftarrow{\mathbf{Y}}$. It does not explicitly take into account that $\mathbf{A}$ and $\mathbf{A}^{-1}$ both contain the same reciprocal eigenvalue pairs $(\lambda, 1/\lambda)$.

Instead of the seemingly natural augmented state approach, a new subspace identification method is proposed below, which estimates a single transfer matrix $\mathbf{A}$ from both $\vec{\mathbf{Y}}$ and $\overleftarrow{\mathbf{Y}}$, by explicitly accounting for the symplectic property of $\mathbf{A}$. It is numerically more efficient and results are more accurate than when two separate identification procedures are averaged. The key element in this new method is to relate the augmented output vector $\mathbf{y}_k$ to the original state vector $\mathbf{x}_k$ without using the augmented state vector $\mathbf{x}_k$, unlike the second equation of (18). In other words, in (5) the output vector will be doubled, i.e., the original vector $\mathbf{y}_k$ will be replaced by $\mathbf{y}_k$, but the original state vector $\mathbf{x}_k$ will remain unchanged. Of course, the matrix $\mathbf{C}$ will not be the same one. This new method is based on the following theorem.

**Theorem 1.** Consider the sequence $\mathbf{x}_k \in \mathbb{C}^{2n}$ recursively generated by

$$ \mathbf{x}_{k+1} = \mathbf{A}\mathbf{x}_k \qquad (20) $$

for $k = 1, 2, \dots, N_o$, from some given $\mathbf{A} \in \mathbb{C}^{2n \times 2n}$ and $\mathbf{x}_1 \in \mathbb{C}^{2n}$. Let $\mathbf{z}_k$ denote the reversed sequence, i.e.,

$$ \mathbf{z}_k = \mathbf{x}_{N_o-k+1}, \quad \text{for } k = 1, 2, \dots, N_o. \qquad (21) $$

If the $2n$ eigenvalues of the matrix $\mathbf{A}$ consist of $n$ reciprocal pairs such as $(\lambda, 1/\lambda)$, then there exists a matrix $\mathbf{\Gamma} \in \mathbb{C}^{2n \times 2n}$ such that

$$ \mathbf{z}_k = \mathbf{\Gamma}\mathbf{x}_k, \quad \text{for } k = 1, 2, \dots, N_o. \qquad (22) $$

A proof of this theorem is presented in Appendix C. To understand its meaning, first remark that $\mathbf{x}_k$ and $\mathbf{z}_k$ are related, because one is the reverse sequence of the other. For example, $\mathbf{z}_1 = \mathbf{x}_{N_o} = \mathbf{A}^{N_o-1}\mathbf{x}_1$ and $\mathbf{z}_2 = \mathbf{x}_{N_o-1} = \mathbf{A}^{N_o-3}\mathbf{x}_2$. In these trivial examples, the matrices $\mathbf{A}^{N_o-1}$ and $\mathbf{A}^{N_o-3}$, which respectively relate $\mathbf{z}_1$ to $\mathbf{x}_1$ and $\mathbf{z}_2$ to $\mathbf{x}_2$, are different. The non-trivialness of the theorem resides in the fact that the matrix $\mathbf{\Gamma}$ relating $\mathbf{z}_k$ to $\mathbf{x}_k$ is independent of $k$. In other words, $\mathbf{\Gamma}$ satisfies simultaneously the $N_o$ vectorial equations (22) corresponding to $k = 1, 2, \dots, N_o$, or more clearly, in matrix form

$$ \begin{bmatrix} \mathbf{z}_1 & \mathbf{z}_2 & \cdots \mathbf{z}_{N_o} \end{bmatrix} = \mathbf{\Gamma} \begin{bmatrix} \mathbf{x}_1 & \mathbf{x}_2 & \cdots \mathbf{x}_{N_o} \end{bmatrix}. $$

This matrix equation totalizes $2nN_o$ scalar equations, but involves $(2n)^2$ unknowns, which are the entries of $\mathbf{\Gamma}$. Since $2nN_o < (2n)^2$ in general, it is thus nontrivial that there exists a $\mathbf{\Gamma}$ satisfying simultaneously all the equations.

Based on this theorem, the symplectic state-space model to be considered for subspace identification is

$$ \left\{ \begin{aligned} \mathbf{x}_{k+1} &= \mathbf{A} \mathbf{x}_k \\ \mathbf{y}_k &= \mathbf{C} \mathbf{x}_k + \mathbf{n}_k \end{aligned} \right. \quad \text{with} \quad \mathbf{C} = \begin{bmatrix} \mathbf{C} \\ \mathbf{C} \mathbf{\Gamma} \end{bmatrix}, \quad \mathbf{y}_k = \left\{ \begin{array}{l} \mathbf{y}_k \\ \mathbf{y}_{N_o-k+1} \end{array} \right\} \quad \forall k \in [1, N_o - 1]. \qquad (23) $$
Then the proposed new method consists in applying the standard subspace identification method described in Section 2.2 to the newly formulated state-space model (23). More specifically, in the Hankel matrix $\mathcal{H}$ in (9), each block $\mathbf{y}_k$ is replaced by $\mathbf{y}_k$, and the procedure described in Section 2.2 is applied to the concatenated dataset (17).

## 2.4. Algorithm of the subspace-based wavenumber identification (SWAN)

The Algorithm 1 summarizes the proposed method for wavenumber identification. Recall that the standard state-space model of (5) can be used for identification either with the measurement collection $\mathbf{Y}$ from (8) or $\mathbf{\tilde{Y}}$ from (14). In this paper, the proposed symplectic state-space model (23) is used to present the numerical and experimental results of the next sections. This symplectic state-space model uses the measurement collection $\mathbf{Y}$, which is the vertical concatenation of $\mathbf{\tilde{Y}}$ and $\mathbf{\tilde{Y}}$ as described in (17). Lastly, the spatial horizon is more often used as $p = 2$ as in the original ESPRIT version [34], to increase the parameter $s$ and thus enhance the estimation of $\mathcal{O}$.

**Algorithm 1 SWAN**

1: **function** SWAN($\mathbf{Y}_s, n_w; p = 2$) 	 $\triangleright \mathbf{Y}_s$ can be chosen as $\mathbf{Y}$ Eq. (17), $\mathbf{\tilde{Y}}$ Eq. (8) or $\mathbf{\tilde{Y}}$ Eq. (14)
2: Build the block-Hankel matrix $\mathcal{H}$ from Eq. (9). 	 $\triangleright p \ge 2$, from Section 2.2.
3: Compute the SVD of $\mathcal{H}$, truncate $\mathbf{U}$ up to $n_w$ to have $\hat{\mathcal{O}}$ from Eq. (12).
4: Estimate $\hat{\mathbf{A}}$ by least squares from Eq. (13).
5: Compute $\hat{\mathbf{\lambda}}$ from the eigendecomposition of $\hat{\mathbf{A}}$.
6: **return** $\hat{\mathbf{\lambda}}$

## 3. An adaptive and general-purpose spectral mask for evanescent wave rejection

In this section, a straightforward strategy for wavenumber selection is derived. The proposed selection strategy acts directly on the identification results and can be used for any wavenumber identification method (therefore, the term general-purpose). In most wavenumber identification methods, the number of estimated waves $n_w$ is a user-defined parameter – the number of waves propagating in the structure is known a priori. A statistical model selection can be performed a priori (e.g., ESTER criterion [31, 43]) or a posteriori (e.g., Akaike information criterion [52]) to identify the parameter $n_w$ automatically. Alternatively, the proposed strategy is based on the physical properties of the identified wavenumbers, which results in their selection among the $n_w$ identified waves – keeping $n_w$ as an user-defined parameter. The selection procedure can be seen as a masking one, since some of the $n_w$ solutions obtained at a specific frequency will be kept and the rest rejected (or masked). It should be stated that the spectral mask is built ad hoc. Still, its application is physically intuitive, and its underlying parameters could be visually inferred from a first glance at the identification results.

<figure><img src="image_4.png" /><figcaption>Figure 3: Two spectral mask examples: (a) Using individual polar equation, see (24), with parameters $\rho(\theta, \phi = \pi, r_1 = 0.95, r_2 = 0, \eta = 90)$, where a notch is created to admit the eigenvalues located in the first Brillouin edge. (b) Using general polar equation, see (25), with parameters $\rho_m(\theta, \phi = \pi, r_1 = 0.95, r_2 = 0, \eta = 90)$, $\rho(\theta, \phi = -\pi/3, r_1 = 0.95, r_2 = 0.5, \eta = 20)$), where an additional notch is created to admit evanescent waves.</figcaption></figure>
A propagative wave is characterized by $|λ| ≈ 1$, otherwise, it is an evanescent wave. In periodic media, a Bragg bandgap occurs when $arg(λ) = π$ – in these regions, $|λ| < 1$. Suppose the propagation constants are plotted in a polar coordinate system. In that case, these two facts prioritize the selection of propagation constants near the unit circle – when they are propagative – or at an angle of $π$, the first Brillouin edge – where Bragg bandgaps exist. A polar curve that fits well this description is depicted in Figure 3(a), where one can see it mostly near the unit circle, but a sharp notch is created at an angle of $π$ to allow the selection of the propagation constants identified in bandgap regions. To mask the estimated values of $λ$ without physical meaning, we can limit their modulus to the one given by this polar curve, proposed in the following equation:

$$ \rho(\theta, \phi, R, r, \eta) = R - (R - r) \left( \frac{|\operatorname{Arg}(e^{i(\theta - \phi + \pi)})|}{\pi} \right)^{\eta} \quad (24) $$

where $Arg(·)$ is the principal argument, i.e., the argument of the complex number in the interval $]-π, π]$, $φ$ is the targeted angle for the notch, $R$ and $r$ are the outer and inner radii, and $η$ is an exponent controlling the width of the notch. For masking purposes, all propagation constants $λ$ should be inside the unit circle: if a certain $λ$ is outside of it ($|λ| > 1$), its reciprocal value $1/λ$ is used to decide whether it should be masked or not.

As a general guideline, the four manually defined parameters $φ, R, r$, and $η$ – which are based on the physical meaning of the identified wave modes – have numerical values around: $φ = π$ for the first Brillouin edge, $R ≈ 1$ for the unit circle, $r ≈ 0.5$ to capture evanescent waves at Brillouin edges, and $η ≈ 90$ for moderately sharp notches. The reference values for $R$ and $φ$ are mostly the same for different measurement setups and data. However, $r$ and $η$ should be adjusted based on the filtered dispersion curve examination, since their values are more sensitive to measurement noise, as it will be shown in Section 5.3.1.

A more general polar equation $ρ_m = ρ_m(θ, ...)$ can be written as a combination of $k$ different masks, as depicted in Figure 3(b) and given by the equation:

$$ ρ_m(θ, ..., ) = \min(ρ_1(θ, ..., ), ..., ρ_k(θ, ..., )). \quad (25) $$

Finally, a piecewise function of $p$ segments in the frequency domain of interest is described by $P(f, θ, ...)$ as:

$$ P(f, \theta, \dots) = \begin{cases} \rho_{m,1}(\theta, \dots) & \forall f \in [f_0, f_1] \\ \rho_{m,2}(\theta, \dots) & \forall f \in [f_1, f_2] \\ \vdots & \\ \rho_{m,p}(\theta, \dots) & \forall f \in [f_{p-1}, f_p]. \end{cases} \quad (26) $$

This wavenumber selection strategy will be applied to the following sections’ numerical and experimental characterization results.

## 4. Numerical validation

This section describes a validation of the proposed approach with a numerical example. The experiment consists of harmonic motion simulations for a periodic waveguide – an aluminum beam with periodic material removal which will be also considered in the experimental section. Results of the identification method for a virtual experiment are presented and compared to those provided by the Prony’s method adapted for Bloch wave identification [10] and INCOME [19]. In the scientific literature, these are the only two available methods to identify wavenumbers in periodic waveguides from multiple measurements per unit cell. These methods are thus considered *multi-output*. Moreover, as the proposed symplectic state-space model (23), INCOME is a method that takes into account the reciprocity between the propagation constants $λ$ and $1/λ$ through the change of variable $μ = λ + 1/λ$. It should be noted that this is not the case for the Prony’s method. Both methods are multimodal – i.e., several physical modes can be identified at a given frequency – but INCOME’s implementation is not straightforward for $n_w > 4$, since a change of variable has to be applied to obtain a palindromic polynomial [19].

A geometric model of the unit cell and its 3D quadratic mesh is first created in Gmsh [53], and the corresponding finite element model is built with the material properties of the unit cell. The unit cell is depicted and described in Figure 4(a). The dispersion curves of the unit cell are then computed by solving a $κ(ω)$ generalized eigenvalue problem within the WFEM framework [54, 55]. These dispersion curves are shown in Figure 4(b) (real part $\Re(κ)$ – upper part – and imaginary part $\Im(κ)$ – lower part). A color scale is superimposed on each curve to show the degree of evanescence of the wave mode by the modulus of the propagation constant – a wave mode is said to be propagative if $|λ| ≈ 1$.
<figure><img src="image_5.png" /><figcaption>Figure 4: (a) A 3D model of the periodic unit cell. Its mechanical properties are: Young's modulus $E = 69$ GPa, Poisson's ratio $\nu = 0.33$, density $\rho = 2700$ kg/m$^3$ and structural loss factor $\eta = 0.005$. Its geometric properties are a length $\Delta = 40$ mm, a width $w = 35$ mm, a thickness $t_r = 1.6$ mm (reduced thickness length is $3/4\Delta$). (b) Dispersion curves for the periodic waveguide (upper: real part of the complex wavenumber $\Re(\kappa)$, lower: imaginary part $\Im(\kappa)$). In-plane motion (shear and compression) get coupled with out-of-plane motion (bending and torsion) between 8 and 11 kHz.</figcaption></figure>

To simulate the data for the identification problem, a periodic waveguide is assembled from the concatenation of $N = 11$ unit cells. The periodic waveguide FE matrices are created from the assembly of the unit cell FE matrices. Finite element and modal analyses use a truncated modal basis that includes the waveguide's first 100 normal modes, with a maximum frequency of 21679 Hz. The frequency response functions are generated for a unit point force $F_0 = 1$ N applied at $x_0 = 25$ mm (along the length), $y_0 = 10$ mm (along the width), $z_0 = 0$ mm (along the height, i.e., on the plane surface). The spatial measurements of mobility (point velocity/injected force at a fixed frequency) are used for identification. Three different levels of noise are added to the spatial measurements to compare the identification method performances with a small number of unit cells – which is often the case in experimental setups. The added white Gaussian noise is generated to provide signal-to-noise ratios (SNR) of 100, 10, and 1, which represent different measurement scenarios, including favorable (100), less favorable (10), and clearly unfavorable (i.e., 1, where the noise power is comparable to the signal power). For the experimental tests carried on this work, the lowest SNR values are already above 100 (vibrometry) and 10 (deflectometry).

The signal-to-noise ratio is the ratio between the power of the noiseless signal to the power of the noise, which for the considered case is expressed as:

$$ \text{SNR} = \frac{P_{\text{signal}}}{P_{\text{noise}}} = \frac{\sum_{i=1}^{N_o \times l} |\mathbf{s}_i|^2}{\sum_{i=1}^{N_o \times l} |\mathbf{n}_i|^2} \quad (27) $$

where $\mathbf{s} = \mathbf{s}(\omega)$ contains the $N_o \times l$ noiseless spatial measurements – e.g., mobility – at the angular frequency $\omega$. The corresponding noise terms are contained in $\mathbf{n} = \mathbf{n}(\omega)$, which is derived from a predefined SNR value.

The symplectic state-space model (23) is employed for identification. The number of outputs per unit cell is $l = 259$, the data length is $N_o = 9$, the number of waves to estimate is $n_w = 4$ and the spatial horizon is $p = 2$. The size of the block-Hankel matrix $\mathcal{H}$ is $518 \times 8$. The first $n_w = 4$ left singular vectors of $\mathcal{H}$ are used as an estimate of the observability matrix $\tilde{\mathcal{O}}$. Accordingly, four $\tilde{\lambda}$ eigenvalues of $\tilde{\mathcal{A}}$ are retained. Before the conversion of propagation constants into complex wavenumbers, the proposed wavenumber selection strategy summarized in (26) is applied to keep the identified wave modes with physical meaning (detailed parameters of this equation will be presented in the next section). Results of this identification and a posteriori selection of wave modes are depicted in Figure 5(c). Results of INCOME and Prony's method are shown in Figures 5(a) and 5(b), respectively.

Computation times per frequency for INCOME, Prony and Subspace methods were 0.12 ms, 0.23 ms and 0.45 ms respectively. In Figure 5, 12800 iterations were considered for a frequency resolution of 1 Hz, which implies that dispersion curves were computed in 1.5 s, 2.9 s, and 5.8 s respectively. These computation times are quite small for the purpose addressed in this work – i.e., experimental characterization – and they are representative of a regular procedure: dispersion curve identification takes less than 10 s.

Some comments can be drawn from Figure 5. First, the three implemented methods give correct estimations
<figure><img src="image_6.png" /><figcaption>Figure 5: Complex wavenumber identification results for the periodic aluminum beam from simulated vibration measurement fields (multi-output measurements with *l* = 259 measurement points per unit cell and *N*<sub>o</sub> = 9 unit cells considered for identification). Results for the proposed method (subspace identification + spectral mask) are compared against two relevant methods for wavenumber identification in periodic waveguides: Prony and INCOME. Identification results are obtained from noiseless simulations (top row) and noisy ones with three different SNR values : SNR=1 (2nd row), SNR=10 (3rd row) and SNR=100 (4th row).</figcaption></figure>

of the imaginary part of the wavenumber within the bandgap region (between 4 and 6.5 kHz, approximately). However, Prony's method does not provide two superimposed branches of the imaginary part of the wavenumber in the bandgap region. This is attributed to the fact that no symplectic considerations are included in this formulation. Regarding the wavenumber's real part, the proposed method gives results similar to Prony and INCOME results. Despite the small number of unit cells, all three methods generally provide good results (recall that *N*<sub>o</sub> = 9). Regarding the behavior of the identification methods in presence of noise, the proposed method produces more clean curves at all the SNR levels, including the noise-free case. Finally, note that the spectral mask was applied to the proposed subspace method, but these simulations' lack of noisy estimates does
not showcase the mask’s potential, which will be addressed quantitatively in terms of false positive rates in the
experimental results section.

**5. Experimental applications**

**5.1. General context for full-field vibration measurements**

This section validates the proposed approach to wavenumber identification in 1D periodic structures through
experimental applications. Full-field measurements are carried out to collect sufficient measurements inside
each unit cell of the periodic waveguide. This strategy enhances the least-squares minimization involved in the
identification methods for periodic waveguides. In the case of homogeneous waveguides, an arbitrary unit cell
can be created with the same goal. Since the proposed approach can be combined with any full-field vibration
measurement, we briefly introduce the two experimental measurement methods that were used.

Over the last three decades, several contactless and potentially full-field vibration measurement techniques have emerged [56]. Significant progress in laser technology first led to the development of integrated contactless measurements using laser Doppler vibrometry (LDV) in the late 1980s. The basic principle is detecting the Doppler frequency shift of the laser beam reflected/scattered by a moving surface [45]. Single-point measurements evolved to scanning and multi-point LDV systems. Considered a gold standard for vibration measurements [57], scanning LDV is used for the first application case addressed in this work, i.e., the periodic aluminum beam.

In parallel, digital cameras paved the way for developing full-field measurement methods. Indeed, the 2D nature of digital sensors – CCD or CMOS – enables capturing spatial information at a typically high frame rate, avoiding scanning operation and allowing practical study of transient phenomena. Among optical methods such as digital image correlation (DIC) [57] or digital holography [56], the optical deflectometry – which has its origins in optical metrology [46] to inspect specular reflecting surfaces – has proved to be a practical and cost-effective solution to time-resolved and dense spatial vibration measurements [58]. In the present case, the deflectometry test uses a high-speed camera, a printed regular grid, and a tested structure with a specular surface to reflect the printed grid. The deflectometry technique has even been extended to infrared measurements [59] to enable testing of structures with non-specularly reflective surfaces. Full-field optical measurement techniques allow measuring all the points simultaneously over the studied structure’s surface, reducing the testing time required by the point-wise sequential measurements in LDV. Deflectometry tests have been carried out for the second application case addressed in this work, a homogeneous beam.

**5.2. Experimental setup**

**5.2.1. Periodic beam**

A 3D scanning laser Doppler vibrometer (S-LDV) has been utilized to collect vibration measurements on a periodic aluminum beam. This beam’s geometrical and mechanical properties were detailed in the previous section. The beam’s boundary conditions were clamped-free, and it was subjected to point impulse excitations using an automatic impact hammer (Maul-Theet vImpact-63 model, positioned at 10 mm along the width and 25 mm along the height). Frequency-domain vibration measurements consist of transfer function spectra from 1 to 12 800 Hz (computed as a mean from 10 hammer impacts), with a frequency resolution of 1 Hz. Regarding the spatial resolution of the measurements, 990 measurement points have been placed on a regular rectangular grid, so that each one of the *N* = 11 unit cells of the beam contains *l* = 90 measurement points per unit cell (9 × 10 measurement points with 4 mm spacing in each direction). Figure 6(a) depicts the experimental setup for this periodic beam. A spatially-averaged transfer function over all measurement points (|*H*₁| mobility, point velocity/injected force) is depicted in Figure 6(b). A comparison with the simulation results (see the previous section) indicates good agreement across most of the frequency range of interest. Some mismatch is observed from 7 kHz onwards, potentially linked to the wave coupling phenomena between in-plane and out-of-plane wave modes.

**5.2.2. Homogeneous beam**

A high-speed camera (Fastcam mini AX200) was used to collect images of the dynamic motion of a homogeneous
aluminum beam. The tested beam has similar dimensions to the periodic one, but does not include a periodic
thickness removal and has a constant thickness of *h* = 3.2 mm. As for the periodic beam, clamped-free boundary
conditions are considered. However, the beam was subjected to a stationary excitation using a vibration shaker
(B&K 4810 model) – white noise excitation over the 0 – 20 kHz frequency range. An impedance head placed
at 10 mm along the width and 25 mm along the height was used to measure the injected force (PCB 288D01
model). A LED spotlight was needed to obtain sufficient brightness in the images at such a high sample rate,
inversely related to the exposure time. The experimental setup for this homogeneous beam is depicted in
Figure 7.
<figure><img src="image_7.png" /><figcaption>Figure 6: (a) The periodic aluminum beam consisting of $N = 11$ unit cells being tested in a clamped-free condition with a 3D-laser Doppler vibrometer and an automatic impact hammer. (b) Average mobility spectrum (point velocity / injected force). A comparison with the results from a numerical model is also shown.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Figure 7: Experimental setup for the deflectometry tests. The high-speed camera captures the images of the regular grid reflected on the beam with mirror-like finishing. The aluminum beam is excited with a shaker, and an impedance head measures the injected force. A LED spotlight is used to increase the brightness of the images captured at a high sample rate.</figcaption></figure>

The deflectometry post-processing necessary to convert dynamic motion images into measurement spectra is briefly described in the following steps (for a detailed explanation, the reader can refer to [56]). Firstly, each grid captured in the images – as shown in Figure 8(b) – is compared to that of a reference image (when the structure is not excited) in a region of interest of the image to retrieve the spatial phase maps that are encoded by the grid – and therefore the beam – dynamics. Secondly, the phase maps are converted to slope fields of the structure with a 2D convolution kernel. The third step involves a spatial integration of the slope fields to recover the displacement fields (the integration constant is chosen such that the displacement of a reference point, as close as possible to the clamping zone, is zero). Finally, frequency-domain response functions are computed using the displacement fields and the injected force measurement.

For the measurements conducted in this study, 174 730 images were collected at 20 000 frames per second (measurement time was 8.74 s). Regarding the spatial resolution of the measurements, displacement fields were recovered as a regular rectangular grid of 12 × 220 measurement – 2 640 points in total. The spacing between points in each direction was 2 mm so that the measurement area was 24 mm × 440 mm. This spatial step of 2 mm was the practical limitation of the set-up, and therefore the highest achievable spatial resolution to capture the full beam with the camera up to the desired frequency. Frequency-domain vibration measurements consist of transfer function spectra from 1 to 4 500 Hz, computed from cross power density estimates with 50% overlapping windows using Welch's periodogram method, and targeting a frequency resolution of 1 Hz. A spatial average of the 2 640 transfer function estimations moduli $|H_1|$ in the form of receptance (point displacement/injected force) is depicted in Figure 8(c), where a comparison with the vibration response results from a numerical model is done, showing acceptable agreement across most of the frequency range of interest. Some modal mismatch is generally observed, potentially linked to the more complex numerical modeling of the impedance head, stinger and cables affecting the dynamic motion of the beam. It is important to note that the average spectrum is noisier than the one obtained with the S-LDV in Figure 6(b). Thus, wavenumber identification is expected to be of lower accuracy for this deflectometry test.
<figure><img src="image_9.png" /><figcaption>Figure 8: (a) The homogeneous aluminum beam with mirror-like finishing being tested in a clamped-free condition with the deflectometry technique. (b) Example of an image captured by the high-speed camera of the regular grid reflected on the aluminum beam. Comparing each image of the dynamic motion with a reference image, the 2D spatial phase maps are retrieved. (c) Average receptance spectrum (point displacement / injected force). A comparison with the results from a numerical model is also shown.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Figure 9: Complex wavenumber identification results for the periodic aluminum beam from multi-output vibration measurements ($l = 90$ measurement points per unit cell, $N_o = 9$ unit cells considered for identification). Results for the proposed method (subspace identification + spectral mask) are compared against two relevant methods for wavenumber identification in periodic waveguides: Prony's method and INCOME.</figcaption></figure>

## 5.3. Experimental results and discussions

The following presents complex wavenumber identification results with the proposed subspace method and spectral mask for both the periodic and the homogeneous (artificially periodic) beams. These results are compared to those provided by Prony's method adapted for Bloch wave identification [10] and INCOME [19], which are the only multi-output and multi-modal identification methods reported in the scientific literature.

### 5.3.1. Periodic beam

Complex dispersion curves for the periodic aluminum beam are illustrated in Figure 9. The identification results were obtained with the symplectic state-space model (23). Impact excitation is on the lower end of the beam (near the clamped end), so that the evanescent wave is propagated over a certain distance, but does not reach to the free end. For the multi-output identification, there are $l = 90$ measurement points per unit cell and $N_o = 9$ unit cells considered. The spatial horizon is $p = 2$ and the number of waves to be identified is $n_w = 4$ ($n'_w = n_w/2 = 2$ for INCOME, which estimates the pair of reciprocal propagation constants under the change of variable $\mu = \lambda + 1/\lambda$).

Some comments can be drawn from Figure 9. Firstly, all three implemented methods yield good results for the real part of the wavenumber $\Re(\kappa)$ (shown in red). Regarding its imaginary part $\Im(\kappa)$, Prony's method presents some difficulties characterizing the bandgap. Using INCOME improves results, where a trend of $\Im(\kappa)$
<figure><img src="image_11.png" /><figcaption>Figure 10: Results of applying a spectral mask for the identified propagation constants. (a) In the pass band region (up to $f_1 = 3.6$ kHz), the spectral mask is simply a circle with radius $R = 0.97$. (b) In the bandgap region (torsional bandgap between $f_1$ and $f_2 = 6.7$ kHz), a small notch is created at an angle $\phi = \pi$ to permit the selection of eigenvalues (propagation constants) located in the first Brillouin edge. (c) In the more intricate coupling region (between $f_2$ and $f_{max} = 12.8$ kHz), a more complex polar equation is used to capture evanescent behavior in the coupling branches: a first notch is located at an angle $\phi = 0$ which captures waves located in the second Brillouin edge, two more notches capture partially the wave coupling phenomena observed in Figure 4 between 8 and 11 kHz.</figcaption></figure>

in the bandgap region can be observed. The proposed subspace method also permits us to identify that trend at approximately the same accuracy level as INCOME. However, for a clear interpretation of $\Im(\kappa)$, the spectral mask proposed in Section 3 has been used. The $\Im(\kappa)$ resolution in the bandgap region is enhanced in this way. Moreover, the spectral mask applied in the bandgap region effectively denoises the $\Re(\kappa)$ results.

The proposed frequency-dependent spectral mask was applied for this case in three different frequency ranges: (i) in the pass band region – where elastic waves are purely propagative –, from $f_0 = 0$ to $f_1 = 3.6$ kHz, (ii) in the torsional bandgap region, between $f_1$ and $f_2 = 6.7$ kHz, and (iii) in the in-plane / out-of-plane coupling region, between $f_2$ and $f_{max} = 12.8$ kHz. A detailed description and results of the masking strategy are shown in Figure 10. After applying the spectral mask, 27% of the identified eigenvalues – the propagation constants – were masked, which can be seen both in Figures 9 and 10. This means that 27% of the identified propagation constants are considered as noise, and therefore, the false positive rate is 27%.

To emphasize the fact that a significant benefit is obtained by including several measurement points per unit cell (multi-output identification), identification results considering just one measurement point per unit cell
<figure><img src="image_12.png" /><figcaption>Figure 11: Complex wavenumber identification results for the periodic aluminum beam from single-output vibration measurements ($l = 1$ measurement point per unit cell, $N_o = 9$ unit cells considered for identification). Results for the proposed method (subspace identification + spectral mask) are compared against two relevant methods for wavenumber identification in periodic waveguides: Prony's method and INCOME.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Figure 12: Complex wavenumber identification results for the homogeneous aluminum beam from multi-output vibration measurements ($l = 240$ measurement points per artificial unit cell, $N_o = 9$ artificial unit cells considered for identification). Results for the proposed method (subspace identification + spectral mask) are compared against two relevant methods for wavenumber identification in periodic waveguides: Prony's method and INCOME.</figcaption></figure>

($l = 1$, single-output identification) are shown in Figure 11. It is clearly visible that identification results for $\Re(\kappa)$ are of lower quality than the results obtained with $l = 90$ in Figure 9. It should be noted that for $l = 1$, the data matrix (9) is Hankel instead of block-Hankel. Thus, the results obtained with the subspace method in Figure 11 are equivalent to those that would have been obtained with the current adaptation of ESPRIT method for wavenumber identification [31]. Note, however, that the spectral mask enhances the interpretation of the results, with a false positive rate of 33%.

### 5.3.2. Homogeneous beam (arbitrary periodicity)

Complex dispersion curves for the homogeneous aluminum beam are illustrated in Figure 12. The identification results were obtained with the symplectic state-space model from (23), since there are only propagating waves to be identified – i.e., reciprocal pairs are present through the whole frequency range. For the multi-output identification, $l = 240$ measurement points per unit cell and $N_o = 9$ artificial unit cells are considered. The spatial horizon is $p = 5$ to improve the least-squares estimation. The number of waves to be identified is $n_w = 4$ ($n'_w = n_w/2 = 2$ for INCOME, which estimates the pair of reciprocal propagation constants under the change of variable $\mu = \lambda + 1/\lambda$). When presenting the INCOME method, the idea of considering artificial unit cells is reported in [19].

Some comments can be drawn from Figure 12. Firstly, all three implemented methods yield fairly good results for the real part of the wavenumber $\Re(\kappa)$ (shown in red), following analytical dispersion curves for the bending and torsional waves. The proposed subspace method gives less noisy estimates, thanks to the symplectic
identification considered in the state-space model (23). The imaginary part $\Im(\kappa)$ is not shown, since the noisy measurements do not provide good estimations of spatial attenuation, which is very low in this aluminum beam. The spectral mask proposed in Section 3 has been used for this homogeneous beam example. It improves the interpretation of $\Re(\kappa)$, and a single polar equation from (24) is used with $R = r = 0.97$, i.e., a circle. For this example, the false positive rate is 77%.

## 6. Conclusions

This paper proposes an original subspace identification method for complex wavenumber estimation in periodic waveguides. The harmonic motion has been formulated through a wave finite element framework, but also in a state-space form. The derived state-space model has been directly used for both standard and symplectic subspace identification. It was proved that the proposed symplectic identification leads to the original state-transition matrix.

Subspace identification schemes are naturally developed in a multi-output fashion – multiple measurements per unit cell in the waveguide. The experimental examples provided in this work illustrate that complex wavenumber identification in periodic waveguides can be enhanced with multiple measurements per unit cell, mainly when bandgaps occur. Moreover, the proposed spectral mask improves the dispersion curves interpretation, especially in the bandgap regions. In that sense, the proposed subspace method and the spectral mask will ultimately lead to a practical assessment of vibration attenuation in periodic waveguides, and enhanced characterization of complex periodic media.

## Appendix A. The output equation involving sensor measurements inside unit cells

The output equation in (5) implies that $y_k$ depends linearly on $x_k$. If $y_k$ contains measurements made inside the $k$-th unit cell (off cell boundaries), then the corresponding components of $y_k$ are not included in the state vector $x_k$ defined in (5). In this case, $y_k$ is not directly related to $x_k$. Before the efficient solution proposed towards the end of this appendix, two obvious solutions based on augmented states are first described, in order to contrast the advantages of the proposed solution.

### Appendix A.1. Augmenting the state with all the inner nodal displacements

The state vector $x_k$ defined in (5) is known as the condensed state vector, which is associated to the condensed state equation formulated in (5). Instead of this condensed version, an augmented state vector, obtained by completing $x_k$ with more components, can include $y_k$ as a sub-vector. In particular, a state vector and the associated state equation can be formulated based on the harmonic equation (2) before the dynamic condensation step summarized in equations (3) and (4). The state vector involving both inner and boundary nodal displacements

$$x_k^{\text{full}} = \begin{Bmatrix} u_k^L \\ u_k^I \end{Bmatrix} \qquad (A.1)$$

satisfies the state equation

$$x_{k+1}^{\text{full}} = A^{\text{full}} x_k^{\text{full}} \qquad (A.2)$$

with a large transfer matrix $A^{\text{full}}$ relating $x_k^{\text{full}}$ to $x_{k+1}^{\text{full}}$. Any measurement point, after finite element discretization, corresponds to a component of $x_k^{\text{full}}$, then the output equation

$$y_k = C^{\text{full}} x_k^{\text{full}} + n_k \qquad (A.3)$$

holds with a matrix $C^{\text{full}}$ filled with 0's and 1's.

The size of this transfer matrix $A^{\text{full}}$ is $(m+n) \times (m+n)$, which is much larger than the size $2n \times 2n$ of the condensed transfer matrix $A$ defined in (5).

### Appendix A.2. Augmenting the state with duplicated state vectors

A second solution is based on the fact that, for a $\Delta$-periodic waveguide, the positions of the unit cell boundaries can be arbitrarily chosen, provided they are separated by the same unit cell length $\Delta$, which is inherent to the waveguide. Consider different configurations of unit cell boundaries indexed by $j = 1, 2, \dots, l$. For the $j$-th configuration, a state vector $x_k^{(j)}$ corresponding to the left boundary nodes of the $k$-th unit cell is defined, in a
way similar to the definition of the state vector $\mathbf{x}_k$ in (5). The associated transfer matrix is denoted by $\mathbf{A}^{(j)}$. Concatenating all the duplicated state vectors and their associated state equations then yields

$$ \mathbf{x}_{k+1}^{\text{dup}} = \mathbf{A}^{\text{dup}} \mathbf{x}_k^{\text{dup}}, \qquad (A.4) $$

with

$$ \mathbf{x}_k^{\text{dup}} = \begin{Bmatrix} \mathbf{x}_k^{(1)} \\ \mathbf{x}_k^{(2)} \\ \vdots \\ \mathbf{x}_k^{(l)} \end{Bmatrix}, \quad \mathbf{A}^{\text{dup}} = \begin{bmatrix} \mathbf{A}^{(1)} & 0 & \cdots & 0 \\ 0 & \mathbf{A}^{(1)} & \cdots & \vdots \\ \vdots & \vdots & \ddots & 0 \\ 0 & \cdots & \cdots & \mathbf{A}^{(l)} \end{bmatrix}. \qquad (A.5) $$

Suppose that the different configurations of unit cell boundaries have been chosen so that each measurement point falls on a boundary. Then the components of $\mathbf{y}_k$ containing measurements are included in the sub-vectors of $\mathbf{x}_k^{\text{dup}}$. Therefore,

$$ \mathbf{y}_k = \mathbf{C}^{\text{dup}} \mathbf{x}_k^{\text{dup}} + \mathbf{n}_k \qquad (A.6) $$

holds with a matrix $\mathbf{C}^{\text{dup}}$ filled with 0's and 1's. Each of the diagonal blocks of $\mathbf{A}^{\text{dup}}$, namely $\mathbf{A}^{(j)}$, has the same size as the condensed transfer matrix $\mathbf{A}$ introduced in (5).

Similar to the previous solution, this one is also based on an augmented state vector and has the same drawback of large state vector and large transfer matrix.

### Appendix A.3. Output equation without augmenting the state

It is in fact possible to keep the condensed state vector $\mathbf{x}_k$ as defined in (5), yet to relate $\mathbf{y}_k$ containing inner measurement points to $\mathbf{x}_k$, with a matrix $\mathbf{C}$ filled with complex numbers, instead of the 0's and 1's in the previous solutions. The key idea is to express the inner nodal displacements $\mathbf{u}_k^I$ as a linear function of $\mathbf{x}_k$. In (3) $\mathbf{u}_k^I$ is already expressed in terms of $\mathbf{u}_k^L$ and $\mathbf{u}_{k+1}^L$. The former is part of $\mathbf{x}_k$ as defined in (5), but not the latter. However, $\mathbf{u}_{k+1}^L$ is part of $\mathbf{x}_{k+1}$, which is equal to $\mathbf{A}\mathbf{x}_k$ according to the state equation in (5). By combining (3) and (5),

$$ \mathbf{u}_k^I = \mathbf{F}^I \mathbf{x}_k \qquad (A.7) $$

with the condensed state vector $\mathbf{x}_k$ as defined in (5), and

$$ \mathbf{F}^I = \left[ \left( \mathbf{D}_{II}^{-1} \mathbf{D}_{IR} \tilde{\mathbf{D}}_{LR}^{-1} \tilde{\mathbf{D}}_{LL} - \mathbf{D}_{II}^{-1} \mathbf{D}_{IL} \right) \left( -\mathbf{D}_{II}^{-1} \mathbf{D}_{IR} \tilde{\mathbf{D}}_{LR}^{-1} \right) \right]. \qquad (A.8) $$

Suppose that $\mathbf{y}_k$ contains not only part of $\mathbf{u}_k^L$, but also part of $\mathbf{u}_k^I$, which consists of the inner nodal displacements. Then

$$ \mathbf{y}_k = \mathbf{G}^L \mathbf{u}_k^L + \mathbf{G}^I \mathbf{u}_k^I + \mathbf{n}_k, \qquad (A.9) $$

with some matrices $\mathbf{G}^L$ and $\mathbf{G}^I$ filled with 0's and 1's selecting the components of $\mathbf{u}_k^L$ and $\mathbf{u}_k^I$ contained in $\mathbf{y}_k$. It then yields

$$ \mathbf{y}_k = \mathbf{G}^L \mathbf{F}^L \mathbf{x}_k + \mathbf{G}^I \mathbf{F}^I \mathbf{x}_k + \mathbf{n}_k, \qquad (A.10) $$

where $\mathbf{F}^L$ is a matrix filled with 0's and 1's selecting first half components of $\mathbf{x}_k$ so that $\mathbf{u}_k^L = \mathbf{F}^L \mathbf{x}_k$. Let

$$ \mathbf{C} = \mathbf{G}^L \mathbf{F}^L + \mathbf{G}^I \mathbf{F}^I. \qquad (A.11) $$

then the output equation

$$ \mathbf{y}_k = \mathbf{C}\mathbf{x}_k + \mathbf{n}_k \qquad (A.12) $$

relates $\mathbf{y}_k$ to the condensed state vector $\mathbf{x}_k$ satisfying the state equation in (5) involving the condensed transfer matrix $\mathbf{A}$.
## Appendix B. Impact of the number of unit cells on identification accuracy

To evaluate the impact of having a large number of unit cells for the identification, a noisy virtual experiment has been simulated. The unit cell considered here is taken from the numerical validation presented in Section 4. Harmonic out-of-plane velocity fields ($l = 259$ points per unit cell) for two periodic beams consisting of $N = 11$ unit cells (short beam with a length of 44 cm) and $N = 500$ unit cells (very long beam with a length of 20 m) have been computed using the WFEM [55]. This has been done in order to emulate the vibration measurements that have been obtained with a laser Doppler vibrometer in Section 5. White Gaussian noise is added to the spatial measurements at each frequency with a signal-to-noise ratio (SNR) of 10.

Subspace-based complex wavenumber identification has been done on the short beam with $N_o = 9$ unit cells. As recommended for short beams, the spatial horizon is chosen as $p = 2$. The number of waves to identify is $n_w = 4$. This configuration has been used for the numerical validation presented in Section 4. Identification results are presented in Figure B.13(a). These results are close to those presented in Section 4, which were based on FEM simulation, while the WFEM is used here in order to simulate a very long beam ($N = 500$). A FEM simulation of the long beam would involve nearly 1.9 million degrees of freedom (the unit cell contains 3786 degrees of freedom).

<figure><img src="image_14.png" /><figcaption>Figure B.13: Complex wavenumber identification results for the periodic aluminum beam from simulated noisy vibration measurement fields (SNR=10). Multi-output measurements are taken with $l = 259$ measurement points per unit cell and $N_o = 9$ unit cells considered for identification. Results for the proposed subspace method are compared in different configurations: (a) Low number of unit cells ($N = 11$): this condition has been tested experimentally. (b) Large number of unit cells ($N = 500$) with the same SNR: identification accuracy is improved for both imaginary and real parts of the wavenumbers.</figcaption></figure>

Wavenumber identification has been done on the very long beam with $N_o = 498$ unit cells. As for short beams, the spatial horizon was maintained in $p = 2$, also the number of waves to identify is $n_w = 4$. Note that the block Hankel matrix has a size of $518 \times 497$ at each frequency, which increases the computation time of the identification method. For the full frequency range, the identified dispersion curves were obtained in 1 m 36 s. Identification results are presented in Figure B.13(b). These results are already more accurate than those of Figure B.13(a), especially in the bandgap region. The two branches of the imaginary part of the wavenumber in the Bragg bandgap region are not superimposed, and this could be explained by the fact that these evanescent waves are not present along the full length of the 20 m waveguide, and therefore, their presence is significantly covered by noise in the middle region of the waveguide.

## Appendix C. Proof of Theorem 1

In the statement of Theorem 1, the matrix $A$ and the vector $x_k$ are respectively of size $2n \times 2n$ and $2n \times 1$. The shorter notation

$$q = 2n \qquad (C.1)$$

will be used in this proof.

The proof will be made by induction. More specifically, it will be shown that

* (Recurrence) if $z_l = \Gamma x_l$ holds for $l \le k - 1$, then so does $z_k = \Gamma x_k$;

* (Initialization) $z_k = \Gamma x_k$ for $k = 1, 2, \dots, q$.
Let us first show the Recurrence part, that is, if there exists a matrix $\Gamma \in \mathbb{C}^{q \times q}$ such that $\mathbf{z}_l = \Gamma \mathbf{x}_l$ holds for $l \le k-1$, then the same matrix $\Gamma$ satisfies also $\mathbf{z}_k = \Gamma \mathbf{x}_k$. This proof considers $k$ in the range $q+1 \le k \le N$, since, as indicated above, $\mathbf{z}_k = \Gamma \mathbf{x}_k$ for $k=1, 2, \dots, q$ will be proved directly in the Initialization part,

The essential step of this proof will show that $\mathbf{x}_k$ is a linear combination of its $q$ previous instances $\mathbf{x}_l$, and so is $\mathbf{z}_k$, i.e.,

$$ \mathbf{x}_k = \sum_{j=0}^{q-1} c_j \mathbf{x}_{k-q+j} \qquad (C.2) $$

$$ \mathbf{z}_k = \sum_{j=0}^{q-1} c_j \mathbf{z}_{k-q+j} \qquad (C.3) $$

with the same coefficients $c_j \in \mathbb{C}$ in the two linear combinations. Then the Recurrence part of the proof will be established by multiplying both sides of (C.2) with $\Gamma$. Before that step, three simple equalities are first evidenced and the Cayley-Hamilton theorem is recalled, as a preparation step.

Since the $2n$ eigenvalues of the matrix $\mathbf{A}$ consist of $n$ pairs such as $\lambda$ and $1/\lambda$, zero cannot be an eigenvalue of $\mathbf{A}$, which is therefore invertible. Then the reversed sequence $\mathbf{z}_k$ satisfies the backward recursion, i.e.,

$$ \mathbf{z}_{k+1} = \mathbf{A}^{-1} \mathbf{z}_k. \qquad (C.4) $$

For two integers $k_1$ and $k_2$ in appropriate ranges, repeatedly applying (20) yields

$$ \mathbf{x}_{k_1+k_2} = \mathbf{A}^{k_1} \mathbf{x}_{k_2}, \qquad (C.5) $$

and similarly, (C.4) leads to

$$ \mathbf{z}_{k_1+k_2} = \mathbf{A}^{-k_1} \mathbf{z}_{k_2}. \qquad (C.6) $$

To prove (C.2) and (C.3), let us recall the Cayley-Hamilton theorem: The characteristic polynomial $p(\lambda)$ of $\mathbf{A}$, which satisfies $p(\lambda) = 0$ for any of its eigenvalues $\lambda$, satisfies also $p(\mathbf{A}) = 0$. The size of the matrix being $q \times q$, its characteristic polynomial is of degree $q$. Then

$$ \mathbf{A}^q = c_0 \mathbf{I}_q + c_1 \mathbf{A} + c_2 \mathbf{A}^2 + \dots + c_{q-1} \mathbf{A}^{q-1} \qquad (C.7) $$

$$ = \sum_{j=0}^{q-1} c_j \mathbf{A}^j \qquad (C.8) $$

for some coefficients $c_j \in \mathbb{C}$.

Because the $2n$ eigenvalues of the matrix $\mathbf{A}$ consist of $n$ pairs such as $\lambda$ and $1/\lambda$, $\mathbf{A}$ and $\mathbf{A}^{-1}$ have the same set of eigenvalues, hence they have the same characteristic polynomial. Therefore,

$$ \mathbf{A}^{-q} = \sum_{j=0}^{q-1} c_j \mathbf{A}^{-j} \qquad (C.9) $$

holds with the same coefficients $c_j$ as in (C.8). This ends the preparation step.

The proof of (C.2) starts from equality (C.5), with $k_1 = q$ and $k_2 = k-q$, yielding

$$ \mathbf{x}_k = \mathbf{A}^q \mathbf{x}_{k-q}. \qquad (C.10) $$

In this equality, substitute $\mathbf{A}^q$ from (C.8), then

$$ \mathbf{x}_k = \left( \sum_{j=0}^{q-1} c_j \mathbf{A}^j \right) \mathbf{x}_{k-q} \qquad (C.11) $$

$$ = \sum_{j=0}^{q-1} c_j \mathbf{A}^j \mathbf{x}_{k-q} \qquad (C.12) $$

$$ = \sum_{j=0}^{q-1} c_j \mathbf{x}_{k-q+j} \qquad (C.13) $$
where the last equality follows (C.5) with $k_1 = j$ and $k_2 = k - q$. This result confirms the equality announced in (C.2). Then (C.3) is proved in the same way, starting from (C.6) with $k_1 = q$ and $k_2 = k - q$, and using (C.9) instead of (C.8).

Now left multiply both sides of (C.2) with $\Gamma$, yielding

$$ \Gamma \mathbf{x}_k = \sum_{j=0}^{q-1} c_j \Gamma \mathbf{x}_{k-q+j}. \qquad (C.14) $$

If $\mathbf{z}_l = \Gamma \mathbf{x}_l$ holds for $l \le k - 1$, then $\mathbf{z}_{k-q+j} = \Gamma \mathbf{x}_{k-q+j}$ at the right hand side of (C.14). Therefore,

$$ \Gamma \mathbf{x}_k = \sum_{j=0}^{q-1} c_j \mathbf{z}_{k-q+j} \qquad (C.15) $$

$$ = \mathbf{z}_k, \qquad (C.16) $$

where the last equality follows (C.3). It is then established at this stage that, if $\mathbf{z}_l = \Gamma \mathbf{x}_l$ holds for $l \le k - 1$, then $\mathbf{z}_k = \Gamma \mathbf{x}_k$.

This result is now completed by the Initialization part, i.e., $\mathbf{z}_l = \Gamma \mathbf{x}_l$ actually holds for $l \in \{1, 2, \dots, q\}$ and for some $\Gamma$, so that the equality is confirmed recursively for all for all $l \in \{1, 2, \dots, N\}$.

Let

$$ \Gamma = \begin{bmatrix} \mathbf{z}_1 & \mathbf{z}_2 & \cdots & \mathbf{z}_q \end{bmatrix} \begin{bmatrix} \mathbf{x}_1 & \mathbf{x}_2 & \cdots & \mathbf{x}_q \end{bmatrix}^{-1}. \qquad (C.17) $$

Rewrite this equation as

$$ \begin{bmatrix} \mathbf{z}_1 & \mathbf{z}_2 & \cdots & \mathbf{z}_q \end{bmatrix} = \Gamma \begin{bmatrix} \mathbf{x}_1 & \mathbf{x}_2 & \cdots & \mathbf{x}_q \end{bmatrix}. \qquad (C.18) $$

This matrix equality implies that the corresponding columns of both sides are equal, hence

$$ \mathbf{z}_1 = \Gamma \mathbf{x}_1 \qquad (C.19) $$

$$ \mathbf{z}_2 = \Gamma \mathbf{x}_2 \qquad (C.20) $$

$$ \vdots \qquad (C.21) $$

$$ \mathbf{z}_q = \Gamma \mathbf{x}_q. \qquad (C.22) $$

It is then confirmed that, with $\Gamma$ chosen as in (C.17), $\mathbf{z}_l = \Gamma \mathbf{x}_l$ indeed holds for $l \in \{1, 2, \dots, q\}$.

By combining the proved Recurrence and Initialization parts, Theorem 1 is then established by induction.

## Acknowledgments

The authors would like to thank the financial support of the Collège Doctoral de Bretagne as part of the outgoing international mobility program for carrying out this work.

## References

[1] S. Dalela, P. Balaji, D. Jena, A review on application of mechanical metamaterials for vibration control, Mechanics of advanced materials and structures 29 (22) (2022) 3237–3262.

[2] J. R. F. Arruda, Manipulating elastic waves aiming at reducing vibration and noise: bandgaps, metamaterials, and topological modes, in: ISMA 2022, 2022, pp. 1–26.

[3] D. Yu, Y. Liu, G. Wang, L. Cai, J. Qiu, Low frequency torsional vibration gaps in the shaft with locally resonant structures, Physics Letters A 348 (3) (2006) 410–415.

[4] X. Wu, L. Sun, S. Zuo, P. Liu, H. Huang, Vibration reduction of car body based on 2D dual-base locally resonant phononic crystal, Applied Acoustics 151 (2019) 1–9.

[5] W. H. Organization, Burden of Disease from Environmental Noise: Quantification of Healthy Life Years Lost in Europe, World Health Organization, Copenhagen, Denmark, 2011.
[6] J. Plisson, A. Pelat, F. Gautier, V. R. Garcia, T. Bourdon, Experimental evidence of absolute bandgaps in phononic crystal pipes, Applied Physics Letters 116 (20) (2020) 201902.

[7] A. Bergamini, M. Miniaci, T. Delpero, D. Tallarico, B. Van Damme, G. Hannema, I. Leibacher, A. Zemp, Tacticity in chiral phononic crystals, Nature communications 10 (1) (2019) 4525.

[8] D. Beli, J. Arruda, M. Ruzzene, Wave propagation in elastic metamaterial beams and plates with inter-connected resonators, International Journal of Solids and Structures 139 (2018) 105-120.

[9] M. Ichchou, J. Berthaut, M. Collet, Multi-mode wave propagation in ribbed plates. Part II: Predictions and comparisons, International Journal of Solids and Structures 45 (5) (2008) 1196-1216.

[10] L. Ribeiro, V. Dal Poggetto, B. Huallpa, J. Arruda, Bloch wavenumber identification of periodic structures using Prony's method, Mechanical Systems and Signal Processing 178 (2022) 109242.

[11] Q. Leclère, F. Ablitzer, C. Pézerat, Practical implementation of the corrected force analysis technique to identify the structural parameter and load distributions, Journal of Sound and Vibration 351 (2015) 106-118.

[12] A. Santoni, S. Schoenwald, B. Van Damme, P. Fausti, Determination of the elastic and stiffness characteristics of cross-laminated timber plates from flexural wave velocity measurements, Journal of Sound and Vibration 400 (2017) 387-401.

[13] T. Brion, P. Fossat, M. Ichchou, O. Bareille, A.-M. Zine, C. Droz, Identification of partial differential equations in structural mechanics theory through k-space analysis and design, Composite Structures 304 (2023) 116297.

[14] P. C. David Alleyne, A two-dimensional Fourier transform method for the measurement of propagating multimode signals, The Journal of the Acoustical Society of America 89 (3) (1991) 1159-1168.

[15] M. Ruzzene, Frequency-wavenumber domain filtering for improved damage visualization, Smart materials and structures 16 (6) (2007) 2116.

[16] T. E. Michaels, J. E. Michaels, M. Ruzzene, Frequency-wavenumber domain analysis of guided wavefields, Ultrasonics 51 (4) (2011) 452-466.

[17] T. J. Plona, B. K. Sinha, S. Kostek, S.-K. Chang, Axisymmetric wave propagation in fluid-loaded cylindrical shells. II: Theory versus experiment, The Journal of the Acoustical Society of America 92 (2) (1992) 1144-1155.

[18] K. Grosh, E. G. Williams, Complex wave-number decomposition of structural vibrations, The Journal of the Acoustical Society of America 93 (2) (1993) 836-848.

[19] R. Boukadia, C. Claeys, C. Droz, M. Ichchou, W. Desmet, E. Deckers, An INverse COnvolution MEthod for wavenumber extraction (INCOME): Formulations and applications, Journal of Sound and Vibration 520 (2022) 116586.

[20] J. Berthaut, M. Ichchou, L. Jezequel, K-space identification of apparent structural behaviour, Journal of Sound and Vibration 280 (3-5) (2005) 1125-1131.

[21] M. N. Ichchou, J. Berthaut, M. Collet, Multi-mode wave propagation in ribbed plates: Part I, wavenumber-space characteristics, International Journal of Solids and Structures 45 (5) (2008) 1179-1195.

[22] B. Van Damme, A. Zemp, Energy distribution and exchange between spatial harmonics in bending wave phononic crystals, Physical Review Applied 10 (1) (2018) 014001.

[23] B. Van Damme, A. Zemp, Measuring dispersion curves for bending waves in beams: a comparison of spatial Fourier transform and inhomogeneous wave correlation, Acta Acustica united with Acustica 104 (2) (2018) 228-234.

[24] L. Van Belle, C. Claeys, E. Deckers, W. Desmet, On the impact of damping on the dispersion curves of a locally resonant metamaterial: Modelling and experimental validation, Journal of Sound and Vibration 409 (2017) 1-23.

[25] J. G. McDaniel, W. S. Shepard Jr, Estimation of structural wave numbers from spatially sparse response measurements, The Journal of the Acoustical Society of America 108 (4) (2000) 1674-1682.
[26] L. Junyi, D. Balint, An inverse method to determine the dispersion curves of periodic structures based on wave superposition, Journal of Sound and Vibration 350 (2015) 41–72.

[27] L. Junyi, V. Ruffini, D. Balint, Measuring the band structures of periodic beams using the wave superposition method, Journal of Sound and Vibration 382 (2016) 158–178.

[28] F. Aloschi, O. S. Bursi, A. Palermo, A. Marzani, A time domain procedure for the identification of periodic structures, International Journal of Mechanical Sciences 269 (2024) 109054.

[29] X. Li, M. Ichchou, A. Zine, C. Droz, N. Bouhaddi, An Algebraic Wavenumber Identification (AWI) technique under stochastic conditions, Mechanical Systems and Signal Processing 188 (2023) 109983.

[30] X. Li, R. Rumpler, H. Mao, T. Brion, M. Ichchou, P. Göransson, Generalized Algebraic K-Space Identification technique for multidimensional signals: Application to wave and energy propagation characterization of curved structures, Mechanical Systems and Signal Processing 225 (2025) 112304.

[31] P. Margerit, A. Lebée, J.-F. Caron, X. Boutillon, High resolution wavenumber analysis (HRWA) for the mechanical characterisation of viscoelastic beams, Journal of Sound and Vibration (2018).

[32] P. Margerit, A. Lebée, J.-F. Caron, K. Ege, X. Boutillon, The high-resolution wavevector analysis for the characterization of the dynamic response of composite plates, Journal of Sound and Vibration 458 (2019) 177–196.

[33] P. Margerit, Extraction of dispersion surfaces: a review of applications, methods and open questions, in: ISMA-USD Noise and Vibration Engineering Conference 2022, 2022, pp. 1–12.

[34] R. Roy, T. Kailath, ESPRIT-estimation of signal parameters via rotational invariance techniques, IEEE Transactions on acoustics, speech, and signal processing 37 (7) (1989) 984–995.

[35] K. Ege, X. Boutillon, B. David, High-resolution modal analysis, Journal of sound and vibration 325 (4-5) (2009) 852–869.

[36] A. L. Swindlehurst, B. Ottersten, R. Roy, T. Kailath, Multiple invariance ESPRIT, IEEE Transactions on signal processing 40 (4) (1992) 867–881.

[37] V. F. Pisarenko, The retrieval of harmonics from a covariance function, Geophysical Journal International 33 (3) (1973) 347–366.

[38] R. Schmidt, Multiple emitter location and signal parameter estimation, IEEE transactions on antennas and propagation 34 (3) (1986) 276–280.

[39] R. Roy, A. Paulraj, T. Kailath, Estimation of signal parameters via rotational invariance techniques- ESPRIT, in: MILCOM 1986-IEEE Military Communications Conference: Communications-Computers: Teamed for the 90's, Vol. 3, IEEE, 1986, pp. 41–6.

[40] Y. Hua, T. K. Sarkar, Matrix pencil method for estimating parameters of exponentially damped/undamped sinusoids in noise, IEEE Transactions on Acoustics, Speech, and Signal Processing 38 (5) (1990) 814–824.

[41] P. Van Overschee, B. De Moor, Subspace algorithms for the stochastic identification problem, Automatica 29 (3) (1993) 649–660.

[42] A. Swindlehurst, R. Roy, B. Ottersten, T. Kailath, System identification via weighted subspace fitting, in: 1992 American Control Conference, IEEE, 1992, pp. 2158–2163.

[43] R. Badeau, B. David, G. Richard, Selecting the modeling order for the ESPRIT high resolution method: an alternative approach, in: 2004 IEEE international conference on acoustics, speech, and signal processing, Vol. 2, IEEE, 2004, pp. ii–1025.

[44] P. Mellinger, M. Döhler, L. Mevel, Variance estimation of modal parameters from output-only and in- put/output subspace-based system identification, Journal of Sound and Vibration 379 (2016) 1–27.

[45] S. Rothberg, M. Allen, P. Castellini, D. Di Maio, J. Dirckx, D. Ewins, B. Halkon, P. Muyshondt, N. Paone, T. Ryan, H. Steger, E. Tomasini, S. Vanlanduit, J. Vignola, An international review of laser doppler vibrometry: Making light work of vibration measurement, Optics and Lasers in Engineering 99 (2017) 11–22, laser Doppler vibrometry.

[46] M. C. Knauer, J. Kaminski, G. Hausler, Phase measuring deflectometry: a new approach to measure specular free-form surfaces, Optical Metrology in Production Engineering 5457 (2004) 366–376.
[47] D. J. Mead, A general theory of harmonic wave propagation in linear periodic systems with multiple coupling, Journal of Sound and Vibration 27 (2) (1973) 235-260.

[48] D. Duhamel, B. Mace, M. Brennan, Finite element analysis of the vibrations of waveguides and periodic structures, Journal of Sound and Vibration 294 (1) (2006) 205-220.

[49] L. Ljung, System Identification: Theory for the User, Prentice-Hall, 1999.

[50] A. Gavilán Rojas, Q. Zhang, C. Droz, A computationally efficient $k(\omega)$-spectral form for partial dispersion analyses within the wave finite element framework, Journal of Sound and Vibration 593 (2024) 118652.

[51] P. V. Overschee, B. D. Moor, Subspace Identification for Linear Systems: Theory – Implementation – Applications, Vol. 254 of Lecture Notes in Control and Information Sciences, Springer, 1996.

[52] H. Akaike, A new look at the statistical model identification, IEEE Transactions on Automatic Control 19 (6) (1974) 716-723.

[53] C. Geuzaine, J.-F. Remacle, Gmsh: A 3-d finite element mesh generator with built-in pre- and post-processing facilities, International Journal for Numerical Methods in Engineering (2009).

[54] C. Droz, J.-P. Lainé, M. Ichchou, G. Inquiété, A reduced formulation for the free-wave propagation analysis in composite structures, Composite Structures 113 (2014) 134-144.

[55] R. F. Boukadia, C. Droz, M. N. Ichchou, W. Desmet, A Bloch wave reduction scheme for ultrafast band diagram and dynamic response computation in periodic structures, Finite Elements in Analysis and Design 148 (2018) 1-12.

[56] P. O'Donoughue, F. Gautier, E. Meteyer, T. Durand-Texte, M. Secail-Geraud, F. Foucart, O. Robin, A. Berry, M. Melon, C. Pézerat, et al., Comparison of three full-field optical measurement techniques applied to vibration analysis, Scientific Reports 13 (1) (2023) 3261.

[57] P. L. Reu, D. P. Rohe, L. D. Jacobs, Comparison of DIC and LDV for practical vibration and modal measurements, Mechanical Systems and Signal Processing 86 (2017) 2-16, full-field, non-contact vibration measurement methods: comparisons and applications.

[58] C. Devivier, F. Pierron, P. Glynne-Jones, M. Hill, Time-resolved full-field imaging of ultrasonic Lamb waves using deflectometry, Experimental Mechanics 56 (3) (2016) 345-357.

[59] O. Robin, P. O'Donoughue, A. Berry, V. Farley, K. Prithipaul, Full field vibration measurements on a cantilever beam under impact using visible and infrared deflectometry, Applied Acoustics 183 (2021) 108294.