Estimation of flexural wavenumbers in the presence of multiple wave components Milena Watanabe Bavaresco 1 Faculty of Engineering and Physical Sciences University of Southampton, Southampton, SO17 1BJ, UK Neil Ferguson 2 Institute of Sound and Vibration Research University of Southampton, Southampton, SO17 1BJ, UK Claus Hessler Ibsen 3 Vestas aircoil, Smed Hansens Vej 13, 6940 Lem, Denmark Atul Bhaskar 4 Department of Mechanical Engineering University of Sheffield, Mappin Street, Sheffield, S1 3JD, UK

ABSTRACT The separation of closely spaced spectral components is a known problem in signal processing, be it

response data either in the spatial or temporal domain. Herein an experimental method based on the

maximisation of the correlation between measured data and an arbitrary plane wave field is

investigated for its capability in distinguishing closely spaced flexural wavenumbers. Limitations

arising from data sampling, windowing and differences in the amplitudes and phases of the two wave

components are considered. These processes are used along with numerical simulations as a basis

to establish the boundaries within which the method successfully detects the components present.

Results show that differences in phases between the wave components have a higher impact on the

method’s capability on resolving the two wavenumbers. Differences in the wave amplitudes are only

limiting when their ratio is comparable to that of the side lobe of the dominant wave in the k-space;

it is inherently related to the process of windowing and its resulting spectral leakage.

1 mwb1f19@soton.ac.uk 2 nsf@isvr.soton.ac.uk 3 ci@vestas-aircoil.com 4 A.Bhaskar@sheffield.ac.uk

1. INTRODUCTION

Identifying wave dispersion properties is a fundamental aspect in many vibro-acoustic applications. Within that framework, the problem of resolution, i.e., being able to discern between closely spaced spectral components, is a recurring one. Among the existing methods to estimate wavenumbers experimentally ([1-5]), the one by McDaniel and Shepard [3] is capable of detecting multiple coexisting wavenumbers considering that the different number of wave types that propagate in a structure is known and given as an input. It relies on numerically optimizing the wave properties of the assumed wave fields to best fit data obtained experimentally.

In a different approach, the empirical mode decomposition (EMD) [6] separates the different frequency components completely empirically and has been investigated for its capability on resolving different frequency components by Rilling and Flandrin [7]. Although it is often referred to in time-frequency conversion, it can be applied to the equivalent space-wavenumber problem.

The method proposed by Ferguson, et al. [8] relies on yet another approach. This method is based on the correlation between measured data and a plane wave with arbitrary wavenumber, from which maximal results are used to indicate the most suitable wavenumbers present between the trial range searched. This method has been used in several applications, such as variability studies by Souza, et al. [9] where it was applied to a metamaterial beam, in damage detection by El Masri, et al. [10] in reinforced concrete beams and in characterisation of coupled beam and plate structures by Thompson, et al. [11]. The latter two cited studies specifically involve applications where wavenumbers of multiple simultaneously travelling waves were estimated. Improvement to the directional resolution of this method has been previously proposed by Thite and Ferguson [12] through the use of singular value decomposition process.

In this paper, the correlation-based method by Ferguson, et al. [8] is investigated for its capability on discerning between closely spaced wavenumbers. The investigation of resolution in this correlation-based method was motivated by previous work of the same authors, where helicoidal grooves break the axisymmetry of a pipe and subsequently the formation of travelling wave components comprising the splitting of the bending modes, resulting in two branches of the flexural wavenumbers which are closely valued and thus difficult to separate experimentally.

Through a correlation method [8], it is possible to observe that there are test setup specifications which can interfere and affect the proper estimation of wavenumbers. These dependencies are studied through numerical simulations in this paper and variations on the wave properties, such as the different relative wave amplitudes and phases, are also taken into account. For this numerical investigation, a similar approach as the one used by Rilling and Flandrin [7] was implemented.

