Benchmarking of methods for the identification of flexural wavenumbers inwooden plates
Federica MORANDI(1) , Olivier ROBIN(2) , Luca BARBARESI(3) , Massimo GARAI(4) , Noureddine ATALLA(5) , Nicolas
QUAEGEBEUR(6) , Patrice MASSON(7)(1) DIN, University of Bologna, Italy, federica.morandi6@unibo.it(2) GAUS, Université de Sherbrooke, Canada, olivier.robin@usherbrooke.ca(3) DIN, University of Bologna, Italy, luca.barbaresi@unibo.it(4) DIN, University of Bologna, Italy, massimo.garai@unibo.it(5) GAUS, Université de Sherbrooke, Canada, noureddine.atalla@usherbrooke.ca(6) GAUS, Université de Sherbrooke, Canada, nicolas.quaegebeur@usherbrooke.ca(7) GAUS, Université de Sherbrooke, Canada, patrice.masson@usherbrooke.ca

Abstract
This paper investigates and compares three methods for the identification of the flexural wavenumbers in thinlaminated wood panels. Laboratory measurements are performed on a plywood panel, an inhomogeneous andorthotropic plate, and the propagation characteristics of flexural waves are characterized along five radial directions. The methods considered for the analysis are the inhomogeneous wave correlation method, the Prony
method and the Time-of-Flight (ToF) method. The first two methods are implemented using the same data: adense line of points obtained with a scanning Laser Doppler Vibrometer (LDV) and a piezoelectric transducer
as an exciter. The ToF method uses a different data set, gathered using four accelerometers and an instrumented
impact hammer. The results show a remarkable match between the first two methods, while the third method
requires a precise time-domain filtering operation to provide accurate results. Nevertheless, the third method
provides a simple and low-cost alternative to the use of the LDV that can be more suited to in situ measurements.
The accuracy limits of each method, especially in terms of targeted frequency bandwidth are discussed.

Keywords: Wave propagation - Dispersion relation - Inhomogeneous wave correlation - Prony analysis

1

INTRODUCTION

The dispersion curve of a propagating wave into a medium relates the group and phase velocities to the frequency. The wavenumber k, inversely proportional to wavelength λ (k = 2π/λ ), depends on the angular frequency ω and the phase velocity c p following the relation k = ω/c p . Accessing the wave propagation characteristics can thus be achieved by performing wavenumber identification or phase speed measurement as a functionof frequency.
Wavenumber measurements in a one- or two-dimensional structure can be used for estimating material’s properties [1, 2], structural loss factor [3, 4], or for damage quantification [5]. Wavenumber estimations can be further
used for calculating a large number of vibroacoustic indicators (sound transmission loss, modal density, radiation efficiency). Wavenumber analysis is also of great interest when highly complex structures like sandwich [6]or ribbed panels [7] are under consideration. Another example is cross-laminated timber [2] which is made ofseveral successive layers, each usually perpendicular to the adjacent one (contrarily to glued-laminated timberthat has laminations all oriented in the same way).
Among the numerous existing wavenumber identification methods, few examples are (1) the iterative method
relying on sparse measurements [8]; (2) the use of laser-ultrasonic technique [9]; (3) the maximum likelihood
method [10]; (4) the spatial Fourier transform, SFT [11]; (5) the inhomogeneous wave correlation, IWC [1,1294

12]; (6) the Prony analysis [13]; (7) the spatial Laplace transform [14] and (8) the Time-of-Flight (ToF) or
phase difference method [15, 16]. Existing methods can be classified according to their frequency range ofapplicability, as was proposed in [17], to their robustness to test conditions, their ability to distinguish angledependent wavenumbers and finally with the final result they provide (real part of the wavenumber, imaginary
part or both). In [18], the SFT and IWC methods were compared in the case of the measurement of dispersion
curves in a beam made of cross-laminated timber and a polyamid beam with periodically varying thickness.
In this paper, the respective performances of IWC, Prony analysis and TOF method are compared in the caseof a wooden plate (an inhomogeneous and orthotropic plywood panel). Each method is briefly introduced inSec. 2, and their pros and cons are then discussed in Sec. 3 in terms of setup, post-processing complexity andresults precision as a function of a target frequency bandwidth.

2

TESTED WAVENUMBER IDENTIFICATION METHODS

2.1 Inhomogeneous Wave Correlation Method
The inhomogeneous wave correlation (IWC) method compares the signal acquired along a mono-dimensional
array with the propagation of an inhomogeneous wave travelling in the same direction õ(x, k) = eikx , k being
an unknown complex wavenumber with values that can be varied between the limits defined by the distance
between the measurement points [1].
The IWC function is defined as:| s̃ ( f , x) · õ∗ (x, k) dx|
R

