Estimating Dispersion Curves from Frequency Response Functions via Vector-Fitting

Estimating Dispersion Curv es from Frequenc y Response Functions via V ector -Fitting Mohammad I. Albakri a,1, ∗ , V ijaya V . N. Sriram Malladi b,1 , Serkan Gugercin d , Pablo A. T arazaga c a Department of Mechanical Engineering, T ennessee T echnological University b Department of Mechanical Engineering - Engineering Mec hanics, Michigan T echnological Univer sity c Department of Mechanical Engineering, V irginia P olytechnic Institute and State University d Department of Mathematics, V ir ginia P olytechnic Institute and State University Abstract Driven by the need for describing and under standing wave pr opagation in structural materials and com- ponents, sever al analytical, numerical, and experimental techniques have been developed to obtain dis- persion curves. Accurate char acterization of the structure (wave guide) under test is needed for analytical and numerical appr oaches. Experimental appr oac hes, on the other hand, r ely on analyzing waveforms as the y pr opagate along the structure . Material inhomogeneity , r eflections fr om boundaries, and the physi- cal dimensions of the structur e under test limit the fr equency r ange over whic h dispersion curves can be measur ed. In this work, a new data-driven modeling appr oach for estimating disper sion curves is de veloped. This appr oach utilizes the r elatively easy-to-measure , steady-state F r equency Response Functions ( FRF s ) to de velop a state-space dynamical model of the structur e under test. The de veloped model is then used to study the transient r esponse of the structur e and estimate its dispersion curves. This paper lays down the foundation of this appr oach and demonstrates its capabilities on a one-dimensional homogeneous beam using numerically calculated FRF s . Both in-plane and out-of-plane FRF s corr esponding , r espectively , to longitudinal (the first symmetric) and flexur al (the first anti-symmetric) wave modes ar e analyzed. The e ff ects of boundary conditions on the performance of this appr oach ar e also addr essed. K e ywords: data dri ven modeling, dispersion curv es, v ector fitting, spectral element method, wa v e propagation 1. Introduction Understanding wa ve propagation in structures is essential for numerous applications such as structural health monitoring, material characterization, stress-state identification [ 1 ], e vent detection and localiza- tion [ 2 – 4 ], vibration suppression, and elastic meta-structures [ 5 ]. K ey to this understanding is the ability to describe the frequency-dependent nature of wav e propagation characteristics gov erned by dispersion curves. Dispersion relations of elastic wa ves are determined by geometric and material characteristics of the wa ve guide, as well as the nature of the deformations induced by the propagating wav e (the wa v e mode) [ 6 ]. Analytical models hav e been dev eloped to deriv e dispersion relations for numerous materials and wa v eguides [ 6 ]. For commonly used homogeneous and heterogeneous materials, dispersion relations ∗ Address all correspondence to this author . Email: malbakri@vt.edu 1 These authors hav e contrib uted equally Pr eprint submitted to J anuary 1, 2020 hav e been documented in v arious handbooks and te xtbooks [ 7 , 8 ]. Numerical methods, such as time and frequency domain spectral element method [ 9 – 12 ], the semi-analytical finite element method [ 13 , 14 ], the wa v e finite element method [ 15 , 16 ], and transfer function and transfer matrix methods [ 17 , 18 ] ha v e also been dev eloped to describe wa ve propagation along uniform and periodically varying wav e guides. For realistic wa ve propagation characteristics to be obtained, accurate description of material and geo- metric characteristics of the wav e guide is required [ 19 ]. Accurate characterization of wa ve guides can be challenging, especially when inhomogeneous materials, such as wood and steel-reinforced concrete, and uncon ventional structures, such as 3D printed and smart structures with integrated sensors and actuators, are being considered. Experimental methods for measuring dispersion relations ha v e also been dev el- oped where a tone b urst is induced in the structure under test and then analyzed as it propag ates along the structure. Howe ver , material inhomogeneity , complex boundary conditions, and the physical dimensions of the structure under test pose challenges to such techniques. In particular , the ability of current tech- niques to measure the low frequency portions of dispersion curv es, which are v ery important for ev ent localization and stress-state identification applications, is limited by structure’ s dimensions, damping, and reflections from structure’ s boundaries. Data-dri ven modeling pro vides an alternati v e approach to create numerical models that capture the dynamic response of the structure under test. Such models rely completely on measured responses, thus, any kno wledge of the underlying physical characteristics of the structure is not required. In this context, the data amounts to a large set of frequency response measurements and the modeling amounts to constructing a low-dimensional rational function to fit this data in an appropriate sense. This is a heavily studied topic with varying approaches based on the e ventual goal; see, e.g., [ 20 – 23 ] for rational approximations that enforce interpolation of the measured data; [ 24 – 28 ] for rational approximations that minimize a least-squares de viation from the measured data; and [ 29 , 30 ] for rational approximations that combine the interpolation and least-squares formulation. Here, the least-squares based V ector Fitting [ 26 ] frame work will be emplo yed. In this work, a ne w data-driv en approach for estimating dispersion curves is dev eloped. Steady-state dynamic responses, in the form of Frequency Response Functions ( FRF s ), are used to build the data- dri ven model of the structure under consideration. The de v eloped model is then used to simulate the transient dynamic response of the structure which, in turn, is analyzed to reconstruct dispersion curves. The capabilities of this approach are demonstrated on a simple one-dimensional structure for which exact dispersion curves are attainable. The work b uilds on a pre vious e ff ort by the authors [ 31 ] and discusses in detail the dev elopment of this approach using numerical experiments. The applicability of this approach to symmetric (longitudinal) and anti-symmetric (fle xural) wav e modes is in v estigated. These wa v e modes are selected as they are the simplest to analyze theoretically and experimentally . Other w a v e modes along with wa ve propagation in two-dimensional structures, such as plates, will be addressed in future studies. This ne w approach is proposed as an alternativ e to the traditional e xperimental techniques, with the potential of overcoming the aforementioned limitations, especially when low-frequenc y portions of the dispersion curves are of interest. Furthermore, the fact that the de veloped approach relies on the steady- state dynamic response to estimate dispersion curves, as opposed to transient wav eform measurements, relaxes sampling rate requirements and impro v es signal-to-noise ratio in testing. This paper is structured as follo ws: Section 2 briefly presents the de velopment of the numerical model used to generate the FRF s to construct data-driv en models. Elementary rod and T imoshenko beam ap- proximate theories, along with the frequency domain spectral element method, are used to dev elop this model, which also pro vides the exact dispersion curves for comparison. Data-driven rational approxi- mation is discussed in Section 3 where the simulated FRF s are used to build a data-driv en, state-space model that describes the dynamics of the beam under test. The layout of the newly de veloped approach is discussed in Section 4 where dispersion curv es of the first anti-symmetric wa v e mode are estimated. The 2 applicability of this approach to estimate dispersion relations for symmetric wav e modes are presented in Section 5. Section 6 discusses the e ff ects of boundary conditions on the performance of this approach. Finally , concluding remarks are presented in Section 7. 2. Mathematical Model and Spectral Element Formulation Consider a long beam excited with a pair of identical piezoelectric actuators, as sho wn in Figure 1 . The piezoelectric actuators can be excited either in-phase to generate pure longitudinal deformations or out-of-phase for pure flexural excitation. This excitation configuration is selected here since it is com- monly used in experimental studies. Due to the structure symmetry around its neutral axis, longitudinal and flexural deformations are completely uncoupled. For simplicity , the piezoelectric actuators are re- placed with two pairs of longitudinal forces acting at the beam’ s upper and lower surfaces, as sho wn in Figure 1 . Readers are referred to [ 32 , 33 ] for detailed formulation of piezoelectric-augmented structures. P ZT 1 P ZT 2 Beam Neutral Axis 𝒙 , 𝒖 𝟎 𝝓 𝒛 , 𝒘 𝟎 𝑭 𝒙𝟏 𝑭 𝒙𝟐 − 𝑭 𝒙𝟏 − 𝑭 𝒙𝟐 𝑭 𝒙𝟏 𝑭 𝒙𝟐 − 𝑭 𝒙𝟏 − 𝑭 𝒙𝟐 Axial Ex cit ation Fle x ural Ex cit ation Figure 1: Schematic of a beam excited with surface mounted piezoelectric actuators (PZT). The mass and sti ff ness of the piezoelectric actuators are assumed to be negligible compared to those of the beam. Follo wing the Elementary rod and the T imoshenk o beam approximate theories, the elasto- dynamic equations of motion can be expressed as ρ A ∂ 2 u 0 ∂ t 2 − ∂ u 0 ∂ x E A ∂ u 0 ∂ x ! = q A , ρ A ∂ 2 w 0 ∂ t 2 − G A ¯ K ∂ ∂ x ∂ w 0 ∂ x − φ ! = P , ρ I ∂ 2 φ ∂ t 2 − E I ∂ 2 φ ∂ x 2 − G A ¯ K ∂ w 0 ∂ x − φ ! = 0 , (1) with the boundary conditions E A ∂ u 0 ∂ x = F x , G A ¯ K ∂ w 0 ∂ x − φ ! = F z , and E I ∂φ ∂ x = M , (2) where u 0 and w 0 denote the longitudinal and lateral displacements of the beam’ s neutral axis, respec- ti vely; φ denotes the angle of rotation of the neutral axis normal vector; ρ , E , and G are, respecti vely , the material’ s volumetric mass density , elasticity modulus, and rigidity modulus; A is the beam’ s cross sectional area; and I is its second moment of area. The T imoshenk o correction factor , ¯ K , is determined by matching the high-frequency wav e speed to that of Rayleigh wa ves [ 9 ]. q ( x , t ) is the externally applied 3 axial body force per unit volume and P ( x , t ) is the externally applied transverse distributed load. Both q and P are non-existent in the current analysis. The externally applied concentrated longitudinal forces, lateral forces, and bending moments are denoted by F x , F z , and M , respectiv ely . The frequenc y domain Spectral Element Method (SEM) is used to solve the elastodynamic equations of motion and to obtain FRF s corresponding to longitudinal and flexural deformations. This method is selected over other numerical methods, such as the FEM, due to its superior accuracy [ 11 ]. This is especially important when the high-frequency dynamic response is of interest. The spectral element formulation is briefly presented in this section. W ith the SEM, all variables appearing in the equations of motion, along with boundary conditions, are first transformed to frequenc y domain using the discrete F ourier transform, and a spectral solution of the follo wing form is assumed u ( x , t ) = 1 ˆ N ˆ N − 1 X n = 0 M X m = 1 r lm A m e − ˙ ı ( k mn x − ω n t ) , (3) where u ( x , t ) = [ u 0 ( x , t ) w 0 ( x , t ) φ ( x , t ) ] T is the response variables vector in time domain; ˆ N is the number of spectral components considered in the discrete Fourier transform, M is the number of wav e- modes contributing to the displacement-field, k mn is the wav e number corresponding to the m th mode at the n th angular frequency ( ω n ), and ˙ ı = √ − 1. Fourier coe ffi cients at ω n are represented in terms of the scaling matrix, r , and the amplitude vector A . The variable l takes the value of 1, 2, or 3 for u 0 , w 0 , and φ , respecti vely . Fourier coe ffi cients corresponding to the coupled variables, w 0 and φ , can be scaled relativ e to either one of them. In this w ork, the coe ffi cient corresponding to w 0 is chosen to be the scaling coe ffi cient, hence, r 2 m = 1 , ∀ m . The uncoupled axial deformations u 0 , on the other hand, can be solved for indepen- dently . Thus, no scaling factor is required for this component. Substituting the spectral solution, ( 3 ) into the gov erning equations ( 1 ) yields  − ρω 2 n + E k 2 mn  A m = 0 , " G A ¯ K k 2 mn − ρ A ω 2 n − i G A ¯ K k mn i G A ¯ K k mn E I k 2 mn − ρ I ω 2 n + G A ¯ K # " 1 r 3 m # A m = 0 , (4) For each frequency of interest, ω n , the characteristic equations are solved for the w a venumber k mn , which results in three wav emodes. These wav emodes represent the propagating first symmetric and first anti- symmetric modes, along with the e v anescent second anti-symmetric mode. The later wa v emode transfers to propagating upon passing through the cut-o ff frequency ω c = p G A ¯ K /ρ I . Group velocity for each propagating mode is then calculated as V Gm = ∂ω n /∂ k mn . These results are used in later sections to assess the performance of the data-driven approach. Once the characteristic equations are solved, the scaling constants, r 3 m , are calculated for each mode m . T o calculate FRF s corresponding to longitudinal and flexural deformations of the structure, the dy- namic sti ff ness matrix, K ( ˙ ıω ) is first e v aluated. This matrix relates the nodal displacements vector , d , and the nodal forces and moments vector , F , in the frequency domain. The dynamic sti ff ness matrix is defined in terms of the shape functions matrix, Ψ ( ˙ ıω ), and the boundary conditions matrix, G ( ˙ ıω ), as follo ws F = K (˙ ıω ) d = Ψ ( ˙ ıω ) − 1 G (˙ ıω ) d , (5) 4 where d , F , Ψ , and G for a two-node spectral finite element are defined in Appendix A. Although spatial discretization is not required for uniform, homogeneous structures, three spectral finite elements are used for the current analysis. This is done to allow the application of the desired excitation configuration described in Figure 1 . Standard assembly procedures are applied to obtain the global dynamic sti ff ness matrix. The beam considered here is a 48-in. long rectangular beam, with a 1 × 1 / 8 − in . 2 cross sectional area. The beam is made of Aluminum, with E = 69 G Pa , G = 26 G Pa and ρ = 2700 kg / m 3 . Piezoelectric wafers are bonded to the beam’ s upper and lower surfaces, 18 . 5 in . from its left end. A total of 23 receptance FRF s are calculated. These are the driving-point FRF (calculated at the right edge of the piezoelectric actuator) along with 22 transfer FRF s calculated o ver a span of 22 in., with 1-in. increments. FRF s are calculated ov er the frequency range of 0 to 50 k H z with 0 . 25 H z resolution. Sev eral boundary conditions are in v estigated in this study . The analysis is first presented for the free-free case. Other boundary conditions are discussed in Section 6 . 3. Data-driven Rational A ppr oximation This section briefly outlines the mathematical framework emplo yed to construct rational approxima- tions (state-space models) from a giv en set of FRF measurements. T o make the notation more clear and the discussion more concise, the methodology for only scalar measurements is presented, i.e., it is as- sumed that FRF measurements come from a single-input single-output dynamical system. Theory and numerical implementations are already established for multi-input multi-output systems, and the reader is referred to [ 26 , 34 ] for those details. Let H (˙ ıω ) ∈ C denote the FRF of the underlying dynamics ev aluated at the frequency ω . Then, the starting point is a set of FRF samples at selected frequencies, denoted by {H (˙ ıω 1 ) , H ( ˙ ıω 2 ) , . . . , H (˙ ıω N ) } where N is a positi ve integer . In this paper where the theoretical framework is established, these FRF s are obtained by sampling the analytical transfer functions deri v ed from Eq. ( 5 ). More specifically , the components of the in v erse of the dynamic-sti ff ness matrix, i.e., K − 1 (˙ ıω ), corresponding to the out-of- plane and in-plane displacements are sampled ov er the frequency range of interest. In practice, these samples will be obtained via experimental measurements. Gi ven the samples, {H (˙ ıω j ) } for j = 1 , 2 , . . . , N , the goal is, now , to construct a degree- r dynamical system described by its degree- r rational function ( FRF ), e H ( s ) = e  ( s I − e  ) − 1 e  , that fits the data in an appropriate sense. The degree- r means that e  ∈ R r × r , e  ∈ R r × 1 , and e  ∈ R 1 × r . In this work, the least-squares ( LS ) measure is used to judge the quality of the data-dri ven model: Find the degree- r rational function e H ( s ) that minimizes N X j = 1     H ( ˙ ıω j ) − e H ( ˙ ıω j )     2 . (6) The V ector Fitting ( VF ) method of Gustavsen et al. [ 26 ] is used (and adopted to our setting) in solving the LS problem ( 6 ). This approach is briefly explained in this section. Let n ( s ) = α 0 + α 1 s + · · · + α r − 1 s r − 1 and d ( s ) = β 0 + β 1 s + · · · + β r − 1 s r − 1 + s r denote, respecti vely , the numerator and denominator of the degree- r rational function e H ( s ), i.e., e H ( s ) = n ( s ) d ( s ) = α 0 + α 1 s + · · · + α r − 1 s r − 1 β 0 + β 1 s + · · · + β r − 1 s r − 1 + s r . The numerator n ( s ) has degree- r − 1 (as opposed to degree- r ) since H ( s ) does not have a direct feed through term; thus, it is a strictly proper rational function. The formulation can be easily modified 5 to allo w a degree- r numerator; but for the problems presented herein, H ( s ) is strictly proper . No w notice that the unknowns appear both in n ( s ) and d ( s ); thus the minimization problem ( 6 ) is a nonlinear LS problem and cannot be solved with direct methods in one step, as is the case for linear LS problems. An iterati v e scheme is needed: the nonlinear LS problem ( 6 ) is solved iterativ ely by solving a sequence of linear LS problems. This is achie ved by linearizing the error function H ( s ) − e H ( s ). Using e H ( s ) = n ( s ) d ( s ) , the LS error in ( 6 ) can be re written as N X j = 1     H ( ˙ ıω j ) − e H ( ˙ ıω j )     2 = N X j = 1    d (˙ ıω j ) H ( ˙ ıω j ) − n (˙ ıω j )    2    d (˙ ıω j )    2 , (7) which is still nonlinear in the v ariables. Starting with the d (0) ( s ) ≡ 1, Sanathanan and Koerner [ 35 ] proposed an iterati ve scheme where at the k th step the LS error N X k = 1       n ( k + 1) ( ξ i ) − d ( k + 1) ( ξ i ) H ( ξ i ) d ( k ) ( ξ i )       2 is minimized by solving for n ( k + 1) ( s ) and d ( k + 1) ( s ) . (8) Note the di ff erence from the earlier problem. The relaxed LS problem in ( 8 ) is linear in the v ariables n ( k + 1) ( s ) and d ( k + 1) ( s ), and can be solved using linear LS solution techniques. Then, the process is repeated using the ne w denominator d ( k + 1) ( s ) until con v er gence is achiev ed. This is the SK iteration to solve the original nonlinear LS problem ( 6 ). The VF method is a further improv ement on the SK iteration and uses a di ff erent parameterization for e H ( s ); specifically , VF uses the barycentric-form for a rational function: At the k th step of the VF iteration, the approximation is parameterized as e H ( k ) ( s ) = ˜ n ( k ) ( s ) ˜ d ( k ) ( s ) = P r j = 1 φ ( k ) j / ( s − λ ( k ) j ) 1 + P r j = 1 ψ ( k ) j / ( s − λ ( k ) j ) , where φ ( k ) j , ψ ( k ) j , λ ( k ) j ∈ C . (9) The adv antage of this formulation is that the poles λ ( k ) j are an arbitrary set of mutually distinct points. Note that λ ( k ) j are not the poles of e H ( k ) ( s ) unless ψ ( k ) j = 0 for j = 1 , 2 , . . . , r . As in the SK iteration, for fixed poles λ ( k ) j , the LS problem for e H ( k ) ( s ) in ( 9 ) becomes linear in the remaining variables φ ( k ) j and ψ ( k ) j . Then, the VF algorithm continues the iteration by updating the poles λ ( k ) j as the zeros of the numerator in Eq. ( 9 ), i.e., as the zeroes of ˜ d ( k ) ( s ) = 1 + P r j = 1 ψ ( k ) j / ( s − λ ( k ) j ). This iteration is run until con v er gence of the poles λ ( k ) j , upon which the denominator ˜ d ( k ) ( s ) → 1 and the final rational approximant is giv en by e H ( s ) = r X j = 1 φ j s − λ j = e  ( s I − e  ) − 1 e  . (10) For further details, we refer the reader to [ 26 , 36 ]. VF has successfully been used in many applications and various modifications hav e been proposed. The pole relocation step has been improv ed [ 37 ], the MIMO case has been e ffi ciently parallelized [ 38 ], it has been combined with quadrature-based sampling [ 27 ], and it has been analyzed in detail in terms of numerical robustness [ 34 ]. The reader is referred to [ 28 ] where a rational Krylov toolbox has been dev eloped to solve similar nonlinear LS problems. The Loe wner framework, introduced by Mayo and Antoulas [ 22 ], is another commonly used approach to construct rational approximants from FRF measurements. In this case, the approximant is sought to be an (approximant) interpolant as opposed to a LS fit. The recently dev eloped Adaptiv e Anderson-Antoulas 6 algorithm [ 29 ] is a hybrid approach where a rational interpolant is constructed to interpolate a subset of the data and to minimize LS error in the rest. As pointed out earlier , in this w ork, the focus is on the LS frame work where the re gular VF framework is used to solv e the data-dri v en modeling problem. 4. Data-driven Model based Dispersion Curv es f or Flexural W aves The proposed approach for estimating dispersion curves utilizes FRF s to generate a single-input multi- output (SIMO), data-driv en, state-space model, e H ( s ), of the structure under consideration. This model is then used to simulate the transient dynamic response of the structure to a gi v en excitation. Simulated responses are then analyzed in the frequency domain to reconstruct dispersion curves. The analysis is first presented for out-of-plane (fle xural) receptance FRF s , corresponding to the first anti-symmetric wa v e mode. The analysis is repeated for the first symmetric wa v e mode, using the in-plane (longitudinal) FRF s , in Section 5 . The process of estimating dispersion curves from FRF s can be split into two main stages: (i) the de velopment of the data-driv en model using VF and (ii) transient response simulation and analysis. The two stages are discussed in detail in the follo wing subsections. 4.1. Data-driven Modeling using VF In this section, a SIMO, data-driv en, state-space model, e H ( s ), is de veloped based on the simulated FRF s for the beam ov er the frequenc y range of 0 . 01 to 50 k H z . FRF s are sampled such that the frequenc y resolution is 0 . 25 H z over the entire frequency range of 50 k H z . Such fine frequency resolution is needed to better capture lo wer-order resonant frequencies, resulting in approximately 200,000 samples. As dis- cussed in detail belo w , an accurate e H ( s ) for the full frequency spectrum requires an approximation order of r > 200. Since we are measuring 23 FRF s , this means that e very VF iteration requires solving a linear least-squares problem with a dense coe ffi cient matrix ha ving more than 2 . 5 × 10 6 ro ws and 4800 columns. Since the coe ffi cient matrix changes at ev ery step, this is not a trivial numerical task and needs careful numerical implementation such as those in [ 34 , 38 ]. Ho we v er , the main issue in trying to fit the FRF s directly on the full frequency spectrum is the sensiti v- ity to the initial pole selection. As discussed in Section 3 , the VF algorithm is initialized by a set of poles. While for modest r values and for a smaller number of measurements, the VF algorithm does a good job of correcting the poles during the iteration, for the scale of the problems and for the complex dynam- ics considered here, a not-well-informed initial pole selection leads to inaccurate rational approximants together with potentially highly ill-conditioned least-squares problems to solve; in other words, for the setting considered in this paper , the quality of the fitted model is more dependent on the choice of starting poles; see [ 26 , 28 , 37 ] for related work. Moreover , when experimentally measured noisy FRF s are used, these issues will be further magnified. Therefore, in order to dev elop a frame w ork that is much better suited for experimentally collected data, a more comprehensive initialization procedure is de veloped for VF . For this purpose, the full frequency range is divided into multiple smaller bands. Since high-quality rational approximations in smaller frequency bands can be obtained via smaller approximation orders, for each frequency band, the initial, b ut now much smaller number of poles for VF , are selected by taking into account the quantity and the location of resonant frequencies. Additionally , due to the free-free boundary conditions, rigid-body modes appear in the FRF s at 0 H z . While the VF algorithm can handle fitting these modes, numerical instabilities arising from the 0 H z pole renders simulation results inaccurate. T o av oid such numerical instabilities, rigid-body modes are eliminated in the current analysis. Once these smaller frequency bands are fitted, the resulting poles are put together to initialize the VF algorithm for the entire frequency range. W ith this much improved initial pole selection, VF con ver ges quickly for the whole frequenc y range and yields an accurate rational approximation, as illustrated in the numerical examples belo w . The follo wing steps further discuss the details of this process. 7 Step 1: VF of FRF Partitioned into Multiple Bands Since for dispersiv e media, such as beams and plates, the frequenc y spacing between resonant frequencies is not uniform, random assignment of initial poles would not be practical. For instance, in 1-dimensional dispersi ve structures, the frequency spacing between consecutiv e natural frequencies increases with mode number . As a result, modal density is higher at lower frequency ranges, i.e., there are more resonant peaks in the 0 − 1 k H z range than, for example, the 2 − 3 k H z range. This becomes e ven more important when dealing with noisy experimental data. As a result, the VF algorithm is initialized using poles that are obtained by fitting FRF s ov er smaller frequenc y bands. When fitting out-of-plane FRF s , the full frequency range is divided into smaller bands in order to achie ve high-quality rational approximations. Seven frequency bands, with comparable modal density , were selected for this purpose. Details of the number of resonant peaks and the corresponding number of poles used to fit FRF in each band is tabulated in T able 1 . Each resonant peak of the FRF is represented by a pair of complex-conjugate poles, where the n th pair of poles are of the form p 2 n − 1 = α + i β , p 2 n = α − i β . Therefore, the number of poles of the fitted model are at least twice the number of resonant frequencies in the frequency band of interest. While the imaginary part ( β ) of these poles indicates the oscillatory behavior , the real part ( α ) represents the damping characteristic of the dynamical system. Resonant peak frequencies are used as the initial v alues for β , whereas the real part is initialized as α = β/ 100, follo wing the recommendation of [ 26 ]. The e ff ect of initial pole placement and model order on the performance of this approach will be further in vestig ated in future studies. T able 1 summarizes fitting results for each frequency band. The quality of the fit is determined using a relativ e L 2 error defined as E rel L 2 = 1 N v u u u u u u u u u u u u u u u t N X i = 1     H ( ˙ ıω i ) − e H ( ˙ ıω i )     2 F N X i = 1 k H ( ˙ ıω i ) k 2 F . (11) where k · k F denotes the Frobenious norm. T able 1 demonstrates that the constructed models fit the data to a high relati ve accurac y . While the minimum number of poles required to fit each frequency band is twice the number of resonant peaks in that band, adding extra poles to account for near by resonant peaks, that hav e strong influence on the FRF s in that band, was found to be necessary . For example, at least 28 complex poles were needed to fit the FRF within the frequency band of 0 - 1 k H z , which only has 13 resonant peaks. A closer examination of the fitted poles shows that the two additional poles correspond to a resonant peak at 1038 H z , which is very close to the current frequency band of interest. T able 1 sho ws the number of resonant peaks in each frequency band along with the number of poles needed to fit FRF s ov er each band. Step 2: VF of the Entir e FRF Once the poles of the state-space models for each frequency band are determined, they are combined together to initialize VF over the entire frequency range. While combining these poles, it is important to av oid duplicating them. The complete FRF with 106 resonant peaks is fitted with 212 poles, exactly the minimum number of poles required to fit all resonant peaks ov er this frequency range. The result is a data-dri ven, state-space model of the form: e H ( s ) = e  ( s I − e  ) − 1 e  , (12) 8 T able 1: Number of resonant peaks and initial poles in partitioned frequency bands Band Frequency range Resonant peaks No. of poles E rel L 2 1 0 . 01 - 1 k H z 13 28 1 . 02 × 10 − 5 2 1 - 5 k H z 18 40 3 . 78 × 10 − 6 3 5 - 10 k H z 13 32 2 . 96 × 10 − 6 4 10 - 20 k H z 21 44 1 . 48 × 10 − 6 5 20 - 30 k H z 16 38 1 . 45 × 10 − 6 6 30 - 40 k H z 13 30 1 . 44 × 10 − 6 7 40 - 50 k H z 12 26 1 . 42 × 10 − 6 where the state matrices are e  ∈ I R 212 × 212 , e  ∈ I R 212 × 1 , and e  ∈ I R 23 × 212 . Figure 2 depicts the magni- tude and phase of the fitted FRF ( e H ( ˙ ıω )) as compared to the out-of-plane component of the FRF matrix ( K − 1 ( ˙ ıω )), deri ved from Eq. ( 5 ) using the SEM, for the frequency band of 30 k H z to 40 k H z . Resonant peaks and anti-resonant valle ys are accurately captured by the fitted model. It should be noted that an in verse weighting scheme is adopted with the least square fit during the VF process. Thus, the VF algo- rithm places more importance on fitting the anti resonances than the resonant peaks. This is particularly important as the intended application of the data-driv en model is to simulate the transient response of the system over the entire frequency range. The relative L 2 error , E rel L 2 , for indi vidual frequency bands is provided in the T able 1 and the error for the complete frequency range is E rel L 2 = 2 . 13 × 10 − 7 . 4.2. T ransient Response Analysis and Dispersion Curves Estimation Once the data-driven, state-space model e H ( s ) in Eq. ( 12 ) is dev eloped, the transient response of the structure under test can be numerically approximated by exciting e H ( s ) with tone-burst input signals. Dispersion curves are then estimated by analyzing the wa v eforms propagating along the structure. Details of the this process are further presented through the follo wing steps: Step 1: T ransient Response Simulations The transient response of the structure under test is simulated using the de veloped data-driv en model. An amplitude-modulated, sine wav e, tone burst is selected as the excitation wav eform. This wa veform is selected as it is found to minimize the separation distance between the measurement location and structure’ s boundaries required for obtaining a reflection-free response [ 39 ]. In this study , the number of cycles is selected to be 2 when the strongly dispersi ve flexural wa ves are being considered and 1 for the weakly dispersi ve longitudinal w a ves. The data-dri ven model is then used to calculate the response of the structure at all 22 points on the beam where FRF s are originally obtained. Figure 3 .a sho ws the simulated response 10 in . away from the e xcitation point as a 2-c ycle, sine wa v e, tone b urst of a central frequency of 5 k H z is applied. As noted in the figure, due to the finite dimensions of the beam, reflections from boundaries is present in the simulated response. Structure’ s boundary conditions, free-free in this case, determine how the incident wa v e is reflected. This information is carried by the FRF s that hav e been used to obtain the data-driv en model. Boundary condition e ff ects on estimating dispersion curves are further discussed in Section 6 . Step 2: Incident W a vef orm Extraction and Processing Once the transient response is simulated at a gi v en location, the incident wa v eform is e xtracted. This is performed by determining the first dominant peak in the simulated response, defining a time-windo w 9 30 32 34 36 38 40 Frequency (kHz) 10 -10 10 -9 10 -8 10 -7 Mag.- Disp./Force (m/N) SEM Vector Fit 30 32 34 36 38 40 Frequency (kHz) 0 50 100 150 Phase (Degrees) Figure 2: Out-of-plane receptance FRF s of the beam obtained by the data-driv en state space model compared to the SEM predictions ov er the frequency range of 30 kH z to 40 k H z . FRF s are calculated 1 in . away from the e xcitation point. containing the incident wa ve, and then applying a strongly decaying exponential functions outside that windo w . The size of this window is determined based on the number of cycles in the excitation signal, excitation frequency , and the location where the response is simulated. Figure 3 .a sho ws the e xtracted incident wav e, labeled Pr ocessed Response as compared to the original simulated response, 10 in . away from the e xcitation point. Figure 3 .b sho ws the processed responses 2 in . and 10 in . a way from the excitation point. The dispersi ve nature of the flexural wav es at this frequency range is e vident in the figure where wa v eform distortion and spread can be noticed as the wa v e propagates along the beam. 10 0 0.5 1 1.5 2 2.5 Time (ms) -1 -0.5 0 0.5 1 Normalized Amplitude Simulated Response Processed Response (a) 0 0.5 1 1.5 Time (ms) -1 -0.5 0 0.5 1 Normalized Amplitude Processed, x = 2 in. Processed, x = 10 in. (b) Figure 3: Fle xural w av eforms showing (a) simulated and processed responses 10 in . aw ay from the excitation point and (b) processed responses 2 in . and 10 in . a way from the excitation point. Responses are simulated with a 2-cycle sine wave tone burst e xcitation signal, with 5 k H z central frequency . Step 3: Dispersion Curves Estimation Gi ven the Pr ocessed Response at locations i and i + 1, the w av enumber , k , corresponding to the w a ve mode under consideration can be calculated as follo ws U i + 1 ( ω ) = U i ( ω ) e − i k ( x i + 1 − x i ) , (13) where U i is the vector of Fourier coe ffi cients of the signal measured at location x i . At the frequency range of interest, the only wa ve mode contributing to the out-of-plane receptance FRF s is the first anti- symmetric mode. Thus, the wa v enumber vector , k , only includes that wa v e mode. Equation ( 13 ) is solv ed for k and the group velocity is then calculated as V G = ∂ω/∂ k . The process is then repeated where the central frequenc y of the excitation wa v e is v aried to sweep the frequenc y range of interest. Equation ( 13 ) requires the response at two distinct points to be known. W ith the response being simulated at 23 points in this study , 253 di ff erent combinations can be obtained, and hence, 253 estimates of the dispersion curves can be calculated. This allows for statistical analysis to be conducted and confidence bands to be defined for the estimated dispersion curve. Such analysis is important when dealing with experimental data where noise and measurement errors can be present. Simulated responses, on the other hand, do not su ff er from such sources of error , thus statistical analysis has not been included in this study . Figure 4 depicts the group velocity of the flexural (the first anti-symmetric) wav e mode calculated using the proposed data-driv en modeling approach as compared to that predicted by the SEM. As depicted in the figure, the results of the data-dri ven modeling approach accord very well with the SEM predictions. It should be noted that at high frequency , greater than 45 k H z , the data-driv en model fails to predict the dispersion curve. This can be ascribed to the fact that the FRF s that have been used to obtain the data- dri ven model o v er the frequency range of 0 − 50 k H z . Since a 2 − cycle sinusoidal w a veform is relati v ely broad-band, high-frequency signals of this form ha ve frequency content that extends well beyond the 50 k H z limit of the model. Increasing the number of cycles in the excitation tone burst reduces the bandwidth of the e xcitation signal. This allo ws dispersion curv es to be correctly estimated at frequencies closer to the measurement ceiling. Ho we v er , the model remains limited by the frequency range over which FRF s ha ve been obtained. 11 0 10 20 30 40 50 Frequency (kHz) 0 0.5 1 1.5 2 2.5 Group Velocity (km/s) Spectral Element Method Data-Driven Approach Figure 4: Comparison of group velocity curves estimated using the proposed data-dri ven approach (red) and that predicted by the SEM (blue) for the first anti-symmetric (flexural) wa v e mode. Close e xamination of Figure 4 also re veals that the very low-frequenc y part of the dispersion curv e is missing. This is due to the large wav elength of flexural wav es at such low frequencies compared to the length of the beam under test. This, in turn, hinders the separation of incident and reflected wa veforms, a prerequisite for Eq.( 13 ) to be applicable. While this is an inescapable limitation for con ventional wav e- propagation-based experimental techniques, the use of the proposed data-driv en approach allo ws for a number of solutions to be applied. One such solution can be to introduce artificial damping to the data- dri ven model in order to attenuate reflected w a ves. This will be the focus of future studies. The results presented in this section highlight the capabilities of the proposed approach where a model deri v ed from the steady-state dynamic response of the structure is emplo yed to accurately capture its transient response and estimate its dispersion curves. The following section extends this discussion to the longitudinal wa v e mode. 5. Data-driven Modeling based Dispersion Curv es f or Longitudinal W aves In this section, the applicability of the proposed approach to calculate the dispersion curves of the first symmetric (longitudinal) wav e mode is in vestig ated. The process presented in this section closely follo ws the one discussed in Section 4 , the main di ff erence is that in-plane receptance FRF s of the beam are utilized to obtain the data-dri ven model. In-plane FRF s are obtained by exciting the upper and lo wer faces of the beam in-phase, as shown in Figure 1 . Since the beam under test is perfectly symmetric with respect to its neutral axis, in-plane and out-of-plane deformation are completely uncoupled. Thus, the aforementioned excitation configuration results in pure longitudinal deformations without any contribution from flexural deformations. Figure 5 sho ws the transfer , longitudinal, receptance FRF , 6 − in a way from the excitation location, as obtained by the SEM. As is the case for the out-of-plane deformations, the analysis is limited to the frequency range of 0 − 50 k H z . Follo wing the procedure outlined in Section 4.1 , VF is utilized to fit the simulated longitudinal FRF s of the beam. In the frequency range of interest, the 24 resonant peaks are fitted with 48 complex-conjugate poles. The nondispersiv e nature of the longitudinal wa v es is e vident from the harmonic nature of the resonant peaks. Gi v en the relativ ely small number of resonant peaks in this frequency range, FRF parti- tioning was unnecessary in this case. T able 2 presents the order of the fitted state-space matrices and the 12 10 20 30 40 50 Frequency (kHz) 10 -10 10 -5 Mag.- Disp./Force (m/N) SEM Vector Fit 10 20 30 40 50 Frequency (kHz) 0 50 100 150 Phase (Degrees) Figure 5: In-plane receptance FRF of the beam obtained by the data-driven state space model compared to the SEM predic- tions. FRF s are calculated 6 in . away from the e xcitation point o ver the frequenc y range of 0 − 50 k H z . corresponding E rel L 2 error between the fitted FRF s and the simulated ones. Additionally , the quality of the fit is illustrated by comparing the fitted FRF to the SEM prediction as sho wn in Figure 5 . The data-driv en state-space model is then used to simulate the transient response of the system. As the longitudinal wa v e mode is nondispersiv e at this frequency range, a 1-cycle sine-w a ve tone burst is used as the e xcitation signal. Incident wav eforms are then extracted from the simulated response. Figure 6 .a sho ws the extracted incident wa v e, labeled Pr ocessed Response as compared to the original simulated response 2 in . aw ay from the excitation point. Figure 6 .b shows the processed responses 2 in . and 12 in . aw ay from the excitation point. The nondispersi ve nature of this wa ve mode is seen in the figure where wa v eforms propagate along the beam without noticeable distortion. While simple time-of-flight calculations can be used to estimate wav e speed of nondispersiv e wav e modes, the frequenc y-domain analysis presented in the pre vious section is follo wed. W ith the frequency- domain analysis, a given w a veform can be used to estimate wa v e speed over a relati v ely wide frequency band, as opposed to a single frequency estimate with the time-of-flight calculations. Thus, statistical techniques can be implemented to define confidence intervals for the estimated dispersion curves, which is important when dealing with experimental measurements. Follo wing Eq. ( 13 ), the processed longitu- dinal wa veforms at locations i and i + 1 are analyzed and the wav enumber , k , corresponding to the first symmetric wa ve mode is obtained for each frequency in the excitation signal. Group velocity is then calculated and the results are shown in Figure 7 . As depicted in the figure, the group velocity calcu- lated using the data-driv en approach accord very well with the SEM predictions. textcolorredAt high 13 T able 2: Fitting results and errors for the v arious data-dri ven models used in this study . Boundary Conditions State Matrices E rel L 2 Free - Free e  ∈ I R 212 × 212 , e  ∈ I R 212 × 1 , e  ∈ I R 23 × 212 2 . 13 × 10 − 7 Clamped - Free e  ∈ I R 214 × 214 , e  ∈ I R 214 × 1 , e  ∈ I R 16 × 214 2 . 75 × 10 − 7 Pinned -Pinned e  ∈ I R 212 × 212 , e  ∈ I R 212 × 1 , e  ∈ I R 16 × 212 3 . 09 × 10 − 7 Free - Free (inplane) e  ∈ I R 48 × 48 , e  ∈ I R 48 × 1 , e  ∈ I R 23 × 48 3 . 12 × 10 − 7 0 0.1 0.2 0.3 Time (ms) -1 -0.5 0 0.5 1 Normalized Amplitude Simulated Response Processed Response (a) 0 0.05 0.1 0.15 0.2 Time (ms) -1 -0.5 0 0.5 1 Normalized Amplitude Processed, x = 2 in. Processed, x = 12 in. (b) Figure 6: Longitudinal w av eforms showing (a) simulated and processed responses 2 in . away from the e xcitation point, and (b) processed responses 2 in . and 12 in . away from the excitation point. Responses are simulated with a 1-cycle sine wav e tone burst e xcitation signal, with 20 k H z central frequency . frequency , > 40 k H z , the proposed approach fails to predict the dispersion curve. This is due to the broad- band nature of the 1 − cycle tone bursts and the limited frequency range over which FRF s are obtained. The very low-frequenc y part of the dispersion curve is also missing since the large wav elength of longi- tudinal wa v es at such lo w frequencies hinders the separation of incident and reflected wa v eforms. 6. On the E ff ects of Boundary Conditions In this section, the e ff ects of boundary conditions on the performance of the proposed approach are in vestigated. While end conditions do not a ff ect wav e propagation characteristics along the wav e guide, they determine ho w wa v es are reflected at the boundaries, which has a profound impact on measured FRF s . The goal of this section is to demonstrate that such changes in FRF s will not a ff ect dispersion curves estimates. For this purpose, two additional combinations of boundary conditions are in vestigated, these are clamped-free and pinned-pinned conditions. out-of-plane receptance FRF s for these boundary conditions are obtained using the SEM. The VF algorithm is then used to create a data-driv en, SIMO model for each set of boundary conditions. T able 2 summarizes fitting results for these boundary condi- tions, along with the free-free conditions discussed in the previous sections. The dimensions of the state matrices and the relative fit error for each case are giv en in the table. Figure 8 depicts out-of-plane re- 14 0 10 20 30 40 50 Frequency (kHz) 4 4.5 5 5.5 6 Group Velocity (km/s) Spectral Element Method Data-Driven Approach Figure 7: Comparison of group velocity curves estimated using the proposed data-dri ven approach (red) and that predicted by the SEM (blue) for the first symmetric (longitudinal) wa ve mode. ceptance FRF s of the beam with clamped-free and pinned-pinned boundary conditions, where simulated and fitted FRF s are compared. The impact of boundary conditions on FRF s is e vident. 10 12 14 16 18 20 Frequency (kHz) 10 -10 10 -9 10 -8 10 -7 10 -6 Mag.- Disp./Force (m/N) CF SEM CF Fit PP SEM PP Fit Figure 8: Out-of-plane receptance FRF s of the beam with clamped-free and pinned-pinned boundaries over the frequency range of 10 k H z to 20 k H z . Simulated and fitted FRF s are compared. FRF s are calculated 1 in . away from the excitation point. Follo wing the procedure outlined in Section 4 , dispersion curv es are estimated using the clamped- free, and the pinned-pinned, data-dri ven models. The results are summarized in Figure 9 . Although the FRF s used to obtain the data-driv en models are a ff ected by changes in boundary conditions, such e ff ects are not reflected on the estimated dispersion curves. Both clamped-free and pinned-pinned models are capable of accurately estimating dispersion curv es, as suggested by the figure. It should be noted that the previously discussed limitations at the upper- and lo wer-ends of the frequency range are applicable to the current cases. While the analysis presented in this section is limited to flexural deformations, the conclusions are equally valid for longitudinal deformations. Observ ations based on the numerical experiments presented in this section highlight the capability of the approach in dealing with v arious boundary conditions. In future studies, the robustness of the proposed approach to uncertainties and noise in the experimental data of structures under real-boundary conditions will be in vestigated. 15 0 10 20 30 40 50 Frequency (kHz) 0 0.5 1 1.5 2 2.5 Group Velocity (km/s) Spectral Element Method Data-Driven Approach (a) 0 10 20 30 40 50 Frequency (kHz) 0 0.5 1 1.5 2 2.5 Group Velocity (km/s) Spectral Element Method Data-Driven Approach (b) Figure 9: Group velocity curves for the first anti-symmetric wave mode estimated using (a) fixed-free, and (b) pinned-pinned data-driv en models. 7. Conclusions In this paper , a new data-dri v en modeling approach for estimating dispersion curv es is presented. The V ector Fitting method is adopted to dev elop a single-input-multi-output, state-space, dynamical model of the beam under test based on receptance FRF s . FRF s corresponding to longitudinal and lateral deforma- tions of the beam are simulated ov er the frequency range of 0 − 50 k H z . The data-dri ven model is able to accurately capture the dynamic behaviour of the system within the frequenc y range of interest with an error of E rel L 2 = 2 . 13 × 10 − 7 for lateral deformations, and E rel L 2 = 3 . 12 × 10 − 7 for longitudinal deformations. The data dri ven model is then used to study wav e propagation along the beam and obtain its dis- persion curves. Group velocity curves corresponding to symmetric and anti-symmetric wa ve modes are estimated with this approach ov er the frequency range of 2 − 40 k H z . The results are found to be in very good agreement with the numerical predictions of the SEM. Dispersion curv es at frequencies lo wer than 2 k H z cannot be estimated directly using this approach. This is due to the large wa v elength at such lo w frequencies compared to the length of the beam under test, which hinders the separation of incident and reflected wa v eforms. While this is a limitation of all current experimental practices, the proposed approach allo ws for potential solutions to be implemented, such as the introduction of artificial damping to attenuate reflected wa v es, which will be addressed in future studies. The e ff ects of boundary conditions on the performance of the dev eloped approach hav e also been studied in this work. Although boundary conditions ha v e a profound impact on the FRF s that are used for the data-driv en models, estimated group velocity curves are found to be una ff ected by such conditions. This allo ws for this approach to be used regardless of boundary conditions of the structure under test, which is crucial for many practical applications. The work presented herein demonstrates the feasibility of estimating dispersion curves based on steady-state FRF s using this ne w data-dri ven modeling approach. The experimental implementation of this technique will in vestigated in future studies, where the issues of model order selection, parameter identification, and the introduction of artificial damping will be in vestigated. Future studies will also ad- dress the scenarios where di ff erent wa v e modes are coupled or simultaneously excited where the impact of coupling at the FRF -le vel on predicted dispersion curv es will be in vestigated. 16 Acknowledgment T arazaga would like to ackno wledge the support provided by the John R. Jones III Faculty Fello wship. The work of Gugercin was supported in parts by NSF through Grant DMS-1522616 and DMS-1720257. The work of Albakri was supported in parts by the Federal Railroad Administration. A ppendix A: Spectral Element Matrices For a tw o-node spectral finite element, vectors d and F in Eq. 5 as defined as follows d = n U 0 ( x 1 ) W 0 ( x 1 ) Φ ( x 1 ) U 0 ( x 2 ) W 0 ( x 2 ) Φ ( x 2 ) o T , F = n − ¯ F x ( x 1 ) − ¯ F z ( x 1 ) − ¯ M ( x 1 ) ¯ F x ( x 2 ) ¯ F z ( x 2 ) ¯ M ( x 2 ) o T , where U 0 , W 0 , Φ , ¯ F x , ¯ F z , and ¯ M are, respectiv ely , the longitudinal displacement, the lateral displacement of the beam’ s neutral axis, the angle of rotation of the neutral axis normal vector , the externally applied axial force, lateral force, and bending moments in the frequency domain. x 1 and x 2 are the coordinates of the left and right nodes of the spectral finite element. The shape functions matrix, Ψ (˙ ıω ), in Eq. 5 , is defined for a two-node spectral finite element as follo ws Ψ (˙ ıω ) =                          ζ 11 ζ 21 0 0 0 0 0 0 ζ 31 ζ 41 ζ 51 ζ 61 0 0 r 33 ζ 31 r 34 ζ 41 r 35 ζ 51 r 36 ζ 61 ζ 12 ζ 22 0 0 0 0 0 0 ζ 32 ζ 42 ζ 52 ζ 62 0 0 r 33 ζ 32 r 34 ζ 42 r 35 ζ 52 r 36 ζ 62                          , where ζ m j = e − ˙ ı k m ( ˙ ıω ) x j with m = 1 , 2 , ..., 6, and j = 1 , 2. The non-zero elements of the boundary conditions matrix, G ( ˙ ıω ) are G 1 s = ˙ ı k sn E Ae − ˙ ı k sn x 1 G 2 t = G A ¯ K ( − ˙ ı k tn + r 3 t ) e − ˙ ı k tn x 1 G 3 t = − ˙ ı E I k tn e − ˙ ı k tn x 1 G 4 s = ˙ ı k tn sn E Ae − ˙ ı k sn x 2 G 5 t = G A ¯ K ( − ˙ ı k tn + r 3 t ) e − ˙ ı k tn x 2 G 6 t = − ˙ ı E I k tn e − ˙ ı k tn x 2 where s = 1 , 2 and t = 3 , 4 , 5 , 6. References [1] M. I. Albakri, V . S. Malladi, P . A. T arazaga, Low-frequenc y acoustoelastic-based stress state char - acterization: Theory and experimental validation, Mechanical Systems and Signal Processing 112 (2018) 417–429. 17 [2] J. Schloemann, V . V . N. S. Malladi, A. G. W oolard, J. M. Hamilton, R. M. Buehrer , P . A. T arazaga, V ibration ev ent localization in an instrumented b uilding, in: Experimental T echniques, Rotating Machinery , and Acoustics, V olume 8, Springer , 2015, pp. 265–271. [3] J. D. Poston, J. Schloemann, R. M. Buehrer , V . V . N. S. Malladi, A. G. W oolard, P . A. T arazaga, T owards indoor localization of pedestrians via smart building vibration sensing, in: Localization and GNSS (ICL-GNSS), 2015 International Conference on, IEEE, 2015, pp. 1–6. [4] A. G. W oolard, P . A. T arazaga, Applications of dispersion compensation for indoor vibration ev ent localization, Journal of V ibration and Control (2017) 1077546317744997. [5] C. Caloz, Metamaterial dispersion engineering concepts and applications, Proceedings of the IEEE 99 (10) (2011) 1711–1719. [6] K. F . Gra ff , W a ve motion in elastic solids, Courier Corporation, 2012. [7] A. Chaigne, Structural acoustics and vibrations, in: Springer Handbook of Acoustics, Springer , 2014, pp. 941–1000. [8] G.-R. Liu, Z. Xi, Elastic wa v es in anisotropic laminates, CRC press, 2001. [9] J. F . Doyle, W a ve propagation in structures: spectral analysis using fast discrete F ourier transforms, Springer , 1997. [10] S. Gopalakrishnan, A. Chakraborty , D. R. Mahapatra, Spectral finite element method: wa v e prop- agation, diagnostics and control in anisotropic and inhomogeneous structures, Springer Science & Business Media, 2007. [11] U. Lee, Spectral element method in structural dynamics, John W iley & Sons, 2009. [12] M. Palacz, Spectral methods for modelling of wa v e propagation in structures in terms of damage detectiona re vie w , Applied Sciences 8 (7) (2018) 1124. [13] A. Marzani, E. V iola, I. Bartoli, F . L. Di Scalea, P . Rizzo, A semi-analytical finite element formula- tion for modeling stress wav e propagation in axisymmetric damped wa ve guides, Journal of Sound and V ibration 318 (3) (2008) 488–505. [14] R. Joseph, L. Li, M. F . Haider , V . Giurgiutiu, Hybrid safe-gmm approach for predictiv e modeling of guided wa v e propagation in layered media, Engineering Structures 193 (2019) 194–206. [15] B. R. Mace, E. Manconi, Modelling wav e propagation in two-dimensional structures using finite element analysis, Journal of Sound and V ibration 318 (4-5) (2008) 884–902. [16] C. Zhou, J. Lain ´ e, M. Ichchou, A. Zine, W av e finite element method based on reduced model for one-dimensional periodic structures, International Journal of Applied Mechanics 7 (02) (2015) 1550018. [17] H. Al Ba’ba’a, M. Nouh, T . Singh, Formation of local resonance band gaps in finite acoustic meta- materials: A closed-form transfer function model, Journal of Sound and V ibration 410 (2017) 429– 446. 18 [18] F . Li, C. Zhang, C. Liu, Acti v e tuning of vibration and wa ve propagation in elastic beams with periodically placed piezoelectric actuator / sensor pairs, Journal of Sound and V ibration 393 (2017) 14–29. [19] S. Bilbao, W ave and scattering methods for numerical simulation, John W iley & Sons, 2004. [20] B. D. O. Anderson, A. C. Antoulas, Rational interpolation and state variable realization, Linear Algebra and Its Applications, Special Issue on Matrix Problems 137 / 138 (1990) 479–509. [21] A. C. Antoulas, S. Lefteriu, A. C. Ionita, A tutorial introduction to the Loewner frame work for model reduction, in: Model Reduction and Approximation, SIAM, Philadelphia, 2017, pp. 335– 376. [22] A. Mayo, A. Antoulas, A frame work for the solution of the generalized realization problem, Linear Algebra and Its Applications 425 (2-3) (2007) 634–662. [23] P . Schulze, B. Unger , C. Beattie, S. Gugercin, Data-dri ven structured realization, Linear Algebra and its Applications 537 (2018) 250–286. [24] J. M. Hokanson, Projected nonlinear least squares for exponential fitting, SIAM Journal on Scien- tific Computing 39 (6) (2017) A3107–A3128. [25] P . Gonnet, R. Pach ´ on, L. N. T refethen, Robust rational interpolation and least-squares, Electronic T ransactions on Numerical Analysis 38 (2011) 146–167. [26] B. Gustavsen, A. Semlyen, Rational approximation of frequency domain responses by vector fitting, IEEE T ransactions on po wer deli very 14 (3) (1999) 1052–1061. [27] Z. Drma ˇ c, S. Gugercin, C. Beattie, Quadrature-based vector fitting for discretized H 2 approxima- tion, SIAM J. Sci. Comp. 37 (2) (2015) A625–A652. [28] M. Berljafa, S. G ¨ uttel, The RKFIT algorithm for nonlinear rational approximation, SIAM Journal on Scientific Computing 39 (5) (2017) 2049–2071. [29] Y . Nakatsukasa, O. S ` ete, L. N. T refethen, The AAA algorithm for rational approximation, SIAM Journal on Scientific Computing 40 (3) (2018) A1494–A1522. [30] Y . Nakatsukasa, L. N. Trefethen, An algorithm for real and complex rational minimax approxima- tion, arXi v preprint arXi v:1908.06001. [31] V . V . S. Malladi, M. I. Albakri, P . A. T arazaga, S. Gugercin, Data-driv en modeling techniques to estimate dispersion relations of structural components, in: ASME 2018 Conference on Smart Materials, Adaptiv e Structures and Intelligent Systems, American Society of Mechanical Engineers, 2018, pp. V001T03A026–V001T03A026. [32] U. Lee, D. Kim, I. Park, Dynamic modeling and analysis of the pzt-bonded composite timoshenko beams: spectral element method, Journal of Sound and V ibration 332 (6) (2013) 1585–1609. [33] M. I. Albakri, P . A. T arazaga, Dynamic analysis of a piezoelectric augmented beam system with adhesi ve bonding layer e ff ects, Journal of Intelligent Material Systems and Structures 28 (2) (2017) 178–194. 19 [34] Z. Drma ˇ c, S. Gugercin, C. Beattie, V ector fitting for matrix-valued rational approximation, SIAM Journal on Scientific Computing 37 (5) (2015) A2346–A2379. [35] C. Sanathanan, J. K oerner , T ransfer function synthesis as a ratio of two complex polynomials, IEEE T rans. Autom. Control 8 (1) (1963) 56–58. [36] Z. Drmac, S. Gugercin, C. Beattie, Quadrature-based vector fitting for discretized h 2 approxima- tion, SIAM Journal on Scientific Computing 37 (2) (2015) A625–A652. [37] B. Gustavsen, Improving the pole relocating properties of vector fitting, IEEE T ransactions on Po wer Deli v ery 21 (3) (2006) 1587–1592. [38] A. Chinea, S. Gri vet-T alocia, On the parallelization of V ector Fitting algorithms, IEEE T ransactions on Components, Packaging and Manuf acturing T echnology 1 (11) (2011) 1761–1773. [39] M. I. Albakri, P . A. T arazaga, A novel acoustoelastic-based technique for stress measurement in structural components, in: Dynamics of Civil Structures, V olume 2, Springer , 2016, pp. 49–56. 20