# Forward/Backward Decomposition for Dispersive Wave Propagation Measurements
## 3 Algorithmic Implementation
Thus, \(\hat{A}e^{-ik_{1}x_{i}}\) in terms of real and imaginary components is  
\[\hat{A}e^{-ik_{1}x_{i}} =\big{(}\texttt{Re}(\hat{A})+i\texttt{Im}(\hat{A})\big{)}\big{(} \cos k_{1}x_{i}-i\sin k_{1}x_{i}\big{)}\] \[=\underbrace{\big{[}\texttt{Re}(\hat{A})\cos k_{1}x_{i}+\texttt{ Im}(\hat{A})\sin k_{1}x_{i}\big{]}}_{\texttt{Re}(\hat{A}e^{-ik_{1}x_{i}})}+i \underbrace{\big{[}\texttt{Im}(\hat{A})\cos k_{1}x_{i}-\texttt{Re}(\hat{A}) \sin k_{1}x_{i}\big{]}}_{\texttt{Im}(\hat{A}e^{-ik_{1}x_{i}})}\]  
and  
\[\hat{C}e^{ik_{1}x_{i}} =\big{(}\texttt{Re}(\hat{C})+i\texttt{Im}(\hat{C})\big{)}\big{(} \cos k_{1}x_{i}+i\sin k_{1}x_{i}\big{)}\] \[=\underbrace{\big{[}\texttt{Re}(\hat{C})\cos k_{1}x_{i}-\texttt{ Im}(\hat{C})\sin k_{1}x_{i}\big{]}}_{\texttt{Re}(\hat{C}^{\hat{A}_{1}x_{i}})}+i \underbrace{\big{[}\texttt{Im}(\hat{C})\cos k_{1}x_{i}+\texttt{Re}(\hat{C}) \sin k_{1}x_{i}\big{]}}_{\texttt{Im}(\hat{C}^{\hat{A}_{1}x_{i}})}.\]  
This allows us to represent our complex-valued (vector) equations as two scalar equations:  
\[\texttt{Re}(\hat{u}(x_{i})) =\texttt{Re}(\hat{A})\cos k_{1}x_{i}+\texttt{Im}(\hat{A})\sin k_ {1}x_{i}+\texttt{Re}(\hat{C})\cos k_{1}x_{i}-\texttt{Im}(\hat{C})\sin k_{1}x_ {i}\] \[\texttt{Im}(\hat{u}(x_{i})) =\texttt{Im}(\hat{A})\cos k_{1}x_{i}-\texttt{Re}(\hat{A})\sin k_ {1}x_{i}+\texttt{Im}(\hat{C})\cos k_{1}x_{i}+\texttt{Re}(\hat{C})\sin k_{1}x_ {i}.\]  
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 \(\hat{A},\hat{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.