The following section provides an overview of the correlation-based method under investigation. A study follows on the effect of the measurement grid chosen for the test setup in the capability of the method to discern between closely spaced wavenumbers. The performance of the method is then assessed for different wave properties in the wave components, i.e. relative wave amplitudes and phases between the two existing wave components are simulated. A performance measure is defined for that purpose. The main goal of this investigation is to provide the future user of the method with information in regards to the conditions under which the method is able to identify the multiple existing wavenumbers prior to performing an experiment.

2. THE CORRELATION METHOD

This method was introduced by Ferguson et al. [8] and comprises in finding, for each frequency, the wavenumber value that provides the maximum absolute value of the complex correlation function:

𝑊(𝑘 𝑡𝑥 , 𝑘 𝑡𝑦 , 𝜔) = ∫ ∫ 𝑤(𝑥, 𝑦, 𝜔) ∞ −∞ 𝑒 −𝑗𝑘 𝑡𝑥 𝑥 𝑒 −𝑗𝑘 𝑡𝑦 𝑦 𝑑𝑥𝑑𝑦 ∞ −∞ , (1)

where the correlation of the measurements 𝑤(𝑥, 𝑦, 𝜔) with the wavefield 𝑒 −𝑖𝑘 𝑡𝑥 𝑥 𝑒 −𝑖𝑘 𝑡𝑦 𝑦 for a 2D - plate is evaluated. This can be assessed for any ( 𝑘 𝑡𝑥 𝑘 𝑡𝑦 ) pair and these are designated as the trial wavenumbers.

Considering discrete measurements locations inside a rectangular data window of the available measurement domain and propagation in only one direction, one can rewrite Equation 1 as:

𝐿

𝑊 ̂ (𝑘 𝑡 , 𝜔) ≈(

𝑁 )∑ 𝑤(𝑥 𝑖 , 𝜔) 𝑁 𝑖=1 𝑒 −𝑗𝑘 𝑡 𝑥 𝑖 , (2)

where 𝐿 corresponds to the length of the measured spatial array and 𝑥 𝑖 the spatial locations of the measurements.

This method can be implemented following the procedure:

1. Define an array 𝑘 𝑡 as the arbitrary wavenumbers that will be used as the search range for the

existing wavenumber(s) (the maximum value of 𝑘 𝑡 for an aliasing-free result is dependent on the spacing between measurements at 𝑥 𝑖 , as described in the next session). 2. Calculate the correlation 𝑊 ̂ (𝑘 𝑡 , 𝜔) for all frequency points of interest and for all measured

𝑤(𝑥 𝑖 , 𝜔) . 3. For each frequency 𝜔 , find the maximum value of the correlation 𝑊 ̂ (𝑘 𝑡 , 𝜔) and normalize

𝑊 ̂ (𝑘 𝑡 , 𝜔) at all 𝑘 𝑡 at that frequency. 4. Extract one or multiple 𝑘 𝑡 that generate maxima in 𝑊 ̂ (𝑘 𝑡 , 𝜔) at a given frequency 𝜔 .

On step 4, as this algorithm is implemented on Matlab, the function findpeaks is used to find the maxima and a threshold for the minimal permissible magnitude of the correlation is defined. This is done to avoid selecting peaks in k -space which are related to windowing and spectral leakage. When selecting the measurement points 𝑥 𝑖 , there are two resulting phenomena that influence the accuracy of the results from this correlation-based method: aliasing and spectral leakage. These influences are analysed next.

2.1. Aliasing and spacing between the measurements Aliasing is strongly related to the sampling spacing 𝛥𝑥 , and for an even distribution of the measurement points, i.e. when the sampling spacing is constant, it dictates the maximum range for the search for the estimated wavenumber as:

𝑘 𝑚𝑎𝑥 = 𝜋

𝛥𝑥 (3)

This is the maximum value of wavenumber that can be estimated through the correlation method for aliased-free results. This is equivalent to what occurs when using the spatial discrete Fourier Transform (SDFT) and has been reported in Ferguson, et al. [8].

