Enhanced Mobile Multiple-Input Multiple-Output Underwater Acoustic Communications

Abstract

This paper focuses on mobile multiple-input multiple-output (MIMO) underwater acoustic communications (UAC) over double-selective channels subject to both intersymbol interference and Doppler scaling effects. Temporal resampling is implemented to effectively convert the Doppler scaling effects to Doppler frequency shifts. Under the assumption that the channels between all the transmitter and receiver pairs experience the same Doppler frequency, a variation of the recently proposed generalization of the sparse learning via iterative minimization (GoSLIM) algorithm, referred to as GoSLIM-V, is employed to estimate the frequency modulated acoustic channels. GoSLIM-V is user parameter free and is easy to use in practical applications. This paper also considers turbo equalization for retrieving the transmitted signal. In particular, this paper reviews the linear minimum mean-squared error (LMMSE) based soft-input soft-output equalizer involved in the turbo equalization scheme and adopts a fast implementation of the equalizer that achieves negligible detection performance degradation compared to its direct implementation counterpart. The effectiveness of the considered MIMO UAC scheme is demonstrated using both simulated data and measurements recently acquired during the MACE10 in-water experiment.

1. Introduction

Achieving reliable underwater acoustic communications (UAC) with high data rate is difficult owing to the unique challenges imposed by the underwater acoustic environment [1, 2]. In typical UAC, the difference in the propagation time between the earliest and latest arrivals could span tens to hundreds of symbol periods [3], which translates into long channel impulse response (CIR) and severe intersymbol interference (ISI) at the receiver side. Moreover, the presence of Doppler effects, owing to the relative motions between the transmitter and receiver platforms and the dynamic underwater acoustic medium, induces temporal scaling (stretching or compression) to the transmitted signals [4]. Doppler-induced scaling effects impair the reliability of UAC, especially in the case of a phase-coherent detection scheme [3]. Furthermore, the scarcely available bandwidth permitted by the acoustic channel imposes an upper bound on the attainable symbol rate [2]. Therefore, the pursuit of high data rate in UAC leverages the multiple-input multiple-output (MIMO) scheme, which offers enhanced reliability and/or increased data rates compared to its single-input counterpart [5–7]. The focus of the present paper is on effective mobile MIMO UAC over double-selective acoustic channels suffering from both ISI and Doppler scaling effects.

Converting the double-selective channel into an ISI channel via temporal resampling is an effective way to tackle mobile UAC difficulties [4]. Although the Doppler scaling effects can be largely mitigated via such a temporal resampling process, the residual Doppler still causes frequency shift on the received measurements. Coherent UAC requires the receiver to acquire knowledge of the underlying channel after temporal resampling via channel estimation [7]. Channel estimation could be conducted either in the training-directed mode, using known training sequences, or in the decision-directed mode, using the detected payload symbols [5, 6]. A preferable tool to characterize a channel subject to both ISI and Doppler frequency shift is the scattering function (SF), which essentially decouples the acoustic channel into a bank of paths that experience different delays and Doppler frequencies [8]. The major concern in SF-based channel estimation is that the problem becomes over parameterized with too many degrees of freedom. It is practically more beneficial to look for a channel model with the smallest number of parameters, but one that still sufficiently reflects the defining characteristics of the acoustic channel of interest. Along this line of thought, it is assumed in [9, 10] that, at each receiver, the channel taps for all the transmitters experience the same Doppler frequency, but different receivers experience different Doppler shifts. The number of unknowns in the frequency dimension, as a consequence, is significantly reduced. Under this assumption, CIRs and the underlying Doppler frequency could be estimated in a separate manner [9] or in a joint manner by employing the generalization of the sparse learning via iterative minimization (GoSLIM) algorithm [10]. It is demonstrated in [10] that GoSLIM outperforms the separate estimation algorithm proposed in [9] in terms of estimation accuracy and robustness against suboptimal training sequences.

In [11], the aforementioned channel model is further simplified by assuming that the channel taps for all the transmitter and receiver pairs experience the same Doppler frequency. As a consequence, the impact of the Doppler frequency shift on the received measurements across all the receivers is taken into account through one unknown common frequency. Accordingly, a variation of GoSLIM, referred to as GoSLIM-V (V stands for variation), is proposed for channel estimation in [11]. Like GoSLIM, GoSLIM-V addresses sparsity through a hierarchical Bayesian model, and because GoSLIM-V is user parameter free, it is easy to use in practical applications. It is demonstrated in [11] that the employment of GoSLIM-V not only reduces the overall complexity in the channel estimation stage but also slightly improves the detection performance compared to its GoSLIM counterpart. Due to this reason, GoSLIM-V is used as the channel estimation algorithm in the present paper.

Following the channel estimation is the design of the detection scheme for extracting the transmitted signals. The channel-induced phase shift should be first compensated out using the Doppler frequency estimate [9]. Such phase compensation task, along with the aforementioned temporal resampling process, effectively converts a double-selective channel subject to both Doppler scaling effects and ISI to an ISI channel, which allows for the employment of various equalization techniques that can effectively combat ISI. We use a linear minimum mean-squared error (LMMSE) based filter for symbol detection. In a MIMO setup, on top of ISI, multiple simultaneously transmitted signals act as interferences to one another. Therefore, interference cancellation scheme also plays a critical role in the overall detection performance. A hard decision based interference cancellation scheme, including vertical BLAST (V-BLAST) [12] and RELAX-BLAST [5], subtracts out the hard decisions of detected signals from the received measurements to aid the detection of the remaining signals. By combining V-BLAST with the cyclic principle of the RELAX algorithm [13], RELAX-BLAST provides superior detection performance over V-BLAST at the cost of slightly increased complexities [5, 6].

The detection performance can be further enhanced by employing a soft interference cancellation scheme, including turbo equalization [14–16]. For a receiver employing turbo equalization, both the equalizer and decoder involved are configured as soft-input soft-output. The detection performance improves as the soft information cycles between the equalizer and decoder. The main drawback of the turbo equalization scheme is the increased computational complexity compared to its hard-decision-based counterparts. To address this problem, we consider a low complexity approximation of soft-input soft-output equalizer [14]. We will show via numerical and experimental examples that the employment of the proposed approximate equalizer enjoys a computational complexity comparable to RELAX-BLAST and provides only slightly degraded detection performance compared to an exactly implemented turbo equalizer.

The rest of the paper is organized as follows. Section 2 presents a system outline. Section 3 describes a model for the acoustic channel subject to both ISI and Doppler scaling effects and reviews the temporal resampling procedure. Section 4 formulates the channel estimation problem in both training-directed and decision-directed modes and then introduces GoSLIM-V as the channel estimation algorithm. Section 5 first formulates the symbol detection problem and then details the LMMSE based soft-input soft-output equalizer and its low complexity approximation. Section 6 presents the simulation results of the turbo equalization scheme, followed by the experimental results obtained from analyzing the MACE10 in-water measured data. The paper is concluded in Section 7.

Notation

Vectors and matrices are denoted by boldface lowercase and uppercase letters, respectively, $\parallel ·\parallel$ denotes the Euclidean norm of a vector, $|·|$ is the modulus and $(·{)}^{*}$ is the complex conjugate of a scalar. $(·{)}^T$ and $(·{)}^H$ denote the transpose and conjugate transpose, respectively, of a matrix or vector, I denotes an identity matrix of appropriate dimension, and $\widehat{\mathbf{x}}$ denotes the estimate of x. ${\mathrm{diag}}_{}(\mathbf{v})$ represents a diagonal matrix in which the elements of v are on the diagonal. $\mathrm{Re}(·)$ and $\mathrm{Im}(·)$ represent the real and the imaginary components of a complex-valued scalar, respectively. $\mathbf{A}\otimes \mathbf{B}$ denotes the Kronecker product of two matrices A and B. Other mathematical symbols are defined after their first appearance.

2. System Outline

Consider an $N\times M$ mobile MIMO UAC system equipped with N transmit transducers and M receive hydrophones. The transmitted payload sequences are divided into multiple blocks, each of which is encoded separately. Figure 1(a) demonstrates the construction of a single payload symbol block (the construction of other blocks follows the same procedure). Denote $a(k)\in \{\mathrm{0,1}\}$ as the kth source bit for $k=\mathrm{1},\dots ,NK$ .  $\{a(k){\}}_{k=\mathrm{1}}^{NK}$ are first fed into a $\mathrm{1}/\mathrm{2}$ rate convolutional encoder with generator polynomials $(\mathrm{1}\hspace{0.5em}\mathrm{0}\hspace{0.5em}\mathrm{0}\hspace{0.5em}\mathrm{1}\hspace{0.5em}\mathrm{1})$ and $(\mathrm{1}\hspace{0.5em}\mathrm{1}\hspace{0.5em}\mathrm{0}\hspace{0.5em}\mathrm{1}\hspace{0.5em}\mathrm{1})$ . The encoded bits $\{b(k){\}}_{k=\mathrm{1}}^{\mathrm{2}NK}$ are then passed to a random interleaver, followed by a quadrature phase-shift keying (QPSK) modulation using Gray code mapping. In Figure 1, interleaver and deinterleaver modules are represented by Π and ${\Pi }^{-\mathrm{1}}$ , respectively. Next, the so-obtained QPSK payload symbols $\{x(k){\}}_{k=\mathrm{1}}^{NK}$ are demultiplexed into the N payload blocks, each consisting of K symbols, across the N transmitters in a round-robin fashion (this is where the “DEMUX” module in Figure 1(a) comes into play). More specifically, in our design, $x(n+(q-\mathrm{1})N)$ corresponds to the qth symbol sent by the nth transmitter, denoted as $x_n(q)$ , for $q=\mathrm{1},\dots ,K$ and $n=\mathrm{1},\dots ,N$ . Accordingly, we denote $c_n(\mathrm{2}q-\mathrm{1})$ and $c_n(\mathrm{2}q)$ as the two consecutive interleaved bits in $\{c(k){\}}_{k=\mathrm{1}}^{\mathrm{2}NK}$ that map to $x_n(q)$ according to the formula given below:

$$\begin{array}{l}{x}_n\left(q\right)=\frac{\mathrm{1}}{\sqrt{\mathrm{2}}}\left[j{\left(-\mathrm{1}\right)}^{c_n\left(\mathrm{2}q-\mathrm{1}\right)}+{\left(-\mathrm{1}\right)}^{c_n\left(\mathrm{2}q\right)}\right],\hfill \\ \phantom{sssssss}q=\mathrm{1},\dots ,K,\hspace{0.17em}\hspace{0.17em}n=\mathrm{1},\dots ,N.\hfill \end{array}$$

(1)

Since $c(k)\in \{\mathrm{0,1}\}$ , the support of $\{x_n(k)\}$ is a 4-element alphabet set $𝒮=\{(\pm j\pm \mathrm{1})/\sqrt{\mathrm{2}}\}$ .

Figure 1

An $N\times M$ MIMO UAC system. (a) Transmitter structure. (b) Receiver structure by employing turbo equalization.

The structure of a receiver employing a turbo equalization scheme is shown in Figure 1(b). The measurements acquired by the M receive hydrophones are first resampled, followed by channel estimation and phase compensation. After phase compensation, the double-selective channel is converted to an ISI channel, and the turbo equalization scheme is employed herein to retrieve the transmitted information. The superior detection performance promised by turbo equalization is mainly due to its mechanism of cycling soft information between the equalizer and the decoder [14]. Accordingly, turbo equalization consists of two key modules, namely, a soft-input soft-output equalizer and a soft-input soft-output decoder [17]. The soft information of a generic bit $a\in \{\mathrm{0,1}\}$ , commonly known as the log-likelihood ratio (LLR), is defined as

