Forward/Backward Decomposition for Dispersive Wave Propagation Measurements

Abstract

Two complicating features commonly found in wave propagation applications include dispersion, i.e., frequency-dependent propagation velocity, and reflections, which introduce coherent noise. In this work, we present a signal processing technique which can be applied in a variety of applications to decompose signals into their forward- and backward-propagating components. The theory is presented, along with algorithmic implementation and experimental validation on a Timoshenko beam. The implications and potential utility of the method are briefly discussed.

1. Introduction

Propagating waves appear in a variety of engineering and science applications, ranging from radar and communications [1] to cosmology [2], seismology [3], acoustics [4], and structural mechanics [5]. In particular, measuring the propagation of vibration waves is popular in structural engineering applications. Contrasting with the widely used modal analysis approach to vibration testing, which is highly dependent on the boundary conditions of a structure or test article, wave propagation in principle is agnostic to boundary conditions: a propagating wave is unaffected by boundaries until it reaches them. Applications in which the boundary conditions are not of interest are therefore often prime candidates for wave propagation-based testing methods. Such applications in structural engineering include damage localization, material characterization, and parameter estimation, which are all examples seen in structural health monitoring (SHM). Very similar techniques can be applied in other settings as well, such as for footstep localization in smart structures [6] or target detection in radar [1]. Generally, a wide variety of different types of waves are used in various applications, ranging from electromagnetic waves to acoustic waves and ultrasonic waves, each having different characteristics. Nonetheless, the mathematical basis for analyzing propagating waves is often similar despite the different applications and wave properties.

One frequent limitation of existing analysis techniques is the requirement to isolate the initial propagating wave signal, known as the incident wave. Often, reflections of the incident wave (from boundaries or other objects in the propagation path) are the primary hindrance in isolating the incident wave. Consequently, wave propagation techniques typically use short burst signals and/or high-frequency signals with short wavelengths to aid in distinguishing the incident wave. Furthermore, the most common analysis techniques are based on leading or trailing edge detection and are only suitable for nondispersive waves which maintain their shape as they propagate (dispersive waves feature frequency-dependent propagation velocities; this causes multi-frequency signals such as pulses to not retain their shape as they propagate). Since some types of propagating waves become dispersive in certain frequency ranges, this imposes a further restriction on the frequencies which can be used in testing. Flexural waves in beams, for example, become dispersive at low frequencies, precluding the possibility of using low frequencies for non-destructive testing with traditional methods. However, for certain applications, such as rail neutral temperature (RNT) and stress measurement in rails [7], it is desirable to use lower frequencies given their enhanced sensitivity to the parameters of interest. Even in ultrasonic applications, there are examples in the literature where dispersion and reflections push researchers to use higher frequencies than they might otherwise use [8]. For nondispersive signals, it is relatively straightforward to compensate for reflections by using time-shifts inversely proportional to the propagation velocity to subtract the reflected signal components. Such methods are not applicable to dispersive waves though, which have frequency-dependent propagation velocities. On the other hand, dispersion compensation in the absence of reflections [7,9] has been addressed in the literature, but to the authors’ knowledge, the issue of compensating for reflections in dispersive waves has until now been unaddressed.

In this work, we present a signal processing technique for separating forward- and backward-propagating signal components in a set of dispersive wave measurements. Of particular interest is the case where the incident wave is not separable from reflections in the time domain in all of the measurements. The technique presented here is based on the Fourier decomposition of the signals into the frequency domain, where different wave modes’ Fourier coefficients superimpose according to simple vector addition. With a sufficient number of measurements, the various vectors corresponding to different individual wave modes can be decomposed at each frequency. After presenting the theory in Section 2, an algorithmic implementation of the technique is presented in Section 3, along with practical considerations to account for random noise in experiments. We then present experimental validation of the proposed approach on a clamped–clamped rectangular beam in Section 4 before providing concluding remarks in Section 5.

2. Theory

The typical experimental setup for wave propagation measurements consists of an actuator and a set of sensors providing measurements ${\mathbf{u}}_i$ at locations ${\mathbf{x}}_i$ as a function of time t, such as in Figure 1. Depending on the sensor type, the measurements ${\mathbf{u}}_i$ may be displacement, velocity, or acceleration measurements.

In this work, we will consider one-dimensional wave propagation, such that the wave measurements take the form $u(x_i,t)$. Furthermore, we consider systems for which a measurement of the response consisting of N time samples over a period of T can be decomposed using spectral analysis using the discrete Fourier transform as

(1)$u(x,t)=\sum _{n=0}^{N-1}{\widehat{u}}_n\left(x\right)e^{i{\omega }_{n}t},$

where ${\widehat{u}}_n$ is the $n\mathrm{th}$ Fourier coefficient at frequency ${\omega }_n={\displaystyle \frac{2\pi n}{T}}$. Generally, the Fourier coefficients ${\widehat{u}}_n$ can be further decomposed as

(2)${\widehat{u}}_n\left(x\right)={\widehat{U}}_n\widehat{G}\left(k_{n}x\right),$