2.2. Spectral leakage and the length of measurement grid The total length of the measurement grid is also shown to have an impact on the correlation method. Although it does not dictate the spectral resolution directly as it does for the SDFT – and thus results are improved in comparison to the SDFT [8] – it represents the length of the windowed signal to be converted to the wavenumber domain and consequently, it leads to leakage affecting the resolution of the method in the k -space.

The rectangular window is the one to produce the narrowest main lobe at the cost of having a higher side lobe amplitude when compared to the other windows. The width of the main lobe can be calculated as

𝐵 𝑘 = 4𝜋

𝐿 𝑀 (4)

and the peak sidelobe amplitude is expected to be -13 dB of that of the main lobe. For a single wave component, these two aspects are represented for an increasing value of the window length 𝐿 in Figure 1, where 𝐿 1 is an arbitrary length.

Figure 1 - Effect of measurement length 𝐿 on k-space: narrowing of the main lobe for longer lengths 𝐿 and the constant amplitude of the side lobe of -13 dB in comparison to the main one.

Because the resolution of the method is closely linked to spectral leakage, the measurement length is also investigated in the numerical simulations presented next. For this analysis, it is considered that the proper spacing between measurement points is being considered to avoid aliasing. 3. NUMERICAL INVESTIGATION FOR DIFFERENT RELATIVE WAVE AMPLITUDES AND PHASES

After having outlined the dependency of the method on the measurement grid chosen for the test setup, the differences between wave amplitudes and phases of two coexisting wave components are analysed. The main goal is to determine the conditions under which the method can discern between the two wave components given these differences in the wave properties. For that purpose, the following signal with two waves, which could represent two propagating waves in an infinite beam, is defined:

𝑤(𝑥, 𝜔) = |𝐴 1 |𝑒 𝑖𝜑 1 𝑒 −𝑖𝑘 1 𝑥 + |𝐴 2 |𝑒 𝑖𝜑 2 𝑒 −𝑖𝑘 2 𝑥 , (5)

with 𝐴 , 𝜑 and 𝑘 being the amplitude, phase and wavenumber with the two subscripts referring to each wave component respectively. The parameters of the second component are defined in relation to that of the first one, i.e. 𝐴 2 = 𝑟 𝐴 𝐴 1 and 𝑘 2 = 𝑘 1 −∆𝑘 , where the following ranges are set 𝑟 𝐴 ∈ ]0.1,1[ , 𝜑 2 ∈ ]0, 𝜋[ and ∆𝑘∈ ]0,15[. The amplitude 𝐴 1 is set to unity and the phase 𝜑 1 to zero. In this way, the properties of the second wave component are considered in relation to those of the first one.

The accuracy of the extracted wavenumbers is evaluated by comparing them to the existing wavenumbers through a quantitative performance metric:

‖𝑘 𝑡1 −𝑘 1 ‖ 2

‖𝑘 𝑡2 −𝑘 2 ‖ 2

1

|𝑘 2 | ) , (6)

𝑣=

2 (

|𝑘 1 | +

where 𝑘 1 and 𝑘 2 are the existing wavenumbers and 𝑘 𝑡1 𝑎𝑛𝑑 𝑘 𝑡2 are the ones estimated through the method. These parameters are sorted so that the highest estimated wavenumber is compared to the

highest existing one, etc. If the method does not find two peaks to allocate 𝑘 𝑡1 and 𝑘 𝑡2 , then both are set to zero and the value of 𝑣 is unity. If the two wavenumbers are correctly estimated, then the numerator is nulled and 𝑣 is zero.

3.1. Waves with distinct wavenumbers and different phases Keeping the two wave components with the same amplitude and varying the phase between them, the method is assessed through the performance measure from Equation 6. The results for a phase difference of 𝜋 radians are shown in Figure 2.

Figure 2 – Performance of the correlation-based method in resolving two wavenumbers for waves of the same amplitude and with a phase difference of π radians as a function of the measurement length 𝐿 and wavenumber spacing ∆𝑘 . Colour map: zero values indicate accurate estimates whereas unity indicates lack of resolution.

The colormap in Figure 2 indicates the regions where the method successfully discerns between the two wavenumbers ( 𝑣= 0 ) and where it does not ( 𝑣= 1 ). Two very distinct areas can be observed, and generally the wavenumbers are badly separated for small measurement grid lengths and for closely spaced wavenumbers. The dashed line in the same Figure represents the criterion ∆𝑘𝐿> ~8.4 introduced by Ferguson, et al. [8] as an approximation of the conditions where the method successfully separates the two wavenumbers. For this difference of phase of 𝜋 radians (Figure 2), this approximation matches the contour colormap, i.e. the numerical results match the criterion of ∆𝑘𝐿> ~8.4 .

a) b) c) d)