$$\begin{array}{cc}& L\left(a\right)=\ln \frac{P\left(a=\mathrm{0}\right)}{P\left(a=\mathrm{1}\right)},\end{array}$$

(2)

where $P(a=\mathrm{0})\in [\mathrm{0,1}]$ represents the probability of a being 0. As shown in Figure 1(b), the multiplexed and deinterleaved version of $\{L_e(c_n(k)){\}}_{n=\mathrm{1},\hspace{0.17em}k=\mathrm{1}}^{N\mathrm{2}K}$ , the a posteriori extrinsic LLR generated by the equalizer, forms the a priori inputs $\{L_e(b(k)){\}}_{k=\mathrm{1}}^{\mathrm{2}NK}$ to the decoder. Conversely, the interleaved and demultiplexed version of $\{L_d(b(k)){\}}_{k=\mathrm{1}}^{\mathrm{2}NK}$ , the a posteriori extrinsic LLR generated by the decoder, serves as a priori information to the equalizer. The subscript e or d reminds us that the LLRs are generated by the equalizer or the decoder, respectively. The soft information is cycled between the equalizer and the decoder multiple times before making hard decisions on the source bits. Note that the interleaver and deinterleaver involved at the transmitter and receiver have the same structure, whereas the “DEMUX” module inside the dashed rectangle in Figure 1(b) is different from that in Figure 1(a) in the sense that the former and the latter demultiplex, respectively, the soft information $\{L_d(c(k)){\}}_{k=\mathrm{1}}^{\mathrm{2}NK}$ and QPSK symbols $\{x(k){\}}_{k=\mathrm{1}}^{NK}$ . Once $\{\tilde{a}(k){\}}_{k=\mathrm{1}}^{NK}$ , the hard decisions on the source bits, are available, we follow the steps in the symbol generation process shown in Figure 1(a): $\{\tilde{a}(k){\}}_{k=\mathrm{1}}^{NK}$ are fed into the convolutional encoder, followed by random interleaving, QPSK mapping, and demultiplexing. This way, an error free decoding ensures a perfect recovery of the transmitted QPSK symbols $\{x_n(k){\}}_{n=\mathrm{1},k=\mathrm{1}}^{NK}$ . The recovered payload symbols will be used in the decision-directed channel estimation stage; see Figure 1(b).

3. Double-Selective Channel with Doppler Scaling Effects

In this section, we start with the modeling of the double-selective channel suffering from both ISI and Doppler scaling effects. Then we describe the temporal resampling procedure to mitigate the Doppler scaling effects. After that, we provide a practical approach to estimate the Doppler scaling factor.

3.1. Channel Model

By adopting a single-carrier communication scheme, at the nth transmitter, the continuous baseband signal $x_n(t)$ (generated by passing the discrete payload symbols $\{x_n(k)\}$ to a pulse shaping filter) and its corresponding frequency modulated signal ${\breve{x}}_n(t)$ are related through

$$\begin{array}{cc}& {\breve{x}}_n\left(t\right)=\mathrm{Re}\left\{x_n\left(t\right)e^{j\mathrm{2}\pi f_{c}t}\right\},\hspace{1em}n=\mathrm{1},\dots ,N,\end{array}$$

(3)

where $f_c$ represents the carrier frequency. For simplicity, the pulse shaping filter, frequency modulation, and real component extraction operation are not shown in Figure 1(a).

Due to multipath effects, the actual transmitted signals ${\{{\breve{x}}_n(t)\}}_{n=\mathrm{1}}^N$ can reach the receive hydrophones via different propagation paths with different delays. Herein, the underlying acoustic channel between each transmitter and receiver pair is characterized by R resolved paths. The rth path between the nth transmitter and mth receiver pair ( $r=\mathrm{1},\dots ,R$ , $n=\mathrm{1},\dots ,N$ , and $M=\mathrm{1},\dots ,M$ ) will affect the transmitted signal ${\breve{x}}_n(t)$ in three aspects, namely, amplitude attenuation, Doppler scaling, and delay, which are denoted, respectively, by three real-valued scalars ${\kappa }_{n,m}(r)$ , ${\alpha }_{n,m}(r)$ , and ${\tau }_{n,m}(r)$ . The signal transmitted via the rth path and acquired by the mth receiver can be written as ${\kappa }_{n,m}(r)·{\breve{x}}_n({\alpha }_{n,m}(r)t-{\tau }_{n,m}(r))$ . By taking into account all of the N transducers and R resolved paths, the received signal at the mth hydrophone can be expressed as (for simplicity, the noise term is omitted for the time being)

$$\begin{array}{l}{\breve{z}}_m\left(t\right)=\sum _{n=\mathrm{1}}^N‍\sum _{r=\mathrm{1}}^R‍{\kappa }_{n,m}\left(r\right)·{\breve{x}}_n\left({\alpha }_{n,m}\left(r\right)t-{\tau }_{n,m}\left(r\right)\right),\hfill \\ \phantom{{\kappa }_{n,m}\left(r\right)·{\breve{x}}_n\left({\alpha }_{n,m}\left(r\right)t1\left(r\right)\right)}m=\mathrm{1},\dots ,M.\hfill \end{array}$$

(4)

We assume that the propagation paths for all the transmitter and receiver pairs experience a common Doppler scaling factor and the resolved paths are synchronized among all the transmitter and receiver pairs, that is, ${\alpha }_{n,m}(r)=\alpha$ and ${\tau }_{n,m}(r)=\tau (r)$ (interested readers are referred to [18] for a detailed treatment of synchronization procedure). By using these assumptions, (4) reduces to

$$\begin{array}{cc}& {\breve{z}}_m\left(t\right)=\sum _{n=\mathrm{1}}^N‍\sum _{r=\mathrm{1}}^R‍{\kappa }_{n,m}\left(r\right)·{\breve{x}}_n\left(\alpha t-\tau \left(r\right)\right),\hspace{1em}m=\mathrm{1},\dots ,M.\end{array}$$

(5)

Substituting (3) into (5) yields

$$\begin{array}{cc}& {\breve{z}}_m\left(t\right)=\mathrm{Re}\left\{\sum _{n=\mathrm{1}}^N‍\sum _{r=\mathrm{1}}^R‍{\kappa }_{n,m}\left(r\right)·x_n\left(\alpha t-\tau \left(r\right)\right)e^{j\mathrm{2}\pi f_c(\alpha t-\tau (r\left)\right)}\right\}.\end{array}$$

(6)

3.2. Temporal Resampling

By resampling the received measurements $\{{\breve{z}}_m(t)\}$ using a factor β, the resampled signal ${\breve{y}}_m(t)$ is given by [4, 19, 20]

$$\begin{array}{cc}& {\breve{y}}_m\left(t\right)={\breve{z}}_m\left(\frac{t}{\beta }\right).\end{array}$$

(7)

Then, the baseband received signal $y_m(t)$ , which is related to ${\breve{y}}_m(t)$ via ${\breve{y}}_m(t)=\mathrm{Re}\{y_m(t)e^{j\mathrm{2}\pi f_{c}t}\}$ , can be expressed as

$$\begin{array}{cc}& y_m\left(t\right)=e^{-\mathrm{2}j\pi \left(\right(\beta -\alpha )/\beta )f_{c}t}\sum _{n=\mathrm{1}}^N‍\sum _{r=\mathrm{1}}^R‍h_{n,m,r}·x_n\left(\frac{\alpha }{\beta }t-\tau \left(r\right)\right),\end{array}$$

(8)

where $h_{n,m,r}\triangleq {\kappa }_{n,m}(r)e^{-j\pi f_c\tau (r)}$ represents the rth channel tap between the nth transmitter and the mth receiver pair, for $r=\mathrm{1},\dots ,R$ , $n=\mathrm{1},\dots ,N$ , and $m=\mathrm{1},\dots ,M$ . It can be readily verified that as long as $\alpha /\beta \approx \mathrm{1}$ , we have $x_n((\alpha /\beta )t-\tau (r))\approx x_n(t-\tau (r))$ . Accordingly, (8) can be approximated as

$$\begin{array}{cc}& y_m\left(t\right)\approx e^{-\mathrm{2}j\pi \left(\right(\beta -\alpha )/\beta )f_{c}t}\sum _{n=\mathrm{1}}^N‍\sum _{r=\mathrm{1}}^R‍h_{n,m,r}·x_n\left(t-\tau \left(r\right)\right).\end{array}$$

(9)

One observes from (9) that effective temporal resampling (meaning $\alpha \approx \beta$ ) converts the Doppler scaling effects to Doppler frequency shifts with the frequency given below:

$$\begin{array}{cc}& f=\left(\frac{\beta -\alpha }{\beta }\right)f_c.\end{array}$$

(10)

Therefore, the determination of the resampling factor β plays a crucial role in the effective mitigation of the Doppler scaling effects.

3.3. Resampling Factor Estimation

We take advantage of the preamble and the postamble of a data packet to estimate β [19] (the structure of a data packet will be discussed in Section 6). By cross-correlating the received signal with the known preamble and postamble, the receiver estimates the time duration of a packet ${\widehat{T}}_{rx}$ [4]. By comparing ${\widehat{T}}_{rx}$ with $T_{tx}$ , the duration of the same packet at the transmitter side, the Doppler scaling factor can be estimated as

$$\begin{array}{cc}& \widehat{\beta }=\frac{{\widehat{T}}_{rx}}{T_{tx}}.\end{array}$$

(11)

Although this method is conceptually simple and easy to implement, its accuracy is sensitive to the signal-to-noise-ratio (SNR).

More accurate Doppler scaling factor estimate can be achieved via channel estimation instead of cross correlation. Based on the CIRs estimated from the two measurement segments in response to the preamble and postamble, the change in the time duration ${\widehat{T}}_d$ imposed on the packet can be inferred from the tap shift of the principal arrivals of these two CIRs. Then the Doppler scaling factor estimate can be computed as

$$\begin{array}{cc}& \widehat{\beta }=\frac{T_{tx}+{\widehat{T}}_d}{T_{tx}}.\end{array}$$

(12)

The $\widehat{\beta }$ obtained using (12) is more robust against the noise contamination than the one from (11). We will show later on in Section 6.2.1 via the MACE10 in-water experimental data that the method in (12) works well in practice.

4. Channel Estimation

Since the $\widehat{\beta }$ obtained using (12) can never be perfectly accurate, after temporal resampling, Doppler frequency shifts (see (10)) still exist, although Doppler scaling effects become negligible. We start below with the problem formulation of channel estimation in both training-directed and decision-directed modes [5, 6]. Then, we propose the GoSLIM-V algorithm for jointly estimating the underlying CIRs and Doppler frequency.

In what follows, ${\{y_m(t)\}}_{m=\mathrm{1}}^M$ and ${\{x_n(t)\}}_{n=\mathrm{1}}^N$ in (9) are represented in discrete-time form. Unless otherwise stated, it is assumed that the channel taps for all the $N\times M$ transmitter-receiver pairs experience the same Doppler frequency f.

4.1. Training-Directed Mode

The initial task of the receiver is to acquire knowledge of the underlying channel between all transmitter and receiver pairs using the training sequences. By adopting the cyclic prefix scheme in [7], the training sequence at the nth transmitter ( $n=\mathrm{1},\dots ,N$ ) is given by