where ${\widehat{U}}_n$ is called the amplitude spectrum and $\widehat{G}\left(k_{n}x\right)$ is called the phase spectrum for a particular wave mode. The phase spectrum $\widehat{G}\left(k_{n}x\right)$ represents a transfer function defining how the wave mode specified by the amplitude spectrum ${\widehat{U}}_n$ manifests at different locations x. The quantity $k_n$ is known as the wavenumber corresponding to the frequency ${\omega }_n$. As written, Equation (2) represents the response containing a single wave mode; in reality, several wave modes may exist concurrently, which can be accounted for by introducing additional modes to the right-hand side of Equation (2) to form a superposition solution.

This representation of a signal is valid for a wide variety of systems which are governed by linear, separable PDEs [5]. In this paper in particular, we consider elementary beam theory for modeling the response of a linear elastic beam with constant cross-sectional and material properties, in which case the transfer function simplifies to $\widehat{G}\left(k_{n}x\right)=e^{-i{k}_{n}x}$; however, the formulation presented herein is also valid for other transfer functions, so long as the transfer function is given by a model. Nonetheless, for common structural elements such as beams (flexural waves) and bars (longitudinal or torsional waves), the simplification $\widehat{G}\left(k_{n}x\right)=e^{-i{k}_{n}x}$ is sufficient. With this, the general solution is given by the so-called spectral solution

(3)$u(x,t)=\widehat{U}{e}^{-i(kx-\omega t)},$

where for notational simplicity we have omitted the summations over n (corresponding to the different frequency components of the signal) and over the various wave modes which may exist, which we will define next.

Since the governing equations are linear, several solutions (wave modes) may exist which can be linearly superimposed to form the total response. In elementary beam theories such as the Euler–Bernoulli or Timoshenko beam theories, four solutions exist. The four wavenumber solutions are obtained by inserting the spectral solution into the governing equations and rearranging for k; in the Timoshenko beam theory, this yields

(4)$\begin{array}{c}\hfill k_1=\pm \sqrt{{\displaystyle \frac{-b+\sqrt{b^2-4ac}}{2a}}}\mathrm{and}{k}_2=\mp \sqrt{{\displaystyle \frac{-b-\sqrt{b^2-4ac}}{2a}}},\end{array}$

with $a=EI,b=F_0-\rho {\omega }^2{\displaystyle \frac{EI}{G\kappa }}-\rho {\omega }^{2}I,\mathrm{and}c={\rho }^2{\omega }^4{\displaystyle \frac{I}{G\kappa }}-\rho A{\omega }^2$. For the low frequency ranges considered in this work, the positive and negative solutions for $k_1$, which are real-valued, represent forward- and backward-propagating waves, respectively, whereas the positive and negative solutions for $k_2$, which are complex-valued, represent forward- and backward-evanescent waves, respectively. At very high frequencies, $k_2$ switches from imaginary to real, converting from an evanescent mode into a propagating mode; however, this switch occurs at frequencies significantly higher than those considered in this work, so we will refer to the components corresponding to $k_2$ as evanescent modes. For a more detailed discussion on the interesting phenomenon of evanescent waves, which represent spatially decaying (rather than time-decaying) oscillations, we again refer to [5].

The quantities represented here are the following: E the elastic modulus of the material, I the cross-sectional moment of inertia (i.e., the second moment of area), $F_0$ the axial loading on the beam, $\rho$ the density of the material, $\omega$ the angular frequency, G the flexural modulus, $\kappa$ the Timoshenko correction factor, and A the cross-sectional area. Note that the wavenumber solution depends on the angular frequency quadratically and not linearly, which results in dispersive wave behavior. The general response according to superposition is therefore

(5)$\begin{array}{cc}\hfill u(x,t)& =\left(\underset{\begin{array}{c}\mathrm{forward}\\ \mathrm{propagating}\end{array}}{\underbrace{\widehat{A}{e}^{-i{k}_{1}x}}}+\underset{\begin{array}{c}\mathrm{forward}\\ \mathrm{evanescent}\end{array}}{\underbrace{\widehat{B}{e}^{-i{k}_{2}x}}}+\underset{\begin{array}{c}\mathrm{backward}\\ \mathrm{propagating}\end{array}}{\underbrace{\widehat{C}{e}^{i{k}_{1}x}}}+\underset{\begin{array}{c}\mathrm{backward}\\ \mathrm{evanescent}\end{array}}{\underbrace{\widehat{D}{e}^{i{k}_{2}x}}}+\underset{\begin{array}{c}\mathrm{higher}-\mathrm{order}\\ \mathrm{modes}\end{array}}{\underbrace{\cdots }}\right)e^{i\omega t}.\hfill \end{array}$

Again, recall that this solution exists for every frequency component, though we have neglected the summation over each frequency ${\omega }_n$ for notational simplicity. Also note how Equation (5) corresponds to Equation (2) with four wave modes. For further reading on the spectral solution approach to wave propagation problems, we refer the reader to [5].

Neglecting noise and higher-order terms from nonlinear responses, each measurement $u(x_i,t)$ of the structure’s response measured by a sensor at a particular location $x_i$ is a snapshot of the general solution for the location $x_i$:

(6)$\begin{array}{cc}\hfill u(x_i,t)& =\left(\widehat{A}{e}^{-i{k}_1{x}_i}+\widehat{B}{e}^{-i{k}_2{x}_i}+\widehat{C}{e}^{i{k}_1{x}_i}+\widehat{D}{e}^{i{k}_2{x}_i}\right)e^{i\omega t}.\hfill \end{array}$

Thus, taking the Fourier transform of a single measurement at location $x_i$ using the Fast Fourier Transform (FFT) algorithm gives (for every frequency) the quantity

(7)$\begin{array}{cc}\hfill \widehat{u}\left(x_i\right)& =\widehat{A}{e}^{-i{k}_1{x}_i}+\widehat{B}{e}^{-i{k}_2{x}_i}+\widehat{C}{e}^{i{k}_1{x}_i}+\widehat{D}{e}^{i{k}_2{x}_i}.\hfill \end{array}$

In practice, one receives a vector of all of the Fourier coefficients when using an FFT algorithm, and the corresponding frequencies can be constructed according to the time discretization parameters.

In Equation (7), the wavenumber solutions $k_1,k_2$ are known based on a model incorporating known material and geometric properties of the waveguide (here, we use the Timoshenko beam theory solutions in Equation (4)). The location of the measurement $x_i$ is also presumed to be measurable. Therefore, for each measurement, Equation (7) contains four unknown quantities: $\widehat{A},\widehat{B},\widehat{C},\widehat{D}\in \mathbb{C}$; these represent the amplitude spectra for the various wave modes which can be present in the structure. To solve for these four waves, four equations are needed, which can be obtained by taking four measurements at distinct locations. The equations given by all of the measurements form a linear system $\mathbf{Ax}=\mathbf{b}$, with 

(8)$\mathbf{A}=\left[\begin{array}{cccc}{e}^{-i{k}_1{x}_1}& e^{-i{k}_2{x}_1}& e^{i{k}_1{x}_1}& e^{i{k}_2{x}_1}\\ e^{-i{k}_1{x}_2}& e^{-i{k}_2{x}_2}& e^{i{k}_1{x}_2}& e^{i{k}_2{x}_2}\\ e^{-i{k}_1{x}_3}& e^{-i{k}_2{x}_3}& e^{i{k}_1{x}_3}& e^{i{k}_2{x}_3}\\ e^{-i{k}_1{x}_4}& e^{-i{k}_2{x}_4}& e^{i{k}_1{x}_4}& e^{i{k}_2{x}_4}\end{array}\right],\mathbf{x}=\left[\begin{array}{c}\widehat{A}\\ \widehat{B}\\ \widehat{C}\\ \widehat{D}\end{array}\right],\mathrm{and}\mathbf{b}=\left[\begin{array}{c}\widehat{u}\left(x_1\right)\\ \widehat{u}\left(x_2\right)\\ \widehat{u}\left(x_3\right)\\ \widehat{u}\left(x_4\right)\end{array}\right].$

A unique solution to this linear system exists if the matrix is nonsingular, which requires the locations $x_i$ and wavenumbers $k_1,k_2$ to be distinct. This system must be solved at every frequency ${\omega }_n$ corresponding to the $n\mathrm{th}$ Fourier coefficient ${\widehat{u}}_n\left(x_i\right)$ in each measurement $u(x_i,t)$. If the measurement locations $x_i$ are too close together (relative to the wavelengths/wavenumbers of interest), the matrix will be ill conditioned and close to singular; in practice, it is therefore necessary to space out the sensors appropriately such that the spatial resolution is sufficient to capture the largest wavelength or lowest frequency of interest.

The solution to the linear system given by Equation (8) consists of the unknown amplitude spectra $\mathbf{x}={\left[\begin{array}{cccc}\widehat{A}& \widehat{B}& \widehat{C}& \widehat{D}\end{array}\right]}^T$. These can be inserted back into the general solution in Equation (6) to reconstruct the full system response. However, the motivation of this work is to separate reflections in the measured signals. Thus, instead of reconstructing the full response, we can reconstruct the individual wave mode responses by simply omitting the others. In particular, the forward-propagating component which contains the incident wave is given by

(9)$\begin{array}{cc}\hfill u(x_i,t)& =\widehat{A}{e}^{-i{k}_1{x}_i}{e}^{i\omega t},\hfill \end{array}$

which is in practice obtained by performing an inverse Fourier transform on the quantity $\widehat{A}{e}^{-i{k}_1{x}_i}$.

3. Algorithmic Implementation

The system of equations in Equation (8) (brought forth by Equations (6) and (7)) forms the theoretical foundation of the method presented in this paper for eliminating reflections from experimental wave propagation data. In this section, we detail how we go about solving this system of equations, along with experimental considerations relating to implementing the technique in a robust manner.

In general, Equation (8) has four unknowns: two propagating waves, and two evanescent waves. While in general all four wave modes may exist, often in practice the two evanescent modes decay quickly in the spatial direction; thus, if sensors are placed sufficiently far from boundaries, it is permissible to only consider the two propagating waves. To be clear, the method still functions with the evanescent modes included, but they often contribute little to the system’s response and are of little interest to engineers compared with the propagating waves; thus, it is practical to neglect them so that there are half as many unknowns to solve for. As a result, each equation in the reduced system has the form