IWC(k) = pR|s̃ ( f , x)|2 dx · |õ (x, k)|2 dx
R(1)

where s̃ ( f , x) is the transfer function measured at different points, and õ∗ (x, k) is the complex conjugate of theinhomogeneous plane wave described in the frequency domain. The integral in Eq. 1 is approximated by asummation for implementation purposes. For each frequency, the IWC function will reach a maximum whenthe measured signal better correlates with the inhomogeneous wave and therefore the estimated wavenumber canbe read at this point.
2.2 Prony Analysis
Prony Analysis (PA) is a parametric method that fits a sum of damped complex exponentials to a signal sampled at equally spaced data points. While Fourier series model a signal with undamped complex exponentials,
this method allows the estimate of the damping coefficients, besides frequency, phase and amplitude [19]. The
Prony’s method was implemented using the built-in Matlab R prony function, that computes the z-transform of
the matrix of the input signals s̃( f , x) as a ratio of polynomials. The numerator and denominator coefficientsare first calculated, the ratio of the two polynomials is then converted into a partial fraction expansion and theresidues and the poles are finally used to compute amplitude, complex wavenumber and damping of propagating
modes (each mode is described by different coefficients in different frequency ranges). The following procedure was followed: (1) the dispersion relations are calculated using Prony’s method; (2) a theoretical modelis used to identify the mode of interest, namely Rindel’s approximation to Mindlin’s theory for acousticallythick plates [15]; (3) the points belonging to the target mode are selected though the definition of a tolerance
range; (4) for each point of the vector k, the point that lies closer to the theoretical model is selected and theassociated coefficient is stored; (5) the same coefficient is used to select wavenumber, amplitude and damping.
2.3 Time-of-Flight (ToF) or phase difference method
When the ToF method is applied in the time domain [2], cross-correlation analysis between impulse responses
measured at a couple of receivers is typically used. The use of the Fourier transform of the cross-correlation [15,16], the Cross-Power Spectral Density (CPSD), is generally preferred and will be used in this work (the imaginary part of the CPSD of signals measured at two locations corresponds to their unwrapped phase difference).

1295

Figure 1. (a) Global view of the measurement setup; (b) Laser Doppler Vibrometer (LDV) measurements; (c)Impact hammer measurements.
Using an impact hammer as a transient source, and two accelerometers spaced by a distance r, the impulse
responses measured at the accelerometers are first windowed in the time domain so as to separate the direct
field contribution from reflections on the panel’s edges. Applying the CPSD and dividing the unwrapped phasedifference by r provides the real part of the bending wavenumber. The Matlab R built-in cpsd function has been
used, implementing Welch’s method. Compared to other methods for the ToF estimation, this method has theadvantage that it can be applied to impulse responses suitably windowed, and does not require the generationof single frequency bursts, which would require the use of a shaker or a piezo as an excitation source [2].

3

MEASUREMENTS

Measurements were performed at the laboratories of the Groupe d’Acoustique de l’Université de Sherbrooke(GAUS), Canada. The test plate is a 5-layer red spruce plywood panel, with dimensions 0.6 m x 1.22 m, 16 mm
thickness and mass density ρ = 480 kg/m3 . The low quality of the panel does not ensure the single layers
having a homogeneous thickness or mechanical parameters from section to section (the respective thicknessesof the layers are approximately 2.5/4/3/4/2.5, in mm). The orthotropic factor, calculated as the ratio of therotational inertia along the two principal directions, is 4.6. The surface layers have a high quality while theinner layers can have missing and non-missing knots. This low homogeneity can be reasonably overcome whenthe measurement mesh spans a dense line, while it can become critical when only a few measurement points
are considered.
A global view of the setup is provided in Fig. 1(a). Measurements were performed along five radial directions
uniformly distributed between 0 and π/2 radians. In the following, ‘0’ angle refers to the direction parallel tothe grain of the outer ply, having the highest rotational inertia, while ‘π/2’ angle will refer to the direction1296