$$\begin{array}{l}{\mathbf{x}}_n=[\underset{L_{\text{CP}}\hspace{0.17em}\hspace{0.17em}\text{prefix}\hspace{0.17em}\hspace{0.17em}\text{symbols}}{\underbrace{x_n\left(P-L_{\text{CP}}+\mathrm{1}\right),\dots ,x_n\left(P\right)}},\hfill \\ \phantom{ssssds}\underset{P\hspace{0.17em}\hspace{0.17em}\text{core}\hspace{0.17em}\hspace{0.17em}\text{training}\hspace{0.17em}\hspace{0.17em}\text{symbols}}{\underbrace{x_n\left(\mathrm{1}\right),x_n\left(\mathrm{2}\right),\dots ,x_n\left(P\right)}}],\hfill \end{array}$$

(13)

where $[x_n(\mathrm{1}),\dots ,x_n(P)]$ is the core training sequence and the leading $L_{\text{CP}}$ symbols form the cyclic prefix. In general, we have $P>L_{\text{CP}}\ge R-\mathrm{1}$ . From an amplifier efficiency point of view, it is practically desirable to use unit modulus (unimodular) training sequences, that is, $|x_n(p)|=\mathrm{1}$ for $n=\mathrm{1},\dots ,N$ and $p=\mathrm{1},\dots ,P$ .

For MIMO UAC over acoustic channels subject to both ISI and Doppler frequency shifts, the measurement vectors can be written as [8, 21]

$$\begin{array}{cc}& {\mathbf{y}}_m=\stackrel{-}{\mathrm{\Lambda }}\sum _{n=\mathrm{1}}^N‍{\mathbf{X}}_n{\mathbf{h}}_{n,m}+{\mathbf{e}}_m,\hspace{1em}m=\mathrm{1},\dots ,M,\end{array}$$

(14)

where

$$\begin{array}{cc}& {\mathbf{y}}_m={[y_m\left(\mathrm{1}\right),\dots ,y_m\left(P\right)]}^T,\hspace{1em}m=\mathrm{1},\dots ,M,\end{array}$$

(15)

contains the P synchronized measured symbols (for instance, $\{y_m(\mathrm{1})\}$ maps to $\{x_n(\mathrm{1})\}$ ) at the mth receiver. ${\mathbf{X}}_n\in {ℂ}^{P\times R}$ is given by

$$\begin{array}{cc}& {\mathbf{X}}_n=\left[\begin{array}{cccc}{x}_n\left(\mathrm{1}\right)& x_n\left(P\right)& \cdots & x_n\left(P-R+\mathrm{2}\right)\\ x_n\left(\mathrm{2}\right)& x_n\left(\mathrm{1}\right)& \cdots & x_n\left(P-R+\mathrm{3}\right)\\ ⋮& ⋮& \ddots & ⋮\\ x_n\left(P\right)& x_n\left(P-\mathrm{1}\right)& \cdots & x_n\left(P-R+\mathrm{1}\right)\end{array}\right],\end{array}$$

(16)

where $n=\mathrm{1},\dots ,N$ and ${\mathbf{e}}_m$ represents additive noise at the mth receiver. In addition,

$$\begin{array}{cc}& {\mathbf{h}}_{n,m}={[h_{n,m,\mathrm{1}},\dots ,h_{n,m,R}]}^T,\end{array}$$

(17)

characterizes the channel of length R between the nth transmitter and the mth receiver, where $n=\mathrm{1},\dots ,N$ and $m=\mathrm{1},\dots ,M$ ( $\{h_{n,m,r}\}$ is defined after (8)). Finally, the so-called Doppler shift matrix $\stackrel{-}{\mathrm{\Lambda }}\in {ℂ}^{P\times P}$ in (14) is constructed as

$$\begin{array}{cc}& \stackrel{-}{\mathrm{\Lambda }}=\mathrm{diag}\left(\left[\mathrm{1},e^{-\mathrm{2}j\pi f{T}_s},\dots ,e^{-\mathrm{2}j\pi f{T}_s(P-\mathrm{1})}\right]\right),\end{array}$$

(18)

for $m=\mathrm{1},\dots ,M$ , where f and $T_s$ represent the Doppler frequency and symbol period, respectively.

The ISI and Doppler shift effects can be viewed separately in (14). More specifically, the term ${\sum }_{n=\mathrm{1}}^N‍{\mathbf{X}}_n{\mathbf{h}}_{n,m}$ indicates the net contribution of N ISI channels, while the impact of the Doppler effects on the measurements comes through $\stackrel{-}{\mathrm{\Lambda }}$ only, which corresponds to the assumption that all the $NMR$ CIR taps involved (recall that we have N transmit transducers, M receive hydrophones, and each transmitter-receiver pair corresponds to an R-tap channel) experience the same Doppler frequency f. The purpose of setting the first diagonal element of $\stackrel{-}{\mathrm{\Lambda }}$ to 1 (see (18)) is to eliminate ambiguities. In our example, relative to $y_m(\mathrm{1})$ , a generic measurement, say $y_m(p)$ , experiences a phase shift of $-f{T}_s(p-\mathrm{1})$ .

We express (14) in a more compact form as

$$\begin{array}{cc}& {\mathbf{y}}_m=\stackrel{-}{\mathrm{\Lambda }}\stackrel{-}{\mathbf{X}}{\mathbf{h}}_m+{\mathbf{e}}_m,\hspace{1em}m=\mathrm{1},\dots ,M,\end{array}$$

(19)

where $\stackrel{-}{\mathbf{X}}=[{\mathbf{X}}_{\mathrm{1}},\dots ,{\mathbf{X}}_N]$ and ${\mathbf{h}}_m=[{\mathbf{h}}_{\mathrm{1},m}^T,\dots ,{\mathbf{h}}_{N,m}^T{]}^T$ . By stacking up the measurements, (19) can be rewritten as

$$\begin{array}{cc}& \mathbf{y}=\mathrm{\Lambda }\mathbf{X}\mathbf{h}+\mathbf{e},\end{array}$$

(20)

where $\mathbf{y}=[{\mathbf{y}}_{\mathrm{1}}^T,{\mathbf{y}}_{\mathrm{2}}^T,\dots ,{\mathbf{y}}_M^T{]}^T$ , $\mathbf{h}=[{\mathbf{h}}_{\mathrm{1}}^T,{\mathbf{h}}_{\mathrm{2}}^T,\dots ,{\mathbf{h}}_M^T{]}^T$ , $\mathbf{e}=[{\mathbf{e}}_{\mathrm{1}}^T,{\mathbf{e}}_{\mathrm{2}}^T,\dots ,{\mathbf{e}}_M^T{]}^T$ ,  $\mathrm{\Lambda }={\mathbf{I}}_{M\times M}\otimes \stackrel{-}{\mathrm{\Lambda }}$ , and $\mathbf{X}={\mathbf{I}}_{M\times M}\otimes \stackrel{-}{\mathbf{X}}$ . Then the training-directed channel estimation reduces to estimating h and f from the measurement vector y and known X. The subject of synthesizing unimodular training sequences, coupled with the employment of the cyclic prefix scheme, to facilitate ISI channel estimation is thoroughly treated in [6]. The shifted PeCAN waveforms [22] are used as the training sequences in the MACE10 in-water experimentations.

4.2. Decision-Directed Mode

The decision-directed channel estimation problem is only a slight twist of its training-directed counterpart. For the former, we use the hard decisions of the previously estimated payload symbols, instead of the training symbols, to estimate the channels; see Figure 1(b). Accordingly, (14) can still be used, where

$$\begin{array}{cc}& {\mathbf{y}}_m={\left[⇳y_m\left(t_i\right),\dots ,y_m\left(t_f\right)\right]}^T,\hspace{1em}m=\mathrm{1},\dots ,M,\end{array}$$

(21)

contains the measurements at the mth receiver belonging to the time index interval $[t_i,t_f]$ , and

$$\begin{array}{cc}& {\mathbf{X}}_n=\left[\begin{array}{cccc}{\tilde{x}}_n\left(t_i\right)& {\tilde{x}}_n\left(t_i-\mathrm{1}\right)& \cdots & {\tilde{x}}_n\left(t_i-R+\mathrm{1}\right)\\ {\tilde{x}}_n\left(t_i+\mathrm{1}\right)& {\tilde{x}}_n\left(t_i\right)& \cdots & {\tilde{x}}_n\left(t_i-R+\mathrm{2}\right)\\ ⋮& ⋮& \ddots & ⋮\\ {\tilde{x}}_n\left(t_f\right)& {\tilde{x}}_n\left(t_f-\mathrm{1}\right)& \cdots & {\tilde{x}}_n\left(t_f-R+\mathrm{1}\right)\end{array}\right],\end{array}$$

(22)

for $n=\mathrm{1},\dots ,N$ , where ${\tilde{x}}_n(t_i-R+\mathrm{1})$ and ${\tilde{x}}_n(t_f)$ represent the hard decisions of the first and last previously estimated symbols (some of them could be the known training symbols), respectively, used for updating the channel. (For notational simplicity, ${\mathbf{X}}_n$ is used in both (16) and (22) to represent two similar but different quantities. The use of ${\mathbf{X}}_n$ , however, should be clear from the context.) The tracking length is represented as $L_{TR}=t_f-t_i+\mathrm{1}$ , that is, the number of rows of ${\mathbf{X}}_n$ . To conform with the matrix dimensions, the Doppler shift matrix $\stackrel{-}{\mathrm{\Lambda }}$ now is $L_{TR}$ by $L_{TR}$ , constructed as $\stackrel{-}{\mathrm{\Lambda }}={\mathrm{diag}}_{}([\mathrm{1},e^{-\mathrm{2}j\pi f{T}_s},\dots ,e^{-\mathrm{2}j\pi f{T}_s(L_{TR}-\mathrm{1})}])$ . Similar to the training-directed mode, the channel estimation problem in the decision-directed mode aims to estimate h and f from the measurement vector y and known X formed from the decision-directed $\{{\mathbf{X}}_n{\}}_{n=\mathrm{1}}^N$ in (22).

4.3. Channel Estimation Algorithm: GoSLIM-V

The channel estimation algorithm, in either training- or decision-directed mode, has the generic form given by (20). We remark that the number of elements in ${\mathbf{y}}_m$ , namely, $d_y$ , might vary in different modes. e in (20) is assumed to contain circularly symmetric independent and identically distributed (i.i.d.) complex-valued Gaussian random variables with zero-mean and variance η, denoted as $\mathbf{e}~𝒞𝒩(\mathbf{0},\eta \mathbf{I})$ . The problem is then to estimate f (contained in Λ) and h given y and X. In UAC systems, the channel h is usually sparse; that is, although it contains $NMR$ unknowns, many of them can be approximated as zero [23]. We present the GoSLIM-V algorithm to solve this channel estimation problem. Note that since h contains the CIRs of all $N\times M$ transmitter-receiver pairs, the GoSLIM-V algorithm will estimate them simultaneously.

Consider the following hierarchical Bayesian model:

$$\begin{array}{cc}& \mathbf{y}\hspace{0.17em}|\hspace{0.17em}\mathbf{h},\mathrm{\Lambda },\eta ~𝒞𝒩\left(\mathrm{\Lambda }\mathbf{X}\mathbf{h},\eta \mathbf{I}\right),\end{array}$$

(23)

$$\begin{array}{cc}& \mathbf{h}\hspace{0.17em}|\hspace{0.17em}\mathbf{p}~𝒞𝒩\left(\mathbf{0},\mathbf{P}\right),\end{array}$$

(24)

