# Subspace-based wavenumber identification in periodic waveguides adapted to full-field vibration measurements

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

► To cite this version:

Alvaro Gavilán Rojas, Qinghua Zhang, Olivier Robin, Christophe Droz. Subspace-based wavenumber identification in periodic waveguides adapted to full-field vibration measurements. IOMAC 2025 - 11th International Operational Modal Analysis Conference, May 2025, Rennes, France. pp.1-7. hal-05204511

HAL Id: hal-05204511

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

Submitted on 30 Dec 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
# Subspace-based wavenumber identification in periodic waveguides adapted to full-field vibration measurements

Alvaro Gavilán Rojas <fnref n="1" />,<fnref n="2" />, Qinghua Zhang <fnref n="1" />, Olivier Robin <fnref n="2" />, Christophe Droz <fnref n="1" />

## ABSTRACT

In the context of wave propagation in periodic media (e.g., composite or architected materials), recent studies have underlined the relevance of subspace identification algorithms to identify wave propagation properties and, thus, to characterize a complex structure experimentally. While such algorithms are usually applied to data sampled in the discrete time domain, as in Operational Modal Analysis (OMA), here the data are collected in the frequency domain from successive periodic unit cells, in the case of 1D-periodic waveguides. Instead of the modal parameters estimated in OMA, such as natural frequencies, the aim here is to estimate the real and imaginary parts of the structural wavenumber (related to the wavelength and the spatial decay, respectively) and the Bloch wave modes (physical wave modes, such as torsional and compressional). These wave-related invariant parameters describe the vibrational behavior of the waveguide, since the displacement field can be represented by wave-mode superposition. The number of measurement points per unit cell is rapidly increasing with the development of full-field vibration measurement techniques. Taking advantage of this data abundance to counterbalance the usually low number of unit cells (a few periods) in the collected data is not of common practice in the community. This paper proposes a subspace identification framework to benefit from multiple points measurements inside each unit cell for periodic waveguides. The approach is illustrated with a simulated periodic beam.

<footnote n="1">¹ Univ. Gustave Eiffel, Inria, I4S team, Rennes, France
{alvaro.gavilan-rojas, qinghua.zhang, christophe.droz}@inria.fr</footnote>

<footnote n="2">² Centre de Recherche Acoustique-Signal-Humain, Université de Sherbrooke, Canada,
{alvaro.gavilan.rojas, olivier.robin}@usherbrooke.ca</footnote>
# 1. INTRODUCTION

Effective vibration control and attenuation of acoustic radiation are essential for managing Noise, Vibration, and Harshness (NVH) in various engineering applications, particularly in the automotive, aerospace, and industrial sectors. To this end, one of the means to control vibration is to act on the underlying elastic waves that propagate in the structure (the waveguide) and build up its vibratory motion. Controlling elastic waves can be effectively achieved by two main strategies: locally resonant materials (LRM) [11] or Bragg scattering (e.g., phononic crystals [9]). This second strategy takes advantage of Bragg bandgaps that occur in periodic structures, that are regions in the frequency domain where waves are strongly attenuated. In the conceptual design stage, the study of the behavior of these periodic structures mainly involves the generation of dispersion curves to examine the evolution of the complex wavenumber through the frequency domain. In the experimental validation stage, identification methods are applied to measured vibration data to extract the underlying wave properties [2, 5, 10].

Within the periodic waveguide context, the propagation constants and its related Bloch modes, representing the free wave motion in the waveguide, are analogous to the modal parameters (natural frequencies, damping ratios, and normal modes) in modal analysis. The propagation constant is expressed as $\lambda = \exp(jk\Delta) \in \mathbb{C}$, where $k \in \mathbb{C}$ is the complex wavenumber and $\Delta \in \mathbb{R}$ is the unit cell length. The real part of the wavenumber $\mathbb{R}(k)$ represents the harmonic motion – which is analogous to the natural frequency in modal analysis. Its imaginary part, $\mathbb{I}(k)$, is related to the spatial attenuation constant – which is analogous to the damping ratio in modal analysis. Finally, the Bloch modes represent the elastic wave modes (e.g., torsional, flexural, compressional). In the modal analysis field, the three modal parameters can be derived from a mechanical model by solving a generalized eigenvalue problem involving the mass, damping and stiffness finite element (FE) matrices of the structure. Similarly, the propagation constants and Bloch modes can be derived from a standard eigenvalue problem involving the transfer matrix [8], which is built from the FE matrices of the unit cell. This analogy suggests subspace methods, which are widely applied in modal analysis, for wave propagation problem. Successful applications of such a method to wavenumber identification in homogeneous structures have been reported in [7].