Figure 3 - Performance of the correlation-based method in resolving two wavenumbers for waves of same amplitude with phase difference of a)

3𝜋

𝜋

𝜋

4 rad, b)

2 rad, c)

4 rad, d) zero as a function of the measurement length 𝐿 and wavenumber spacing ∆𝑘 .

3𝜋

The colour maps of Figure 3 show the performance of the method for phase differences of

4 radians to zero, along with the analytical prediction in the dashed line for reference. An interesting phenomenon can be observed as the area of colour map where the method resolves the two wavenumbers increases with a reduction of the phase difference between the wave components. To visualize what occurs in the k-spectrum for a point in that domain, Figure 4 illustrates the correlation magnitudes at one specific point, where 𝐿= 1 m and ∆𝑘= 7.5 rad/m. These are considered at the

[wpe av Lim

3𝜋

same phases of

4 radians to zero as in Figure 3. The vertical lines represent the specified existing wavenumbers (dashed black line) and the estimated ones (continuous blue line). In Figure 4.a and Figure 4.b the peaks corresponding to the two existing wavenumbers merge into one, thus the correlation method is able to estimate only one wavenumber at a mean value between the two existing ones. In Figure 4.c and Figure 4.d, in contrast, the changes in the phase difference results in a k- spectrum with two visible peaks, two wavenumbers are thus estimated.

a) b) c) d)

Figure 4 – Scaled correlation magnitudes for two waves present of the same amplitude with phase difference of a)

3𝜋

𝜋

𝜋

4 rad, b)

2 rad, c)

4 rad, d) zero. The measurement length 𝐿= 1 m and wavenumber spacing ∆𝑘= 7.5 rad/m. Vertical lines indicate the existing wavenumbers (dashed black lines) and estimations (blue continuous lines).

Referring back to Figure 3, the case where the two wave components have the same phase is the most advantageous for the method, as a wider combination of ∆𝑘𝐿 can be solved. Graphically, the product of ∆𝑘𝐿> ~4.2 is sufficient for a proper estimation of the two existing wavenumbers under those conditions.

3.2. Waves with distinct wavenumbers and different amplitudes To analyse the effect of the relative amplitude between the two wave components only the two extreme phase differences from the previous section are considered, i.e. 𝜑 2 = 0 and 𝜑 2 = 𝜋 radians.

The plot in Figure 3.d is representative of the results for 𝜑 2 = 𝜋 radians for amplitude ratios 𝑟 𝐴 ranging from 1 (same amplitude) down to 0.3. At this lower limit, because the amplitude of the correlation 𝑊 ̂ (𝑘 𝑡2 ) at the smallest existing wavenumber is similar to that of the side lobes of the largest wavenumber, the method fails in estimating the two wavenumber components. The same occurs for the case where 𝜑 2 = 0 when the amplitude ratio is of 0.4 and below. 4. CONCLUSIONS