where (23) follows directly from the assumption $\mathbf{e}~𝒞𝒩(\mathbf{0},\eta \mathbf{I})$ . Let $p_{n,m,r}$ be the variance of $h_{n,m,r}$ for $n=\mathrm{1},\dots ,N$ ,  $m=\mathrm{1},\dots ,M$ , and  $r=\mathrm{1},\dots ,R$ , and define ${\mathbf{p}}_{n,m}=[p_{n,m,\mathrm{1}},p_{n,m,\mathrm{2}},\dots ,p_{n,m,R}{]}^T$ ,   ${\mathbf{p}}_m=[{\mathbf{p}}_{\mathrm{1},m}^T,{\mathbf{p}}_{\mathrm{2},m}^T,\dots ,{\mathbf{p}}_{N,m}^T{]}^T$ , and $\mathbf{p}=[{\mathbf{p}}_{\mathrm{1}}^T,{\mathbf{p}}_{\mathrm{2}}^T,\dots ,{\mathbf{p}}_M^T{]}^T$ . Then the covariance matrix P in (24) is constructed as $\mathbf{P}=\mathrm{diag}(\mathbf{p})$ .

Furthermore, by considering a flat prior on f, η, and $\{p_{n,m,r}{\}}_{n=\mathrm{1},\hspace{0.17em}m=\mathrm{1},\hspace{0.17em}r=\mathrm{1}}^{NMR}$ , the channel vector h, Doppler frequency f, the covariance matrix P (or more precisely, its diagonal elements p), and the noise power η can be estimated based on the MAP criterion

$$\begin{array}{cc}& \underset{\mathbf{h},\mathbf{p},\eta ,f}{\max }\hspace{0.17em}p\left(\mathbf{h},\mathbf{p},\eta ,f\hspace{0.17em}|\hspace{0.17em}\mathbf{y}\right)=\underset{\mathbf{h},\mathbf{p},\eta ,f}{\max }\hspace{0.17em}p\left(\mathbf{y}\hspace{0.17em}|\hspace{0.17em}\mathbf{h},\eta ,f\right)p\left(\mathbf{h}\hspace{0.17em}|\hspace{0.17em}\mathbf{p}\right).\end{array}$$

(25)

By combining (23), (24), and (25), and by taking the negative logarithm of the cost function, the optimization problem formulated in (25) becomes

$$\begin{array}{l}\underset{\mathbf{h},\mathbf{p},\eta ,f}{\min }(d_{y}M\log \eta +\frac{{\parallel \mathbf{y}-\mathrm{\Lambda }\mathbf{X}\mathbf{h}\parallel }^{\mathrm{2}}}{\eta }\hfill \\ \phantom{bbb000}+\sum _{n=\mathrm{1}}^N‍\sum _{m=\mathrm{1}}^M‍\sum _{r=\mathrm{1}}^R‍\log p_{n,m,r}+\sum _{n=\mathrm{1}}^N‍\sum _{m=\mathrm{1}}^M‍\sum _{r=\mathrm{1}}^R‍\frac{{\left|h_{n,m,r}\right|}^{\mathrm{2}}}{p_{n,m,r}}),\hfill \end{array}$$

(26)

which can be solved using an alternating approach; at each iteration, one of the parameters h, p, η, and f is updated while keeping the other three fixed. In this way, the single difficult joint optimization problem is divided into 4 simpler separate subproblems. GoSLIM-V keeps iterating until a predefined iteration number is reached. Under mild conditions, the cyclic optimization scheme guarantees that the GoSLIM-V algorithm converges, at least to a local minimum of (26) [24].

The 5 steps of the GoSLIM-V algorithm at the tth iteration are outlined below.

• (1)

Given ${\mathbf{h}}^{(t-\mathrm{1})}$ , the optimal ${\mathbf{P}}^{(t)}$ that minimizes the cost function in (26) is given by

$$\begin{array}{cc}& p_{n,m,r}^{\left(t\right)}={\left|h_{n,m,r}^{(t-\mathrm{1})}\right|}^{\mathrm{2}},\end{array}$$

(27)

for $n=\mathrm{1},\dots ,N$ ,    $m=\mathrm{1},\dots ,M$ , and $r=\mathrm{1},\dots ,R$ . For better numerical stability, we set $p_{n,m,r}^{(t)}$ (or equivalently $h_{n,m,r}^{(t)}$ ) to zero if $p_{n,m,r}^{(t)}<\mathrm{1}{\mathrm{0}}^{-\mathrm{15}}$ .

(2)

Once ${\mathbf{P}}^{(t)}$ is available, the CIR is updated as

$$\begin{array}{cc}& {\mathbf{h}}^{\left(t\right)}={\left[{\mathbf{X}}^H\mathbf{X}+{\eta }^{(t-\mathrm{1})}{\left({\mathbf{P}}^{\left(t\right)}\right)}^{-\mathrm{1}}\right]}^{-\mathrm{1}}{\left({\mathrm{\Lambda }}^{(t-\mathrm{1})}\mathbf{X}\right)}^H\mathbf{y}.\end{array}$$

(28)

While inverting ${\mathbf{P}}^{(t)}$ , its zero diagonal entries are removed, and the associated columns in X are discarded.

(3)

Next, using the most recently obtained $\{{\mathbf{h}}_m^{(t)}\}$ in (28), we estimate the Doppler frequency f. For ease of exposition, we denote $z_m^{(t)}(i)=y_m^{*}(i){\breve{x}}_m^{(t)}(i)$ , where $y_m(i)$ and ${\breve{x}}_m^{(t)}(i)$ represent, respectively, the ith element of the measurement vector ${\mathbf{y}}_m$ and ${\breve{\mathbf{x}}}_m^{(t)}$ with ${\breve{\mathbf{x}}}_m^{(t)}=\stackrel{-}{\mathbf{X}}{\mathbf{h}}_m^{(t)}$ , $i=\mathrm{1},\dots ,d_y$ . It is easy to verify that

$$\begin{array}{l}{\parallel \mathbf{y}-\mathrm{\Lambda }\mathbf{X}{\mathbf{h}}^{\left(t\right)}\parallel }^{\mathrm{2}}=\text{const}\hfill \\ \phantom{=}-\mathrm{2}\mathrm{Re}\left[\sum _{i=\mathrm{1}}^{d_y}‍\left(\sum _{m=\mathrm{1}}^M‍z_m^{\left(t\right)}\left(i\right)\right)e^{-\mathrm{2}j\pi f{T}_s(i-\mathrm{1})}\right].\hfill \end{array}$$

(29)

Since the constant term in (29) is not a function of f, minimizing the cost function in (26) is equivalent to solving

$$\begin{array}{cc}& f^{\left(t\right)}=\arg \underset{f}{\max }\mathrm{Re}\left[\sum _{i=\mathrm{1}}^{d_y}‍\left(\sum _{m=\mathrm{1}}^M‍z_m^{\left(t\right)}\left(i\right)\right)e^{-\mathrm{2}j\pi f{T}_s(i-\mathrm{1})}\right].\end{array}$$

(30)

Since the summation term within the parenthesis above is nothing but the discrete-time Fourier transform (DTFT) of the sequence $\{{\sum }_{m=\mathrm{1}}^M‍z_m^{(t)}(i){\}}_{i=\mathrm{1}}^{d_y}$ evaluated at frequency f, $f^{(t)}$ is obtained as the location of the dominant peak of the real part of the DTFT.

(4)

Using the ${\mathbf{h}}^{(t)}$ and ${\mathrm{\Lambda }}^{(t)}$ most recently obtained via (28) and (30), respectively, we finally estimate the noise power as

$$\begin{array}{cc}& {\eta }^{\left(t\right)}=\frac{\mathrm{1}}{d_{y}M}{\parallel \mathbf{y}-{\mathrm{\Lambda }}^{\left(t\right)}\mathbf{X}{\mathbf{h}}^{\left(t\right)}\parallel }^{\mathrm{2}}.\end{array}$$

(31)

(5)

Set $t=t+\mathrm{1}$ . Go back to Step 1 if t is less than a predefined iteration number, or terminate otherwise.

In the training-directed mode, the channel characteristics in general are not available a priori. In our examples, ${\mathbf{h}}^{(\mathrm{0})}$ is initialized using the standard matched filter, $f^{(\mathrm{0})}$ is initialized as 0, and the noise power ${\eta }^{(\mathrm{0})}$ is initialized with a small positive number, for instance, $\mathrm{1}{\mathrm{0}}^{-\mathrm{10}}$ . Our empirical experience suggests that the GoSLIM-V algorithm does not provide significant performance improvements after 15 iterations or less.

5. Symbol Detection

In this section, we proceed to study the detection of the payload symbols given the estimates of CIRs and Doppler frequency f obtained by GoSLIM-V. The detection task is achieved via two steps: phase compensation followed by turbo equalization. As shown in Figure 1(b), turbo equalization consists of an equalizer and a decoder, both configurated as soft-input soft-output. The decoder is conventionally implemented by the Max-Log-MAP algorithm [17], and our focus herein is on the soft-input soft-output equalizer. We first formulate the symbol detection problem and then describe the phase compensation procedure. After that, we elaborate the LMMSE-based turbo equalization design and discuss its low complexity approximation.

5.1. Problem Formulation

Treating the transmitted symbols as the unknowns and the CIRs and Doppler frequency as known, the measurement vector in (14) can be expressed as [8, 21]

$$\begin{array}{cc}& {\mathbf{y}}_m\left(k\right)=\widehat{\mathrm{\Lambda }}\left(k\right)\sum _{n=\mathrm{1}}^N‍{\widehat{\mathbf{H}}}_{n,m}{\mathbf{x}}_n\left(k\right)+{\mathbf{e}}_m,\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}m=\mathrm{1},\dots ,M,\end{array}$$

(32)

where the estimated CIR matrix ${\widehat{\mathbf{H}}}_{n,m}\in {ℂ}^{\hspace{0.17em}\hspace{0.17em}R\times (\mathrm{2}R-\mathrm{1})}$ is given by

$$\begin{array}{cc}& {\widehat{\mathbf{H}}}_{n,m}=\left[\begin{array}{ccccc}{\widehat{h}}_{n,m,R}& \cdots & {\widehat{h}}_{n,m,\mathrm{1}}& & \mathbf{0}\\ ⋮& \ddots & & \ddots & ⋮\\ \mathbf{0}& & {\widehat{h}}_{n,m,R}& \cdots & {\widehat{h}}_{n,m,\mathrm{1}}\end{array}\right],\end{array}$$

(33)

for $n=\mathrm{1},\dots ,N$ and $m=\mathrm{1},\dots ,M$ . The entry ${\widehat{h}}_{n,m,r}$ here represents the estimate of $h_{n,m,r}$ in (17) given by GoSLIM-V at the conclusion of the iteration. Also,

$$\begin{array}{l}{\mathbf{x}}_n\left(k\right)={[x_n\left(k-R+\mathrm{1}\right),\dots ,x_n\left(k\right),\dots ,x_n\left(k+R-\mathrm{1}\right)]}^T,\hfill \\ \phantom{xxxxxxxxxxxxxxxxxxxxxxxxxxxxxx}n=\mathrm{1},\dots ,N,\hfill \end{array}$$

(34)

$$\begin{array}{l}{\mathbf{y}}_m\left(k\right)={[y_m\left(k\right),\dots ,y_m\left(k+R-\mathrm{1}\right)]}^T,\hspace{1em}m=\mathrm{1},\dots ,M.\hfill \end{array}$$

(35)