$\begin{array}{cc}\hfill \widehat{u}\left(x_i\right)& =\widehat{A}{e}^{-i{k}_1{x}_i}+\widehat{C}{e}^{i{k}_1{x}_i},i=1,2,\cdots ,\mathit{number}\mathit{of}\mathit{sensors},\hfill \end{array}$

where we have essentially imposed equations setting $\widehat{B}$ and $\widehat{D}$ to zero based on physical intuition for our experimental setup. In theory then, only two measurement locations are necessary to provide the two equations to solve for the two unknown amplitude spectra.

However, since noise and experimental uncertainty are common in test data, it is often desirable to use more measurements than are strictly required and solve in a least-squares sense. The theoretical formulation presented in Section 2 easily accommodates additional measurements at different locations as extra equations which overdetermine the system. Then, instead of solving the system $\mathbf{Ax}=\mathbf{b}$ exactly, we may solve the least-squares problem ${min}_{\mathbf{x}}\left|\right|\mathbf{Ax}-{\mathbf{b}\left|\right|}_2$. While our system is linear, future work is interested in an extension of the algorithm involving a nonlinear system of equations. We therefore employ Matlab’s built-in lsqnonlin() function, which uses an iterative optimization algorithm also suitable for nonlinear systems of equations, to solve the least-squares problem.

It is important to note a few details about the linear system that is to be solved. The first key detail is that the system is complex-valued, since the Fourier coefficients live in the frequency domain and are complex. The complex unknowns can be thought of as vectors in ${\mathbb{R}}^2$, which may be represented in either rectangular or polar coordinates, i.e., in terms of real and imaginary components or amplitude and phase in the complex plane. Attempting to solve this system in polar coordinates introduces solution uniqueness difficulties relating to the wrap-around at $2\pi$, which manifests as an ill-behaved non-convex objective function. Solving in rectangular coordinates eliminates this convergence issue.

Recall the general form for one of the measurement equations in our system of equations:

$\begin{array}{cc}\hfill \widehat{u}\left(x_i\right)& =\widehat{A}{e}^{-i{k}_1{x}_i}+\widehat{C}{e}^{i{k}_1{x}_i},i=1,2,\cdots ,\mathit{number}\mathit{of}\mathit{sensors}.\hfill \end{array}$

To represent this in rectangular coordinates, we employ Euler’s formula, which relates the complex exponential to rectangular coordinates: $e^{-ikx}=coskx-isinkx$. Additionally, the unknowns $\widehat{A},\widehat{C}\in \mathbb{C}$ can be represented as $\widehat{A}=\mathtt{Re}\left(\widehat{A}\right)+i\mathtt{Im}\left(\widehat{A}\right)$ and $\widehat{C}=\mathtt{Re}\left(\widehat{C}\right)+i\mathtt{Im}\left(\widehat{C}\right)$. Thus, $\widehat{A}{e}^{-i{k}_1{x}_i}$ in terms of real and imaginary components is

$\begin{array}{cc}\hfill \widehat{A}{e}^{-i{k}_1{x}_i}& =\left(\mathtt{Re}\left(\widehat{A}\right)+i\mathtt{Im}\left(\widehat{A}\right)\right)\left(cos{k}_1{x}_i-isin{k}_1{x}_i\right)\hfill \\ \hfill & =\underset{\mathtt{Re}\left(\widehat{A}{e}^{-i{k}_1{x}_i}\right)}{\underbrace{[\mathtt{Re}\left(\widehat{A}\right)cos{k}_1{x}_i+\mathtt{Im}\left(\widehat{A}\right)sin{k}_1{x}_i]}}+i\underset{\mathtt{Im}\left(\widehat{A}{e}^{-i{k}_1{x}_i}\right)}{\underbrace{[\mathtt{Im}\left(\widehat{A}\right)cos{k}_1{x}_i-\mathtt{Re}\left(\widehat{A}\right)sin{k}_1{x}_i]}}\hfill \end{array}$

and

$\begin{array}{cc}\hfill \widehat{C}{e}^{i{k}_1{x}_i}& =\left(\mathtt{Re}\left(\widehat{C}\right)+i\mathtt{Im}\left(\widehat{C}\right)\right)\left(cos{k}_1{x}_i+isin{k}_1{x}_i\right)\hfill \\ \hfill & =\underset{\mathtt{Re}\left(\widehat{C}{e}^{i{k}_1{x}_i}\right)}{\underbrace{[\mathtt{Re}\left(\widehat{C}\right)cos{k}_1{x}_i-\mathtt{Im}\left(\widehat{C}\right)sin{k}_1{x}_i]}}+i\underset{\mathtt{Im}\left(\widehat{C}{e}^{i{k}_1{x}_i}\right)}{\underbrace{[\mathtt{Im}\left(\widehat{C}\right)cos{k}_1{x}_i+\mathtt{Re}\left(\widehat{C}\right)sin{k}_1{x}_i]}}.\hfill \end{array}$

This allows us to represent our complex-valued (vector) equations as two scalar equations:

$\begin{array}{cc}\hfill \mathtt{Re}\left(\widehat{u}\left(x_i\right)\right)& =\mathtt{Re}\left(\widehat{A}\right)cos{k}_1{x}_i+\mathtt{Im}\left(\widehat{A}\right)sin{k}_1{x}_i+\mathtt{Re}\left(\widehat{C}\right)cos{k}_1{x}_i-\mathtt{Im}\left(\widehat{C}\right)sin{k}_1{x}_i\hfill \\ \hfill \mathtt{Im}\left(\widehat{u}\left(x_i\right)\right)& =\mathtt{Im}\left(\widehat{A}\right)cos{k}_1{x}_i-\mathtt{Re}\left(\widehat{A}\right)sin{k}_1{x}_i+\mathtt{Im}\left(\widehat{C}\right)cos{k}_1{x}_i+\mathtt{Re}\left(\widehat{C}\right)sin{k}_1{x}_i.\hfill \end{array}$

Subtracting the right-hand side from the left gives an error, which is fed into lsqnonlin() as the objective function to minimize.

The next important implementation detail is the use of bounds and an initial guess with the lsqnonlin() iterative solver. Since the system of equations must be solved at every frequency for the corresponding amplitude spectrum solution, we must loop over the frequency values. The first frequency component corresponds to the DC offset, which we set to zero to have a zero-mean signal. All of the subsequent frequencies instead use the previous frequency’s solution as the initial guess; in practice, the amplitude and phase for the signals used tend to be continuous, so using the previous solution as an initial guess starts the iteration close to the current solution.

Finally, an optional argument for the lsqnonlin() function is an upper and lower bound for each solution component; while not strictly necessary, since the magnitude of the amplitude spectrum is known from the input, we include solution bounds to restrict the solution space for good measure. A lower bound of zero helps to ensure that there is no confusion between forward- and backward-propagating components. For the upper bound, the magnitude of the amplitude spectrum of the input is normalized to the largest component in the measurement; this just provides a conservative upper bound on the energy that is possible in the signal based on the known input. This bound is applied to both the real and imaginary components in the solution variables.

After solving the nonlinear system of equations at every frequency, the solution is the amplitude spectra $\widehat{A},\widehat{C}$ corresponding to the forward- and backward-propagating waves. The pseudocode for the entire process for decomposing a set of wave propagation measurements into forward- and backward-propagating components is provided in Algorithm 1. Note that the algorithm does not require any pre-processing of the data in the form of windowing the response or filtering the response; preliminary investigations suggest that windowing the response imposes undesirable artifacts on the results. If filtering of certain frequencies is desired, they can simply be skipped in the solution process and set to zero, or a standard filter may be applied to the data before or after decomposing the signal into forward and backward components.

Algorithm 1: Forward/backward decomposition algorithm for decomposing dispersive signals in the presence of reflections.

[Figure omitted. See PDF]

4. Experimental Validation

To validate the performance of the wave decomposition technique on real-world data, an experiment was conducted using a rectangular aluminum beam with clamped ends. The test rig is a beam tensioning apparatus which has previously been used in [7]. The objective of the experiment was to validate the performance of the method with respect to its ability to identify reflections in the presence of other sources of random noise.

4.1. Experimental Setup

The test article, seen in Figure 2, is an aluminum 6061 rectangular beam with dimensions and material properties shown in Table 1.

The experimental apparatus, shown in Figure 2, was developed by Albakri et al. [7]. It consists of a steel beam with two end brackets featuring a tensioning mechanism. In this work, no tension was applied, so the test article featured simple clamped–clamped boundary conditions. The actuator used to provide excitation pulses was a Macro Fiber Composite (MFC) 2814 adhered near the center of the beam. By virtue of being adhered to the side of the beam, displaced from the neutral axis, and expanding/contracting primarily in the longitudinal direction of the beam, the MFC provides a bending moment that has been found to primarily excite the flexural propagating modes of the beam. The actuator also may excite some longitudinal (bar) modes and torsional modes due to any asymmetries in the placement of the MFC or the alignment of the clamps at the boundary, but these modes are small in amplitude and can be ignored.

A two-cycle sine wave with center frequency $\Omega$ and amplitude modulated by a Hanning window is used as the excitation signal. The time domain representation, along with the magnitude of the amplitude spectrum in the frequency domain, can be seen in Figure 3.

Center frequencies of 500 Hz and 1500 Hz were selected to demonstrate the algorithm’s performance at two different frequency ranges with different amounts of dispersion. The excitation signal was generated in MATLAB and output through an NI 9263 analog voltage output DAQ card. The NI DAQ used to send the excitation signal was housed in an NI cDAQ-9184 four-slot chassis, which was connected via Ethernet to a Windows 10 laptop with an AMD Ryzen 9 5900HS CPU, 16 GB Ram, and running MATLAB 2021a. A Trek 50/750 high-voltage amplifier with gain set to 15 V/V was used to amplify the signal to an amplitude of 150 Volts (zero to peak) for the actuator. This excitation level was found to provide a sufficiently strong excitation while avoiding pushing the actuator into its nonlinear regime.

A PSV-400 laser Doppler vibrometer was used to obtain non-contact vibration measurements to avoid mass-loading the structure with accelerometers. Forty scan points with a spacing of approximately 2 cm were used, providing measurements on an 80 cm region of the beam between the actuator and one boundary. A PSV-A-420 geometry scan unit was used to accurately measure the locations of the scan points on the beam, and the data acquisition was triggered by the excitation signal sent to the actuator.

