1   Introduction

Underwater acoustic channels are double-selective, multipath, and time-varying (Stojanovic and Preisig, 2009; Song, 2016). The multipath effect causes severe inter-symbol interference (ISI), making high-rate underwater acoustic communications (UACs) challenging. To eliminate ISI, channel equalization must compensate for the frequency selectivity caused by multipath. The commonly used equalizers in UAC are decision feedback equalizers (DFEs), which include direct adaptive DFE (Yang, 2004; Roy et al., 2007) and channel estimation-based DFE (CE-DFE) (Wang et al., 2012; Pajovic and Preisig, 2015). DFE is a compromise between complexity and performance.

Compared with DFE, maximum likelihood sequence estimation (MLSE) based on the Viterbi algorithm (VA) is the optimal equalization scheme (Forney, 1972; Viterbi, 2006; Barhumi and Moonen, 2009). However, applying this scheme in practical scenarios is difficult because the complexity of VA increases exponentially with the channel length. Thus, some promising algorithms have been proposed to reduce the complexity of MLSE in wireless communication. For example, the partitioned VA divides the higher-order modulated signal constellation graph into several subsets to reduce the number of states (Eyuboglu and Qureshi, 1988; Miller et al., 2001; Turner and Taylor, 2009). However, this approach only reduces the modulation order of the signal, and the complexity remains exponential with the channel length. The delayed decision feedback equalization was proposed to effectively reduce the number of states per moment and eliminate ISI in the tail by feedback (Duel-Hallen and Heegard, 1989; Wu and Sun, 1998; Gerstacker and Schober, 2002). However, this approach can lead to error propagation due to incorrect decision symbols. Another algorithm called the iterative MLSE detector has a bit error rate (BER) performance close to VA, and its complexity is quadratic in sequence and channel lengths (Myburgh and Olivier, 2008). However, the algorithm is more complex than VA when the channel length is short or the data sequence is long. To further reduce complexity, some algorithms have been proposed to reduce the survival paths. Sphere decoding (SD) was first applied to multiple-input multiple-output (MIMO) wireless communication systems to resolve inter-channel interference (Wang, 2011; Niu, 2014; Kurniawan et al., 2019). During the state trellis search, SD discards the surviving paths with a small likelihood and keeps only those with a large likelihood. SD can significantly reduce computational complexity, especially in an environment with a high signal-to-noise ratio (SNR) (Vikalo et al., 2006). Notably, computational complexity is greatly affected by the size of the sphere radius, which is not specified. If the chosen radius is too large, a large number of surviving paths will be retained. However, there will be no surviving paths if the chosen radius is too small.

Therefore, to make MLSE applicable to UAC, this paper proposes a two-step complexity reduction approach. First, a minimum mean square error (MMSE)-based channel shortening reduces the number of channel taps while maintaining a high level of complexity. The state of each moment is then reduced by selecting a few consecutive taps with the main concentrated energy for SD. To improve BER performance, we employ an iterative method to eliminate ISI caused by the remaining channel taps (Qu et al., 2010). The initial SD radius is provided by recursive least square (RLS-DFE). The radius is updated during the iterative process, and the radius of one iteration is determined by the hard decision output of the previous iteration.

The remainder of this study is organized as follows. Section 2 presents the system model and the proposed channel shortening and iterative SD in detail. Sections 3 and 4 discuss the simulation and experimental results of the proposed algorithm, respectively. Finally, Section 5 presents the conclusion.

2   Model and complexity-reduced MLSE methods

2.1   System model

For a single-input–single-output system model with a time-invariant channel, the received signal y can be expressed as follows:

where x represents the transmit signal with length N, H represents the Toeplitz matrix of the channel vector h with length L and dimension (N+L − 1) ×N, n is the additive Gaussian white noise with length N+L−1, and y is the received signal with length N+L − 1. Moreover, the channel matrix H can be expressed as follows:

The passband signal received from the hydrophone is first filtered using a bandpass filter to remove out-of-band noise and then processed using a down-mixing and matching filter. After the carrier frequency offset (CFO) estimation and compensation, the received baseband symbols are obtained by resampling. The baseband processing diagram is shown in Figure 1. The number of original equivalent baseband channel taps is reduced by channel shortening. After shortening, consecutive taps with the main concentrated energy are then selected for MLSE along with SD. The residual ISI is iteratively eliminated based on the previous hard decision output. Several iterations are required to produce a final demodulation output that represents the desired transmit sequence.