characterized by the lowest rotational inertia, perpendicular to the first. The panel was kept in place using a
supporting frame with three clamps, and the panel’s edges were covered with a damping paste to minimize
reflections.
Two different measurement setups were used. Prony’s method and the IWC method are implemented on a setof data measured along lines. The acquisition of the out-of-plane velocity of the plate was made using a 3D
Polytec Laser Doppler Vibrometer. To enhance the Signal-to-Noise Ratio (SNR), small strips of reflective tapewere attached to the panel. Flexural waves were generated using a rectangular piezoelectric actuators embedded
in a resin mould with magnetic clamps [20], that helped the correct positioning on the plate (see Fig. 1(b)).
Measurements were performed in the frequency domain (resolution of 6.25 Hz, 70 averages) using a pseudo
random excitation. The measurement points were equally spaced (∼3 mm) spanning an overall length of 30
cm, corresponding to 101 measurement points.
A parallel set of measurements was conducted using an instrumented hammer for the generation of pulses andfour accelerometers as receivers, which were fixed on the plate at a distance of 0, 0.025, 0.125 and 0.30 mfrom the origin of the measurement array used for the laser vibrometer (Fig. 1(c)). The sampling frequency
was set to the highest available (65,536 Hz) so as to obtain the smallest time resolution. The real part of thewavenumber was then calculated using the ToF method previously described. Particular attention was paid tothe choice of the suitable time window that should be used to isolate the direct field from reflected waves.
Two different time windows were tested: (A) a Blackmann-Harris window with a length corresponding to thespan in which the signal keeps the same sign; (B) a force-exponential window (a flat portion spans the timelag between the beginning of the impulse response and the first change of sign occurring after the first local
maximum; a steep exponential then damps the signal).

4

RESULTS AND DISCUSSION

Figure 2 reports the results obtained using the three methods in the 300-10,000 Hz frequency range. It should be
noted that the Prony analysis revealed that only the A0 mode propagates in the frequency range 0-20 kHz. The
lowest frequency considered in the analysis is 300 Hz, due to the limitations arising from the finite dimensionsof the measurement mesh that leads to largely biased results (see zoomed sub-picture in the 0-500 Hz frequency
range in Figure 2 for the π/4 angle). The results obtained with the IWC and the Prony’s method give consistent
results over the whole considered frequency range. Besides a spurious dip displayed for the IWC method inalong the π/8 direction around the 2,000 Hz frequency, the two methods provide similar dispersion curves.
For the ToF method and for each radial direction, six dispersion curves are calculated that relate to the different
distances that can be obtained using the four measurement points. For each combination, the two time windows
were used, and their effect was monitored on all the data available. Therefore, a total of twelve dispersion
curves are obtained per each radial direction. In order to provide a first comparison among all data, the resultsof Prony’s method were taken as a reference and, for each dispersion curve, the relative error between theanalyzed curve and Prony’s reference curve was calculated and averaged in frequency. This operation was
repeated for each direction of propagation and the results were sorted out so as to identify the combination offactors (distance, window) that returned the lowest error in a selected frequency range. It is also noted thatin this preliminary study, only a single measurement is considered for each method even if several consecutive
tests were performed. The repeatability of each measurements will be studied in a subsequent study as well as
a comparison of mean relative error between Prony’s and IWC methods.
The results are reported in Table 1. The cells filled in green represent the lowest value of the average absolute
error, while the yellow cells contain values that are ≤ 5%. The miniature impact hammer is designed to havea linear response up to 20,000 Hz, but in the presented case, most of the spectral energy lies below 5,000 Hz.
For this reason, the metrics discussed above are evaluated over two different frequency ranges and the last one,300-5,000 Hz, gives consistent results for all the analyzed distances. The largest deviation from the reference
results are found for the the measurements performed along the 0 radians direction, see Table 1. Nevertheless,
it appears that provided a large separation is used between the two sensors, the phase difference method can

1297

Real wavenumber (rad/m)

150

0

100

50
Prony
IWC
ToF

0
300

1000

150

Real wavenumber (rad/m)

π/8

2000

3000

4000

5000

6000

7000

8000

9000

10000

200
100

100

0
0

250

500

50

0
300

1000

2000

3000

4000

5000

6000

7000

8000

9000

10000

2000

3000

4000

5000

6000

7000

8000

9000

10000

2000

3000

4000

5000

6000

7000

8000

9000

10000

150

Real wavenumber (rad/m)

π/4

100

50

0
300

Real wavenumber (rad/m)

150

1000

3π/8

100

50

0
300

1000

150

Real wavenumber (rad/m)

π/2

100

50
Prony
IWC
ToF

0
300

1000

2000

3000

4000

5000

6000

7000

8000

9000

10000

Frequency (Hz)Figure 2. Comparison between the results along the five considered angular directions obtained using Prony
Analysis, the IWC method and the ToF method (a separation distance of 0.275 m and a force-exponential time
window were used in this case).
1298