The PSV-400 sampling frequency for data acquisition was adjusted according to the center frequency of the excitation signal in order to obtain at least 20 samples per waveform, which is sufficient to capture the frequency content of interest without over- or under-sampling. Since the structure contains little damping, it continues ringing after the initial wave and subsequent reflections pass, so the best total sampling time is not obvious. In the experiments presented here, the sampling time was arbitrarily selected to be 10 times the duration of the excitation pulse, and the response was not windowed or manipulated in any way. In order to minimize the effects of random noise, measurements were time-averaged by sending pulses several times and averaging the resulting responses. This assumes the response is repeatable, which is reasonable in this case since the excitation is identical and synced for each run. The laser vibrometer measurements were averaged 3 times for a $\sqrt{3}$ factor improvement in the signal-to-noise ratio.

4.2. Experimental Results and Discussion

We begin this section by showing some representative examples of the raw measurements. The 1500 Hz test case is shown on the left of Figure 4, and the 500 Hz test case is shown on the right. In each plot, the dark line shows the signal at the first measurement location, whereas the light gray signal is the measurement at the 40th scan point approximately 80 cm away. In both cases, the dispersion effect is apparent, as the toneburst changes shape as it propagates. In the 500 Hz case, the measurement at sensor location 1 exhibits relatively low dispersion, and the incident two-cycle Hanning windowed toneburst is even identifiable in the first sensor measurement. This is because the sensor spacing is relatively small compared to the long wavelength for the low-frequency waves; thus, the wave has had very little time to disperse. In contrast, although it technically is less-dispersive, the high-frequency case exhibits more dispersion already in the first measurement (note the qualitative difference in the pulse shape compared with the original toneburst). This is because, relative to the wavelengths of the frequency components of the signal, the first sensor’s location is already relatively far from the actuator. As a result, the wave has already had time to disperse, even if the incident wave is still identifiable. In both cases, the incident wave is not identifiable in the gray signal representing the last measurement location due to reflections from the boundary of the beam. For techniques which rely on incident wave isolation, these measurements would not be acceptable.

The results of executing the forward/backward decomposition proposed in this paper are shown in Figure 5, Figure 6, Figure 7 and Figure 8. Recall that the process performed by Algorithm 1 involves solving an inverse vector addition problem at each frequency to identify the forward- and backward-propagating components in each of the signals from the 40 scan points. Consider first the high-frequency case, in which we can clearly identify the incident wave and reflected waves in the unprocessed signals. The unprocessed signal from Figure 4 is plotted as the dotted gray line in Figure 5; the purple line is the forward-propagating component identified by the algorithm, whereas the yellow line is the backward-propagating component identified by the algorithm. Finally, the solid black line is the reconstructed signal obtained by combining the identified forward- and backward-propagating signals; ideally, the reconstructed signal should match the measured signal.

The algorithm successfully decomposes the measurements. First, we see that the incident pulse between 0 and 0.002 s is completely categorized as forward-propagating as desired. However, we also see other forward-propagating components, which we will proceed to identify. Consider that the actuator actually sends a pulse symmetrically forward along the structure and backward along the structure; hence, while we only measure the region in front of the actuator and only consider one of these pulses as the incident wave, we should expect to see two reflections (one from each boundary) of these two initial waves. At the location of sensor 1 depicted in Figure 5, which is near the center of the beam, these two reflections arrive at almost the same time (around 0.004 s), so in the time domain measurement it is not obvious that there are actually two signals present. Based on the measurements at other sensor locations though, the algorithm is able to identify and decompose these two reflections accordingly, with one being forward- and the other backward-propagating. At the location of sensor 40 depicted in Figure 6, which is near the boundary, we again see that what initially appear to be individual signal packets are actually various pulses superimposed and traveling in opposite directions. Consider the first pulse, between 0.002 and 0.004 s. As the forward-propagating incident wave reflects from the boundary, it will immediately turn into a backward-propagating reflection and return over the sensor near the boundary. Hence, in this first group, we see the incident wave and a first reflection. Then, at around 0.006 to 0.008 s, we see the same phenomenon, but with the initially backward-propagating incident wave. This wave has at this point traveled backwards toward the rear boundary, and then forward past the sensors where it is measured for the first time, and then reflected back from the front boundary, at which time it is registered for a second time.

An important note here is that every time a reflection is generated, it changes between forward- and backward-propagating, but the algorithm groups all of the signals into these two categories. Hence, the algorithm does not truly isolate the incident wave from all other signals; rather, it isolates forward-propagating components from backward-propagating ones. Further processing involving information about boundary conditions would be required to truly isolate the incident wave from reflections, which is outside the scope of this work.

Next, consider the low-frequency case in Figure 7 and Figure 8. Here, the incident and reflected waves are not clearly identifiable in the unprocessed time domain measurements plotted as the dotted gray lines. We emphasize that this undesirable characteristic typically prompts users to increase the excitation frequency in order to obtain separation between the incident wave and subsequent reflections, and that the motivation of this work is to provide an alternative. As in Figure 5, we can observe that the algorithm correctly identifies the incident wave at the location of sensor 1 shown in Figure 7, where it is mostly uncontaminated, whereas at the location of sensor 40 depicted in Figure 6 we see that there is a reflection from the end of the beam, which immediately contaminates the signal.