Figure  1  Baseband processing diagram

Download: Full-Size Img

2.2   Channel-shortening methods

2.2.1   Passive time reversal

Sea experiments show that passive time reversal (PTR) can effectively suppress the multipath of UAC channels (Rouseff et al., 2001). Before transmitting the message sequence x, a probe sequence needs to be transmitted for channel estimation. Assuming that the estimated channel impulse response (CIR) is $\hat{\boldsymbol{h}}_l$, and the real CIR is h_l, the received signal y′ after PTR can be expressed as follows:

where Q denotes the autocorrelation function of the channel $\hat{\boldsymbol{h}}_l$, and z denotes the noise interference. When CIR has excellent performance in autocorrelation, Q can be approximated as a delta function.

2.2.2   MMSE equalizer

Although the PTR processing is relatively simple, it is ineffective in some UAC channels with poor autocorrelation. To further improve channel autocorrelation, the superposition of multiple hydrophone signals is used to increase the energy ratio between the main and side lobes of the autocorrelation function Q (Rouseff et al., 2001). However, this method requires more receiving elements, leading to high hardware costs. The MMSE equalizer is another channel-shortening method (Melsa et al., 1996; Martin et al., 2004). It takes noise into account and is more generalizable than the PTR method. Figure 2 shows the MMSE shortening equalizer model.

Figure  2  MMSE shortening equalizer model

Download: Full-Size Img

In Figure 2, w denotes the coefficient vector of the shortening equalizer of length t, and b denotes the coefficient vector of the target impulse response (TIR) of length v+1, which is longer than t (Melsa et al., 1996). The main idea of this algorithm is to minimize the mean square error (MSE) between the outputs after the equivalent channel r and after the TIR d. Assume that the convolutional output of channel h and equalizer w is h_eff, the length of h_eff is (L + t−1), and b is the subvector of h_eff starting from the given index s (0 ≤ s ≤ L + t − v − 2). The optimal index s_opt will be explained at the end of this subsection. h_wall represents the remaining taps in h_eff with a length of (L+t–v−2). b and h_wall can be expressed as follows:

Using Eqs. (4) and (5), the energy of h_wall and b can be expressed as follows:

The following MSE and optimization problem are derived according to previous studies (Martin et al., 2004):

where R_n represents the autocorrelation matrix of the noise. Eq. (10) gives the definition of R_n. Here, E() denotes mathematical expectation.

The interpretation of the above optimization problem is to minimize the energy of h_wall for a certain energy of TIR considering noise impact.

The optimal w_opt satisfies as follows:

where λ is the smallest generalized eigenvalue of matrix A + R_n with respect to matrix B.

In a previous study (Melsa et al., 1996), the analytical solution of w_opt is outlined and expressed as follows:

where l_min is the eigenvector corresponding to the smallest eigenvalue of matrix $(\sqrt{\boldsymbol{B}})^{-1}\left(\boldsymbol{A}+\boldsymbol{R}_n\right)\left(\sqrt{\boldsymbol{B}}^{\mathrm{H}}\right)^{-1}$. The solution of w_opt involves the eigenvalue decomposition and inversion of two t-dimensional square matrices, multiplication operation between two t-dimensional square matrices, and multiplication operation between two t-dimensional square matrices and t-dimensional vectors. Therefore, the total number of complex multiplications required to calculate $\boldsymbol{w}_{\mathrm{opt}} \text { is } \frac{14}{3} t^3+2 t^2$.

The final optimal solution of the MMSE equalizer w_{opt_end} and the optimal value of index s_opt can be expressed as follows:

2.3   Iterative SD

Although the complexity is reduced after channel shortening, it still grows exponentially with the number of taps. To further reduce the complexity, we used SD to reduce the number of states.

The traditional MLSE is to maximize the probability.

where X is the set of all possible transmitted sequences. Assume that the noise conforms to the complex Gaussian distribution of zero mean, Eq. (14) can be expressed as follows:

As shown in Eqs. (14) and (15), the state S(k) at moment k is determined using the nearest L inputs.