The variable k represents the time index corresponding to the payload symbols of current interest. Although ${\mathbf{y}}_m$ represents different portions of the received signal in (15), (21), and (35), its use should be clear from the context. Per the discussions following (18), once $\widehat{f}$ is available, the estimated Doppler shift matrix $\widehat{\mathrm{\Lambda }}(k)$ in (32) can be constructed as

$$\begin{array}{cc}& \widehat{\mathrm{\Lambda }}\left(k\right)=\mathrm{diag}\left(\left[e^{-\mathrm{2}j\pi \widehat{f}{T}_s(k-\mathrm{1})},\dots ,e^{-\mathrm{2}j\pi \widehat{f}{T}_s(k+R-\mathrm{2})}\right]\right).\end{array}$$

(36)

When detecting symbols, we use the estimates $\{{\widehat{\mathbf{h}}}_{n,m}\}$ and $\widehat{f}$ obtained from the previous channel update and we treat $\{{\widehat{\mathbf{H}}}_{n,m}\}$ and $\widehat{\mathrm{\Lambda }}(k)$ in (32) as known.

5.2. Phase Compensation

Stacking up all the measurements, (32) can be written as

$$\begin{array}{cc}& \left[\begin{array}{c}{\mathbf{y}}_{\mathrm{1}}\left(k\right)\\ ⋮\\ {\mathbf{y}}_M\left(k\right)\end{array}\right]=\stackrel{°}{\mathrm{\Lambda }}\left(k\right)\sum _{n=\mathrm{1}}^N‍\left[\begin{array}{c}{\widehat{\mathbf{H}}}_{n,\mathrm{1}}\\ ⋮\\ {\widehat{\mathbf{H}}}_{n,M}\end{array}\right]{\mathbf{x}}_n\left(k\right)+\left[\begin{array}{c}{\mathbf{e}}_{\mathrm{1}}\\ ⋮\\ {\mathbf{e}}_M\end{array}\right],\end{array}$$

(37)

or equivalently as

$$\begin{array}{cc}& \mathbf{y}\left(k\right)=\stackrel{°}{\mathrm{\Lambda }}\left(k\right)\sum _{n=\mathrm{1}}^N‍{\widehat{\mathbf{H}}}_n{\mathbf{x}}_n\left(k\right)+\mathbf{e}=\stackrel{°}{\mathrm{\Lambda }}\left(k\right)\widehat{\mathbf{H}}\mathbf{x}\left(k\right)+\mathbf{e},\end{array}$$

(38)

where $\mathbf{y}(k)$ and $\mathbf{e}\in {ℂ}^{\hspace{0.17em}\hspace{0.17em}MR\times \mathrm{1}}$ , $\stackrel{°}{\mathrm{\Lambda }}(k)={\mathbf{I}}_{M\times M}\otimes \widehat{\mathrm{\Lambda }}(k)$ , $\{{\widehat{\mathbf{H}}}_n{\}}_{n=\mathrm{1}}^N\in {ℂ}^{\hspace{0.17em}\hspace{0.17em}MR\times (\mathrm{2}R-\mathrm{1})}$ , $\widehat{\mathbf{H}}=[{\widehat{\mathbf{H}}}_{\mathrm{1}},\dots ,{\widehat{\mathbf{H}}}_N]$ , and $\mathbf{x}(k)=[{\mathbf{x}}_{\mathrm{1}}^T(k),\dots ,{\mathbf{x}}_N^T(k){]}^T$ . The phase compensation task is simply achieved by multiplying $[\stackrel{°}{\mathrm{\Lambda }}(k){]}^H$ to both sides of (38), yielding

$$\begin{array}{cc}& \stackrel{°}{\mathbf{y}}\left(k\right)=\widehat{\mathbf{H}}\mathbf{x}\left(k\right)+\stackrel{°}{\mathbf{e}},\end{array}$$

(39)

where $\stackrel{°}{\mathbf{y}}(k)=[\stackrel{°}{\mathrm{\Lambda }}(k){]}^H\mathbf{y}(k)$ and $\stackrel{°}{\mathbf{e}}=[\stackrel{°}{\mathrm{\Lambda }}(k){]}^H\mathbf{e}$ . Given $\mathbf{e}~𝒞𝒩(\mathbf{0},\eta \mathbf{I})$ , $\stackrel{°}{\mathbf{e}}$ still has the distribution of $𝒞𝒩(\mathbf{0},\eta \mathbf{I})$ since $[\stackrel{°}{\mathrm{\Lambda }}(k){]}^H$ is unitary.

Phase compensation, along with the aforementioned temporal resampling process, effectively converts the original double-selective channel to an ISI channel. Given the phase-compensated measurement vector $\stackrel{°}{\mathbf{y}}(k)$ , the estimated CIR matrix $\widehat{\mathbf{H}}$ , and $\{L_d(c_n(k)){\}}_{n=\mathrm{1},\hspace{0.17em}k=\mathrm{1}}^{N\mathrm{2}K}$ , we consider using an LMMSE-based soft-input soft-output equalizer to compute the $a\hspace{0.17em}posteriori\hspace{0.17em}extrinsic$ information of the transmitted signal.

5.3. LMMSE-Based Soft-Input Soft-Output Equalizer

As shown in Figure 2, an LMMSE-based soft-input soft-output equalizer can be functionally divided into four modules. The $a\hspace{0.17em}priori$ LLR preprocessor calculates the mean and the variance of each QPSK payload symbol $x_n(k)$ , denoted as ${\stackrel{-}{x}}_n(k)$ and $v_n(k)$ , respectively, from the $a\hspace{0.17em}posteriori$    $extrinsic$ information $L_d(c_n(\mathrm{2}k-\mathrm{1}))$ and $L_d(c_n(\mathrm{2}k))$ generated by the decoder for $n=\mathrm{1},\dots ,N$ and $k=\mathrm{1},\dots ,K$ . Next, the transmitted symbol $x_n(k)$ is estimated via LMMSE filtering given $\stackrel{°}{\mathbf{y}}(k)$ and $\widehat{\mathbf{H}}$ in (39), along with $\{{\stackrel{-}{x}}_n(k)\}$ and $\{v_n(k)\}$ . Specifically, as demonstrated in Figure 2, the LMMSE filter is applied to the residual signal generated by subtracting out the so-called soft interferences from the phase-compensated measurements. The soft interferences characterize the contributions of all the payload symbols except $x_n(k)$ , the one of the current interest, in terms of soft information. Based on the symbol estimates ${\widehat{x}}_n(k)$ , the $a\hspace{0.17em}posteriori$ LLR generator provides the $extrinsic$ LLR outputs $L_e(c_n(\mathrm{2}k-\mathrm{1}))$ and $L_e(c_n(\mathrm{2}k))$ ( $n=\mathrm{1},\dots ,N$ , $k=\mathrm{1},\dots ,K$ ), which will be fed into the soft-input soft-output decoder as $a\hspace{0.17em}priori$ LLR; see Figure 1(b). In the following, these modules will be elaborated further.

Figure 2

The structure of the LMMSE-based soft-input soft-output equalizer.

5.3.1. A Priori LLR Preprocessor

In this task, we calculate ${\stackrel{-}{x}}_n(k)$ and $v_n(k)$ from $L_d(c_n(\mathrm{2}k-\mathrm{1}))$ and $L_d(c_n(\mathrm{2}k))$ . According to the definitions, the mean and variance of $x_n(k)$ are given by [14]

$$\begin{array}{l}{\stackrel{-}{x}}_n\left(k\right)=\stackrel{\mathrm{4}}{\underset{i=\mathrm{1}}{\sum ‍}}{\alpha }_i·P\left(x_n\left(k\right)={\alpha }_i\right),\hfill \end{array}$$

(40)

$$\begin{array}{l}{v}_n\left(k\right)=\stackrel{\mathrm{4}}{\underset{i=\mathrm{1}}{\sum ‍}}{\left|{\alpha }_i\right|}^{\mathrm{2}}·P\left(x_n\left(k\right)={\alpha }_i\right)-{\left|{\stackrel{-}{x}}_n\left(k\right)\right|}^{\mathrm{2}}\hfill \\ =\mathrm{1}-{\left|{\stackrel{-}{x}}_n\left(k\right)\right|}^{\mathrm{2}},\hfill \end{array}$$

(41)

where $\{{\alpha }_i{\}}_{i=\mathrm{1}}^{\mathrm{4}}$ denote the four QPSK constellation points of 𝒮 and $|{\alpha }_i|=\mathrm{1}$ for $i=\mathrm{1},\dots ,\mathrm{4}$ ; see the definition of 𝒮 after (1). One can see from (41) that $v_n(k)$ depends on ${\stackrel{-}{x}}_n(k)$ , and the evaluation of ${\stackrel{-}{x}}_n(k)$ in (40) requires $P(x_n(k)={\alpha }_i)$ for $i=\mathrm{1},\dots ,\mathrm{4}$ . Since the interleaved bits $\{c(k)\}$ can be reasonably assumed to be independent of each other due to the employment of the random interleaver (see Figure 1(a)), $P(x_n(k)={\alpha }_i)$ is determined as the product of the probabilities of the two interleaved bits that map to ${\alpha }_i$ . For instance, $c_n(\mathrm{2}k-\mathrm{1})=\mathrm{0}$ and $c_n(\mathrm{2}k)=\mathrm{1}$ map to $x_n(k)=(-\mathrm{1}+j)/\sqrt{\mathrm{2}}$ according to (1), and, therefore, $P(x_n(k)=(-\mathrm{1}+j)/\sqrt{\mathrm{2}})=P(c_n(\mathrm{2}k-\mathrm{1})=\mathrm{0})·P(c_n(\mathrm{2}k)=\mathrm{1})$ , where for a generic interleaved bit $c_n(q)$ we have

$$\begin{array}{cc}& P\left(c_n\left(q\right)=\mathrm{0}\right)=\frac{e^{L_d\left(c_n\left(q\right)\right)}}{\mathrm{1}+e^{L_d\left(c_n\left(q\right)\right)}},\\ P\left(c_n\left(q\right)=\mathrm{1}\right)\hspace{0.17em}\hspace{0.17em}=\frac{\mathrm{1}}{\mathrm{1}+e^{L_d\left(c_n\right(q\left)\right)}},\end{array}$$

(42)

for $n=\mathrm{1},\dots ,N$ and $q=\mathrm{1},\dots ,\mathrm{2}K$ . Equation (42) follows from (2) and $P(c_n(q)=\mathrm{0})+P(c_n(q)=\mathrm{1})=\mathrm{1}$ .

Plugging (42) into (40) gives

$$\begin{array}{l}{\stackrel{-}{x}}_n\left(k\right)=\frac{\mathrm{1}}{\sqrt{\mathrm{2}}}\{j\tanh \left[\frac{L_d\left(c_n\left(\mathrm{2}k-\mathrm{1}\right)\right)}{\mathrm{2}}\right]\hfill \\ \phantom{ssssssg5sss1sg}+\tanh \left[\frac{L_d\left(c_n\left(\mathrm{2}k\right)\right)}{\mathrm{2}}\right]\},\hspace{1em}\hfill \\ n=\mathrm{1},\dots ,N,\hspace{0.17em}\hspace{0.17em}k=\mathrm{1},\dots ,K,\hfill \end{array}$$

(43)

which, combined with (41), yields $v_n(k)$ .

5.3.2. LMMSE Filtering

Depending on whether the a priori LLR information is incorporated or not, two types of LMMSE filters are studied in the following.

In the Absence of A Priori Knowledge