In the low-frequency case, because the wavelengths are so long, the separation in the time domain between “incident waves” and “reflected waves” even in the decomposed signals is diminished. Nonetheless, by going into the frequency domain, the algorithm is able to identify the portion of the signal corresponding to forward-propagating waves and backward-propagating waves.

Based on the intuition of the expected propagating signals, the algorithm results match our expectations. To more quantitatively assess the performance of the algorithm, we can consider the error between the reconstructed signals obtained by combining the identified forward- and backward-propagating components and the measured signals. We consider the relative error, defined as the norm of the absolute error between the reconstructed signal and the measured signal divided by the norm of the measured signal at a particular scan point located at $x_i$, as follows:

$\begin{array}{c}\hfill E_{\mathrm{rel}}\left(x_i\right)={\displaystyle \frac{\left|\right|u_{\mathrm{reconstructed}}(x_i,t)-u_{\mathrm{measured}}(x_i,t){\left|\right|}_2}{\left|\right|u_{\mathrm{measured}}(x_i,t){\left|\right|}_2}}.\end{array}$

Figure 9 depicts the relative errors for both the high- and low-frequency cases at all 40 scan points. The average relative error for the high-frequency case is $17.37\%$, whereas the average relative error for the low-frequency case is $17.71\%$, indicating that the identified signal components account for more than $80\%$ of the measurement data. The relative error metric is very sensitive to errors in phase, so note that even errors on the order of $20\%$ to $30\%$ appear to give decent reconstructions, as seen for example in Figure 5 and Figure 6. It is also noteworthy that the errors trend upwards towards the boundaries of the scan region; this is potentially due to the necessary truncation of the sampling time of the response, which results in some signals having entered the scan region and appearing in the outermost measurements but not having traversed the entire scan region or reaching the interior scan locations. The average error in the center of the scan region, for example, from points 15 to 30, is lower than average, by about $10\%$.

In observing the relative errors at the various scan points, we can see that the 40th scan point for the low-frequency case exhibits a high error compared with most of the other points. In Figure 8, we can see that the measured signal is quite small in amplitude, and similarly the reconstructed signal is small due to apparent destructive interference between the forward- and backward-propagating components. This causes the relative error to be quite high, since the norm in the denominator is relatively small. In Figure 10, we look instead at the 39th scan point for the low-frequency case, which has a relative reconstruction error of $16.75\%$. The reconstructed signal in Figure 10 is more similar to the measured signal, as suggested by the error metric.

5. Conclusions and Future Work

We have presented a method for decomposing wave propagation measurements into their forward- and backward-propagating components in a one-dimensional linear waveguide. The general theory and algorithm based on Fourier analysis were presented, with specific algorithmic implementation details for the case of dispersive vibration waves propagating in a Timoshenko beam. While existing methods in the literature typically require isolating incident waves from reflected waves, this requirement can be relaxed to identifying forward- and backward-propagating components using the proposed approach. This permits the use of lower-frequency signals, longer pulses with more localized frequency content, and measurement locations closer to reflection sources. The method was validated on experimental data from a clamped–clamped Timoshenko beam, for which the algorithm successfully identified the forward- and backward-propagating waves.

Direct applications for this work include non-destructive testing and SHM. By including the wavenumber $k_1$ as an unknown, the system of Equation (8) can be solved for the dispersion relation in addition to the forward- and backward-propagating waves, enabling characterization of material and geometry properties, such as stress. While the algorithm currently satisfactorily identifies the incident wave and its primary reflections, the current time sampling criteria are a known area of improvement. It is possible that introducing some artificial damping into the measurements to give more weight to the incident wave could improve the results. It is also possible that the technique could be built upon with additional information in order to detect reflection sources (i.e., damage localization) or further isolate the true incident wave from other forward-propagating components, although this has not been investigated.

Author Contributions

Conceptualization, N.A.C. and P.T.; methodology, N.A.C. and P.T.; software, N.A.C.; validation, N.A.C.; formal analysis, N.A.C.; investigation, N.A.C.; resources, P.T.; data curation, N.A.C.; writing—original draft preparation, N.A.C.; writing—review and editing, N.A.C. and P.T.; visualization, N.A.C.; supervision, P.T.; project administration, P.T.; funding acquisition, P.T. All authors have read and agreed to the published version of the manuscript.

Data Availability Statement

The original data presented in the study are openly available in the repository at https://doi.org/10.18738/T8/JQSOVI (accessed on 11 November 2024).

Acknowledgments

The authors acknowledge that this work was performed during their time at Virginia Polytechnic Institute and State University in the Mechanical Engineering and Biomedical Engineering and Mechanics Departments.

Conflicts of Interest

The authors declare no conflicts of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript; or in the decision to publish the results.

Abbreviations

The following abbreviations are used in this manuscript:

SHM	Structural Health Monitoring
RNT	Rail Neutral Temperature
FFT	Fast Fourier Transform
MFC	Macro Fiber Composite
PSV	Polytec Scanning Vibrometer