Table 1. Mean value of the relative error (%) calculated between results obtained with ToF analysis and Prony’s
method in two frequency ranges: 300-10,000 Hz and 300-5,000 Hz. Separation distances are sorted by increasing values: 0.025, 0.100, 0.125, 0.175, 0.275 and 0.300 m. The green cells indicate lowest frequency-averaged
relative errors values, while yellow cells the values the lie within a 5% tolerance.
Direction(radians)

Window

Distance 1
0.025 mDistance 20.100 mDistance 30.125 mDistance 40.175 mDistance 50.275 mDistance 60.300 m

12
13
7
6
3
2
2
2
15
17

8
7
2
1
5
6
5
3
7
6

6
5
2
1
7
7
6
3
7
7

23
22
3
3
3
2
2
2
9
7

12
10
2
1
1
1
2
1
1
1

9
8
2
1
1
1
2
2
1
1

300-10,000 Hz
0
π/8
π/4
3 π/8
π/2

A
B
A
B
A
B
A
B
A
B

31
21
8
6
24
26
11
12
17
26

13
14
11
10
12
15
15
7
41
44

12
14
9
9
14
17
14
7
36
40
300-5,000 Hz

0
π/8
π/4
3 π/8
π/2

A
B
A
B
A
B
A
B
A
B

45
35
10
9
6
12
6
7
5
7

11
14
4
4
4
4
5
5
13
10

15
17
4
4
5
5
5
5
11
9

provide very similar results to Prony’s method with reduced instrumentation.
The dispersion relations calculated considering separation distance no. 5 and the time window type B are
plotted in Fig. 2 together with the results estimated with Prony analysis and the IWC method. Though thebest fit between the three methods is clearly found up to 5,000 Hz, the results are shown in a wider frequency
range, from 300 to 10,000 Hz. Note that a zoomed plot is provided in Figure 3b in the ‘0’ radians direction
and in the 2,500-5,000 Hz frequency range, in which it appears that some localized differences can be foundbetween the three methods. When identified wavenumbers are plotted as a function of angle and frequency, theanisotropy of the plywood panel is clearly visible (see Figure 3a). The larger propagation speed (i.e. smaller
wavenumber values) might partly explain the large errors seen in this direction for the ToF.
Concerning computational time, the Prony analysis and the IWC method were implemented on the same set ofdata, and the calculation times were respectively 5 s and 366 s. Prony analysis is thus way much faster anddirectly provides information concerning all propagating modes in the structure, estimating k as a complex value.
For the IWC, estimating the imaginary part of the wavenumber can be troublesome especially for materials withsmall structural loss factor [1]. The computational time required for the computation of the CPSD is based ona different set of signals, and therefore cannot be directly compared to the two other methods, but the timerequired to evaluate the 6 dispersion relations (derived from the combination of source and receivers) is 0.8 s.

1299

Figure 3. Comparison between wavenumber identification results along the five considered angular directions
obtained using the three tested methods; (a) polar representation of the results at 2,500 and 5,000 Hz and (b)zoomed plot on the 2,500-5,000 Hz frequency range of the dispersion curves for the ’0’ radians direction.

5

CONCLUSIONS

The analysis presented in this work aims at evaluating three different methods for the identification of theflexural wavenumber in a plywood panel, and under which conditions they could be considered consistent. The
analysis was carried out through measurements performed in five radial directions. Two different measurements
setup were implemented: the first, using a piezoelectric actuator and a Laser Doppler Vibrometer scanning anarray of 101 points, and processing the data using the Inhomogeneous Wave correlation method and Prony
analysis. The second, using an impact hammer and four measurement points, and computing the ToF from thecalculation of the cross power spectral density of the signals, previously time-windowed.
Results show that the inhomogeneous wave correlation method and Prony analysis are generally in good agreement over the whole frequency range of analysis, and for all radial directions considered. The phase difference
method was analyzed through the comparison of the average error from Prony’s curve, that was taken as areference. For the ToF, the best results were obtained using the force-exponential window (even if no firm difference was found between the two suggested time windows) and using the largest separation distance between
accelerometers. In this case, an average calculation of the relative error in the frequency range 300-5,000 Hz
provides a value of 1% for most of the radial directions.
In the light of this preliminary analysis and if only the real part of the wavenumber is sought, the ToF method
appears to be a good alternative to Prony’s method provided that large separation distance is used between
sensors. If the complex wavenumber has to be identified, one should preferably choose Prony’s method because
it will lead to a clear reduction of computation time compared with IWC method. Future work will address an
exhaustive comparison between the three methods, including repeatability and reproducibility.