In conventional MLSE, the number of states at moment k is M^L, and M is the modulation size. For example, in quadrature phase-shift keying (QPSK) modulation, M = 4. SD reduces the number of states because it only searches for lattice points within the radius (Figure 3), greatly reducing the computational complexity, especially at high SNR. The SD initial radius is provided by DFE. The radius for one following iteration is determined by the hard decision output of the previous iteration. Notably, the energy of the UAC channel is concentrated on a few dominant taps, which are selected for MLSE with SD. ISI caused by weak taps is eliminated iteratively.

Figure  3  Searching scheme of SD

Download: Full-Size Img

The iterative MLSE with SD can be summarized as follows:

Step 1: Select a few dominant consecutive taps for MLSE, and the cost function can be expressed as follows:

where L_b and L_f are the first and second selected tap positions, respectively, and L_b to L_f are consecutive.

Step 2: The MLSE combined with SD is used, and the radius r_k at moment k can be expressed as follows:

where the initial $\hat{\boldsymbol{x}}_{\text {prev }}$ is the DFE result, and $\hat{\boldsymbol{x}}_{\text {prev }}$ thereafter is the hard output of the previous SD. Furthermore, SD removes the path outside the radius at each moment, which further reduces the number of states. The method used here is a combination of VA and SD. For VA, there are a total of $M^{L_f-L_b+1}$ paths at moment $k$. When the radius is calculated, the convolutional output of all possible states and channels under the moment $k$ is calculated based on the radius $r_k$, and these paths are deleted if the sum of squared errors with the actual output is greater than the radius $r_k$. Finally, at moment $k$, the retained paths are less than $M^{L_f-L_b+1}$ paths. Therefore, at the moment $k+1$, the number of states is less than $M^{L_f-L_b+1}$, and we continue to use SD to delete the paths outside the radius $r_{k+1}$. The final retained paths are less than the moment $k$. At the last moment, we use the SD method of orthogonal triangular decomposition to find the path with the smallest metric. In summary, the SD used at moment $\mathrm{k}$ reduces the search paths at moment $k+1$. It can also avoid the complexity of VA searching $M^{L_f-L_b+1}$ paths at each moment. Figure 4 shows the process of combining SD and VA to reduce the complexity. The empty nodes indicate that the radius of the state is greater than the sphere radius, and the corresponding path is deleted.

Figure  4  SD combined with VA to reduce the complexity with four states

Download: Full-Size Img

The current demodulation output $\hat{\boldsymbol{x}}_{\mathrm{cur}}$ is the sequence when Eq. (15) is taken to be the minimum value.

Step 3: Update the receive vector y by eliminating ISI brought by the remaining weak taps. The ith element of the updated vector $\tilde{\boldsymbol{y}}$ is expressed as follows:

Step 4: MLSE combined with SD is used to demodulate the input sequence. Update $\hat{\boldsymbol{x}}_{\text {prev }}, \mathrm{MSE}_{\text {all }}$, and radius $r_k$ at moment $k$ as follows:

The demodulation outputs are as follows:

Step 5: Repeat Step 3 and Step 4.

2.4   Performance of SD combined with channel shortening versus SD

To apply the maximum likelihood method to UACs, we use the channel shortening and SD algorithms described in Subsections 2.2 and 2.3. Using only SD in severe multipath UACs would still have high complexity. Meanwhile, the SD combined with the channel-shortening method used in this paper further reduces the computational complexity and makes it possible to implement the near-maximum likelihood SD algorithm in hardware, which has good practicality. The two-step complexity reduction of SD combined with channel shortening has not been observed in previous studies. Moreover, this paper combines channel shortening with SD for the first time for MLSE at lower complexity.

Figures 5 and 6 show the BER performance and computational complexity of SD versus SD with channel shortening, respectively. It can be seen that the BER curves of SD and SD with channel shortening almost overlap; however, the computational complexity of SD with channel shortening decreases by a factor of 1 to 2.

Figure  5  BERs of SD with channel shortening and SD with a channel length of 7

Download: Full-Size Img

Figure  6  Average calculation time of SD with channel shortening and SD with a channel length of 7

Download: Full-Size Img

3   Complexity-reduced MLSE simulation results

In this section, simulations are performed to verify the performance of channel shortening and SD in MLSE.