The equalizer is performed in the absence of a priori knowledge at the very first iteration before using the decoder. This scenario amounts to setting $L_d(c_n(k))=\mathrm{0}$ for $n=\mathrm{1},\dots ,N$ and $k=\mathrm{1},\dots ,\mathrm{2}K$ , which implies that ${\stackrel{-}{x}}_n(k)=\mathrm{0}$ and $v_n(k)=\mathrm{1}$ according to (43) and (41) for $n=\mathrm{1},\dots ,N$ and $k=\mathrm{1},\dots ,K$ . In this case, the LMMSE filter coefficient vector, denoted as ${\mathbf{f}}_n$ , is given by [5, 25]

$$\begin{array}{cc}& {\mathbf{f}}_n={\left(\widehat{\mathbf{H}}{\widehat{\mathbf{H}}}^H+\widehat{\eta }\mathbf{I}\right)}^{-\mathrm{1}}{\mathbf{s}}_n.\end{array}$$

(44)

Here, $\widehat{\eta }$ represents the noise power estimate given by GoSLIM-V at the conclusion of the iteration, ${\mathbf{s}}_n=[{\widehat{\mathbf{h}}}_{n,\mathrm{1}}^T,\dots ,{\widehat{\mathbf{h}}}_{n,M}^T{]}^T$ denotes the steering vector corresponding to $x_n(k)$ in (38) for $n=\mathrm{1},\dots ,N$ , and ${\widehat{\mathbf{h}}}_{n,m}$ is the estimate of ${\mathbf{h}}_{n,m}$ defined in (17). An estimate of $x_n(k)$ is obtained by applying ${\mathbf{f}}_n$ to the phase-compensated measurement vector obtained in (39) as

$$\begin{array}{cc}& {\widehat{x}}_n\left(k\right)={\mathbf{f}}_n^H\stackrel{°}{\mathbf{y}}\left(k\right),\hspace{1em}n=\mathrm{1},\dots ,N.\end{array}$$

(45)

In the Presence of A Priori Knowledge

In this case, the LMMSE estimate of $x_n(k)$ is given by [14]

$$\begin{array}{cc}& {\widehat{x}}_n\left(k\right)={\stackrel{-}{x}}_n\left(k\right)+v_n\left(k\right){\mathbf{f}}_n^{\prime }{\left(k\right)}^H\left[\stackrel{°}{\mathbf{y}}\left(k\right)-\widehat{\mathbf{H}}E\left(\mathbf{x}\left(k\right)\right)\right],\end{array}$$

(46)

where

$$\begin{array}{cc}& {\mathbf{f}}_n^{\prime }\left(k\right)={\left[\widehat{\mathbf{H}}\mathbf{V}\left(k\right){\widehat{\mathbf{H}}}^H+\widehat{\eta }\mathbf{I}\right]}^{-\mathrm{1}}{\mathbf{s}}_n\end{array}$$

(47)

represents the LMMSE filter coefficient vector. In (46), each component of $E(\mathbf{x}(k))$ is the expected value of the corresponding component of $\mathbf{x}(k)$ calculated in (40). In (47), the covariance matrix $\mathbf{V}(k)=\mathrm{diag}([{\mathbf{v}}_{\mathrm{1}}(k{)}^T,\dots ,{\mathbf{v}}_N(k{)}^T])$ , and

$$\begin{array}{l}{\mathbf{v}}_n\left(k\right)={[v_n\left(k-R+\mathrm{1}\right),\dots ,v_n\left(k\right),\dots ,v_n(k+R-\mathrm{1}\left)\right]}^T,\hfill \\ \phantom{[v_n\left(k-R+\mathrm{1}\right),\dots \mathrm{5},v_n\left(k\right),v_n(k+R-\mathrm{1})]}n=\mathrm{1},\dots ,N.\hfill \end{array}$$

(48)

Each component of ${\mathbf{v}}_n(k)$ is obtained according to (41).

Equation (46) suggests that the estimation of $x_n(k)$ depends on its own $extrinsic$ LLR information $L_d(c_n(\mathrm{2}k-\mathrm{1}))$ and $L_d(c_n(\mathrm{2}k))$ , whose impact on ${\widehat{x}}_n(k)$ comes through ${\stackrel{-}{x}}_n(k)$ and $v_n(k)$ . From the belief propagation theory point of view, the generation of $extrinsic$ information of a payload symbol needs to avoid such dependency [26]. To achieve this goal, we modify (46) as

$$\begin{array}{l}{\widehat{x}}_n\left(k\right)={\mathbf{s}}_n^H{\left[\widehat{\mathbf{H}}\mathbf{V}\left(k\right){\widehat{\mathbf{H}}}^H+\widehat{\eta }\mathbf{I}+\left(\mathrm{1}-v_n\left(k\right)\right){\mathbf{s}}_n{\mathbf{s}}_n^H\right]}^{-\mathrm{1}}\hfill \\ \phantom{=}\times \left[\stackrel{°}{\mathbf{y}}\left(k\right)-\widehat{\mathbf{H}}E\left(\mathbf{x}\left(k\right)\right)+{\stackrel{-}{x}}_n\left(k\right){\mathbf{s}}_n\right].\hfill \end{array}$$

(49)

Compared to (46), the presence of the two additional terms in (49), namely, $(\mathrm{1}-v_n(k)){\mathbf{s}}_n{\mathbf{s}}_n^H$ and ${\stackrel{-}{x}}_n(k){\mathbf{s}}_n$ , resembles a scenario of ${\stackrel{-}{x}}_n(k)=\mathrm{0}$ and $v_n(k)=\mathrm{1}$ (or equivalently, $L_d(c_n(\mathrm{2}k-\mathrm{1}))=L_d(c_n(\mathrm{2}k))=\mathrm{0}$ ) in (46), as if $x_n(k)$ is estimated without incorporating its own LLR information.

Define

$$\begin{array}{cc}& {\mathbf{f}}_n^{\mathrm{\prime \prime }}\left(k\right)={\left[\widehat{\mathbf{H}}\mathbf{V}\left(k\right){\widehat{\mathbf{H}}}^H+\widehat{\eta }\mathbf{I}+\left(\mathrm{1}-v_n\left(k\right)\right){\mathbf{s}}_n{\mathbf{s}}_n^H\right]}^{-\mathrm{1}}{\mathbf{s}}_n,\end{array}$$

(50)

$$\begin{array}{cc}& {\dot{\mathbf{y}}}_n\left(k\right)=\stackrel{°}{\mathbf{y}}\left(k\right)-\left[\widehat{\mathbf{H}}E\left(\mathbf{x}\left(k\right)\right)-{\stackrel{-}{x}}_n\left(k\right){\mathbf{s}}_n\right].\end{array}$$

(51)

Then, (49) can be rewritten as

$$\begin{array}{cc}& {\widehat{x}}_n\left(k\right)={\mathbf{f}}_n^{\mathrm{\prime \prime }}{\left(k\right)}^H{\dot{\mathbf{y}}}_n\left(k\right).\end{array}$$

(52)

In (51), the terms within the square brackets correspond to the output of the soft interference generator in Figure 2. To get ${\widehat{x}}_n(k)$ , the LMMSE filter coefficient vector in (50) is applied to the residual measurement vector ${\dot{\mathbf{y}}}_n(k)$ . Note that (52) includes (45) as a special case when no a priori knowledge is available, that is, $L_d(c_n(k))=\mathrm{0}$ for $n=\mathrm{1},\dots ,N$ and $k=\mathrm{1},\dots ,\mathrm{2}K$ .

5.3.3. A  Posteriori LLR Generator

This task calculates the $extrinsic$ LLR $L_e(c_n(\mathrm{2}k-\mathrm{1}))$ and $L_e(c_n(\mathrm{2}k))$ from the symbol estimates ${\widehat{x}}_n(k)$ obtained in (45) or (52) for $n=\mathrm{1},\dots ,N$ and $k=\mathrm{1},\dots ,K$ .

We assume that, given $x_n(k)={\alpha }_i$ , ${\widehat{x}}_n(k)$ is a circularly symmetric i.i.d. complex-valued Gaussian random process, that is, $P({\widehat{x}}_n(k)\hspace{0.17em}|\hspace{0.17em}{x}_n(k)={\alpha }_i)~𝒞𝒩({\mu }_i,{\sigma }^{\mathrm{2}})$ , where the mean ${\mu }_i$ and variance ${\sigma }^{\mathrm{2}}$ are calculated, respectively, as ${\mu }_i={\alpha }_i{\mathbf{f}}_n^{\mathrm{\prime \prime }}(k{)}^H{\mathbf{s}}_n$ and ${\sigma }^{\mathrm{2}}={\mathbf{f}}_n^{\mathrm{\prime \prime }}(k{)}^H{\mathbf{s}}_n-{\mathbf{f}}_n^{\mathrm{\prime \prime }}(k{)}^H{\mathbf{s}}_n{\mathbf{s}}_n{\mathbf{f}}_n^{\mathrm{\prime \prime }}(k)$ [14]. Under this assumption, the output LLR of the two consecutive bits mapping to $x_n(k)$ is calculated as [14]

$$\begin{array}{cc}& L_e\left(c_n\left(\mathrm{2}k-\mathrm{1}\right)\right)=\frac{\sqrt{\mathrm{8}}\mathrm{Im}\left({\mathbf{f}}_n^{\mathrm{\prime \prime }}{\left(k\right)}^H{\dot{\mathbf{y}}}_n\left(k\right)\right)}{\mathrm{1}-{\mathbf{s}}_n^H{\mathbf{f}}_n^{\mathrm{\prime \prime }}\left(k\right)},\\ L_e\left(c_n\left(\mathrm{2}k\right)\right)=\frac{\sqrt{\mathrm{8}}\mathrm{Re}\left({\mathbf{f}}_n^{\mathrm{\prime \prime }}{\left(k\right)}^H{\dot{\mathbf{y}}}_n\left(k\right)\right)}{\mathrm{1}-{\mathbf{s}}_n^H{\mathbf{f}}_n^{\mathrm{\prime \prime }}\left(k\right)},\end{array}$$

(53)

for $n=\mathrm{1},\dots ,N$ and $k=\mathrm{1},\dots ,K$ .

Let ${\mathbf{R}}^{\prime }(k)=\widehat{\mathbf{H}}\mathbf{V}(k){\widehat{\mathbf{H}}}^H+\widehat{\eta }\mathbf{I}$ and ${\mathbf{R}}_n^{\mathrm{\prime \prime }}(k)=\widehat{\mathbf{H}}\mathbf{V}(k){\widehat{\mathbf{H}}}^H+\widehat{\eta }\mathbf{I}+(\mathrm{1}-v_n(k)){\mathbf{s}}_n{\mathbf{s}}_n^H$ . Then ${\mathbf{f}}_n^{\prime }(k)$ in (47) and ${\mathbf{f}}_n^{\mathrm{\prime \prime }}(k)$ in (50) can be rewritten as ${\mathbf{f}}_n^{\prime }(k)={\mathbf{R}}^{\prime }(k{)}^{-\mathrm{1}}{\mathbf{s}}_n$ and ${\mathbf{f}}_n^{\mathrm{\prime \prime }}(k)={\mathbf{R}}_n^{\mathrm{\prime \prime }}(k{)}^{-\mathrm{1}}{\mathbf{s}}_n$ , respectively. One observes that the derivation of $\{{\mathbf{f}}_n^{\mathrm{\prime \prime }}(k)\}$ requires the inversion of ${\mathbf{R}}_n^{\mathrm{\prime \prime }}(k)$ for each transmitter at each time index, whereas the computation of $\{{\mathbf{f}}_n^{\prime }(k)\}$ needs to invert ${\mathbf{R}}^{\prime }(k)$ at each time index. Consequently, by following (47) and (50) directly, the computational complexity of calculating $\{{\mathbf{f}}_n^{\mathrm{\prime \prime }}(k)\}$ is approximately N times more expensive than obtaining $\{{\mathbf{f}}_n^{\prime }(k)\}$ .