REFERENCES[1] Berthaut J., Ichchou M. N. and Jezequel L. K-space identification of apparent structural behaviour. Journalof Sound and Vibration 280, 1125–1131 (2005).
[2] Santoni A., Schoenwald S., Van Damme B. and Fausti P. Determination of the elastic and stiffness characteristics of cross-laminated timber plates from flexural wave velocity measurements, Journal of Sound andVibration 400, 387–401 (2017).

1300

[3] Rak M., Ichchou M.N. and Holnicki-Szulc J. Identification of structural loss factor from spatially distributed
measurements on beams with viscoelastic layer. Journal of Sound and Vibration, 310(4–5), 801–811 (2008).
[4] Cherif R., Chazot J-D and Atalla N. Damping loss factor estimation of two-dimensional orthotropic structures from a displacement field measurement. Journal of Sound and Vibration 356, 61–71, 2015.
[5] Mesnil O., Leckey C. and Ruzzene M. Instantaneous and local wavenumber estimations for damage quantification in composites, Structural Health Monitoring (2014).
[6] Sargianis J. and Suhr J. Core material effect on wave number and vibrational damping characteristics incarbon fiber sandwich composites, Composites Science and Technology 72 (13), 1493–1499 (2012).
[7] Ichchou M.N., Berthaut J. and Collet M. Multi-mode wave propagation in ribbed plates: Part I,
wavenumber-space characteristics, International Journal of Solids and Structures 45 (5), 1179–1195 (2008).
[8] McDaniel J.G. and Shepard W.S. Estimation of structural wavenumbers from spatially sparse response measurements. The Journal of the Acoustical Society of America, 108, 1674–1682 (2000).
[9] Kley M., Valle C., Jacobs L.J., Qu J. and Jarzynski J. Development of dispersion curves for two-layered
cylinders using laser ultrasonics. The Journal of the Acoustical Society of America 106 (2), 582–588 (1999).
[10] Halkyard C.R. Maximum likelihood estimation of flexural wavenumbers in lightly damped plates. Journalof Sound and Vibration, 300 (1-2), 217–240 (2007).
[11] Alleyne D. and Cawley P. A two-dimensional Fourier transform method for the measurement of propagating multimode signals. The Journal of the Acoustical Society of America, 89 (3), 1159–1168 (1991).
[12] Cherif R. and Atalla N. Experimental investigation of the accuracy of a vibroacoustic model for sandwichcomposite panels. The Journal of the Acoustical Society of America 137 (3), 1541–1550, (2015).
[13] Grosh K. and Williams E.G. Complex wave-number decomposition of structural vibrations. The Journal of
the Acoustical Society of America, 3 (2), 836–848 (1993).
[14] Geslain A., Raetz S., Hiraiwa M., Abi Ghanem M., Wallen S.P., Khanolkar A., Boechler N., LaurentJ., Prada C., Duclos A., Leclaire P., and Groby J.-P. Spatial Laplace transform for complex wavenumber
recovery and its application to the analysis of attenuation in acoustic systems. Journal of Applied Physics
120, 135107 (2016).
[15] Rindel J.H. Dispersion and absorption of structure-borne sound in acoustically thick plates. Applied Acoustics 41, 97–111 (1994).
[16] Roelens I., Nuytten F., Bosmans I. and Vermeir G. In Situ measurement of the stiffness properties ofbuilding elements, Applied Acoustics 52(3/4), 289–309 (1997).
[17] Margerit P., Lebée A., Caron J.-F. and Boutillon X. High Resolution Wavenmber Analysis (HRWA) for themechanical characterisation of viscoelastic beams, Journal of Sound and Vibration 433, 198–211 (2018).
[18] Van Damme B. and Zemp A. Measuring dispersion curves for bending waves in beams: a comparison of
Spatial Fourier Transform and Inhomogeneous Wave Correlation. Acta Acustica united with Acustica 104,
228–234 (2018).
[19] Castaniè F. (Editor). Digital Spectral Analysis - parametric, non-parametric and advanced methods. Wiley,
Hoboken, NJ, USA (2011).
[20] Quaegebeur N., Bouslama N., Bilodeau M., Guitel R., Masson P., Maslouhi A., and Micheau P. Guided
wave scattering by geometrical change or damage: Application to characterization of fatigue crack andmachined notch. Ultrasonics 73, 187–195 (2017).

1301