Forward/Backward Decomposition for Dispersive Wave Propagation Measurements

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.

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 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 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 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 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 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$: 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 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 

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

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 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: 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 and This allows us to represent our complex-valued (vector) equations as two scalar equations: 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.

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.

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.