Footnotes

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Figures and Table

Figure 1. Diagram showing sensors along waveguide. The orange rectangle represents an actuator, and the black circles represent the sensors.

Figure 2. Experimental setup showing the view of the test article from the laser Doppler vibrometer. The scan region is labeled, along with stars highlighting the first and last scan points. The close-up pictures in call-outs depict, from left to right, the MFC actuator, two of the four strain gages, and the tensioning mechanism, which is present at both ends of the beam.

Figure 3. Time and frequency domain representations of the excitation signal.

Figure 4. Beam responses to the two-cycle toneburst with measurements at locations 1 and 40. (a) 1500 Hz center frequency toneburst. (b) 500 Hz center frequency toneburst.

Figure 5. Decomposed high-frequency (1500 Hz) beam responses showing the forward- and backward-propagating wave components for sensor 1. The relative error between the reconstructed signal and the measured signal is [Forumla omitted. See PDF.].

Figure 6. Decomposed high-frequency (1500 Hz) beam responses showing the forward- and backward-propagating wave components for sensor 40. The relative error between the reconstructed signal and the measured signal is [Forumla omitted. See PDF.].

Figure 7. Decomposed low-frequency (500 Hz) beam responses showing the forward- and backward-propagating wave components for sensor 1. The relative error between the reconstructed signal and the measured signal is [Forumla omitted. See PDF.].

Figure 8. Decomposed low-frequency (500 Hz) beam responses showing the forward- and backward-propagating wave components for sensor 40. The relative error between the reconstructed signal and the measured signal is [Forumla omitted. See PDF.].

Figure 9. Relative errors for the signal reconstructions at all 40 scan points. The average relative error for the high-frequency case (1500 Hz) is [Forumla omitted. See PDF.], and the average relative error for the low-frequency case (500 Hz) is [Forumla omitted. See PDF.].

Figure 10. Decomposed low-frequency (500 Hz) beam responses showing the forward- and backward-propagating wave components for sensor 39. The relative error between the reconstructed signal and the measured signal is [Forumla omitted. See PDF.].

Table 1

Test article material and geometry.

Property	Symbol	Value
length	L	1.75 m
width	W	3.81 cm
thickness	T	3.175 mm
cross-sectional area	A	$W·T$
cross-sectional moment of inertia	I	${\displaystyle \frac{1}{12}}W·T^3$
density	$\rho$	2700 kg/m3
elastic modulus	E	69 GPa
Poisson ratio	$\nu$	0.33
flexural modulus	G	${\displaystyle \frac{E}{2(1+\nu )}}$
Timoshenko correction factor	$\kappa$	${\displaystyle \frac{{\pi }^2}{12}}$
axial force	$F_0$	0 N

References

1. Meeks, M.L. Radar propagation at low altitudes. NASA STI/Recon Tech. Rep. A; 1982; 83, 24897.

2. Ezquiaga, J.M.; Zumalacárregui, M. Dark Energy in Light of Multi-Messenger Gravitational-Wave Astronomy. Front. Astron. Space Sci.; 2018; 5, 44. [DOI: https://dx.doi.org/10.3389/fspas.2018.00044]

3. Park, C.B.; Miller, R.D.; Xia, J. Imaging dispersion curves of surface waves on multi-channel record. Seg Tech. Program Expand. Abstr.; 1999; 17, pp. 1377-1380. [DOI: https://dx.doi.org/10.1190/1.1820161]

4. Albert, D.G.; Liu, L.; Moran, M.L. Time reversal processing for source location in an urban environment. J. Acoust. Soc. Am.; 2005; 118, pp. 616-619. [DOI: https://dx.doi.org/10.1121/1.1925849]

5. Doyle, J.F. Wave Propagation in Structures: Spectral Analysis Using Fast Discrete Fourier Transforms; 2nd ed. Mechanical Engineering Series; Springer: New York, NY, USA, 1997.

6. Poston, J.D.; Buehrer, R.M.; Tarazaga, P.A. Indoor footstep localization from structural dynamics instrumentation. Mech. Syst. Signal Process.; 2017; 88, pp. 224-239. [DOI: https://dx.doi.org/10.1016/j.ymssp.2016.11.023]

7. Albakri, M.I.; Malladi, V.S.; Tarazaga, P.A. Low-frequency acoustoelastic-based stress state characterization: Theory and experimental validation. Mech. Syst. Signal Process.; 2018; 112, pp. 417-429. [DOI: https://dx.doi.org/10.1016/j.ymssp.2018.04.011]

8. Shamshirsaz, M.; Bakhtiari-Nejad, F.; Khelghatdoost, M.; Asadi, S. Analytical and experimental analysis of Lamb wave generation in piezoelectrically driven Timoshenko beam. J. Intell. Mater. Syst. Struct.; 2015; 26, pp. 2314-2321. [DOI: https://dx.doi.org/10.1177/1045389X14554134]

9. Wilcox, P. A rapid signal processing technique to remove the effect of dispersion from guided wave signals. IEEE Trans. Ultrason. Ferroelectr. Freq. Control; 2003; 50, pp. 419-427. [DOI: https://dx.doi.org/10.1109/TUFFC.2003.1197965] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/12744398]