Since ${\mathbf{R}}_n^{\mathrm{\prime \prime }}(k)={\mathbf{R}}^{\prime }(k)+(\mathrm{1}-v_n(k)){\mathbf{s}}_n{\mathbf{s}}_n^H$ , the use of the matrix inversion lemma gives

$$\begin{array}{cc}& {\mathbf{R}}_n^{\mathrm{\prime \prime }}{\left(k\right)}^{-\mathrm{1}}={\mathbf{R}}^{\prime }{\left(k\right)}^{-\mathrm{1}}-\frac{\left(\mathrm{1}-v_n\left(k\right)\right){\mathbf{R}}^{\prime }{\left(k\right)}^{-\mathrm{1}}{\mathbf{s}}_n{\mathbf{s}}_n^H{\mathbf{R}}^{\prime }{\left(k\right)}^{-\mathrm{1}}}{\mathrm{1}+\left(\mathrm{1}-v_n\left(k\right)\right){\mathbf{s}}_n^H{\mathbf{R}}^{\prime }{\left(k\right)}^{-\mathrm{1}}{\mathbf{s}}_n}.\end{array}$$

(54)

Right multiplying ${\mathbf{s}}_n$ on both sides of (54) yields

$$\begin{array}{cc}& {\mathbf{f}}_n^{\mathrm{\prime \prime }}\left(k\right)=\frac{{\mathbf{f}}_n^{\prime }\left(k\right)}{\mathrm{1}+\left(\mathrm{1}-v_n\left(k\right)\right){\mathbf{s}}_n^H{\mathbf{f}}_n^{\prime }\left(k\right)},\end{array}$$

(55)

which, combined with (53), follows

$$\begin{array}{cc}& L_e\left(c_n\left(\mathrm{2}k-\mathrm{1}\right)\right)=\frac{\sqrt{\mathrm{8}}\mathrm{Im}\left({\mathbf{f}}_n^{\prime }\left(k\right){\dot{\mathbf{y}}}_n\left(k\right)\right)}{\mathrm{1}-v_n\left(k\right){\mathbf{s}}_n^H{\mathbf{f}}_n^{\prime }\left(k\right)},\\ L_e\left(c_n\left(\mathrm{2}k\right)\right)=\frac{\sqrt{\mathrm{8}}\mathrm{Re}\left({\mathbf{f}}_n^{\prime }\left(k\right){\dot{\mathbf{y}}}_n\left(k\right)\right)}{\mathrm{1}-v_n\left(k\right){\mathbf{s}}_n^H{\mathbf{f}}_n^{\prime }\left(k\right)}.\end{array}$$

(56)

Complexitywise, the LLR calculation formula in (56) is preferable over (53) since, as we just remarked, it is more efficient to calculate $\{{\mathbf{f}}_n^{\prime }(k)\}$ than $\{{\mathbf{f}}_n^{\mathrm{\prime \prime }}(k)\}$ . Due to this reason, LLR is calculated according to (56) in our numerical and experimental examples provided later on.

5.4. Low-Complexity Approximate LMMSE Filtering

Although the calculation of a posteriori LLR according to (56) is more efficient than (53), it still constitutes the major computational bottleneck in turbo equalization mainly because $\{{\mathbf{f}}_n^{\prime }(k)\}$ needs to be calculated at each time index. To further reduce the computational complexity, we consider a low-complexity approximate LMMSE filter whose coefficient vector is given by [14]

$$\begin{array}{cc}& {\mathbf{f}}_n^{\prime }={\left(\widehat{\mathbf{H}}\stackrel{-}{\mathbf{V}}{\widehat{\mathbf{H}}}^H+\widehat{\eta }\mathbf{I}\right)}^{-\mathrm{1}}{\mathbf{s}}_n,\end{array}$$

(57)

where $\stackrel{-}{\mathbf{V}}=(\mathrm{1}/K){\sum }_{k=\mathrm{1}}^K‍\mathbf{V}(k)$ . Since $\{{\mathbf{f}}_n^{\prime }\}$ is constant for each transmitter over one payload block (hence the time index k is dropped in (57)), the overall complexity of calculating $\{{\mathbf{f}}_n^{\prime }\}$ is approximately K times faster than deriving $\{{\mathbf{f}}_n^{\prime }(k)\}$ according to (47). Substituting ${\mathbf{f}}_n^{\prime }(k)$ in (56) with ${\mathbf{f}}_n^{\prime }$ yields

$$\begin{array}{cc}& L_e\left(c_n\left(\mathrm{2}k-\mathrm{1}\right)\right)=\frac{\sqrt{\mathrm{8}}\mathrm{Im}\left({\mathbf{f}}_n^{\prime }{\dot{\mathbf{y}}}_n\left(k\right)\right)}{\mathrm{1}-v_n\left(k\right){\mathbf{s}}_n^H{\mathbf{f}}_n^{\prime }},\\ L_e\left(c_n\left(\mathrm{2}k\right)\right)=\frac{\sqrt{\mathrm{8}}\mathrm{Re}\left({\mathbf{f}}_n^{\prime }{\dot{\mathbf{y}}}_n\left(k\right)\right)}{\mathrm{1}-v_n\left(k\right){\mathbf{s}}_n^H{\mathbf{f}}_n^{\prime }}.\end{array}$$

(58)

We hereafter refer to the turbo equalization scheme that calculates the a posteriori extrinsic information according to (56) and (58) as exact-LMMSE turbo and approximate-LMMSE turbo, respectively.

Note that matrix inversion is an indispensable stage in calculating the LMMSE filter coefficients in (44), (47), and (57). To expedite the calculation, we can make use of the conjugate gradient (CG) method and fast Fourier transform (FFT) operations, as elaborated in [27]. Although [27] focuses on the efficient calculation of the LMMSE filter coefficients in the form of (44), the extension to a more general scenario in (47) or (57) is straightforward. In the present paper, both exact-LMMSE turbo and approximate-LMMSE turbo are implemented using the FFT-based CG method.

6. Numerical and Experimental Results

6.1. Numerical Results

Consider transmitting four payload blocks simultaneously over time-invariant ISI channels using a MIMO UAC system equipped with $N=\mathrm{4}$ transmitters and $M=\mathrm{12}$ receivers. Block length is fixed at $K=\mathrm{250}$ . The four payload blocks across the $N=\mathrm{4}$ transmitters are constructed from a randomly generated binary source sequence of length $NK=\mathrm{1000}$ according to the procedure detailed in Section 2. We simulate $N\times M=\mathrm{48}$ frequency-selective channels involved in the MIMO UAC system. To resemble practical UAC scenarios, these simulated CIRs are estimated from MACE10 in-water experimental data and each CIR has $R=\mathrm{50}$ taps. CIRs have been normalized to 1, that is, ${\parallel {\mathbf{h}}_{n,m}\parallel }^{\mathrm{2}}=\mathrm{1}$ for $n=\mathrm{1},\dots ,N$ and $m=\mathrm{1},\dots ,M$ . The received data samples are then constructed according to (14). Since Doppler effects are not considered in this example, $\stackrel{-}{\mathrm{\Lambda }}=\mathbf{I}$ . The noise vector $\{{\mathbf{e}}_m\}$ is assumed to contain circularly symmetric i.i.d. complex-valued Gaussian random variables with zero-mean and variance ${\sigma }^{\mathrm{2}}$ . The simulation of ISI channels, combined with the assumption that each receiver has perfect knowledge on the channel characteristics $\{{\mathbf{h}}_{n,m}\}$ , suggests that we can bypass the temporal resampling, channel estimation, and phase compensation modules in Figure 1(b) and apply exact-LMMSE turbo, approximate-LMMSE turbo, and RELAX-BLAST directly to the received measurements. Figures 3(a) and 3(b) show the average coded bit error rate (BER) given by exact-LMMSE turbo and approximate-LMMSE turbo, respectively, along with the RELAX-BLAST performance at different SNRs, where SNR is defined as $\mathrm{1}/{\sigma }^{\mathrm{2}}$ . Each point is averaged over 500 Monte-Carlo trials. The binary source sequence and the noise pattern vary from one trial to another. The curve labeled as “No Iteration” is obtained by employing the equalizer and the decoder only once, that is, the feedback loop is yet to be formed. In addition, the average coded BER given by RELAX-BLAST is obtained after three iterations. We can see from Figure 3 that both types of turbo equalization schemes effectively reduce the coded BER as the iteration proceeds and significantly outperform RELAX-BLAST, and exact-LMMSE turbo provides only slightly better detection performance than approximate-LMMSE turbo. Complexity wise, the average time required to finish one trial is 18.64 s, 0.49 s, and 0.19 s on an ordinary workstation (Intel Xeon E5506 processor 2.13 GHz, 12 GB RAM, Windows 7 64-bit, and MATLAB R2010b) for exact-LMMSE turbo, approximate-LMMSE turbo, and RELAX-BLAST, respectively. Consequently, approximate-LMMSE turbo is preferred over its exact-LMMSE turbo counterpart since the former provides almost the same detection performance as the latter but with a computational complexity on the same order as RELAX-BLAST.

Figure 3

(a) Coded BER performance by using exact-LMMSE turbo along with RELAX-BLAST performance. (b) Coded BER performance by using approximate-LMMSE turbo along with RELAX-BLAST performance. Each point is averaged over 500 Monte Carlo trials. In this simulation, $N=\mathrm{4}$ , $M=\mathrm{12}$ , $R=\mathrm{50}$ , and $K=\mathrm{250}$ .

6.2. MACE10 In-Water Experimental Results

6.2.1. Experiment Specifics

The MACE10 in-water experiment was conducted by the Woods Hole Oceanographic Institution (WHOI) off the coast of Martha's Vineyard, MA, USA, in June 2010. A source array consisting of 4 transducers was vertically deployed at a depth of 80 m and towed by a vessel. At the receiver side, a 12-element hydrophone array was mounted on a buoy. The vessel moved from the minimum range of 500 m away from the receiving array outbound to the maximum range of 4000 m and then inbound back to the minimum range. The carrier frequency, sampling frequency, and symbol rate employed in the MACE10 experiment were 13 kHz, 39.0625 kHz, and 3.90625 kHz, respectively. By transmitting $N=\mathrm{4}$ sequences simultaneously and incorporating the measurements acquired from all of the $M=\mathrm{12}$ receiver elements for analysis, we established a $\mathrm{4}\times \mathrm{12}$ MIMO UAC system.

The structure of a transmitted data package is shown in Figure 4. Each package consists of 4 packets. The first packet conveys a grayscale Gator mascot and the subsequent 3 packets are combined from a colored mascot. The RGB components of the colored image were transmitted in the 2nd, 3rd, and 4th packets, respectively. Each pixel of the Gator grayscale image is represented by 5 bits, corresponding to 32 different intensities (e.g., pure white and pure dark pixels are represented by 11111 and 00000, resp.). The 64-pixel by 100-pixel grayscale mascot image, as a consequence, is represented by a total of 32 k source bits. Accordingly, a colored mascot image is represented by 96 k bits. The contrast of the grayscale image, as well as the hue of the colored image, has been carefully adjusted so that the image carries approximately equal numbers of 1s and 0s.

Figure 4

