Estimating Experimental Dispersion Curves from Steady-State Frequency Response Measurements
V. V. N. Sriram Malladi, Mohammad I. Albakri, Manu Krishnan, Serkan Gugercin, Pablo A. Tarazaga
EE stimating E xperimental D ispersion C urves from S teady -S tate F requency R esponse M easurements A P reprint
V. V. N. Sriram Malladi ∗† Department of Mechanical Engineering -Engineering MechanicsMichigan Technological University,Houghton, Michigan, 49931 [email protected]
Mohammad I. Albakri † Department of Mechanical EngineeringTennessee Technological University,Cookeville, Tennessee, 38505 [email protected]
Manu Krishnan
Department of Mechanical EngineeringVirginia TechBlacksburg, Virginia, 24060 [email protected]
Serkan Gugercin
Department of MathematicsVirginia TechBlacksburg, Virginia, 24060 [email protected]
Pablo A. Tarazaga
Department of Mechanical EngineeringVirginia TechBlacksburg, Virginia, 24060 [email protected]
January 5, 2021 A bstract Dispersion curves characterize the frequency dependence of the phase and the group velocities ofpropagating elastic waves. Many analytical and numerical techniques produce dispersion curves fromphysics-based models. However, it is often challenging to accurately model engineering structureswith intricate geometric features and inhomogeneous material properties. For such cases, thispaper proposes a novel method to estimate group velocities from experimental data-driven models.Experimental frequency response functions (FRFs) are used to develop data-driven models, which arethen used to estimate dispersion curves. The advantages of this approach over other traditionally usedtransient techniques stem from the need to conduct only steady-state experiments. In comparison,transient experiments often need a higher-sampling rate for wave-propagation applications and aremore susceptible to noise.The vector-fitting ( VF ) algorithm is adopted to develop data-driven models from experimental in-planeand out-of-plane FRFs of a one-dimensional structure. The quality of the corresponding data-drivenestimates is evaluated using an analytical Timoshenko beam as a baseline. The data-driven model(using the out-of-plane FRFs) estimates the anti-symmetric ( A ) group velocity with a maximumerror of 4% over a 40 kHz frequency band. In contrast, group velocities estimated from transientexperiments resulted in a maximum error of 6% over the same frequency band. K eywords data-driven models · dispersion curves · least-squares · vector-fitting algorithm · longitudinal and flexuralmodels Numerous applications incorporate guided waves as essential components. Many structural health monitoring algorithms,for instance, rely on guided-waves to detect and localize structural damages [13, 42]. These techniques depend onaccurately discerning variations in guided-wave features [33]. However, inhomogeneity in material properties and ∗ Address all correspondence to this author. † These two authors have contributed equally. a r X i v : . [ phy s i c s . d a t a - a n ] J a n preprint - J anuary
5, 2021complex geometric features makes it challenging to determine their actual attributes. As a result, there is limitedapplicability of these techniques to in-homogeneous systems such as reinforced-concrete floors and non-uniformcomposites.Over the years, several approaches have been developed to experimentally determine dispersion curves and capturethe frequency-dependent nature of wave numbers associated with di ff erent wave modes. These approaches can bebroadly categorized based on the nature of testing as time-domain transient approaches or frequency-domain steady-state approaches. Traditionally, most of the dispersion estimation relies on the time domain analysis of transientmeasurements. Popular techniques include time-of-flight analysis [48], cross-correlation techniques [48], wavelet-based methods [31], and model-based approaches [27]. Depending on the structure under test, there are advantages anddisadvantages of using time-domain techniques. The advantages extend from the need for minimal spatial information.For instance, estimating the time taken by a guided waveform to travel between two know spatial locations is su ffi cientto establish its phase and group velocity. However, experimental noise and contamination of measured signals fromreflections hamper the reliability of the transient approaches.Frequency-domain approaches, on the other hand, such as the spatial Fourier transform [45], the Bayesian estima-tion [44], and the in-homogeneous wave correlation [9] rely on the steady-state dynamic response in the form ofFrequency Response Functions (FRFs) and Operating Deflection Shapes (ODSs). The advantages of steady-stateapproaches stem from the relative ease of measuring FRFs and ODSs as opposed to propagating waveforms [11, 14, 35].In time-domain tests, accurate measurement of a propagating wave necessitates a sampling rate of at least 10-15times the maximum frequency content of the waveform. In contrast, the need for strict temporal resolution is relaxedwith steady-state dynamical measurements where a sampling rate of 2 - 2.5 times the maximum frequency of interestis su ffi cient. Generally, repeatability is higher with steady-state measurements, and averaged FRFs have superiorsignal-to-noise ratio compared to transient time-domain measurements. However, successful implementation of theaforementioned frequency-domain approaches requires fine spatial resolution to avoid spatial aliasing. As a result, thesteady-state response needs to be measured at a large number of points, which adds to the cost and complexity of theunderlying experiments.In a recent e ff ort by the authors [1], a novel data-driven-based approach that brings together the advantages of boththe time-domain and frequency-domain testing has been developed. This approach uses FRFs to develop a data-driven dynamical model for the structure under test. Dispersion curves are then determined using the developedmodel by applying transient time-domain techniques to numerically simulated waveforms. Thus, the advantages ofsteady-state experimentation is combined with the relaxed spatial resolution requirements of the transient approaches.While there is a plethora of approaches for constructing data-driven rational approximations from measured data (see,e.g., [4–6, 32, 37] for interpolatory methods; [7, 15, 26, 29, 30] for least-squares based approaches; and [10, 22, 39, 40]for rational approximations that combine the interpolatory and least-squares frameworks), the least-squares basedvector-fitting ( VF ) method [26] was the method of choice in [1] and the ability of the FRF-based data-driven models tosimulate the transient dynamics was established. Through numerical experiments, the capabilities of this new approachhas been demonstrated where accurate estimates of dispersion curves corresponding to longitudinal and flexural wavemodes in one-dimensional structures were obtained.This paper extends this investigation by tailoring the VF algorithm to experimentally measured data. Other rationaldata-driven modeling frameworks such as those listed above can be employed here as well; however this work usesthe VF algorithm to experimentally validate and test the authors’ earlier work [1], which also employed VF . Unlikenumerically generated counterparts, experimentally measured FRFs present several challenges to data-driven modeling.This includes noise, multi-mode excitation (resulting in mixed axial, flexural, and torsional FRFs), and non-idealboundary conditions. The current work addresses practical issues in fitting such noisy data and provides a VF frameworkto build data-driven models from experimentally measured FRFs. The impact of such artifacts on data-driven models,their ability to capture the dynamic system’s response, and resulting predictions of dispersion curves is investigated.The current paper is structured as follows: Section 2 describes the experimental setup, and the test procedure followedto generate in-plane and out-of-plane FRF. A brief overview of the nonlinear least-squares problem and the iterativevector fitting techniques is presented in Section 3. Implementation of the VF technique to develop out-of-plane andin-plane data-driven models is discussed in Section 4. The performance of these models in estimating dispersion curvesis the topic of discussion in Section 5. Towards the end, an overview of the results is summarized along with finalremarks for future work. The experimental set up used to measure in-plane, and out-of-plane FRFs is depicted in Figure 1. A 48 in . longaluminum beam with a rectangular cross-section of (1 in . × / in . ) has been selected for this study. Free-free boundary2 preprint - J anuary
5, 2021conditions are approximated by suspending the beam under test with two bungee cords. Two piezoceramic wafers(PZT-5H) of dimensions L p × W p (1 in . × / in . ) are bonded to the upper and lower faces of the beam, 18 . in . from itsleft end. The poling direction of the piezoceramics (PZTs), along with the two possible electrical configurations, areshown in the figure. Configuration 1 is designed to amplify the net bending moments generated by actuating the PZTswhile cancelling in-plane forces. This allows for the measuring out-of-plane FRFs with minimal contribution fromin-plane deformations. Configuration 2, on the other hand, is designed to neutralize bending moments and amplifyin-plane forces to facilitate the measurement of in-plane FRFs. However, perfect decoupling between in-plane andout-of-plane deformations is not possible due to PZTs mounting imperfections and misalignment. Thus, a combinationof vibrational responses corresponding to longitudinal, flexural, and torsional modes always appear in experimentalmeasurements. An in-depth experimental modal analysis would di ff erentiate these modes of vibration. The vector-fittingalgorithm can conveniently include the dynamics corresponding to all of these vibrational modes in the data-drivenmodels. Aluminum beam3D Scanning Laser Doppler VibrometerPolytec controllerand DAQ Piezoelectric patch attached with beam
PZT Patch 1PZT Patch 2
Trek amplifierPZD350A-2-L
Beam side view
Configuration 2Configuration 1
Figure 1: Experimental setup measuring the in-plane and the out-of-plane vibrations of the aluminum beam whenexcited with two piezoceramicsA 3D Scanning Laser Doppler Vibrometer (3D-SLDV) is used to measure the beam’s dynamic response when excitedwith PZTs under both configurations. The 3D-SLDV has three laser vibrometers that simultaneously measure thevelocity response at every scanning point on the beam. The Polytec software then analyzes the measurements andtransforms the data to compute the velocity components along the three perpendicular axes. As a result, it is possible tomeasure both in-plane and out-of-plane responses of the beam. 3D-SLDV also triangulates the geometric coordinatesof the points as it measures the response along the beam. This provides accurate information about the coordinates anddistance between measurement points, which is later used for estimating dispersion curves. It is important to note thatthere is an uncertainty of 0 . mm while measuring the coordinates of the measurement points, which will be reflected ingroup velocity calculations. In the present set of experiments, the velocity measurements are recorded at 19 equallyspaced scan points along a straight line of about 18 in , starting 2 in away from the PZTs.Using the experimental setup shown in Figure 1, two sets of experiments have been conducted. The first set targeted thesteady-state dynamic response of the beam in the form of FRFs. For this purpose, the beam is excited with a swept-sinesignal covering the frequency range of 0 − kHz . These FRFs are then used to construct the data-driven model ofthe beam. A sampling frequency of 128 kS / s is used for this experiment, and the measured FRFs have a frequencyresolution of 0 . Hz . At each point, the mean FRF and the mean coherence estimates are determined for a set of 10repeated trials. The second set of experiments targeted the transient dynamic response of the beam. For this purpose,the beam is excited with a 2 − or 4 − cycle, amplitude-modulated, tone burst with a center frequency ranging from 1 kHz to 50 kHz . This set of experiments aims to provide the data necessary to directly measure the beam’s dispersion curves,3 preprint - J anuary
5, 2021which is used for comparison purposes. A sampling rate of 1 . MS / s is used for this set of experiments, which is tentimes the rate used for measuring FRFs in the first set. For both sets of experiments, in-plane and out-of-plane responsesare excited and measured using the configurations mentioned earlier. Excitation signals are generated using the Polytecfunction generator and a Trek (PZD350A-2-L) power amplifier, as shown in the figure. During the experiments, the Polytec data acquisition system provides the frequency response of the beam at n out = v ( t ) = [ v ( t ) , ..., v n out ( t )] T ) to a voltage input f ( t ) applied to the PZT. Let V ( s ) and F ( s ) denote the Laplacetransforms of v ( t ) and f ( t ), respectively. Then, the corresponding frequency response function H ( s ) satisfies V ( s ) = H ( s ) F ( s ). Our experimental set-up samples this frequency response function H ( s ) at multiple frequencies s j = ˙ ıω j for j = , , . . . , N . In our experiments, N = / single output (SISO) FRF at afrequency s = ˙ ıω is then calculated using H estimator and / or H estimator [19] given by, H (˙ ıω ) = G v f (˙ ıω ) G f f (˙ ıω ) and H (˙ ıω ) = G vv (˙ ıω ) G f v (˙ ıω ) , (1)where the cross-power spectral densities and auto-power spectral densities are defined as G v f (˙ ıω ) = E (cid:2) V (˙ ıω ) F (˙ ıω ) ∗ (cid:3) , G vv (˙ ıω ) = E (cid:2) V (˙ ıω ) V (˙ ıω ) ∗ (cid:3) and G f f (˙ ıω ) = E (cid:2) F (˙ ıω ) F (˙ ıω ) ∗ (cid:3) , respectively; and E (cid:2) . (cid:3) denotes the estimator functionand {·} ∗ is the complex-conjugate of a the matrix. In the absence of noise, H (˙ ıω )and H (˙ ıω ) result in identical FRFs.However, in the presence of noise in either the input or the output measurements, both methods approximate theinherent dynamics of the structure. Generally, the H (˙ ıω ) formulation is sensitive to noise on the input signal f ( t )and underestimates H (˙ ıω ). On the contrary, the H (˙ ıω ) formulation is sensitive to noise on the output signal v ( t ) andoverestimate the true FRF. Therefore, the quality of the measured FRFs can be determined by the coherence function γ (˙ ıω ) defined as the ratio of the two formulations, γ (˙ ıω ) = H (˙ ıω ) H (˙ ıω ) = |G v f (˙ ıω ) | G vv (˙ ıω ) G f f (˙ ıω ) . (2)Typically, coherence values lie in the range of [0 − / multi-output (MIMO) systems as in the present case,the FRF formulation can be written as H (˙ ıω ) = G v f (˙ ıω ) G − f f (˙ ıω ) . (3)The output cross-spectral density can be determined using two methods G vv (˙ ıω ) = H (˙ ıω ) G f f (˙ ıω ) H (˙ ıω ) ∗ and ˜ G vv (˙ ıω ) = H (˙ ıω ) G f v (˙ ıω ) . (4)Although these two formulation are equivalent theoretically, experimental noise results in numerical discrepancies. Thecoherence function of a MIMO system is defined as the correlation coe ffi cient between each output ( v i ( t )) to the input( f ( t )). Therefore, the multiple coherence function γ v i : f (˙ ıω ) is the ratio of the diagonal elements of ˜ G vv and G vv , i.e., γ v i : f (˙ ıω ) = ˜ G v i v i (˙ ıω ) G v i v i (˙ ıω ) . (5)This work utilizes the H estimates of the FRFs to develop Single-Input Multi-Output (SIMO) data-driven models. Forsimplicity, the subscript is dropped and the experimental FRFs are denoted by H (˙ ıω ) throughout the paper. Therefore,the experimental FRF data are H (˙ ıω j ) ∈ C × for j = , , . . . , . (6)Given the FRF data in (6), the objective is to construct a low-order rational approximation (a reduced-order frequencyresponse function) (cid:101) H ( s ) that approximates the data, and thus the true frequency response function H ( s ), in anappropriate measure. This measure will be clarified below. Once (cid:101) H ( s ) is obtained, it could be used to study theunderlying dynamics and determine the system parameters such as the poles and the residues of the system. Generally,most of the modal-identification techniques in the frequency domain start by assuming a rational function to fit theFRFs. 4 preprint - J anuary
5, 2021A rational function could be represented in many equivalent forms. To simplify the presentation, for time being, assumethat H ( s ) is a SISO system; thus so is the synthesized FRF , (cid:101) H ( s ). Let (cid:101) H ( s ) have order r . Then, it could be simplywritten as (cid:101) H ( s ) = n ( s ) d ( s ) = α + α s + · · · + α r − s r − β + β s + · · · + β r − s r − + s r , (7)where the numerator n ( s ) and denominator d ( s ) are polynomial of degree r − r , respectively,. The numeratorcoe ffi cients { α i } and the denominator coe ffi cients { β i } fully describe it. Since, in our structural dynamics, the underlyingfrequency response function is strictly proper, i.e., the order of the denominator is strictly greater than that of thenumerator, (cid:101) H ( s ) in (7) is chosen to have the same form. This is not a restriction and the general case can be easilyhandled, but is not relevant to our set-up.The unknowns coe ffi cients { α k } and { β k } are typically estimated through the least-squares ( LS ) approach. In other words,given the FRF data in (6), the goal is to construct (cid:101) H ( s ) such the LS error J = N (cid:88) j = w j (cid:12)(cid:12)(cid:12)(cid:12) H (˙ ıω j ) − (cid:101) H (˙ ıω j ) (cid:12)(cid:12)(cid:12)(cid:12) (8)is minimized where w j ’s are the weights corresponding to the frequency ω j . This LS problem is nonlinear due to itsnonlinear dependence on the denominator coe ffi cients. The first step in most techniques popular in the modal-testingcommunity [3] is to linearize this nonlinear problem via H ( s ) − (cid:101) H ( s ) = H ( s ) − n ( s ) d ( s ) (cid:32) d ( s ) H ( s ) − n ( s ) . Then, instead of the original nonlinear LS error in (8), one minimizes, as done in, e.g., [20], the linearized LS error J lin = N (cid:88) j = w j (cid:12)(cid:12)(cid:12) d (˙ ıω j ) H (˙ ıω j ) − n (˙ ıω j ) (cid:12)(cid:12)(cid:12) . (9)Similarly, other modal identification techniques linearize the rational-polynomial to apply a LS procedure to estimatecoe ffi cients [19, 24, 36]. Allenmang et.al. [2, 3] have developed a unified matrix polynomial framework to identifythe similarities and di ff erences between multiple modal identification techniques available in the literature. Theseapproaches convert the nonlinear least-square problem (8) to a linear one (9) and the solution is obtained, without anyiteration, in one step by solving the resulting linear problem. In this work, to better handle the non-linearity, an iterativeprocedure will be employed: The nonlinear least-problem will be converted to solving a sequence of linear least-squareproblems.Using (cid:101) H ( s ) = n ( s ) / d ( s ), the cost function J in (8) can be written as J = N (cid:88) j = w j (cid:12)(cid:12)(cid:12)(cid:12) H (˙ ıω j ) − (cid:101) H (˙ ıω j ) (cid:12)(cid:12)(cid:12)(cid:12) = N (cid:88) j = w j (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) H (˙ ıω j ) d (˙ ıω j ) − n (˙ ıω j ) d (˙ ıω j ) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) . (10)The linearization approach mentioned above amounts to ignoring the denominator in this expression to reach a linear LS problem.Instead, Sanathanan and Koerner [43] proposed a technique in which the denominator is incorporated into the problem,leading to an iterative algorithm solving a sequence of linearized LS problem. In this approach, as opposed to ignoringthe denominator in (10), one replaces it with the denominator from the previous step.Let (cid:101) H ( k ) ( s ) = n ( k ) ( s ) / d ( k ) ( s ) denote the synthesised / identified model at the k th step of the Sanathanan and Koerner ( SK )iteration. Starting with the initial guess d (0) ( s ) ≡
1, the SK iteration, at the k th step, minimizes the linear LS error N (cid:88) j = w j (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) n ( k + (˙ ıω j ) − d ( k + (˙ ıω j ) H ( ω j ) d ( k ) ( ω j ) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) (11)to solve for n ( k + ( s ) and d ( k + ( s ) and then set (cid:101) H ( k + ( s ) = n ( k + ( s ) / d ( k + ( s ). Note that the LS problem in (11) is linearin the variables n ( k + ( s ) and d ( k + ( s ). The iteration continues until a desired convergence tolerance is achieved. Readersare referred to [46, 47], where the SK iteration in conjunction with a Gauss-Newton based optimization approach [41] isemployed to fit MIMO FRFs in the modal-testing community.5 preprint - J anuary
5, 2021Writing (cid:101) H ( s ) as a ratio of two polynomials is only one way of representing a rational function. A mathematicallyequivalent but numerically more stable representation is the so-called barycentric form [8] where (cid:101) H ( s ) is represented asa ratio of two rational functions. The Vector-Fitting ( VF ) algorithm [26] precisely does this; it performs the SK iterationusing the barycentric form. As in the case of the SK iteration, VF is also an iterative procedure. At the k th step, (cid:101) H ( s ) isparameterized as (cid:101) H ( k ) ( s ) = ˜ n ( k ) ( s )˜ d ( k ) ( s ) = (cid:80) ri = φ ( k ) i / ( s − λ ( k ) i )1 + (cid:80) ri = ψ ( k ) i / ( s − λ ( k ) i ) , (12)where φ ( k ) i , ψ ( k ) i , λ ( k ) i ∈ C and { λ ( k ) i } are an arbitrary set of mutually distinct points. Therefore, in VF , the model (cid:101) H ( s )is parameterized by φ ( k ) j , ψ ( k ) j , and λ ( k ) j . The SK iteration, in this new barycentric form, can be considered as runningthe iteration (11) for fixed { λ ( k ) i } . VF takes advantage of the fact that { λ ( k ) i } are free and can be updated in every stepof the iteration. It does so by choosing the zeros of the numerator ˜ d ( k ) ( s ) = + (cid:80) rj = ψ ( k ) j / ( s − λ ( k ) j ) in (12) as the nextset of { λ ( k ) i } . Similar to the SK iteration, in every step of VF , a linear LS problem is solved to update the parameters φ ( k ) i , ψ ( k ) i , and λ ( k ) i ∈ C . The iteration and updates continue until convergence, upon which, due to the updating scheme,˜ d ( k ) ( s ) →
1, yielding the final form of the approximant (cid:101) H ( s ) = r (cid:88) k = φ k s − λ k , (13)in the pole-residue form: the sum of ‘ r ’ modal components, i.e., the poles ( { λ i } ) and residues ( { φ i } ). Therefore, eventhough λ ( k ) i ’s are not the poles of (cid:101) H ( s ) in (12) initially, upon convergence, they do become the poles of (cid:101) H ( s ). Due tothis connection, we will call { λ ( k ) i } as the poles; but this should be understood so in the context of convergence. A briefremark on the convergence of VF is in order: Even though one can construct examples where VF does not converge [34],when { λ i } are initialized appropriately, it usually converges quickly. Proper initialization becomes even more importantif the order r grows modestly. This is indeed the case in our setting and will be discussed in more detail below. Forfurther details on VF , see [12, 15, 16, 23, 25, 26] and the references therein. In the modal-testing community, a recenttechnique, called Maximum Likelihood estimation of a Modal Model [17, 18], uses the pole-residue parametrization ofthe transfer function as in (13), i.e., takes the residues φ k and poles λ k as the variables, and employs the Levenberg-Marquardt [38] scheme to fit these parameters to FRF samples. Either using the pole-residue or the barycentric form,one can indeed employ a variety of well-established optimization techniques for the resulting nonlinear LS problems;see for example [21, 28, 29] and the references therein.The formulation above assumed, for simplicity, a SISO model. In our experimental set-up, we have a SIMO systemwith n out =
19 outputs. The theoretical analysis above directly extends to our SIMO setting. For the data in (6), the LS error in (8) is now given by J = N (cid:88) j = w j (cid:13)(cid:13)(cid:13)(cid:13) H (˙ ıω j ) − (cid:101) H (˙ ıω j ) (cid:13)(cid:13)(cid:13)(cid:13) , (14)where (cid:107) · (cid:107) denotes the 2-norm. Then, the only revision in the barycentric form in (12) is that φ i is now a column vectorof size 19 ×
1. Therefore, upon the convergence of VF , we have (cid:101) H ( s ) = r (cid:88) j = φ j s − λ j where φ i ∈ C × and λ i ∈ C for i = , . . . , r . (15)It is usually preferred to have a state-space formulation for frequency-response function. Therefore, we will represent(15) equivalently in a real state-space form, i.e„ (cid:101) H ( s ) = r (cid:88) j = φ j s − λ j = (cid:101) (cid:131) ( sI − (cid:101) (cid:129) ) − (cid:101) (cid:130) where (cid:101) (cid:129) ∈ R r × r , (cid:101) (cid:131) ∈ R × r , and (cid:101) (cid:130) ∈ R r . (16)Note that the realness of the state-space realization in (16) is guaranteed by implicitly including the conjugate data, i.e., {H ( − ˙ ıω j ) } , in the FRF data. 6 preprint - J anuary
5, 2021
This section investigates the performance of two SIMO data-driven models. While one simulates the flexural response,the other simulates the beam’s in-plane response over a frequency band of 0 − kHz . There is a significant number ofnatural frequencies in the frequency range of interest. Notably, there are more out-of-plane natural frequencies for thecurrent 1D-structure. Consequently, it is challenging to work with out-of-plane FRFs. A multi-step modeling procedureis employed to fit the experimental out-of-plane data-set. On the other hand, it is feasible to vector fit all in-plane FRFsin a single setting. In the current work, a data-driven model is developed over a frequency bandwidth of [0 − kHz ]. Typically, naturalfrequencies of dispersive structures are unevenly distributed over such a broad frequency bandwidth; thus, modaldensity (number of poles in a frequency interval) is higher at lower frequencies than at higher frequency bandwidths.Modal density also plays a role in fitting FRFs, leading to over / underfitting parts of the frequency bandwidth as theexperimental FRFs are sampled uniformly. Vector fitting a broad bandwidth of 50 kHz requires many poles, resultingin a large model order r . These factors make it challenging to uniformly fit experimental measurements with 208400samples for each of the 19 FRFs. Consequently, while fitting such big-datasets, the choice of initial / starting polessignificantly impacts the convergence of the VF algorithm and the quality of the final rational approximations. Therefore,a modified approach is considered for pooling a rich choice of starting poles.First, the total frequency bandwidth is divided into multiple bands, and the VF algorithm is run independently for eachsub-band. This approach has several advantages. As the size of the FRF data-set on individual sub-bands is drasticallysmaller, a low-order model can appropriately fit the FRFs, which, in turn, significantly simplifies the pole-initializationstep for each band. The model order (the number of poles in each band) is determined (approximated) based on thenumber of peaks / phase-change in each sub-band. Additionally, di ff erent weighting functions can be applied to fitresonant peaks and anti-resonant valleys uniformly. Once each sub-bands are fitted via VF , all the fitted poles areconcatenated to form a larger-set of poles. Care is taken to avoid duplication of poles in the concatenated set. This finalset of combined poles is used to initialize the VF algorithm to fit the entire frequency bandwidth. Due to this improvedpole-initialization step, the VF algorithm converges faster and yields high-fidelity rational approximants for the full-dataset, i.e., for the full frequency band.The quality of the fitted model is determined by comparing synthesized FRFs with experimental FRFs. In the presentwork, two types of relative errors are defined for evaluating the quality of the fit, namely E L = N (cid:118)(cid:117)(cid:117)(cid:117)(cid:117)(cid:117)(cid:117)(cid:117)(cid:117)(cid:117)(cid:117)(cid:117)(cid:117)(cid:117)(cid:117)(cid:117)(cid:117)(cid:117)(cid:116) N (cid:88) j = (cid:13)(cid:13)(cid:13)(cid:13) H (˙ ıω j ) − (cid:101) H (˙ ıω j ) (cid:13)(cid:13)(cid:13)(cid:13) N (cid:88) j = (cid:13)(cid:13)(cid:13) H (˙ ıω j ) (cid:13)(cid:13)(cid:13) and E wL = N (cid:118)(cid:117)(cid:117)(cid:117)(cid:117)(cid:117)(cid:117)(cid:117)(cid:117)(cid:117)(cid:117)(cid:117)(cid:117)(cid:117)(cid:117)(cid:117)(cid:117)(cid:117)(cid:116) N (cid:88) j = w j (cid:13)(cid:13)(cid:13)(cid:13) H (˙ ıω j ) − (cid:101) H (˙ ıω j ) (cid:13)(cid:13)(cid:13)(cid:13) N (cid:88) j = w j (cid:13)(cid:13)(cid:13) H (˙ ıω j ) (cid:13)(cid:13)(cid:13) . (17)While E L is the relative error between experimental FRF - H (˙ ıω j ) and the synthesized FRF - (cid:101) H (˙ ıω j ), the error metric E wL also considers the role of the weighting scheme ( w j ) in fitting the FRFs. The initial focus is to fit the out-of-plane FRFs over shorter frequency bands and then move to the full frequency band.In the present case, eight sub-bands are identified such that the number of peaks in each sub-band does not exceed 30.Table 1 lists the details of these sub-bands. This approach facilitates the handling of large data sets. Fitting individualsub-bands yields a richer distribution of starting poles, enabling faster convergence of the VF algorithm in fitting thecomplete FRF.Peak-frequencies in each sub-band are used to initialize the VF algorithm and a unit weighting function is used to obtainan unbiased estimation of poles in each sub-band. The result is eight data-driven models that span the entire frequencyband. Table 1 summarizes the results of fitting FRFs over individual frequency sub-bands including model order andindividual sub-band fitting errors ( E L ). Figure 2 compares the magnitude and the phase of the experimental FRF withthe fitted FRF. The small deviation between the two FRFs highlights the quality of the fit across the frequency range.7 preprint - J anuary
5, 2021
10 12 14 16 18 20
Frequency [kHz] -4 -2 M ag . ( V e l ./ V o l ) [ m / ( s V ) ] Exp. Data Vector Fit Deviation
10 12 14 16 18 20
Frequency [kHz] P ha s e [ D eg r ee s ] Figure 2: Comparison of the magnitude and phase plots of the experimental and the estimated FRF (at point 13) overthe frequency sub-band of 10 − kHz . The VF estimate of the FRF is generated for the model using unit weights.Table 1: Partition of the complete frequency bandwidth into multiple sections. A unit weighting function is used todetermine poles of sub-bands, therefore E wL = E L .Group Frequency range Number of poles E L error1 0 - 1 kHz
32 8 . × − kHz
58 1 . × − kHz
44 1 . × − kHz
56 7 . × − kHz
54 1 . × − kHz
48 5 . × − kHz
40 9 . × − kHz
62 1 . × − An ideal structural dynamics test would produce well-spaced flexural peaks in the FRFs. However, in practice, it ischallenging to isolate flexural vibrations over a broad frequency band. Often, undesired artifacts show up in the FRFmeasurements. While experimental irregularities induce some artifacts, the others are due to the limits of the experiment.For instance, it is interesting to examine the number of poles in each sub-band. The three frequency sub-bands withinthe 10 − kHz frequency band require 56, 54, and 40 poles, which broadly follows the analytical perspective of a1D-beam, i.e., the frequency spacing between consecutive natural frequencies increases. However, the 40 − kHz sub-band significantly deviates from this trend, as it needs 102 poles. 1D flexural systems do not display such suddenchange in characteristics. Therefore, it is valid to assume that the structure deviates from behaving as a 1-dimensional8 preprint - J anuary
5, 2021beam above 40 kHz , which is an important observation. In later sections, the dispersion curves further validate andestablish this observation.
Iteratively solving the nonlinear LS problem with large data-sets is computationally challenging. Therefore, startingpoles that initialize the VF algorithm are picked from the pool of poles fitting FRFs in each sub-band. The selectedsubset of 368 poles spans the entire 50 kHz band. These complex poles are chosen to correspond to the peak frequenciesand other artifacts observed in the experimental FRFs. This approach avoids redundancy of starting poles, therebyavoiding over-fitting. As a result, VF algorithm yields state-space matrices of dimensions (cid:129) ∈ IR × , (cid:130) ∈ IR × ,and (cid:131) ∈ IR × .Ordinary least-squares problem formulation reduces the L error, i.e., squared errors at all data points. Such a formulationinadvertently focuses on accurately capturing high-valued data-points. As a result, the relative error in estimatingnumerically smaller data-points is considerably higher. The VF algorithm also su ff ers from similar weaknesses and,therefore, better captures FRF characteristics at resonant frequencies than at anti-resonant frequencies. A weightedleast-squares problem is one way to address this limitation. Besides, the weighting function’s choice plays a vitalrole in estimating both resonances (peaks) and anti-resonances (valleys) of the FRFs. While accurate pole estimationcorrelates with peak characteristics, fitting valleys correspond with good residue (or zeros) estimation. In this work,four weighting choices are considered. Unit weight ( w j =
1) and weak inverse FRF ( w j (˙ ıω ) = / (cid:112) | H (˙ ıω j ) | ) are two ofthe weighing options, which are also traditionally used in VF. Having a weak inverse FRF as the weighting functionenhances the data-driven model’s ability to capture numerically smaller values near anti-resonances. The remaining twooptions are based on experimental signal-to-noise ratio (SNR) as described next.During the experiments, SNR reduces at anti-resonances. If the noise floor is higher than the actual FRF value,the measured input and the output signals are incoherent. As a result, the value of the coherence function γ (˙ ıω )drops, indicating low confidence in the measured response. This experimental attribute is dominant at anti-resonances,especially in the lower frequency ranges, as represented by the coherence plots in Figure 3.The remaining two weighting functions make use of this characteristic. Coherence function γ (˙ ıω j ) and a hybridfunction ( w j (˙ ıω j ) = γ (˙ ıω j ) / (cid:112) | H (˙ ıω j ) | ) are the other two weighting functions considered in this study. Table 2summarizes the performance of the resulting data-driven models, and Figure 4 displays the quality of the data-drivenmodel, with the hybrid weighting function, at one of the 19 points (location 13). The resulting relatives errors forthe four weighting cases are similar. However, comparing the FRF estimates of all four data-driven models at lowerfrequencies shows the variation in fitted models. Figure 3 compares the out-of-plane FRF estimates at two frequencysub-bands. At higher frequencies, visually, there is little di ff erence between the four FRF estimates. However, atlow-frequency bands, the data-driven models based on the two inverse weighting options perform better in capturingthe FRF characteristics at the anti-resonances. While the weak inverse FRF weighting function appears to over-fit noiseat anti-resonances, the hybrid weighting function avoids this problem. Section 5 uses these out-of-plane data-drivenmodels to estimate dispersion curves for the flexural (the first anti-symmetric) wave mode.Table 2: Performance of VF algorithm in fitting weighted experimental FRFsWeighting function ( w j ) E L E wL Unit weight w j = . × − . × − Weak inv. FRF w j = (cid:112) H (˙ ıω j ) 3 . × − . × − Coherence w j = γ (˙ ıω j ) 2 . × − . × − Hybrid weight w j = γ (˙ ıω j ) (cid:112) H (˙ ıω j ) 3 . × − . × − This section discusses the development of in-plane data-driven models. Working with in-plane FRFs presents a di ff erentset of challenges. In the 1D structure under test, the measured in-plane response is an amalgam of longitudinal modesand flexural / torsional models. Figure 5 presents an in-plane experimental FRF. Although longitudinal modes dominatethe FRF, the contribution of other vibrations is significant. Limitations in testing are the major contributor to the9 preprint - J anuary
5, 2021 -4 -2 M ag . ( V e l ./ V o l ) [ m / ( s V ) ] E x p . C oh . ( )
30 31 32 33 34 35
Frequency [kHz] -4 -2 M ag . ( V e l ./ V o l ) [ m / ( s V ) ] E x p . C oh . ( ) Exp. Unit wt. Weak inv. FRF Coh. wt. Hybrid wt.
Figure 3: Comparison of the magnitude points of the experimental and the estimated FRFs (at point 6) over twofrequency sub-bands of 0 . − . kHz and 30 − kHz . The VF estimate of the FRF is generated for the model for allfour weighting functions. The magnitude of the experimental FRF is plotted in black. The experimental coherence isplotted in orange.non-dominant peaks in the experimental FRFs. One way to manage such data is to fit all peaks, including dominantand non-dominant resonant peaks. This approach results in a state-space model of 416 poles. The resultant state-spacematrices have the dimensions (cid:101) (cid:129) ∈ IR × , (cid:101) (cid:130) ∈ IR × , and (cid:101) (cid:131) ∈ IR × . Figure 5 presents the experimental FRFalong with the fitted FRF. This plot shows that the data-driven FRF accurately captures the dynamics of the longitudinalmodes along with artifacts from flexural / torsional vibrations. The error estimate for the 416 pole data-driven model is2 . × − . The developed data-driven, state-space model of the beam under test is used to study wave-propagation characteristicsalong the beam and calculate its dispersion curves. Figure 6 briefly summarizes the main steps involved in this process.The discussion of the first anti-symmetric (flexural) wave mode is presented first, followed by the first symmetric(longitudinal) wave mode. A 2-cycle, amplitude-modulated, sine wave, tone burst is used as the input to the beam’s data-driven model. The response to this excitation at all locations where FRFs are experimentally measured is then calculated.Incident waveforms at those locations are extracted from the simulated response, and the estimated wavenumber, ˆ k ,corresponding to the wave mode of interest, is then calculated using the following equation: U i + ( ω ) = U i ( ω ) e − ˙ ı ˆk ( x i + − x i ) , (18)10 preprint - J anuary
5, 2021
Frequency [kHz] -4 -2 M ag . ( V e l ./ V o l ) [ m / ( s V ) ] Exp. Data Vector Fit
Figure 4: Comparison of the magnitude of the experimental and estimated FRFs (at point 13) over the entire frequencyband. The VF estimate of the FRF is generated for the model using hybrid weights.
Frequency [kHz] -4 -2 M ag . - V e l ./ V o l t age [ m / s / V ] Exp. Data Vector Fit
Figure 5: Comparison of the magnitude of the experimental and estimated in-plane FRFs (at point 13) over the entirefrequency band.where U i is the vector of Fourier coe ffi cients of the signal obtained at location x i . Group velocity is then calculated as V G = ∂ω/∂ ˆ k . Dispersion curves are constructed one frequency band at a time by varying the central frequency of theexcitation signal and repeating the aforementioned process.A minimum of two measurement locations are required to get an estimate of the dispersion curve. However, FRFs at 19locations have been obtained and used to build the SIMO state-space model. This allows for responses to be simulatedat those 19 locations, and thus, a total of 171 estimates can be obtained. A subset of these combinations is used in thisstudy. Measurement locations closest to the beam’s edge are found to be strongly contaminated with reflections atlow-excitation frequencies. Thus, these locations are discarded in the analysis. Furthermore, measurement locations arepaired such that the distance between them is at least 5 − in . , so that the impact of measurement location uncertainty ondispersion curves estimates is minimized. Figure 7 shows the results obtained for the flexural wave mode based on thedata-driven models developed using the out-of-plane FRFs.In Figure 7, the data-driven estimates of the dispersion curves are presented along with the analytical dispersioncurve.The analytical curve follows the Timoshenko beam’s characteristic equation given by EI ρ A k − IA (cid:18) + EG κ (cid:19) k ω − ω + ρ IGA κ ω = , (19)11 preprint - J anuary
5, 2021Figure 6: A flow chart of the Data-driven-model-based dispersion relations algorithm. In the figure, F (0 , t ) is externallyapplied external force. u i ( x i , t ) is simulated response at location x i , and U i ( x i , f ) is the vector of Fourier coe ffi cientscorresponding to u ( x i , t ) .where E (70 . × Pa ) is the Young’s modulus of the beam, κ is the Timoshenko shear coe ffi cient, k is the theoreticalradial wave number, ρ (2650 kg / m ) is the volume density, and G ( = E / + ν )) is the shear modulus and ν (0 . ω = kc p gives the analytical dispersion relationship of the beam, in terms of thephase velocity c p . Then, the theoretical group velocity ( c g ) is given by c g = d ω dk . Figure 7 compares the experimental(data-driven) estimates of group-velocity ( V G ) against numerical estimates ( c g ).The red-colored curves plotted in Figure 7 are the mean estimate of the group-velocity of flexural waves from data-drivenmodels. The 95% confidence intervals are also provided for these estimates. The error between the numerical estimates(blue curve) and the data-driven estimates (red curve) is presented via the green curves. The four data-driven modelsresult in four group-velocity curves and four error curves. All four data-driven estimates are close to the analytical curvewith a maximum error of 6%. In particular, the data-driven models weighted with the weak inverse FRF and the hybridweight are moderately better at estimating dispersion curves at lower frequencies, with a maximum error under 4%.However, all models appear to deviate from the theoretical value in the frequency range over 40 kHz . This observationin this frequency range agrees with the behavior of the out-of-plane FRFs. Previous notion that the structure’s responseover 40 kHz deviates from a one-dimensional beam is supported by this result. In summary, among all four models, thedata-driven model using hybrid weights has the best estimate for the dispersion curve.For comparison purposes, the frequency domain analysis, descried by (18), is applied to the experimentally measuredwaveforms. These waveforms are obtained with the transient response measurement experiment discussed in Section2. Experimentally measured dispersion curves of the beam under test, corresponding to the first anti-symmetric wavemode, are then obtained. The results are depicted in Figure 8. As one can notice, dispersion curves obtained using theproposed data-driven approach (in Figure 7) and the conventional experimental transient-response-based approach (asseen in Figure 8), are in very good agreement. Experimentally measured dispersion curves are also found to agree wellwith the analytical predictions. With the transient-response-based approach, the high-frequency portion of the dispersioncurves has not been measured accurately as reflections from the beam sides strongly contaminated the response. This isthe same frequency range where plate-like modes are observed in the steady-state FRFs. Reflections from the beam endsare found to hinder the measurement of the low-frequency portion of the dispersion curve, which is also the case withthe proposed data-driven approach. While being an inevitable limitation with the conventional transient-response-basedapproach, the proposed data-driven approach allow for a number of solutions to be implemented to compensate forreflections.As for the first symmetric (longitudinal) wave mode, an alternative estimate of the dispersion curve is obtained usingthe well-known relationship between resonance frequencies and wave speed. For a free-free beam undergoing axialdeformations, the elementary rod theory predicts the n th natural frequency of the beam, ω n , to be ω n = n π cL , n = , , , ... (20)where c is the wave speed of the first symmetric wave mode and L is the length of the beam under test. The classical rodtheory predicts the fist symmetric wave mode to be nondispersive, thus wave speed is frequency independent. Basedon the natural frequencies identified over the frequency range of 0 − kHz , 24 estimates of the wave speed can be12 preprint - J anuary
5, 2021
Frequency [kHz] G r oup V e l o c i t y [ k m / s ] Analytical ApproachData-driven, Mean EstimateData-driven, 95% CI
10 20 30 40
Frequency [kHz] E rr o r % (a) Unit weight Frequency [kHz] G r oup V e l o c i t y [ k m / s ] Analytical ApproachData-driven, Mean EstimateData-driven, 95% CI
10 20 30 40
Frequency [kHz] E rr o r % (b) Weak inverse FRF weight Frequency [kHz] G r oup V e l o c i t y [ k m / s ] Analytical ApproachData-driven, Mean EstimateData-driven, 95% CI
10 20 30 40
Frequency [kHz] E rr o r % (c) Coherence weight Frequency [kHz] G r oup V e l o c i t y [ k m / s ] Analytical ApproachData-driven, Mean EstimateData-driven, 95% CI
10 20 30 40
Frequency [kHz] E rr o r % (d) Hybrid weightFigure 7: Comparison out-of-plane dispersion curves based on four data-driven models (a) Unit weight, (b) Weakinverse FRF weight, (c) Coherence weight, and (d) Hybrid weight.obtained. The results are depicted in Figure 9. While a direct wave-speed-estimate based on resonant frequenciesprovide accurate results, this technique works only when boundary conditions are well-defined. Uncertainty in boundaryconditions will reflect on the predicted wave speed. The current work experimentally investigates the performance of a novel data-driven approach in estimating dispersioncurves from steady-state FRFs. Experimental data of a free-free beam is adopted as a use-case scenario for this study. Atest setup with a Polytec SLDV provides the in-plane and the out-of-plane FRFs for developing data-driven models.13 preprint - J anuary
5, 2021
Frequency [kHz] G r oup V e l o c i t y [ k m / s ] Analytical ApproachExp., Transient Approach
10 20 30 40
Frequency [kHz] E rr o r % Figure 8: Comparison between experimentally measured dispersion curves and analytical dispersion curve. The greencurve displays the error between the two curves.
Frequency [kHz] G r oup V e l o c i t y [ k m / s ] Analytical ApproachData-driven, Mean EstimateData-driven, Rod-theoryData-driven, 95% CI
Frequency [kHz] E rr o r % Figure 9: Comparison of data-driven dispersion estimates, dispersion estimates based on rod-theory and analyticaldispersion curve. The green curve displays the error between the data-driven dispersion estimate and the analyticalvalue.Vector fitting the FRFs resulted in two data-driven models: an out-of-plane model and an in-plane model. Thesemodels are qualitatively assessed by comparing the fit at the poles and zeros of the FRFs. Four weighting functionsare considered for improving the quality of the fitted models. Finally, when the FRF data is weighted with a hybridweighting function, the resulting out-of-plane model with 368 poles has a relative fitting error in the order of O (10 − ).In contrast, the final in-plane model has a similar relative error of order O (10 − ) with 416 poles.In the next step, the transient response of burst-tone inputs is simulated using the data-driven models. Analyzing thesimulated transient response produced the out-of-plane estimates of the phase and group velocities. These data-drivenestimates are in good agreement with analytical and experimental (transient) velocities. The out-of-plane data-drivenmodels weighted with the weak inverse FRF and the hybrid weight estimate dispersion curves with a maximum error ofunder 4%. Additionally, the data-driven estimates are closer to the analytical values than the transient-experimentalestimates. The high-frequency portion of the dispersion curves has not been measured accurately as reflections from thebeam sides strongly contaminated the response. This is the same frequency range where plate-like modes are observedin the steady-state FRFs. Reflections from the beam ends are found to hinder the measurement of the low-frequency14 preprint - J anuary
5, 2021portion of the dispersion curve, which is also the case with the proposed data-driven approach. While being an inevitablelimitation with the conventional transient-response-based approach, the proposed data-driven approach allows for manysolutions to be implemented, such as the introduction of artificial damping. This will be further investigated in a futurework.
Acknowledgment
Tarazaga would like to acknowledge the support provided by the John R. Jones III Faculty Fellowship. The work ofGugercin was supported in parts by NSF through Grant DMS-1819110 and DMS-1720257.
References [1] M. I. Albakri, V. V. S. Malladi, S. Gugercin, and P. A. Tarazaga. Estimating dispersion curves from frequencyresponse functions via vector-fitting.
Mechanical Systems and Signal Processing , 140:106597, 2020.[2] R. Allemang and D. Brown. A complete review of the complex mode indicator function (cmif) with applications.In
Proceedings, international conference on noise and vibration engineering (ISMA), Katholieke UniversiteitLeuven, Belgium , volume 38, pages 36–44, 2006.[3] R. J. Allemang and D. Brown. A unified matrix polynomial approach to modal identification.
Journal of soundand Vibration , 211(3):301–322, 1998.[4] B. D. O. Anderson and A. C. Antoulas. Rational interpolation and state variable realization.
Linear Algebra andIts Applications, Special Issue on Matrix Problems , 137 / Interpolatory methods for model reduction . Computational Scienceand Engineering 21. SIAM, Philadelphia, 2020.[6] A. C. Antoulas, S. Lefteriu, and A. C. Ionita. A tutorial introduction to the Loewner framework for modelreduction. In
Model Reduction and Approximation , pages 335–376. SIAM, Philadelphia, 2017.[7] M. Berljafa and S. Güttel. The RKFIT algorithm for nonlinear rational approximation.
SIAM J. Sci. Comput. ,39(5):2049–2071, 2017.[8] J.-P. Berrut and L. N. Trefethen. Barycentric Lagrange interpolation.
SIAM review , 46(3):501–517, 2004.[9] J. Berthaut, M. Ichchou, and L. Jezequel. K-space identification of apparent structural behaviour.
Journal ofSound and Vibration , 280(3-5):1125–1131, 2005.[10] A. Carracedo Rodriguez and S. Gugercin. The p-AAA algorithm for data driven modeling of parametric dynamicalsystems. arXiv e-print arXiv:2003.06536, 2020.[11] A. Chaigne. Structural acoustics and vibrations. In
Springer Handbook of Acoustics , pages 941–1000. Springer,2014.[12] A. Chinea and S. Grivet-Talocia. On the parallelization of Vector Fitting algorithms.
IEEE Transactions onComponents, Packaging and Manufacturing Technology , 1(11):1761–1773, 2011.[13] F. Ciampa and M. Meo. A new algorithm for acoustic emission localization and flexural group velocity deter-mination in anisotropic structures.
Composites Part A: Applied Science and Manufacturing , 41(12):1777–1786,2010.[14] J. F. Doyle.
Wave propagation in structures . Springer, 1989.[15] Z. Drmaˇc, S. Gugercin, and C. Beattie. Quadrature-based vector fitting for discretized H approximation. SIAM J.Sci. Comp. , 37(2):A625–A652, 2015.[16] Z. Drmaˇc, S. Gugercin, and C. Beattie. Vector fitting for matrix-valued rational approximation.
SIAM Journal onScientific Computing , 37(5):A2346–A2379, 2015.[17] M. El-Kafafy, B. Peeters, T. Geluk, and P. Guillaume. The mlmm modal parameter estimation method: A newfeature to maximize modal model robustness.
Mechanical Systems and Signal Processing , 120:465–485, 2019.[18] M. El-Kafafy, B. Peeters, and P. Guillaume. Incorporating real modes and reciprocity constraints in the iterativemlmm modal parameter estimation method. In
Proceedings of ISMA2016 international conference on noiseand vibration engineering and USD2016 international conference on uncertainty in structural dynamics , pages2933–2946, 2016.[19] D. J. Ewins.
Modal testing: theory and practice , volume 15. Research studies press Letchworth, 1984.15 preprint - J anuary
5, 2021[20] D. Formenti and M. Richardson. Parameter estimation from frequency response measurements using rationalfraction polynomials (twenty years of progress). In
Proceedings of International Modal Analysis Conference XX .Citeseer, 2002.[21] G. H. Golub and V. Pereyra. The di ff erentiation of pseudo-inverses and nonlinear least squares problems whosevariables separate. SIAM Journal on Numerical Analysis , 10(2):413–432, 1973.[22] I. V. Gosea and S. Güttel. Algorithms for the rational approximation of matrix-valued functions. arXiv e-printarXiv:2003.06410, 2020.[23] S. Grivet-Talocia and B. Gustavsen.
Passive macromodeling: Theory and applications , volume 239. John Wiley& Sons, 2015.[24] P. Guillaume, P. Verboven, S. Vanlanduit, H. Van Der Auweraer, and B. Peeters. A poly-reference implementationof the least-squares complex frequency-domain estimator. In
Proceedings of IMAC , volume 21, pages 183–192. AConference & Exposition on Structural Dynamics, Society for Experimental . . . , 2003.[25] B. Gustavsen. Improving the pole relocating properties of vector fitting.
IEEE Transactions on Power Delivery ,21(3):1587–1592, 2006.[26] B. Gustavsen and A. Semlyen. Rational approximation of frequency domain responses by vector fitting.
IEEETransactions on power delivery , 14(3):1052–1061, 1999.[27] J. S. Hall and J. E. Michaels. A model-based approach to dispersion and parameter estimation for ultrasonicguided waves.
The Journal of the Acoustical Society of America , 127(2):920–930, 2010.[28] J. Hokanson.
Numerically stable and statistically e ffi cient algorithms for large scale exponential fitting . PhDthesis, Rice University, 2013.[29] J. M. Hokanson. Projected nonlinear least squares for exponential fitting. SIAM J. Sci. Comput. , 39(6):A3107–A3128, 2017.[30] J. M. Hokanson and C. C. Magruder. Least squares rational approximation. arXiv e-print arXiv:1811.12590, 2018.[31] H. Jeong and Y.-S. Jang. Wavelet analysis of plate wave propagation in composite laminates.
Composite Structures ,49(4):443–450, 2000.[32] D. Karachalios, I. V. Gosea, and A. C. Antoulas. The loewner framework for system identification and reduction.In
Model Reduction Handbook: Volume I: System-and Data-Driven Methods and Algorithms . De Gruyter, 2020.[33] J.-R. Lee, S. Y. Chong, H. Jeong, and C.-W. Kong. A time-of-flight mapping method for laser ultrasound guidedin a pipe and its application to wall thinning visualization.
NDT & E International , 44(8):680–691, 2011.[34] S. Lefteriu and A. C. Antoulas. On the convergence of the vector-fitting algorithm.
IEEE transactions onmicrowave theory and techniques , 61(4):1435–1443, 2013.[35] G.-R. Liu and Z. Xi.
Elastic waves in anisotropic laminates . CRC press, 2001.[36] N. M. M. Maia and J. M. M. e Silva.
Theoretical and experimental modal analysis . Research Studies Press, 1997.[37] A. Mayo and A. Antoulas. A framework for the solution of the generalized realization problem.
Linear Algebraand Its Applications , 425(2-3):634–662, 2007.[38] J. J. Moré. The Levenberg-Marquardt algorithm: implementation and theory. In
Numerical analysis , pages105–116. Springer, 1978.[39] Y. Nakatsukasa, O. Sète, and L. N. Trefethen. The AAA algorithm for rational approximation.
SIAM Journal onScientific Computing , 40(3):A1494–A1522, 2018.[40] Y. Nakatsukasa and L. N. Trefethen. An algorithm for real and complex rational minimax approximation.
SIAMJournal on Scientific Computing , 42(5):A3157–A3179, 2020.[41] J. Nocedal and S. Wright.
Numerical optimization . Springer Science & Business Media, 2006.[42] S. Rosalie, M. Vaughan, A. Bremner, and W. Chiu. Variation in the group velocity of lamb waves as a tool for thedetection of delamination in glare aluminium plate-like structures.
Composite Structures , 66(1-4):77–86, 2004.[43] C. Sanathanan and J. Koerner. Transfer function synthesis as a ratio of two complex polynomials.
IEEE Trans.Autom. Control , 8(1):56–58, 1963.[44] M. R. Souza, D. Beli, N. S. Ferguson, J. R. d. F. Arruda, and A. T. Fabro. A bayesian approach for wavenum-ber identification of metamaterial beams possessing variability.
Mechanical Systems and Signal Processing ,135:106437, 2020. 16 preprint - J anuary
5, 2021[45] B. Van Damme and A. Zemp. Measuring dispersion curves for bending waves in beams: a comparison of spatialfourier transform and inhomogeneous wave correlation.
Acta Acustica united with Acustica , 104(2):228–234,2018.[46] J. Vayssettes and G. Mercère. New developments for experimental modal analysis of aircraft structures. In
MATEC Web of Conferences , volume 20, page 01001. EDP Sciences, 2015.[47] J. Vayssettes, P. Vacher, and G. Mercère. An iterative algorithm for modal analysis based on structured matrixfractions.
IFAC Proceedings Volumes , 45(16):518–523, 2012.[48] S. M. Ziola and M. R. Gorman. Source location in thin plates using cross-correlation.