# Forward/Backward Decomposition for Dispersive Wave Propagation Measurements
## 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  
\[\hat{u}(x_{i})=\hat{A}e^{-ik_{1}x_{i}}+\hat{C}e^{ik_{1}x_{i}},\quad i=1,2, \ldots,\text{ }\text{\emph{number of sensors}},\]  
where we have essentially imposed equations setting \(\hat{B}\) and \(\hat{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{A}\mathbf{x}=\mathbf{b}\) exactly, we may solve the least-squares problem \(\min_{\mathbf{x}}||\mathbf{A}\mathbf{x}-\mathbf{b}||_{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:  
\[\hat{u}(x_{i})=\hat{A}e^{-ik_{1}x_{i}}+\hat{C}e^{ik_{1}x_{i}},\quad i=1,2, \ldots,\text{ }\text{\emph{number of sensors}}.\]To represent this in rectangular coordinates, we employ [PERSON]'s formula, which relates the complex exponential to rectangular coordinates: \(e^{-ikx}=\cos kx-i\sin kx\). Additionally, the unknowns \(\hat{A},\hat{C}\in\mathbb{C}\) can be represented as \(\hat{A}=\texttt{Re}(\hat{A})+i\texttt{Im}(\hat{A})\) and \(\hat{C}=\texttt{Re}(\hat{C})+i\texttt{Im}(\hat{C})\).