The structure of the package used in the MACE10 experiment.

As shown in Figure 4, each packet is constructed as follows: time-marking sequences are placed at the beginning of each packet to facilitate the temporal resampling procedure; two guard intervals, each containing 500 silent symbols, are placed, respectively, before and after the segments containing the payload symbols and training sequences. The payload symbols contain the information of the Gator mascot image. We herein elaborate how to generate the 1st packet from the grayscale Gator mascot image (the packet generation for each of the RGB components of the colored image follows the same procedure). Specifically, the 32 k source bits are first interleaved so that the bits fed into the convolutional encoder module have an equal chance of being 0 or 1; see Figure 5. The so-obtained 32 k interleaved source bits are then divided into 32 groups, each containing 1 k bits. The bits in the $i\text{th}$ group ( $i=\mathrm{1},\dots ,\mathrm{32}$ ) will be used to construct the ith payload symbol block across the 4 transmitted sequences, and the construction procedure follows Figure 1(a). Note that in Figure 5, the depth of the interleavers ${\Pi }^{\prime }$ and Π is 32 k and 2 k, respectively. Figure 5 illustrates a scenario with $i=\mathrm{1}$ . The shifted PeCAN training sequences with length $P=\mathrm{512}$ , in conjunction with $L_{\text{CP}}=\mathrm{99}$ cyclic prefix symbols, form the training section, which is located between the 16th and 17th payload blocks. This MIMO UAC design leads to a net coded data rate of 11.7 kbps. The data package was transmitted periodically and recorded by the receiver array. A total of 120 epochs were available and they are referred to as “E001”−“E120,” respectively.

Figure 5

The structure of the transmitted symbols for the $\mathrm{4}\times \mathrm{12}$ MIMO BLAST scheme used in MACE10.

To estimate the Doppler scaling factor, we treat the time-marking sequences at the beginning of a packet as its preamble and those at the beginning of the subsequent packet as the postamble. Take the 2nd packet of epoch “E002,” for example. For the channel between the 1st transmitter and the 1st receiver, the superimposed modulus of the CIRs obtained by GoSLIM-V from the preamble and postamble is shown in Figure 6. The indexes of the principal arrivals for the preamble and postamble are 12 and 21, respectively. Hence, the time duration change imposed on the packet is ${\widehat{T}}_d=(\mathrm{21}-\mathrm{12})T_s$ , where $T_s$ is the symbol period defined after (18). Then the Doppler scaling factor $\widehat{\beta }$ can be estimated according to (12).

Figure 6

The superimposed modulus of the CIRs obtained from the preamble and postamble, respectively.

To assess the performance of the resampling process, the CIR and Doppler frequency evolutions obtained by GoSLIM-V before we resample the 2nd packet of epoch “E002” are shown in Figures 7(a) and 7(c), respectively. In comparison, Figures 7(b) and 7(d) demonstrate the corresponding CIR and Doppler frequency evolutions obtained after resampling the packet, respectively. We can see from Figure 7 that the temporal resampling procedure successfully reduces the Doppler scaling effects to Doppler frequency shifts. The relative speed between the transmitter and the receiver arrays can be estimated as $\widehat{v}=(\widehat{\beta }-\mathrm{1})·c$ , using a common underwater sound speed of $c=\mathrm{1500}$  m/s. It is interesting to look at Figure 8 where the vessel speed estimated during the resampling stage is plotted on top of the GPS reference information provided by WHOI (the GPS device was equipped on the moving vessel). The good agreement between these two curves verifies the effectiveness of the resampling procedure we employ. The analysis presented hereafter is based on the resampled measurements.

Figure 7

(a) CIR evolution of epoch “E002” before resampling. (b) CIR evolution of epoch “E002” after resampling. (c) Doppler frequency evolution of epoch “E002” before resampling. (d) Doppler frequency evolution of epoch “E002” after resampling.

Figure 8

The relative speed between the transmitter and receiver array given by GPS and estimated during the temporal resampling stage (courtesy of Milica Stojanovic's group).

6.2.2. Performance Evaluation

By fixing the tracking length at $L_{TR}=\mathrm{450}$ , the channel tracking starts with training-directed channel estimation using GoSLIM-V. Then we perform phase compensation on the received measurements, followed by the detection of the first $K=\mathrm{250}$ payload symbols contained in the 17th payload block for each transmitted sequence using RELAX-BLAST, exact-LMMSE turbo, and approximate-LMMSE turbo; see Figure 5. Next, the channels are updated in the decision-directed mode using $L_{TR}=\mathrm{450}$ symbols (containing the previously detected payload symbols, as well as a portion of the training sequence as well). With the updated CIRs and Doppler frequency, after phase compensation, the subsequent $K=\mathrm{250}$ payload symbols contained in the 18th block are estimated using the same symbol detection scheme. This process continues until all of the 16 payload blocks to the right-hand side of the training sequences are detected. This same tracking scheme can be applied in a reverse manner to the detection of the 16 payload blocks ahead of the training sequences.

We deem a packet to be successfully detected if the resulting coded BER is less than 0.1. After analyzing a total of 480 packets available, Table 1 summarizes the successfully detected packet percentage, the zero BER packet percentage, the coded BER averaged over the successful packets, and the time ratio of the time consumed to process a packet on the workstation specified in Section 6.1 to $T_{tx}=\mathrm{2.741}\hspace{0.17em}\text{s}$ ( $T_{tx}$ is defined in (12)) obtained using exact-LMMSE turbo, approximate-LMMSE turbo, and the RELAX-BLAST scheme, respectively. The results are obtained by applying 3 iterations for all of the three types of detection schemes considered. One observes from Table 1 that (1) BER-wise, both exact-LMMSE turbo and approximate-LMMSE turbo outperform RELAX-BLAST significantly, (2) compared to exact-LMMSE turbo, approximate-LMMSE turbo greatly reduces the computational time at the cost of slight BER performance degradation, and (3) compared to RELAX-BLAST, approximate-LMMSE turbo improves the BER performance by two orders of magnitude without significantly increasing the computational complexities. These observations are in line with those made from the numerical examples in Section 6.1. Moreover, we analyze epoch “E054” that leads to perfect recovery of both the grayscale and colored mascots (see Figures 9(b) and 9(d)) using either exact-LMMSE Turbo or approximate-LMMSE turbo. In comparison, the grayscale and colored mascots recovered from epoch “E054” using RELAX-BLAST are shown in Figures 9(a) and 9(c), respectively, with the corresponding coded BERs being $\mathrm{1.8}\times \mathrm{1}{\mathrm{0}}^{-\mathrm{1}}$ and $\mathrm{1.5}\times \mathrm{1}{\mathrm{0}}^{-\mathrm{1}}$ . We note that the turbo equalization schemes are highly effective.

Table 1

A summary of the performance of the three detection schemes (3 iterations applied).

	Successful packet percentage (%)	Zero BER packet percentage (%)	Average coded BER	Time ratio
Exact-LMMSE turbo	100	76.7	$\mathrm{9.2}\times \mathrm{1}{\mathrm{0}}^{-\mathrm{5}}$	488.1
Approximate-LMMSE turbo	100	74.4	$\mathrm{2.1}\times \mathrm{1}{\mathrm{0}}^{-\mathrm{4}}$	17.2
RELAX-BLAST	82.5	4.8	$\mathrm{1.6}\times \mathrm{1}{\mathrm{0}}^{-\mathrm{2}}$	16.9

Figure 9

(a) Grayscale mascot recovered from epoch “E054” using RELAX-BLAST. (b) Grayscale mascot recovered from epoch “E054” using turbo equalization. (c) Colored mascot recovered from epoch “E054” using RELAX-BLAST. (d) Colored mascot recovered from epoch “E054” using turbo equalization.

To further illustrate the detection performance of turbo equalization, Table 2 shows the coded BER averaged over all of the 480 packets at different iteration numbers obtained by exact-LMMSE turbo and approximate-LMMSE turbo. We can see from Table 2 that the coded BER improves with iteration. Empirical experience indicates that the detection performance for both types of turbo equalization converges after three iterations. Next, we choose one payload block and denote $\{L_d(a(k)){\}}_{k=\mathrm{1}}^{\mathrm{1000}}$ as the LLR soft information of the corresponding 1 k source bits $\{a(k){\}}_{k=\mathrm{1}}^{\mathrm{1000}}$ generated by the Max-Log-MAP decoder. Figures 10(a)–10(d) and 10(e)–10(h) show $\{L_d(a(k)){\}}_{k=\mathrm{1}}^{\mathrm{1000}}$ obtained by exact-LMMSE turbo and approximate-LMMSE turbo, respectively, at different iteration numbers. $\{\tilde{a}(k)\}$ in Figure 1(b) are the hard decisions determined from $\{L_d(a(k))\}$ . Specifically, if $L_d(a(k))>\mathrm{0}$ then $\tilde{a}(k)=\mathrm{0}$ ; otherwise, $\tilde{a}(k)=\mathrm{1}$ (see (2)). In Figure 10, the circles indicate bit errors. We can see from Figure 10 that the LLRs of source bits are moving away from zero as the iteration proceeds (the first iteration has the most significant impact), which suggests that with the help of cycling soft information, the decoder is more and more confident about the corresponding source bits being 0 or 1.

Table 2

The average coded BER obtained by exact-LMMSE turbo and approximate-LMMSE turbo, respectively.

	No iteration	1 iteration	2 iterations	3 iterations
Exact-LMMSE turbo	$\mathrm{2.7}\times \mathrm{1}{\mathrm{0}}^{-\mathrm{1}}$	$\mathrm{8.1}\times \mathrm{1}{\mathrm{0}}^{-\mathrm{4}}$	$\mathrm{1.3}\times \mathrm{1}{\mathrm{0}}^{-\mathrm{4}}$	$\mathrm{9.2}\times \mathrm{1}{\mathrm{0}}^{-\mathrm{5}}$
Approximate-LMMSE turbo	$\mathrm{2.7}\times \mathrm{1}{\mathrm{0}}^{-\mathrm{1}}$	$\mathrm{2.2}\times \mathrm{1}{\mathrm{0}}^{-\mathrm{3}}$	$\mathrm{3.2}\times \mathrm{1}{\mathrm{0}}^{-\mathrm{4}}$	$\mathrm{2.1}\times \mathrm{1}{\mathrm{0}}^{-\mathrm{4}}$

Figure 10

The LLR soft information about the source bits at the output of the decoder. ((a)–(d)) are obtained by exact-LMMSE turbo from no iteration to 3 iterations, respectively. ((e)–(h)) are obtained by approximate-LMMSE turbo from no iteration to 3 iterations, respectively.

7. Conclusions

For double-selective channels encountered in mobile MIMO UAC, we have demonstrated via the MACE10 in-water experimental data analysis that it is reasonable to assume a common Doppler scaling factor imposed on the propagation paths among all the transmitter and receiver pairs when the Doppler effects are mainly induced by the relative motions between the transmitter and receiver arrays. Temporal resampling has been used to effectively convert the Doppler scaling effects to Doppler frequency shifts. A data-adaptive sparse channel estimation algorithm, referred to as the GoSLIM-V algorithm, is used to estimate the underlying CIRs and Doppler frequency in a joint manner. For symbol detection, we have investigated the turbo equalization schemes implemented by the LMMSE-based soft-input soft-output equalizer as well as its low complexity approximation. The latter provides only slightly degraded detection performance but at a significantly lower computational complexity compared to the former and is thus preferred. The effectiveness of the considered approaches has been verified using both numerical and the MACE10 in-water experimental results.