To verify the performance of SD in reducing complexity, we choose a channel length of 5 that follows the Rayleigh fading distribution. The symbols are modulated by QPSK, and the channel is estimated using least squares. As classical equalization algorithms in UAC, RLS-DFE, improved proportionate normalized least mean squares DFE (IPNLMS-DFE), and CE-DFE are a compromise between complexity and BER performance (Stojanovic et al., 1993; Qin et al., 2020), which are used in this study for comparison. Since the channel length is only 5 and the computational complexity is acceptable for VA, channel shortening is not performed. The simulation result only verifies the performance of SD in reducing complexity. Figures 7 and 8 show the BER and normalized average calculation time comparison of the five algorithms, RLS-DFE, IPNLMS-DFE, CE-DFE, VA, and SD. When the SNR is greater than 16 dB, the BER performance of SD is almost the same as that of VA. The SNR required for SD is 3, 5, and 2 dB less than that of RLS-DFE, IPNLMS-DFE, and CE-DFE, respectively, at 10⁻³ BER. When the SNR is greater than 19 dB, the complexity of SD is close to that of RLS-DFE, IPNLMS-DFE, and CE-DFE.

Figure  7  BERs of RLS-DFE, IPNLMS-DFE, CE-DFE, VA and SD with a channel length of 5

Download: Full-Size Img

Figure  8  Normalized average calculation time of RLS-DFE, IPNLMS-DFE, CE-DFE, VA, and SD with a channel length of 5

Download: Full-Size Img

To validate the performance of channel shortening and SD in reducing complexity and BER, we choose a channel length of 7 that follows the Rayleigh fading distribution. The longer channel length makes VA computationally challenging; therefore, the channel is also shortened using the MMSE equalizer. After channel shortening, the channel taps are changed from 7 to 6. Figures 9 and 10 show the BER performance and normalized average calculation time comparison of RLS-DFE, IPNLMS-DFE, CE-DFE, VA combined with channel shortening, and SD combined with channel shortening without iteration. The results show that the BER performance of SD is almost the same as that of VA combined with channel shortening when the SNR is greater than 22 dB. The SNR required for SD combined with channel shortening is 1.5, 1.5, and 3.5 dB less than that of RLS-DFE, CE-DFE, and IPNLMS-DFE, respectively, at 10⁻⁴ BER. As the SNR increases above 30 dB, channel-shortening-based SD is twice that of RLS-DFE and one-tenth that of channel-shortening-based VA.

Figure  9  BERs of RLS-DFE, IPNLMS-DFE, CE-DFE, VA with channel shortening and SD with channel shortening and channel length of 7

Download: Full-Size Img

Figure  10  Normalized average calculation time of RLS-DFE, IPNLMS-DFE, CE-DFE, VA with channel shortening, and SD with channel shortening and channel length of 7

Download: Full-Size Img

To verify the performance of channel shortening and iterative SD in a real sea experiment channel, we choose the at-sea measured channel for simulation. Figure 11 shows the CIR for coherent communication in Zhairuoshan Island, Zhoushan, with a carrier frequency of 15 kHz and a bandwidth of 4 kHz. The multipath delay is approximately 6 ms; therefore, the number of channel taps in the baseband is 24. The channel is normalized by energy, and each channel tap is multiplied by a Gaussian distributed variable. After shortening by the MMSE equalizer, the channel lengths are changed from 19 to 10, and the energy of the multipath is mainly concentrated in seven consecutive taps (Figure 12). These seven consecutive channel taps are selected for iterative detection. Meanwhile, SD is used for MLSE to reduce computational complexity and make the algorithm implementable in hardware. Figure 13 shows the BER performance of RLS-DFE, IPNLMS-DFE, CE-DFE, and iterative SD. In the first iteration, SD performs worse than RLS-DFE and CE-DFE when the SNR is greater than 12 dB because there is a larger ISI not eliminated. When the number of iterations increases to 3 and 5, the BER performance of iterative SD is better than that of RLS-DFE and CE-DFE. Iterative SD requires 2 dB less SNR than RLS-DFE and CE-DFE at 10⁻⁴ BER.

Figure  11  CIR of sea experiment channel

Download: Full-Size Img

Figure  12  Original channel versus shortened channel

Download: Full-Size Img

Figure  13  BERs of RLS-DFE, IPNLMS-DFE, CE-DFE, and Iterative SD in an at-sea measured channel