In this paper, data-driven subspace identification method is used to identify complex wavenumbers in periodic structures. This method takes advantage of the multiple measurement points, which can be obtained from a full-field measurement technique (e.g., laser-Doppler vibrometry), to compensate relatively few unit cells in most experiments to characterize periodic waveguides. Identification results for a virtual measurement example illustrate the benefit of this method for efficient estimation of the spatial attenuation in bandgaps.

# 2. SUBSPACE-BASED WAVENUMBER IDENTIFICATION

<figure><img src="image_2.png" /><figcaption>**Figure 1:** Unit cell and its partitioned degrees of freedom (left), and the periodic waveguide formed by assembling the unit cells (right).</figcaption></figure>

To introduce the approach of this paper for wavenumber identification from sensor measurements, we reformulate the wave propagation equations and their link to measurements in a state-space form. In practice, there is no need to know the details of these equations. Only their existence is assumed. The model of a periodic waveguide is based on the wave finite element method [8], hereafter referred to as WFEM, in which the dynamic motion equations of a periodic waveguide are expressed in terms of those of the smallest repetitive element, i.e., the $\Delta$-periodic unit cell. The time-harmonic motion equation at
the angular frequency $\omega$ of a periodic unit cell from a finite waveguide comprised of $N$ unit cells, after a finite element discretization, is written as:

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

where $\mathbb{D}(\omega) \in \mathbb{C}^{n_d \times n_d}$ represents the *dynamic stiffness matrix*, dependent on the angular frequency $\omega \in \mathbb{R}$ and formed with the classical finite element mass $\mathbf{M} \in \mathbb{R}^{n_d \times n_d}$, damping $\mathbf{C} \in \mathbb{R}^{n_d \times n_d}$ and stiffness $\mathbf{K} \in \mathbb{R}^{n_d \times n_d}$ matrices, $\mathbf{u}_k \in \mathbb{C}^{n_d}$ is the vector of nodal displacements and $\mathbf{f}_k \in \mathbb{C}^{n_d}$ is the vector of nodal forces, $n_d$ being the number of degrees of freedom of the finite element discretization. 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, as depicted in Figure 1. The corresponding unit cell vector and matrices are partitioned in the same manner.

When $N > 2$ unit cells are assembled in the periodic direction, as represented in Figure 1, the interface between the $k$-th and $(k + 1)$-th unit cells are coincident. Using exponents “$R$” and “$L$” to indicate respectively the nodes at the right and the left interfaces of each unit cell, we can write the displacement continuity 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}$. Moreover, if no external forces act in the inner region of the unit cell ($\mathbf{f}_k^I = \mathbf{0}$), the inner displacement vector of the $k$-th unit cell can be written – or condensed – in terms of the interface displacement vectors as $\mathbf{u}_k^I = -\mathbb{D}_{II}^{-1}(\mathbb{D}_{IL}\mathbf{u}_k^L + \mathbb{D}_{IR}\mathbf{u}_{k+1}^L)$. This gives us the *condensed equation of motion*:

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

Finally, by reordering (2) we obtain a state-space model in the following form:

$$ \left\{ \begin{matrix} \mathbf{x}_{k+1} = \mathbf{A}\mathbf{x}_k \\ \mathbf{y}_k = \mathbf{C}\mathbf{x}_k + \mathbf{n}_k \end{matrix} \right. \quad \text{with} \quad \mathbf{x}_k = \left\{ \begin{matrix} \mathbf{u}_k^L \\ \mathbf{f}_k^L \end{matrix} \right\} \quad (3) $$

where the state vector $\mathbf{x}_k \in \mathbb{C}^{2n}$ contains the nodal displacements and forces. The output vector $\mathbf{y}_k \in \mathbb{C}^l$ corresponds to the $l$ $\Delta$-periodic measurement points inside the $k$-th unit cell of harmonic displacements – i.e., the observations of the system in the frequency domain –, $\mathbf{C} \in \mathbb{C}^{l \times 2n}$ is the output matrix, and $\mathbf{n}_k \in \mathbb{C}^l$ is the measurement noise, which is assumed to be a zero-mean, stationary, white noise sequence. The matrix $\mathbf{A} \in \mathbb{C}^{2n \times 2n}$ is known as *transfer matrix* in the WFEM community. The eigenvalues of $\mathbf{A}$ are the propagation constants, which are related to the wavenumbers of interest in the considered context.

It is known that this state-space model is linked to subspace identification methods [6], which allow to estimate the eigenvalues of the matrix $\mathbf{A}$ from measurements corresponding to $y_k$. In this approach, the first step is to build a block-Hankel matrix $\mathcal{H} \in \mathbb{C}^{(pl) \times (N-p)}$ from the sequence of measurement vectors $\mathbf{y}_k \forall k = [1, N]$, as follows,

$$ \mathcal{H} = \begin{bmatrix} \mathbf{y}_1 & \mathbf{y}_2 & \cdots & \mathbf{y}_{N-p+1} \\ \mathbf{y}_2 & \mathbf{y}_3 & \cdots & \mathbf{y}_{N-p+2} \\ \vdots & \vdots & \ddots & \vdots \\ \mathbf{y}_p & \mathbf{y}_{p+1} & \cdots & \mathbf{y}_N \end{bmatrix}, \quad (4) $$
where $p \in \mathbb{N}$ denotes the spatial horizon of the observed data. This block-Hankel matrix can be decomposed as (see, e.g., [6]):

$$ \mathcal{H} = \mathcal{O}\mathcal{X} + \mathcal{N}, \quad \text{with} \quad \mathcal{O} = \begin{bmatrix} \mathbf{C} \\ \mathbf{C}\mathbf{A} \\ \vdots \\ \mathbf{C}\mathbf{A}^{p-1} \end{bmatrix} \quad \text{and} \quad \mathcal{X} = [\mathbf{x}_1 \ \mathbf{x}_2 \ \dots \ \mathbf{x}_{N-p+1}] \quad (5) $$

where $\mathcal{N}$ is the block-Hankel matrix of the measurement noise. The matrix $\mathcal{O}$ can be estimated from the singular value decomposition (SVD) of $\mathcal{H}$ as:

$$ \mathcal{H} = \mathbf{U}\mathbf{S}\mathbf{V}^{\mathrm{H}} \quad \text{and} \quad \hat{\mathcal{O}} = \mathbf{U}. \quad (6) $$

Let $\hat{\mathcal{O}}_\uparrow$ and $\hat{\mathcal{O}}_\downarrow$ be $\hat{\mathcal{O}}$ without the last block-row, and without the first block-row respectively (note that $p \ge 2$ is required). An estimate – up to some similarity transformation – of the state-transition matrix $\hat{\mathbf{A}}$ can be obtained from the relationship $\hat{\mathcal{O}}_\downarrow \hat{\mathbf{A}} \approx \hat{\mathcal{O}}_\downarrow$ by solving for $\hat{\mathbf{A}}$ in the least-squares sense:

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

The eigendecomposition of $\hat{\mathbf{A}} = \hat{\mathbf{T}}\hat{\mathbf{\Lambda}}\hat{\mathbf{T}}^{-1}$ provides the diagonal matrix of estimated propagation constants $\hat{\mathbf{\Lambda}}$, which can be converted to an estimated wavenumber matrix by using the relationship $\hat{\mathbf{\Lambda}} = \exp(j\hat{\mathbf{k}}\Delta)$.

### 3. NUMERICAL VALIDATION

To validate the proposed approach to complex wavenumber estimation, we present the results of the identification method for a virtual experiment consisting of harmonic motion simulations for a periodic waveguide – an aluminum beam with periodic material removal.

<figure><img src="image_3.png" /><figcaption>**Figure 2:** Periodic unit cell (left) and the corresponding periodic waveguide (right) including $N = 11$ unit cells. The unit cell geometric properties are a length $\Delta = 40$ mm, a width $w = 35$ mm, a thickness $t = 3.2$ mm, reduced thickness $t_r = 1.6$ mm (reduced thickness length is $3/4 \Delta$). The considered 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$.</figcaption></figure>

Firstly, a finite element model of the unit cell is built [1] by defining the geometry to create a 3D quadratic mesh [4] and the material properties to generate the FE matrices $\mathbf{M}$, $\mathbf{C}$ and $\mathbf{K}$. The unit cell is depicted and described in Figure 2.

To identify the complex wavenumbers of the unit cell, the dispersion curves are computed by solving a $k(\omega)$ generalized eigenvalue problem within the WFEM framework [3] (with $\omega = 2\pi f$ the angular frequency for an harmonic motion). These dispersion curves are shown in Figure 3 (real part $\Re(k)$ – upper part – and imaginary part $\Im(k)$ – lower part). A color scale superimposed on each curve indicates
<figure><img src="image_4.png" /><figcaption>**Figure 3:** Dispersion curves for the periodic waveguide (upper part: real part of the complex wavenumber $\Re(k)$– lower part: imaginary part $\Im(k)$. In-plane motion (compression and shear) get coupled with out-of-plane motion (bending and torsion) between 8 and 10 kHz.</figcaption></figure>

the degree of propagation of the wave mode by the modulus of the propagation constant – a wave mode is said to be propagative if $|\lambda| \approx 1$.

Aiming to generate the data for the identification problem, a periodic waveguide is assembled from the concatenation of $N = 11$ unit cells and depicted in Figure 2. The periodic waveguide FE matrices are created from the assembly of the unit cell FE matrices. Finite element and modal analysis uses a truncated modal basis comprised of the first 100 normal modes of the waveguide, with a maximum frequency of 21679 Hz. Frequency response functions are generated for a unit point load 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). Figure 4 depicts three vibration measurement fields to better understand how elastic waves exhibit a space-harmonic behavior that varies through the frequency domain.

<figure><img src="image_5.png" /><figcaption>**Figure 4:** Three simulated vibration measurement fields for the periodic waveguide in the frequency domain (mobility response functions) at 1.5 kHz (bottom), 5 kHz (middle) and 11 kHz (top).</figcaption></figure>
A spatial averaged mobility spectrum (point velocity / injected force) is illustrated in Figure 5, where the vibration amplitude level can be seen through the frequency range of interest [0, 12.8] kHz

<figure><img src="image_6.png" /><figcaption>**Figure 5:** Average mobility spectrum for the periodic waveguide. The torsional waves between 4.0 and 6.4 kHz are not observable from the vibration measurements, but the flexural modes continue to be propagative in this region, and they can be observed through the three vibration bending modes in this region.</figcaption></figure>

For the identification part, the number of outputs per unit cell is $l = 259$, which are the out-of-plane displacement degrees of freedom of the unit cell nodes after discretization (see Figure 1), the number of unit cells (the data length) is $N = 11$, and the spatial horizon is $p = 4$, so that the size of the block-Hankel matrix $\mathcal{H}$ is $1036 \times 7$. All the left singular vectors of $\mathcal{H}$ are used as an estimate of the observability matrix $\hat{\mathcal{O}}$. Accordingly, all the $\hat{\lambda}$ eigenvalues – the estimated propagation constants – of the estimated state-transition matrix $\hat{\mathcal{A}}$ are retained. Prior to the wavenumber estimation from $\hat{\lambda}$, a filtering criterion based on the modulus of $\hat{\lambda}$ is applied, which helps to find the wave modes with a physical meaning (4 wave modes evolving through the frequency domain among the 7 poles of $\hat{\mathcal{A}}$).

<figure><img src="image_7.png" /><figcaption>**Figure 6:** Identified wavenumbers for the periodic waveguide. Only out-of-plane modes are captured by the out-of-plane virtual vibration measurements (velocities along the z-axis).</figcaption></figure>

The results of the identification procedure are shown in Figure 6. The identified complex wavenumbers match properly the WFEM simulated dispersion curves. In the torsional bandgap region (4 to 6.4 kHz), the imaginary part of the wavenumber is well captured, which is usually difficult to achieve. In the region where in-plane modes are coupled with out-of-plane modes (8 to 10 kHz, the identification performance
is lower, which is due to the complex wave conversion phenomena and also to the decrease of the vi-
bration amplitude, which can be seen for the same region in the average mobility (see Figure 5). These
identification results identification have been obtained for a noiseless virtual measurement data, to let the
identification performance be affected by number of unit cells only.

4. CONCLUSIONS

In this paper, we presented an application of the subspace identification methods for complex wavenumber estimation in periodic waveguides. The periodic waveguide harmonic motion has been presented in its mechanical system formalism – finite element equations – but also in a state-space form, with a subspace identification perspective. Complex wavenumber identification in periodic waveguides can become a complex task, especially when bandgaps occur (i.e., limited or highly modified propagation). Through a numerical validation, the presented method proved to be adequate for the identification of the real but also imaginary parts of the wavenumber in these regions, which ultimately leads to an useful method for the assessment of vibration attenuation in a periodic waveguide.

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] J. Bezanson, A. Edelman, S. Karpinski, and V. B. Shah. Julia: A fresh approach to numerical computing. *SIAM review*, 59(1):65–98, 2017.

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

[3] R. F. Boukadia, C. Droz, M. N. Ichchou, and 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:1–12, 2018.

[4] C. Geuzaine and 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.

[5] X. Li, M. Ichchou, A. Zine, C. Droz, and N. Bouhaddi. An algebraic wavenumber identification (AWI) technique under stochastic conditions. *Mechanical Systems and Signal Processing*, 2023.

[6] L. Ljung. *System Identification: Theory for the User*. Prentice-Hall, 1987.

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

[8] D. J. Mead. A general theory of harmonic wave propagation in linear periodic systems with multiple coupling. *Journal of Sound and Vibration*, 27(2):235–260, Mar. 1973.

[9] J. Plisson, A. Pelat, F. Gautier, V. R. Garcia, and T. Bourdon. Experimental evidence of absolute bandgaps in phononic crystal pipes. *Applied Physics Letters*, 116(20):201902, 05 2020.

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

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