In this paper, the main parameters affecting the capability of the correlation-based method in resolving between two closely valued wavenumbers were investigated. By numerically simulating two waves with different relative phases and amplitudes, it was observed that the criterion of ∆𝑘𝐿> ~8.4 is very accurate in determining the region where the method is reliable for a phase difference of 𝜋 rad. For smaller phase differences, this criterion becomes conservative, and the method is capable of discerning between the wavenumber components at even smaller products of ∆𝑘𝐿> ~4.2 .

Different amplitude ratios for the wave components did not significantly affect the success of the method, except for the threshold when one of the amplitudes is approximately 0.4 or smaller than that of the other component. Below this threshold, the method does not accurately estimate the wavenumbers and leakage in the k -spectrum is found to be the main source of the problem.

Improvements to this performance can be attempted, from which the work by Thite and Ferguson [12] can serve as basis. A comparison of the performance of this method with some of the previously discussed ones might also be beneficial to aid in the selection of the method most appropriate for a given application. These subjects are suggested as future work.

5. ACKNOWLEDGEMENTS

The authors gratefully acknowledge the funding received from the European Union’s Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement No 765636. 6. REFERENCES

1. Muggleton, J. M., Brennan, M. J., Linford, P. W.: Axisymmetric wave propagation in fluid-filled

pipes: wavenumber measurements in in vacuo and buried pipes, Journal of Sound and Vibration, vol. 270, no. 1-2, pp. 171-190, (2004). 2. Berthaut, J., Ichchou, M., Jézéquel, L.: Wavenumbers identification in two-dimensional

structures: Comparison between two new methods, in Proceedings of the Tenth International Congress on Sound and Vibration, pp. 1155-1162, (2003). 3. McDaniel, J., Shepard, W.: Estimation of structural wave numbers from spatially sparse response

measurements, The Journal of the Acoustical Society of America, vol. 108, pp. 1674-82, (2000). 4. Boukadia, R., F., Claeys, C., Droz, C., Ichchou, M., Desmet, W., Deckers, E.: An inverse

COnvolution MEthod for wavenumber extraction (INCOME): Formulations and applications, Journal of Sound and Vibration, p. 116586, (2021). 5. Van Damme, B., 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, vol. 104, pp. 228-234, (2018). 6. Huang, N. E., Shen, Z., Long, S. R., Wu, M. C., Shih, H. H., Zheng, Q., Yen, N.-C., Tung, C. C.,

Liu, H. H.: The empirical mode decomposition and the Hilbert spectrum for nonlinear and non- stationary time series analysis, Proceedings of the Royal Society of London. Series A: Mathematical, Physical and Engineering Sciences, vol. 454, no. 1971, pp. 903-995, (1998). 7. Rilling, G., Flandrin, P.: One or Two Frequencies? The Empirical Mode Decomposition

Answers, IEEE Transactions on Signal Processing, vol. 56, no. 1, pp. 85-95, (2008). 8. Ferguson, N. S., Halkyard, C. R., Mace, B. R., Heron, K. H.: The estimation of wavenumbers in

two-dimensional structures, in Proceedings of Isma 2002: International Conference on Noise and Vibration Engineering, Vols 1-5, pp. 799-806, (2002). 9. Souza, M. R., Beli, D., Ferguson, N. S., de F. Arruda, J. R., Fabro, A. T.: A Bayesian approach

for wavenumber identification of metamaterial beams possessing variability, Mechanical Systems and Signal Processing, vol. 135, p. 106437, (2020). 10. El Masri, E., Ferguson, N., Waters, T.: Wave propagation and scattering in reinforced concrete

beams, The Journal of the Acoustical Society of America, vol. 146, pp. 3283-3294, (2019). 11. Thompson, D. J., Ferguson, N. S., Yoo, J. W., Rohlfing, J.: Structural waveguide behaviour of a

beam–plate system, Journal of Sound and Vibration, vol. 318, no. 1, pp. 206-226, (2008). 12. Thite, A., Ferguson, N.: Wavenumber Estimation: Further Study of the Correlation Technique

and Use of SVD to Improve Propagation Direction Resolution, (2004).