Download: Full-Size Img

4   Sea experiment result

This section presents the sea experimental results for Zhairuoshan Island, Zhoushan, Zhejiang Province, on September 27, 2022. The transducer at the transmitting ship was submerged at a depth of 5 m, and the hydrophone at the receiving ship was submerged at a depth of 3 m. The distance between the transmitter and receiver was 2 000 m.

To increase the data rate in UAC, we used two bands, 8–12 kHz and 13–17 kHz, for data transmission simultaneously. Thus, the symbol rate was 8 000 symbols/s, and the sampling rate f_s of the passband signal was 96 kHz. The modulation scheme was QPSK. Pulse shaping was achieved using a root-rise cosine-matched filter with a roll-off factor of 0.2. The channel was estimated using orthogonal matching pursuit (Li and Preisig, 2007).

Figure 14 shows the structure of the data blocks in the sea experiment. We used a linear frequency modulation signal at the beginning of the data block for data synchronization. The gap between blocks was used to prevent inter-block interference, and the 511-M-sequence was used to estimate and compensate for the CFO of the packet. The m-sequence also served as a training sequence in RLS-DFE. Meanwhile, the m-sequence was used for channel estimation in channel shortening and iterative SD. A total of 3 840 information symbols follow the training symbols.

Figure  14  Block structure of transmitted data

Download: Full-Size Img

Figure 15 shows the CIR of 8–12 and 13–17 kHz in this experiment. The CIRs in Figures 15(a)–15(c) are uniformly normalized using the maximum amplitude of the 8–12 kHz band channel. There were approximately 15 taps in the 8–12 kHz channel, with the main energy concentrated on six continuous taps. Meanwhile, there were approximately 24 taps in the 13–17 kHz channel, with the main energy concentrated on 13 continuous taps. A faster fluctuation rate and stronger frequency selectivity were observed in the higher-band channel (Figures 15(a) and 15(b)). As shown in Figure 15(c), the dominant channel taps were reduced to eight after channel shortening by the MMSE equalizer.

Figure  15  CIRs in two frequency bands

Download: Full-Size Img

For the 8–12 kHz band, the complexity of seven taps used for iterative SD is acceptable; thus, channel shortening is not applied here. Table 1 presents the BER performance of RLS-DFE, IPNLMS-DFE, CE-DFE, and iterative SD in the 8–12 kHz band. At the initial iteration, the BER performance of iterative SD is worse than that of RLS-DFE because the ISI caused by weak channel taps is not eliminated. After five iterations, the average BER of SD is 0.0037, which is approximately half of the BER of RLS-DFE. Moreover, the average computation time per iteration of SD is approximately 2.5 times that of RLS-DFE. The experimental results demonstrate that iterative SD performed better than RLS-DFE, IPNLMS-DFE, and CE-DFE in the 8–12 kHz band.

Table 2 shows the BER performance of RLS-DFE, IPNLMS-DFE, CE-DFE, and iterative SD in the 13–17 kHz band. As shown in Figure 15(b), the BER performance in the 13–17 kHz band is worse than that in the 8–12 kHz band because the multipath in the 13–17 kHz band is more severe, and the channel varies faster. Furthermore, the higher-band channel exhibited a lower SNR because of the increased attenuation of high-frequency signals. Moreover, the average computation time per iteration of SD is approximately 6.2 times that of RLS-DFE. As presented in Table 2, the average BER performance of iterative MLSE with channel shortening and SD is slightly better than that of RLS-DFE, IPNLMS-DFE, and CE-DFE.

5   Conclusions

In this study, we propose an improved equalization scheme for MLSE based on the high complexity of MLSE in severe multipath channels to improve the BER performance. An MMSE equalizer is used for channel shortening. Several consecutive dominant channel taps are then selected for MLSE, and SD is adopted to reduce the number of surviving paths. Iterations are applied to the decoding process to eliminate ISI caused by weak channel taps. The simulation and experimental results demonstrate that the complexity of the proposed algorithm is greatly reduced compared with that of the traditional VA, and the BER performance is improved compared with that of RLS-DFE. The sea experiment shows that the BER performance is reduced from 0.007 6 to 0.003 7 in the 8–12 kHz band and from 0.151 6 to 0.114 5 in the 13–17 kHz band after five iterations.