Channel Estimation and Data Detection in DS-Spread Channels: A Unified Framework, Novel Algorithms, and Waveform Comparison

TL;DR

This paper introduces a unified framework for channel estimation and data detection in DS-spread channels, develops novel algorithms for underwater acoustic communications, and compares multiple waveforms under realistic conditions.

Contribution

It presents a comprehensive framework, new variational Bayesian algorithms for channel estimation, and a detailed comparison of waveforms in DS-spread channels.

Findings

ODSS outperforms other waveforms with low-complexity equalizer

VSSD algorithm equalizes performance among waveforms in coded scenarios

Waveform differences diminish with advanced detection algorithms

Abstract

We present a unified receiver processing framework for communication over delay-scale (DS)-spread channels that arise in underwater acoustic (UWA) communications that addresses both channel estimation (CE) and data detection for different modulation waveforms, namely OFDM, OTFS, OCDM, and ODSS, through a common input--output relation. Using this framework, we conduct a fair and comprehensive comparative study of these waveforms under DS-spread UWA channels and similar receiver complexities. We also develop a novel iterative variational Bayesian (VB) off-grid CE algorithm to estimate the delay and scale parameters of the channel paths, via two approaches: a first-order approximation scheme (FVB) and a second-order approximation scheme (SVB). We propose a low-complexity variational soft symbol detection (VSSD) algorithm that outputs soft symbols and log-likelihood ratios for the data bits, and a data-aided iterative CE and data detection (ICED) scheme that utilizes detected data symbols as \emph{virtual} pilots to further improve the CE and data detection accuracy. Our numerical results reveal the efficacy of the proposed algorithms for CE and data detection. In terms of relative performance of different waveforms, in uncoded communications, (a) with a low-complexity subcarrier-by-subcarrier equalizer, ODSS offers the best performance, followed by OCDM and OTFS, while OFDM performs the worst, and (b) with the VSSD algorithm, OTFS, OCDM, and ODSS perform similarly, and they outperform OFDM. With coded communications, interestingly, all waveforms offer nearly the same BER when the VSSD receiver is employed. Hence, we conclude that when the receiver complexity is constrained, waveform choice matters, especially under harsh channel conditions, whereas with more sophisticated receiver algorithms, these differences disappear.

Equations

$${\mathbf{X}}_{\mathrm{TF}}={\mathbf{F}}_M{\mathbf{X}}_{\mathrm{DD}}{\mathbf{F}}_N^H,$$

$$s_{\mathrm{T}}(t)=\sum _{m=0}^{M-1}\sum _{n=0}^{N-1}{X}_{\mathrm{TF}}[m,n]g_{\mathrm{tx}}(t-nT)e^{j2\pi f_m(t-nT)},$$

$${\mathbf{s}}_{\mathrm{T}}=\mathrm{vec}({\tilde{\mathbf{G}}}_{\mathrm{T}}{\mathbf{F}}_M^H{\mathbf{X}}_{\mathrm{TF}})=({\mathbf{F}}_N^H\otimes {\tilde{\mathbf{G}}}_{\mathrm{T}})\mathbf{x}={\mathbf{G}}_{\mathrm{T}}\mathbf{x},$$

$${\mathbf{s}}_{\mathrm{F}}=({\mathbf{I}}_N^H\otimes {\tilde{\mathbf{G}}}_{\mathrm{T}}{\mathbf{F}}_M^H)\mathbf{x}={\mathbf{G}}_{\mathrm{F}}\mathbf{x},$$

$$s_{\mathrm{C}}(t)=\sum _{m=0}^{M-1}\sum _{n=0}^{N-1}{X}_{\mathrm{C}}[m,n]{\psi }_m(t-nT)e^{j2\pi f_{c}t},$$

$${\mathbf{X}}_{\mathrm{SD}}={\mathbf{T}}_D\mathbf{x},$$

$$\begin{gathered} s_{\mathrm{D}}(t)=\sum _{m=0}^{M-1}\sum _{n=0}^{N(m)-1}{X}_{\mathrm{SD}}[n,m]q^{\frac{m}{2}}{g}_{\mathrm{tx}}\left(q^m\left(t-\frac{n}{q^{m}W}\right)\right) \\ \times e^{j2\pi f_0^D{q}^m\left(t-\frac{n}{q^{m}W}\right)}, \end{gathered}$$

$$r(t)=\sum _{p=1}^P{h}_p\sqrt{{\alpha }_p}s({\alpha }_p(t-{\tau }_p))+w(t),$$

$$\begin{gathered} r(t)=\sum _{p=1}^P{h}_p\sqrt{{\alpha }_p}{\mathcal{F}}^{-1}\left(\mathcal{F}\left(s({\alpha }_p(t-{\tau }_p))\right)\right)+w(t) \\ =\sum _{p=1}^P{h}_p\sqrt{{\alpha }_p}{\mathcal{F}}^{-1}\left(\frac{1}{{\alpha }_p}{e}^{-j2\pi f{\tau }_p}\int s(t)e^{-j2\pi \frac{f}{{\alpha }_p}t}dt\right)+w(t), \end{gathered}$$

$$\mathbf{r}={\mathbf{H}}^t\mathbf{s}+{\mathbf{w}}^t,$$

$$\mathbf{y}={\mathbf{G}}^H\mathbf{r}={\mathbf{G}}^H{\mathbf{H}}^t\mathbf{Gx}+{\mathbf{G}}^H{\mathbf{w}}^t=\mathbf{Hx}+\mathbf{w},$$

$$\mathbf{H}\triangleq {\mathbf{G}}^H{\mathbf{H}}^t\mathbf{G}.$$

$${\mathbf{y}}_{\mathrm{p}}={\mathbf{A}}_{\mathrm{p},\alpha }(\mathbf{\tau },\mathbf{\alpha })\mathbf{h}+{\mathbf{w}}_{\mathrm{p}},$$

$${\mathbf{a}}_{\mathrm{p},\alpha }({\tau }_p,{\alpha }_p)=\sqrt{{\alpha }_p}{\mathbf{G}}_{\mathrm{p}}^H{\mathbf{F}}_{M_{\mathrm{p}}}^H{\mathbf{\Gamma }}_{\mathrm{p},p}{\mathbf{F}}_{M_{\mathrm{p}}}^{\cdot \frac{1}{{\alpha }_p}}{\mathbf{G}}_{\mathrm{p}}{\mathbf{x}}_{\mathrm{p}},$$

$${\mathbf{y}}_{\mathrm{p}}={\mathbf{A}}_{\mathrm{p},\omega }(\bar{\mathbf{\tau }},\bar{\mathbf{\omega }})\bar{\mathbf{h}}+{\mathbf{w}}_{\mathrm{p}},$$

$${\mathbf{A}}_{\mathrm{p},\omega }(\tilde{\mathbf{\tau }},\tilde{\mathbf{\omega }})={\tilde{\mathbf{A}}}_{\mathrm{p},\omega }^{(0)}+{\mathbf{B}}_{\mathrm{p}}\text{ diag}[{\mathbf{\beta }}_{\tilde{\mathbf{\tau }}}]+{\mathbf{C}}_{\mathrm{p}}\text{ diag}[{\mathbf{\beta }}_{\tilde{\mathbf{\omega }}}],$$

$${\tilde{\mathbf{\tau }}}^{(j+1)}$$

$${\tilde{\mathbf{\omega }}}^{(j+1)}$$

$${\tilde{\tau }}_l^{(j+1)}={\tilde{\tau }}_l^{(j)}-\frac{f_{{\tilde{\tau }}_l}^{\prime }({\widetilde{{\tau }_l}}^{(j)})}{f_{{\tilde{\tau }}_l}^{\prime \prime }({\widetilde{{\tau }_l}}^{(j)})},\text{ and}$$

$${\tilde{\omega }}_l^{(j+1)}={\tilde{\omega }}_l^{(j)}-\frac{f_{{\tilde{\omega }}_l}^{\prime }({\widetilde{{\omega }_l}}^{(j)})}{f_{{\tilde{\omega }}_l}^{\prime \prime }({\widetilde{{\omega }_l}}^{(j)})},$$

$$p(\bar{\mathbf{h}}|\mathbf{\delta })=\prod _{l=0}^{N_{\tau }{M}_{\alpha }-1}p({\bar{h}}_l|{\delta }_l)=\prod _{l=0}^{N_{\tau }{M}_{\alpha }-1}\mathcal{C}\mathcal{N}({\bar{h}}_l;0,{\delta }_l^{-1}),$$

$$p(\mathbf{\delta })\triangleq p(\mathbf{\delta };{ϵ}_1,{ϵ}_2)=\prod _{l=0}^{N_{\tau }{M}_{\alpha }-1}\Gamma ({\delta }_l;{ϵ}_1,{ϵ}_2).$$

$$\ln q^{\ast }(z_k)\propto {\mathbb{E}}_{u\ne k}[\ln p(\mathbf{z},{\mathbf{y}}_{\mathrm{p}})],$$

$$p(\mathbf{z},{\mathbf{y}}_{\mathrm{p}})=p(\bar{\mathbf{h}},\mathbf{\delta },\gamma ,{\mathbf{y}}_{\mathrm{p}})=p({\mathbf{y}}_{\mathrm{p}}|\mathbf{h};\gamma )p(\mathbf{h}|\mathbf{\delta })p(\mathbf{\delta })p(\gamma ),$$

$$\ln q^{(j+1)}(\bar{\mathbf{h}})$$

$$\ln q^{(j+1)}(\mathbf{\delta })$$

$$\ln q^{(j+1)}(\gamma )$$

$$\ln q^{(j+1)}(\bar{\mathbf{h}})\propto$$

$$\propto$$

$${\mathbf{\mu }}_{\bar{\mathbf{h}}}^{(j+1)}={\hat{\gamma }}^{(j)}{\mathbf{\Sigma }}_{\bar{\mathbf{h}}}^{(j+1)}{\bar{\mathbf{A}}}_{\mathrm{p},\omega }^{(j{)}^H}{\mathbf{y}}_{\mathrm{p}},\text{ and}$$

$${\mathbf{\Sigma }}_{\bar{\mathbf{h}}}^{(j+1)}=({\hat{\gamma }}^{(j)}{\bar{\mathbf{A}}}_{\mathrm{p},\omega }^{(j{)}^H}{\bar{\mathbf{A}}}_{\mathrm{p},\omega }^{(j)}+{\hat{\mathbf{\Delta }}}^{(j)}{)}^{-1},$$

$$\ln q^{(j+1)}(\mathbf{\delta })\propto$$

$$\propto$$

$${\hat{\delta }}_l^{(j+1)}=⟨{\delta }_l{⟩}_{q^{(j+1)}({\delta }_l)}=\frac{{ϵ}_1+1}{({ϵ}_2+⟨|{\bar{h}}_l{|}^2{⟩}_{q^{(j+1)}(\bar{\mathbf{h}})})},$$

Full text

We also develop a novel iterative variational Bayesian (VB) off-grid CE algorithm to estimate the delay and scale parameters of the channel paths, via two approaches: a first-order approximation scheme (FVB) and a second-order approximation scheme (SVB). We propose a low-complexity variational soft symbol detection (VSSD) algorithm that outputs soft symbols and log-likelihood ratios for the data bits, and a data-aided iterative CE and data detection (ICED) scheme that utilizes detected data symbols as virtual pilots to further improve the CE and data detection accuracy.

Our numerical results reveal the efficacy of the proposed algorithms for CE and data detection. In terms of relative performance of different waveforms, in uncoded communications, (a) with a low-complexity subcarrier-by-subcarrier equalizer, ODSS offers the best performance, followed by OCDM and OTFS, while OFDM performs the worst, and (b) with the VSSD algorithm, OTFS, OCDM, and ODSS perform similarly, and they outperform OFDM. With coded communications, interestingly, all waveforms offer nearly the same BER when the VSSD receiver is employed. Hence, we conclude that when the receiver complexity is constrained, waveform choice matters, especially under harsh channel conditions, whereas with more sophisticated receiver algorithms, these differences disappear.

I Introduction

In the recent literature, a variety of waveforms with associated transceiver architectures have been proposed as candidates for 6G cellular, WiFi, and underwater acoustic communications, with each tailored to specific channel characteristics such as delay spread, delay-Doppler (DD) spread, and delay-scale (DS) spread channels [3, 4, 5, 6, 7, 8]. These studies typically consider a waveform-dependent transceiver structure, different channel models, receive-processing techniques, etc, leaving the following core question unanswered: do these performance differences arise from an inherent difference in the way the waveforms interact with the underlying channel, or do they mainly arise from the differences in the receiver processing (channel estimation and data detection) techniques used while evaluating the performance? Addressing this question requires a consistent receiver design that can be applied across all waveforms, and this motivates a second goal of this work: to develop a unified receiver processing framework that is agnostic to the underlying waveform and encompasses channel estimation, equalization, and soft-symbol detection. Developing such a framework is necessary for comparing different waveforms on a level playing field and thereby answering the core question mentioned above.

In this paper, we present a comprehensive comparative study of the performance of several key waveforms, considering both estimated channel state information at the receiver (ECSIR) and perfect channel state information at the receiver (PCSIR), in the context of delay-scale spread channels that arise in underwater acoustic (UWA) communications. We focus on UWA channels because these channels present the harshest conditions for communications: multipath with large delay spread leading to frequency-selectivity and a spread in the time-scaling of the multipath components due to severe Doppler effects, leading to time-selectivity.

Most radio frequency (RF) wireless communication channels satisfy the narrowband criteria,111The narrowband criteria are: 1. The signaling bandwidth is much smaller than its center frequency; 2. The ratio of the relative speed of the transmitter and receiver to the signal velocity in the medium is much smaller than the time-bandwidth product of the signal. and the relative motion between source, receiver, or scatterers results in a constant frequency shift, also known as Doppler shift, of the signals. One of the prominent ways to model these narrowband channels is the DD-spread representation [3], where each tap of the multipath channel is characterized by a distinct complex gain, propagation delay, and Doppler shift. However, the linear time-varying channels present in UWA communications do not satisfy the narrowband assumptions. Then, the effect of Doppler is to time-scale the transmitted signal [4]. As a result, these wideband channels are better characterized by a DS-spread representation [9], where each tap of the multipath channel is characterized by a distinct complex gain, propagation delay, and time-scale parameter.

Time-scale effects in UWA channels pose significant challenges for data detection, as conventional receivers, such as subcarrier-by-subcarrier processing in orthogonal frequency division multiplexing (OFDM) systems, tend to perform poorly. Orthogonal time frequency space (OTFS) [6] and orthogonal chirp division multiplexing (OCDM) [7] have been introduced in the literature to address the limitations of OFDM in DD-spread channels. Orthogonal delay-scale space (ODSS) [8] is a recently proposed waveform specifically designed for DS-spread channels. Given these advancements, it is pertinent to examine the relative performance of these waveforms under different channel estimation (CE) and data detection methods.

In [10], a detailed study of the challenges associated with OFDM in UWA communications is presented. The BER of OCDM compared to OFDM in shallow water acoustic channels is evaluated in [11, 12], utilizing a channel decoder that processes soft symbols derived from a minimum mean square error (MMSE) equalizer. A comparison of OFDM, OTFS and OCDM in UWA uncoded communications under both PCSIR and ECSIR using MMSE equalizers is presented in [13, 14]. The main conclusion was that, in uncoded communications and with MMSE receivers, OTFS and OCDM outperform OFDM, but the performance gap reduces in the presence of channel estimation errors. Further, [15] proposes a low-complexity MMSE turbo equalization technique for UWA OTFS systems. In [16], a learned denoising-based sparse adaptive channel estimation method is presented in the DD domain, and [17] addresses the peak-to-average power ratio issue by proposing a DD domain MMSE turbo equalizer tailored for single-carrier UWA communications in rapidly time-varying channels. Joint carrier frequency offset estimation gridless channel estimation in OFDM systems is addressed in [18]. In [19], a parametric bilinear generalized approximate message passing-based algorithm to jointly estimate channels and data under nonuniform Doppler shifts in single-carrier systems is developed. However, all these works assume that the UWA channel is narrowband, which does not fully capture its characteristics. UWA channels are inherently wideband and are more accurately modeled as DS-spread channels, as discussed in [4, 9]. Moreover, the relative performance of different waveforms also depends on the receiver processing (e.g., the channel estimation algorithm, channel coding, etc) employed. In this work, we consider all these aspects and present a comprehensive comparative study.

The estimation of DS-spread channels has been actively studied in the literature. In [20], the time-scale is estimated by analyzing the peaks of the matched filter outputs, which is extended to multicarrier systems and validated through experimental studies in [21]. Both works use a linear frequency modulation (LFM) signal as preamble and postamble. Alternatively, [22] and [23, 24] propose the use of OFDM and OCDM, respectively, as a preamble. These works assume a common Doppler scaling factor across all propagation paths, an assumption that does not hold in practical DS-spread UWA channels, where time-scales are path-dependent [25, 26].

In the DS domain (also in the DD domain), the channel response exhibits sparsity due to the small number of propagation paths. This sparsity structure can be exploited to estimate the channel with low pilot/training overhead, by using sparse signal recovery (SSR) algorithms. This is done by considering a grid of points in the DS domain and constructing a large dictionary matrix whose columns capture the effect of a path existing on each grid point on the received samples. Then, one can construct an underdetermined set of linear measurements and solve them using SSR techniques. In [25, 26], a matching pursuit-based SSR algorithm is developed to estimate the DS-spread channel parameters. We developed a variational Bayesian (VB) approach for DS-spread channel estimation algorithm in [1] and extended it to joint channel estimation and data detection in [2]. These methods, including our past work, rely on a fixed dictionary matrix based on finite grid points, assuming on-grid channel parameters (i.e., the channel parameters fall on the grid points determined by the sampling rate at the receiver), which is often unrealistic. In practical DS-spread channels, the multipath delay and scale parameters need not lie on the grid points. This off-grid nature introduces significant challenges in estimating the channel parameters, as the fixed dictionary matrix fails to accurately represent the channel response.

Off-grid DD-spread CE has received much attention in the recent literature [27, 28, 29, 30]. These works employ a first-order Taylor series-based approximation of the dictionary matrix, which works well for DD-spread channels with relatively small delay and Doppler spreads. However, the estimation of off-grid parameters is more challenging in UWA channels with large delay and scale spreads arising from the slow propagation speed of acoustic waves ($1500$ m/s). This motivates the need for developing algorithms to accurately estimate the off-grid delay and scale parameters of UWA channels. To this end, our main contributions in this paper are as follows:

We develop a unified model for the end-to-end communication system that allows us to compare multiple waveforms on a level-playing field (see Sec. II.) Further, we present a unified receiver processing framework that includes both CE and data detection. 2. 2.

Based on the unified framework, we propose a novel two-step iterative off-grid DS-spread CE scheme (see Sec. III.) In the first step, an SSR problem is formulated using a dictionary matrix constructed from a fixed set of grid points in the DS-domain, and the sparse channel parameters are estimated using a VB-based technique (see Sec. IV.) In the second step, the grid points are updated using either a first-order approximation of the dictionary matrix (see Sec. III-A) or a second-order approximation of an objective function (see Sec. III-B), leading to two novel off-grid CE approaches: FVB and SVB. In particular, while prior work has applied a second-order Newtonized step to the orthogonal matching pursuit algorithm for off-grid frequency estimation [31], to the best of our knowledge, no work has integrated second-order approximation into a variational Bayesian framework for delay and scale parameter estimation. This work is the first to do so, and we derive novel update rules for DS grid refinement using Newton iterations. 3. 3.

We benchmark the normalized mean square error (NMSE) of the CE algorithms by deriving a sparsity-aware Cramér-Rao lower bound (CRLB) (see Sec. V.) 4. 4.

Within the unified framework, we design a low complexity VB-based VSSD algorithm (see Sec. VI), which computes the posterior probabilities of the data bits (i.e., the soft symbol estimates). The soft symbols are used to compute the LLRs, which are passed to the channel decoder in coded communications. 5. 5.

We propose a data-aided ICED technique (see Sec. VII) to enhance both CE and data detection. The process begins with an initial estimation of the off-grid DS-spread channel using a short preamble (pilot) through the FVB or SVB algorithm, followed by data symbol detection using either the MMSE or VSSD equalizer. The detected data symbols are then used as virtual pilots to refine the CE. This iterative process of alternating between CE and data detection improves the overall performance. 6. 6.

Finally, using the unified framework, we conduct a fair and comprehensive comparative evaluation of the waveforms in terms of CE and data detection performance under similar channel conditions and receiver complexities. Specifically, the NMSE is evaluated for all the waveforms under ECSIR. For uncoded communications, we compare the BER of all the waveforms under PCSIR and benchmark it against their BER under ECSIR. We also compare the BER under ECSIR with coded communications.

Our simulation results (in Sec. VIII) show that the NMSE in CE of the algorithms developed asymptotically approach the CRLB, when the channel parameters take on-grid values. For practical off-grid channel parameters, the off-grid algorithms outperform the state-of-the-art algorithms. Furthermore, the proposed VSSD equalizer outperforms the conventional MMSE equalizer by $2.2$ dB and $5.4$ dB in uncoded and coded communications, respectively, at a BER of $10^{-2}$. The ICED technique offers a substantial further improvement in performance, with the BER of ICED-based ECSIR degrading only marginally compared to VSSD-based PCSIR.

In terms of the relative performance of the different waveforms, in uncoded communications under PCSIR, with a $1$-tap equalizer,222A $1$-tap equalizer performs subcarrier-by-subcarrier equalization, followed by data detection. ODSS delivers the best performance, followed by OCDM, OTFS, and OFDM. With the VSSD equalizer, ODSS, OCDM, and OTFS exhibit the same performance but outperform OFDM by a large margin. In coded communications, under ECSIR, with MMSE-based LLRs, the performance gap between OFDM and other waveforms significantly reduces. Finally, with VSSD-based LLRs, the BER gap between all the waveforms disappears completely.

Thus, we conclude that different waveforms offer differentiated performances when the receiver complexity is constrained and the channels are harsh. With high-quality soft-symbol estimates and channel coding, or when the channel conditions are more benign, all waveforms perform equally well.

Notation: Matrices and vectors are denoted by bold uppercase and bold lowercase letters, respectively. $(\cdot)^{T}$, $(\cdot)^{H}$, $(\cdot)^{*},\mathrm{Tr(\cdot)}$, $\|\cdot\|$, $\|\cdot\|_{\mathrm{F}}$, $\otimes$, $\odot$, $\lfloor\cdot\rfloor$ and $\mathrm{vec}(\cdot)$ represent the transpose, Hermitian, complex conjugate, trace, $\ell_{2}$ norm, Frobenius norm, Kronecker product, Hadamard product, floor, and vectorization operations, respectively. $\mathbb{E}[\cdot]$ denotes the expectation operator, while $\langle\cdot\rangle_{q(\cdot)}$ represents expectation with respect to the distribution $q(\cdot)$. If the operand is a matrix, $\mathrm{diag}[\cdot]$ denotes a vector containing the diagonal elements of the matrix; if the operand is a vector, it represents a square diagonal matrix with the elements of the vector along the main diagonal. $\mathcal{CN}(\boldsymbol{\mu},\boldsymbol{\Sigma})$ denotes a circularly symmetric complex Gaussian distribution with mean $\boldsymbol{\mu}$ and covariance matrix $\boldsymbol{\Sigma}$, $\Gamma(\lambda_{1},\lambda_{2})$ denotes Gamma distrbution with shape parameter $\lambda_{1}$ and rate parameter $\lambda_{2}$, $\mathcal{U}(\varphi_{1},\varphi_{2})$ denotes uniform distribution with support $[\varphi_{1},\varphi_{2}]$. $\eta^{\boldsymbol{\zeta}}$,$\frac{\partial\boldsymbol{\zeta}}{\partial\eta}$, and $\frac{\partial^{2}\boldsymbol{\zeta}}{\partial\eta^{2}}$ denote element-wise power, first-derivative, and second-derivative of the vector $\boldsymbol{\zeta}$ with respect to a scalar $\eta$, respectively. We use $\mathbf{F}^{\cdot\beta}$ to denote an element-wise power of the entries of $\mathbf{F}$.

II System Model

We consider data transmission in frames, as shown in Fig. 1. Each frame comprises a preamble period of duration $T_{\mathrm{p}}$, a guard interval (with no transmission) of duration $T_{0}$, and a data period of duration $T_{\mathrm{d}}$ seconds. Since the channel we consider is DS-spread, we need to use a guard interval to limit the interference between the preamble and data, and thereby facilitate CE [22]. The data period is used to transmit data symbols that modulate an appropriate waveform. The transmitter and receiver architectures for all the waveforms considered in this work can be depicted using the single block diagram shown in Fig. 2. The system bandwidth is $B$ Hz.

II-A Transmitted Data Signal

The transmitter block, denoted by the matrix $\mathbf{G}$, consists of a waveform-dependent transform followed by a modulator that maps the data symbols $\mathbf{x}$ to a sequence of samples $\mathbf{s}$ to be transmitted over the DS-spread channel. The receiver operation, denoted by the matrix $\mathbf{G}^{H}$, comprises a corresponding demodulator followed by an inverse transform. In the following subsections, we explicitly specify the matrix $\mathbf{G}$ for the OTFS, OFDM, OCDM, and ODSS waveforms.

II-A1 OTFS Waveform

We consider an $M\times N$ OTFS system with subcarrier spacing $\Delta f=B/M$ and total symbol duration of $T_{s}=NT$, where $T=\frac{1}{\Delta f}$. The transmitted DD-domain data symbol matrix $\mathbf{X}_{\mathrm{DD}}\in\mathbb{Q}^{M\times N}$ is converted to the time-frequency (TF) domain matrix $\mathbf{X}_{\mathrm{TF}}\in\mathbb{C}^{M\times N}$ using the inverse symplectic finite Fourier transform as [32]

[TABLE]

where $\mathbb{Q}$ denotes the complex $Q-$QAM constellation, and $\mathbf{F}_{M}\in\mathbb{C}^{M\times M}$ and $\mathbf{F}_{N}\in\mathbb{C}^{N\times N}$ are unitary DFT matrices. The Heisenberg transform converts the TF domain data symbols $X_{\mathrm{TF}}[m,n]$ into the continuous-time domain (TD) passband transmitted signal $s_{\mathrm{T}}(t)$ using a time-limited ($t\in[0,T]$) pulse-shaping filter $g_{\mathrm{tx}}(t)$ as

[TABLE]

where the center frequency of $m^{\mathrm{th}}$ subcarrier $f_{m}=f_{0}+m\Delta f$ with $f_{0}$ be the center frequency of the lowest subcarrier. In discrete-time, we can write (1) in a matrix-vector form as

[TABLE]

where ${\widetilde{\mathbf{G}}}_{\mathrm{T}}=\mathrm{diag}[\widetilde{g}_{\mathrm{T}}[0],\widetilde{g}_{\mathrm{T}}[1],...,\widetilde{g}_{\mathrm{T}}[M-1]]\in\mathbb{C}^{M\times M}$ with $\widetilde{g}_{\mathrm{T}}[m]=g_{\mathrm{tx}}(\frac{mT}{M})e^{j2\pi f_{0}\frac{mT}{M}}$ for $m=0,1,...,M-1$, and $\mathbf{x}=\mathrm{vec}(\mathbf{X}_{\mathrm{DD}})$. Thus, the transmitter matrix for OTFS is $\mathbf{G}_{\mathrm{T}}\triangleq\mathbf{F}_{N}^{H}\otimes\widetilde{\mathbf{G}}_{\mathrm{T}}\in\mathbb{C}^{MN\times MN}$.

In OFDM, we mount the symbols in the TF-domain [3], and the discrete-time transmitted signal is

[TABLE]

where $\mathbf{x}=\mathrm{vec}(\mathbf{X}_{\mathrm{TF}})$, and $\mathbf{G}_{\mathrm{F}}\triangleq\mathbf{I}_{N}^{H}\otimes\widetilde{\mathbf{G}}_{\mathrm{T}}\mathbf{F}_{M}^{H}\in\mathbb{C}^{MN\times MN}$.

II-A2 OCDM Waveform

For an OCDM system with $M$ subcarriers occupying a subcarrier bandwidth $\Delta f$ within the total bandwidth $B=M\Delta f$, the core orthogonal chirp basis for $m^{\mathrm{th}}$ subcarrier can be expressed as $\psi_{m}(t)=e^{j\frac{\pi}{4}}e^{-j\pi\frac{M}{T^{2}}(t-m\frac{T}{M})^{2}},t\in[0,T]$, considering $M$ to be even [7]. Let $N$ time domain OCDM symbols be transmitted over a total duration of $T_{s}=NT$. The continuous-time passband transmitted signal is

[TABLE]

where $X_{\mathrm{C}}[m,n]$ denotes the data symbol mounted on the $m^{\mathrm{th}}$ subcarrier in the $n^{\mathrm{th}}$ symbol, and $f_{c}$ denotes the center frequency. Similar to OTFS, the discrete-time transmitted signal can be written in matrix-vector form as $\mathbf{s}_{\mathrm{C}}=\mathbf{G}_{\mathrm{C}}\mathbf{x}$, where $\mathbf{x}=\mathrm{vec}(\mathbf{X}_{\mathrm{C}}\in\mathbb{Q}^{M\times N})\in\mathbb{C}^{MN\times 1}$. The transmitter matrix for OCDM is $\mathbf{G}_{\mathrm{C}}\triangleq\mathbf{I}_{N}\otimes\widetilde{\mathbf{G}}_{\mathrm{C}}\boldsymbol{\Psi}\in\mathbb{C}^{MN\times MN}$, where $\widetilde{\mathbf{G}}_{\mathrm{C}}=\mathrm{diag}[1,e^{j2\pi f_{c}\frac{T}{M}},...,e^{j2\pi f_{c}(M-1)\frac{T}{M}}]\in\mathbb{C}^{M\times M}$, and the matrix $\boldsymbol{\Psi}\in\mathbb{C}^{M\times M}$ is constructed using the chirp basis, with $(m^{\prime},m)$-th entry $e^{j\frac{\pi}{4}}e^{-j\pi\frac{1}{M}(m^{\prime}-m)^{2}},\,m,m^{\prime}=0,1,...,M-1$.

II-A3 ODSS Waveform

For an ODSS system with $M$ subcarriers, we first choose a constant $q$ and a base subcarrier-width $W$ such that $B=\sum_{m=0}^{M-1}q^{m}W$ [8]. The $m$th subcarrier occupies a bandwidth $q^{m}W$ and a time duration $\frac{1}{q^{m}W}$. Thus, $N(m)=\lfloor q^{m}WT_{s}\rfloor$ symbols can be accommodated on the $m$th subcarrier in a time duration of $T_{s}$. Here, the data symbols are multiplexed onto the 2D Mellin-Fourier (MF) domain, which is mapped into the DS-domain using the ODSS transform matrix $\mathbf{T}_{D}\in\mathbb{C}^{M_{\mathrm{tot}}\times M_{\mathrm{tot}}}$ [8, Eq. (39)] as

[TABLE]

where $M_{\mathrm{tot}}=\sum_{m=0}^{M-1}N(m)$. The vectors $\mathbf{x}\in\mathbb{Q}^{M_{\mathrm{tot}}}$ and $\mathbf{X}_{\mathrm{SD}}\in\mathbb{C}^{M_{\mathrm{tot}}}$ are vectorized versions of the data symbol matrix in the MF- and DS-domains, respectively. The ODSS modulator converts the DS-domain data symbols, $X_{\mathrm{SD}}[n,m]$, to a continuous TD transmitted passband signal, $s_{\mathrm{D}}(t)$, using a time-limited pulse shaping filter, $g_{\mathrm{tx}}(t),t\in[0,\frac{1}{W}]$, as

[TABLE]

where $f_{0}^{D}$ is the center frequency of the lowest ODSS subcarrier. Similar to OTFS, (3) and (4) can be written in matrix-vector form as $\mathbf{s}_{\mathrm{D}}=\mathbf{G}_{\mathrm{D}}\mathbf{x}$, where ODSS transmitter matrix $\mathbf{G}_{\mathrm{D}}\in\mathbb{C}^{M_{\mathrm{tot}}\times M_{\mathrm{tot}}}$ is constructed as in [5, Eq. (16)].

II-B The DS-Spread Channel Model

The received passband signal after propagation through a $P$-path DS-spread channel is given by [4]

[TABLE]

where $w(t)$ is the complex additive white Gaussian noise (AWGN), and the tuple $(h_{p},\tau_{p},\alpha_{p})$ contains the complex channel gain, delay, and time-scale parameters associated with the $p^{\mathrm{th}}$ scattering path. We assume that the channel delays and time-scales are bounded, i.e., $\tau_{p}\in[0,\tau_{\text{max}}]$ and $\alpha_{p}\in[\frac{1}{\alpha_{\text{max}}},\alpha_{\text{max}}]$, where $\tau_{\max}\geq 0$ and $\alpha_{\max}\geq 1$ are the channel delay and scale spreads, respectively. Furthermore, in practical DS-spread channels [9], the delay values are uniformly distributed over the range $[0,\tau_{\text{max}}]$, while the scale values are uniformly distributed in the logarithmic scale, i.e., $\ln\alpha_{p}$ is uniformly distributed over $[-\ln\alpha_{\max},\ln\alpha_{\max}]$. From (5), the received signal can be reformulated as

[TABLE]

where $\mathcal{F}$ and $\mathcal{F}^{-1}$ are the Fourier and inverse-Fourier transform operations, respectively. Using (6), the sampled vector of the received TD signal, after discarding the cyclic prefix (CP), can be derived (similar to [33, Eq. (16)]) as:

[TABLE]

where $\mathbf{H}^{t}\triangleq\sum_{p=1}^{P}h_{p}\sqrt{\alpha_{p}}\mathbf{F}_{M_{\mathrm{d}}}^{H}\mathbf{\Gamma}_{p}\mathbf{F}_{M_{\mathrm{d}}}^{\cdot\frac{1}{\alpha_{p}}}\in\mathbb{C}^{M_{\mathrm{d}}\times M_{\mathrm{d}}}$ is the TD channel matrix and $\mathbf{w}^{t}\in\mathbb{C}^{M_{\mathrm{d}}}$ is the TD AWGN vector. The matrices $\mathbf{\Gamma}_{p}=\frac{1}{\alpha_{p}}\mathrm{diag}[e^{-j2\pi\mathbf{f}\tau_{p}}]\in\mathbb{C}^{M_{\mathrm{d}}\times M_{\mathrm{d}}}$, and $\mathbf{F}_{M_{d}}=\frac{1}{\sqrt{M_{\mathrm{d}}}}e^{-j2\pi\mathbf{f}\mathbf{t}^{T}}$, where $\mathbf{f}=f_{L}+\frac{B}{M_{\mathrm{d}}}[0,...,M_{\mathrm{d}}-1]^{T}$ with lower end frequency $f_{L}$, and $\mathbf{t}=\frac{T_{s}}{M_{\mathrm{d}}}[0,...,M_{\mathrm{d}}-1]^{T}$.

In the sequel, the total number of transmitted data symbols is denoted by $M_{\mathrm{d}}$, which is equal to $MN$ for OTFS, OFDM and OCDM and $M_{\mathrm{tot}}$ for ODSS.

II-C Demodulated Data Symbols

At the receiver, we invert all the operations performed at the transmitter side to get the data symbols in the corresponding domain for all the waveforms. Considering the received pulse shaping filter is the same as the transmitted pulse shaping filter, we get the demodulated data symbol vector for all the waveforms in a unified form as

[TABLE]

where $\mathbf{G}\in\mathbb{C}^{M_{\mathrm{d}}\times M_{\mathrm{d}}}$ is the transmitter matrix corresponding to a waveform; it is denoted by $\mathbf{G}_{\mathrm{T}}$, $\mathbf{G}_{\mathrm{F}}$, $\mathbf{G}_{\mathrm{C}}$, and $\mathbf{G}_{\mathrm{D}}$ for OTFS, OFDM, OCDM, and ODSS, respectively. Also, $\mathbf{w}\triangleq\mathbf{G}^{H}\mathbf{w}^{t}\in\mathbb{C}^{M_{\mathrm{d}}}$ is the AWGN, with $\mathbf{w}\sim\mathcal{CN}(\mathbf{0},\sigma_{\mathrm{d}}^{2}\mathbf{I}_{M_{\mathrm{d}}})$. Further, $\mathbf{H}\in\mathbb{C}^{M_{\mathrm{d}}\times M_{\mathrm{d}}}$ is the effective channel matrix with

[TABLE]

Note that (9) represents the channel in a domain that depends on the modulation waveform: $\mathbf{H}$ is in the TF, DD, and MF domains for OFDM (OCDM), OTFS, and ODSS, respectively.

II-D Preamble-Based CE Problem

We use a preamble signal of short duration $T_{\mathrm{p}}$ as pilots to estimate the DS-spread channel. The preamble consists of $M_{\mathrm{p}}$ known symbols mounted on a waveform. To facilitate CE, we can rewrite the received pilot symbols using (8) as

[TABLE]

where $\mathbf{h}=[h_{1},...,h_{P}]^{T}\in\mathbb{C}^{P},\boldsymbol{\tau}=[\tau_{1},...,\tau_{P}]^{T}\in\mathbb{R}^{P}$, and $\boldsymbol{\alpha}=[\alpha_{1},...,\alpha_{P}]^{T}\in\mathbb{R}^{P}$ contain the true channel gains, delays and scale values, respectively. The matrix $\mathbf{A}_{\mathrm{p},\alpha}(\boldsymbol{\tau},\boldsymbol{\alpha})\in\mathbb{C}^{M_{\mathrm{p}}\times P}$ is constructed using (9) as $\mathbf{A}_{\mathrm{p},\alpha}(\boldsymbol{\tau},\boldsymbol{\alpha})=[\mathbf{a}_{\mathrm{p},\alpha}(\tau_{1},\alpha_{1}),\ldots,\mathbf{a}_{\mathrm{p},\alpha}(\tau_{P},\alpha_{P})],$ where

[TABLE]

for $p=1,2,\ldots,P$, where $\mathbf{G}_{\mathrm{p}}\in\mathbb{C}^{M_{\mathrm{p}}\times M_{\mathrm{p}}}$ is the waveform-dependent transmitter matrix used in the preamble, $\mathbf{\Gamma}_{\mathrm{p},p}\triangleq\frac{1}{\alpha_{p}}\mathrm{diag}[1,e^{-j2\pi\frac{B}{M_{\mathrm{p}}}\tau_{p}},\ldots,e^{-j2\pi\frac{(M_{\mathrm{p}}-1)B}{M_{\mathrm{p}}}\tau_{p}}]\in\mathbb{C}^{M_{\mathrm{p}}\times M_{\mathrm{p}}}$, $\mathbf{x}_{\mathrm{p}}\in\mathbb{Q}^{M_{\mathrm{p}}}$ is the pilot symbol vector, and $\mathbf{w}_{\mathrm{p}}\sim\mathcal{CN}(\mathbf{0},\sigma_{\mathrm{p}}^{2}\mathbf{I}_{M_{\mathrm{p}}})$ is the AWGN.

The advantage of reformulating the received pilot symbols as (10) is that it represents the UWA DS-spread channel in the domain in which it exhibits sparsity, namely, the DS domain. With (10) in hand, we can now use sparsity-promoting techniques, such as VB, to estimate the channel regardless of which waveform is used to generate the preamble signal. Specifically, to estimate the DS-spread channel, we need to determine the number of channel paths and each channel path’s delay, scale, and gain parameters. However, from (11), we see that the relationship between the DS channel parameters $(\tau_{p},\alpha_{p})_{p=1}^{P}$ and measurement matrix $\mathbf{A}_{\mathrm{p},\alpha}(\boldsymbol{\tau},\boldsymbol{\alpha})$ is highly non-linear, which makes jointly estimating them from the pilot measurements in (10) challenging. We address this in the following section.

III Off-Grid DS-Spread CE Framework

We adopt an iterative approach to estimate the channel parameters $(h_{p},\tau_{p},\alpha_{p})_{p=1}^{P}$. We first form a coarse sampling grid of the DS parameters and compute the corresponding measurement matrix using (11). Then, each iteration involves the following two steps: (a) Estimation stage: we apply the VB technique to estimate the delay and scale parameters on the sampling grid along with the corresponding channel gain values. (b) Refinement stage: we refine the delay and scale parameters and update the corresponding grid points as well as the dictionary matrix.

Recall that the continuous-valued delay and scale parameters span $[0,\tau_{\max}]$ and $[\frac{1}{\alpha_{\max}},\alpha_{\max}]$, respectively. A commonly adopted approach to estimating these parameters is to form a finite-sized grid in the DS domain. Then, considering these grid points as candidate channel path parameters, we look for the sparsest possible subset of the grid points and corresponding path gains so that (10) is satisfied. To elaborate, let $N_{\tau}$ and $M_{\alpha}$ denote the size of the sampling grid along the delay and scale axes, respectively. Let $\bar{\boldsymbol{\tau}}=[\bar{\tau}_{0},\ldots,\bar{\tau}_{N_{\tau}M_{\alpha}-1}]^{T}\in\mathbb{R}^{N_{\tau}M_{\alpha}}$ represent the grid points along the delay axis, with $\bar{\tau}_{0},\ldots,\bar{\tau}_{N_{\tau}M_{\alpha}-1}\in[0,\tau_{\max}]$, and $\bar{\boldsymbol{\alpha}}=[\bar{\alpha}_{0},\ldots,\bar{\alpha}_{N_{\tau}M_{\alpha}-1}]^{T}\in\mathbb{R}^{N_{\tau}M_{\alpha}}$ represent the grid points along the scale axis, with $\bar{\alpha}_{0},\ldots,\bar{\alpha}_{N_{\tau}M_{\alpha}-1}\in[1/\alpha_{\max},\alpha_{\max}]$. In correspondence with the respective distributions of the delay and scale values [9], we sample the delay axis linearly and the scale axis geometrically. Due to this, instead of considering a grid in the scale axis, we consider an equivalent log-scale version, $\bar{\boldsymbol{\omega}}=[\bar{\omega}_{0},...,\bar{\omega}_{N_{\tau}M_{\alpha}-1}]^{T}\in\mathbb{R}^{N_{\alpha}M_{\alpha}}$, with $\bar{\omega}_{i}\triangleq\ln\bar{\alpha}_{i}/\ln q_{\alpha}$, where $q_{\alpha}$ is a constant (that converts the logarithm from base $e$ to base $q_{\alpha}$).

In addition, note that the matrix $\mathbf{A}_{\mathrm{p},\alpha}(\boldsymbol{\tau},\boldsymbol{\alpha})$ in (10) can also be written as a function of $\boldsymbol{\omega}$ as $\mathbf{A}_{\mathrm{p},\omega}(\boldsymbol{\tau},\boldsymbol{\omega})=\mathbf{A}_{\mathrm{p},\alpha}(\boldsymbol{\tau},q_{\alpha}^{\boldsymbol{\omega}})$. Thus, in the estimation stage, we use the pilot symbols to formulate an SSR problem (on the sampling grid) as

[TABLE]

where the bar above the variable (e.g., $\bar{\boldsymbol{\tau}}$) is used to denote variable parameters corresponding to the dictionary matrix; variables without the bar (e.g., $\boldsymbol{\tau}$) represent the true channel parameters. Also, $\mathbf{A}_{\mathrm{p},\omega}(\bar{\boldsymbol{\tau}},\bar{\boldsymbol{\omega}})\in\mathbb{C}^{M_{\mathrm{p}}\times N_{\tau}M_{\alpha}}$ denotes the matrix $\mathbf{A}_{\mathrm{p},\alpha}({\boldsymbol{\tau}},{\boldsymbol{\alpha}})$ in (10) evaluated using the grid-based delay parameters $\bar{\boldsymbol{\tau}}$ and scale parameters $q_{\alpha}^{\bar{\boldsymbol{\omega}}}$ instead of $\boldsymbol{\tau}$ and $\boldsymbol{\alpha}$, respectively. As mentioned above, initially, we consider a uniformly spaced sampling grid along the delay axis with the delay resolution $r_{\tau}=\frac{\tau_{\max}}{N_{\tau}}$ and spanning the interval $[0,\tau_{\max}]$. For sampling the scale axis, we consider $M_{\alpha}$ to be odd (the extension to the case of even $M_{\alpha}$ is straightforward) and choose $q_{\alpha}$ such that $\alpha_{\max}=q_{\alpha}^{\frac{M_{\alpha}-1}{2}}$, to obtain an equally spaced sampling grid along the log-scale axis with resolution $r_{\omega}=1$. Thus, we get the initial delay grid vector $\bar{\boldsymbol{\tau}}^{(0)}=[\bar{\tau}_{0}^{(0)},...,\bar{\tau}_{N_{\tau}M_{\alpha}-1}^{(0)}]\in\mathbb{R}^{N_{\tau}M_{\alpha}}$ and the log-scale grid vector $\bar{\boldsymbol{\omega}}^{(0)}=[\bar{\omega}_{0}^{(0)},...,\bar{\omega}_{N_{\tau}M_{\alpha}-1}^{(0)}]\in\mathbb{R}^{N_{\tau}M_{\alpha}}$, where $\bar{\tau}_{n^{\prime}M_{\alpha}+m^{\prime}}^{(0)}=n^{\prime}r_{\tau}$ and $\bar{\omega}_{n^{\prime}M_{\alpha}+m^{\prime}}^{(0)}=\left(-\frac{M_{\alpha}-1}{2}+m^{\prime}\right)r_{\omega}$, $n^{\prime}\in\{0,...,N_{\tau}-1\}$ and $m^{\prime}\in\{0,...,M_{\alpha}-1\}$.

Note that $\bar{\mathbf{h}}\in\mathbb{C}^{N_{\tau}M_{\alpha}}$ is likely to be close to a sparse vector since the channel contains only $P\ll N_{\tau}M_{\alpha}$ physical paths, and only components of $\bar{\mathbf{h}}$ that correspond to these paths will be nonzero. Hence, SSR methods can be applied to solve for $\bar{\mathbf{h}}$ using (12). However, the true delay and scale values may not lie on the grid values used to construct the dictionary matrix $\mathbf{A}_{\mathrm{p},\omega}(\bar{\boldsymbol{\tau}},\bar{\boldsymbol{\omega}})$, leading to the estimated $\bar{\mathbf{h}}$ being only approximately sparse. This mismatch between the true channel parameters and the grid values can be alleviated by performing dictionary refinement to update the dictionary matrix, where we update the grid by estimating the off-grid delay and scale parameters of the DS-spread channel.

Our proposed off-grid DS-spread CE technique works as follows. First, starting from the dictionary matrix $\mathbf{A}_{\mathrm{p},\omega}(\bar{\boldsymbol{\tau}}^{(0)},\bar{\boldsymbol{\omega}}^{(0)})$, we estimate the sparse channel vector $\bar{\mathbf{h}}$ using a VB-based technique (see Sec. IV). Next, in the dictionary refinement step, we use the estimated $\bar{\mathbf{h}}$ to update the grid points $\bar{\boldsymbol{\tau}}$ and $\bar{\boldsymbol{\omega}}$. We update $\bar{\boldsymbol{\tau}}$ and $\bar{\boldsymbol{\omega}}$ using either a first-order approximation of the basis of the dictionary matrix or a second-order approximation of an objective function and use the updated grid to recompute the dictionary matrix. Note that, to reduce complexity, once $\bar{\mathbf{h}}$ is estimated, we can truncate entries of $\bar{\mathbf{h}}$ whose magnitude is below a threshold to zero and remove them from $\bar{\mathbf{h}}$, resulting in a truncated vector $\widetilde{\mathbf{h}}\in\mathbb{C}^{\hat{P}}$, where $\hat{P}$ denotes estimated number of paths. Correspondingly, we truncate $\bar{\boldsymbol{\tau}}$ and $\bar{\boldsymbol{\omega}}$, and denote the truncated versions by $\widetilde{\boldsymbol{\tau}}\in\mathbb{R}^{\hat{P}}$, $\widetilde{\boldsymbol{\omega}}\in\mathbb{R}^{\hat{P}}$, respectively. Next, with the refined dictionary in hand, we re-estimate $\bar{\mathbf{h}}$, and repeat the process iteratively. For clarity of presentation, before describing the overall recipe, we first briefly discuss the dictionary refinement procedure in the following two subsections.

III-A First-Order Approximation (FOA) of the Basis

Due to the finite (and possibly coarse) sampling grid, the true delay and scale parameters need not lie on grid points. Hence, we seek to refine the grid points progressively so that the true channel parameters eventually align with the grid. We accomplish this by estimating the off-grid components of the parameters and using them to recompute the dictionary matrix. Now, given the $\widetilde{\boldsymbol{\tau}}$ and $\widetilde{\boldsymbol{\omega}}$ returned by the VB algorithm (after the thresholding step), let $\boldsymbol{\beta_{\widetilde{\tau}}}\in\mathbb{R}^{\hat{P}}$ and $\boldsymbol{\beta_{\widetilde{\omega}}}\in\mathbb{R}^{\hat{P}}$ denote the additive correction to be applied to the delay and Doppler grid point values, respectively. Then, using a first-order linear approximation of $\mathbf{A}_{\mathrm{p},\omega}(\widetilde{\boldsymbol{\tau}},\widetilde{\boldsymbol{\omega}})\in\mathbb{C}^{M_{\mathrm{p}}\times\hat{P}}$, we can write

[TABLE]

where the matrices $\widetilde{\mathbf{A}}_{\mathrm{p},\omega}^{(0)}=\mathbf{A}_{\mathrm{p},\omega}(\widetilde{\boldsymbol{\tau}}^{(0)},\widetilde{\boldsymbol{\omega}}^{(0)})\in\mathbb{C}^{M_{\mathrm{p}}\times\hat{P}}$, $\mathbf{B}_{\mathrm{p}}=[\mathbf{b}_{\mathrm{p}}(\widetilde{\tau}_{0}^{(0)},\widetilde{\omega}_{0}^{(0)}),\ldots,\mathbf{b}_{\mathrm{p}}(\widetilde{\tau}_{\hat{P}-1}^{(0)},\widetilde{\omega}_{\hat{P}-1}^{(0)})]\in\mathbb{C}^{M_{\mathrm{p}}\times\hat{P}}$ with $\mathbf{b}_{\mathrm{p}}(\tau,\omega)=\frac{\partial{\mathbf{a}_{\mathrm{p},\omega}(\tau,\omega)}}{\partial{\tau}}\in\mathbb{C}^{M_{\mathrm{p}}}$. Also, the matrix $\mathbf{C}_{\mathrm{p}}=[\mathbf{c}_{\mathrm{p}}(\widetilde{\tau}_{0}^{(0)},\widetilde{\omega}_{0}^{(0)}),\ldots,\mathbf{c}_{\mathrm{p}}(\widetilde{\tau}_{\hat{P}-1}^{(0)},\widetilde{\omega}_{\hat{P}-1}^{(0)})]\in\mathbb{C}^{M_{\mathrm{p}}\times\hat{P}}$ with $\mathbf{c}_{\mathrm{p}}(\tau,\omega)=\frac{\partial{\mathbf{a}_{\mathrm{p},\omega}(\tau,\omega)}}{\partial{\omega}}\in\mathbb{C}^{M_{\mathrm{p}}}$, and $\mathbf{a}_{\mathrm{p},\omega}(\tau,\omega)$ is a column of $\mathbf{A}_{\mathrm{p},\omega}(\boldsymbol{\tau},\boldsymbol{\omega})$ computed using parameters $\tau$ and $\omega$. The vectors $\boldsymbol{\beta_{\widetilde{\tau}}}=[\beta_{\widetilde{\tau}_{0}},\ldots,\beta_{\widetilde{\tau}_{\hat{P}-1}}]^{T}$ and $\boldsymbol{\beta_{\widetilde{\omega}}}=[\beta_{\widetilde{\omega}_{0}},\ldots,\beta_{\widetilde{\omega}_{\hat{P}-1}}]^{T}$ represent the vectors containing the off-grid components. We will use (13) along with $\mathbf{y}_{\mathrm{p}}$ in (12) to obtain $\boldsymbol{\beta}_{\widetilde{\tau}}^{(j)}$ and $\boldsymbol{\beta}_{\widetilde{\omega}}^{(j)}$, the estimates of $\boldsymbol{\beta_{\widetilde{\tau}}}$ and $\boldsymbol{\beta_{\widetilde{\omega}}}$, respectively, in the $j$th iteration (see Sec. IV-A1 for explicit expressions.) Then, we update the $\hat{P}$ delay and log-scale grid points at the $(j+1)$th iteration as

[TABLE]

III-B Second-Order Approximation (SOA) of Objective Function

Based on (12), we define the objective function to be minimized as $f(\widetilde{\boldsymbol{\tau}},\widetilde{\boldsymbol{\omega}})=\|\mathbf{y}_{\mathrm{p}}-\mathbf{A}_{\mathrm{p},\omega}(\widetilde{\boldsymbol{\tau}},\widetilde{\boldsymbol{\omega}})\widetilde{\mathbf{h}}\|^{2}$. We update $\widetilde{\boldsymbol{\tau}}$ and $\widetilde{\boldsymbol{\omega}}$ element-wise, in a coordinate-descent manner. Let $\widetilde{\boldsymbol{\tau}}_{\neq l}\in\mathbb{R}^{(\hat{P}-1)}$ and $\widetilde{\boldsymbol{\omega}}_{\neq l}\in\mathbb{R}^{(\hat{P}-1)}$ denote the vector $\widetilde{\boldsymbol{\tau}}$ and $\widetilde{\boldsymbol{\omega}}$, respectively, after deleting the $l$th element. We define the objective function for the $l$th delay parameter at the $(j+1)$th update as $f_{\widetilde{\tau}_{l}}(\tau)=\|\mathbf{y}_{\mathrm{p}}-\mathbf{A}_{\mathrm{p},\omega}(\widetilde{\boldsymbol{\tau}}^{(j)}_{\neq l},\widetilde{\boldsymbol{\omega}}^{(j)}_{\neq l})\widetilde{\mathbf{h}}_{\neq l}-{\mathbf{a}}_{\mathrm{p},\omega}(\tau,\widetilde{\omega}_{l}^{(j)})\widetilde{h}_{l}\|^{2}$, and we minimize it with respect to $\tau$ to get the $(j+1)$th update of $\widetilde{\tau}_{l}$. Using a second-order Newton’s iteration to optimize the element-wise objective functions, we obtain the updates for $(\widetilde{\boldsymbol{\tau}}^{(j+1)},\widetilde{\boldsymbol{\omega}}^{(j+1)})$ as

[TABLE]

respectively, where f_{\widetilde{\tau}_{l}}^{\prime}({\widetilde{\tau}_{l}}^{(j)})=\frac{\partial f_{\widetilde{\tau}_{l}}(\tau)}{\partial\tau}\Big{|}_{\tau=\widetilde{\tau}_{l}^{(j)}}, f_{\widetilde{\tau}_{l}}^{\prime\prime}(\widetilde{\tau_{l}}^{(j)})=\frac{\partial^{2}f_{\widetilde{\tau}_{l}}(\tau)}{\partial\tau^{2}}\Big{|}_{\tau=\widetilde{\tau}_{l}^{(j)}}, f_{\widetilde{\omega}_{l}}^{\prime}(\widetilde{\omega_{l}}^{(j)})=\frac{\partial f_{\widetilde{\omega}_{l}}(\omega)}{\partial\omega}\Big{|}_{\omega=\widetilde{\omega}^{(j)}_{l}}, and f_{\widetilde{\omega}_{l}}^{\prime\prime}({\widetilde{\omega}_{l}}^{(j)})=\frac{\partial^{2}f_{\widetilde{\omega}_{l}}(\omega)}{\partial\omega^{2}}\Big{|}_{\omega=\widetilde{\omega}_{l}^{(j)}}; $f_{\widetilde{\omega}_{l}}(\omega)$ is defined similar to $f_{\widetilde{\tau}_{l}}(\tau)$. We provide explicit expressions for the above updates in Sec. IV-A2.

Regardless of the update methods, the off-grid delay and log-scale components are bounded between $[-\frac{r_{\tau}}{2},\frac{r_{\tau}}{2}]$ and $[-\frac{r_{\omega}}{2},\frac{r_{\omega}}{2}]$, respectively. In the FOA method, the updates to the $\hat{P}$ grid point values are performed using (14) and (15), while in the SOA method, they are carried out using (16) and (17). These refined grid point values are then used to update the $\hat{P}$ columns of the dictionary matrix, resulting in the new dictionary $\widetilde{\mathbf{A}}_{\mathrm{p},\omega}^{(j+1)}=\mathbf{A}_{\mathrm{p},\omega}(\widetilde{\boldsymbol{\tau}}^{(j+1)},\widetilde{\boldsymbol{\omega}}^{(j+1)})$. The remaining $(N_{\tau}M_{\alpha}-\hat{P})$ columns of the dictionary matrix are not altered and we get the $(j+1)$th update of the entire dictionary matrix, $\bar{\mathbf{A}}_{\mathrm{p},\omega}^{(j+1)}=\mathbf{A}_{\mathrm{p},\omega}(\bar{\boldsymbol{\tau}}^{(j+1)},\bar{\boldsymbol{\omega}}^{(j+1)})$.

The following section explains the details of the proposed VB-based off-grid DS-spread CE algorithm.

IV VB-Based Off-grid DS-Spread CE

Recall that the vector $\bar{\mathbf{h}}$ in (12) is sparse because the number of paths, $N_{p}$, is much smaller than the dimension of $\bar{\mathbf{h}}$ which equals the number of grid points, $N_{\tau}M_{\alpha}$. A Bayesian approach to estimating $\bar{\mathbf{h}}$ from (12) involves assigning a prior distribution to $\bar{\mathbf{h}}$ that promotes sparsity while simplifying the computation of the posterior. A commonly used method employs a two-stage hierarchical prior [34]. In the first stage, a complex Gaussian distribution is imposed to the entries of $\bar{\mathbf{h}}$:

[TABLE]

where $\boldsymbol{\delta}=[\delta_{0},...,\delta_{N_{\tau}M_{\alpha}-1}]\in\mathbb{R}_{+}^{N_{\tau}M_{\alpha}}$ is the precision vector, and $\delta_{l}$ is the precision parameter in the distribution of $\bar{h}_{l}$. In the second stage, $\delta_{l}$ are assumed to follow an independent and identically distributed (i.i.d.) Gamma distribution:

[TABLE]

Also, the additive noise is Gaussian distributed, and the inverse of its variance is modeled as $\gamma\sim\Gamma(\epsilon_{3},\epsilon_{4})$, which eliminates the need to know the noise variance beforehand. Choosing $\epsilon_{1},\ldots,\epsilon_{4}$ to be small numbers (of the order $10^{-6}$) renders the prior non-informative and results in sparse solutions. However, the performance is not sensitive to the precise values chosen; any small numbers will do. Let $\mathbf{z}\triangleq\{\bar{\mathbf{h}},\boldsymbol{\delta},\gamma\}$ denote the latent variables. The MMSE estimate of $\mathbf{z}$ is the mean of the conditional posterior distribution $p(\mathbf{z}|\mathbf{y}_{\mathrm{p}})$. Since the exact posterior is hard to compute, we approximate it using a factorized family of distributions $q(\mathbf{z})$ i.e., $q(\mathbf{z})=q(\bar{\mathbf{h}})q(\boldsymbol{\delta})q(\gamma)$.333Note that, for brevity, we use the same notation $q(\cdot)$ for all factors and distinguish them using the argument itself. The optimal marginal distribution for the $k$th element ($k\in\{1,2,3\}$) of $\mathbf{z}$ is obtained by minimizing the Kullback-Leibler divergence, $\text{KL}(p(\mathbf{z}|\mathbf{y}_{\mathrm{p}})||q(\mathbf{z}))$, and is given by [35]:

[TABLE]

where $\mathbb{E}_{u\neq k}[\cdot]$ denotes the expectation with respect to all factorized distributions of $\mathbf{z}$ except $q(\mathbf{z}_{k})$, and $\propto$ denotes equality up to an additive normalization constant. Using Bayes’ theorem, the joint distribution $p(\mathbf{z},\mathbf{y}_{\mathrm{p}})$ can be expressed as

[TABLE]

The approximate posterior distribution is found through an iterative update of the marginals as follows:

[TABLE]

where $q^{(j)}(\cdot)$ denotes the marginals in the $j$th iteration. Detailed derivations for the marginal updates are presented below.

Update of $q(\bar{\mathbf{h}})$: From (18) and (21), (22) is derived as

[TABLE]

where $\widehat{\gamma}^{(j)}=\langle\gamma\rangle_{q^{(j)}(\gamma)}$ and $\widehat{\boldsymbol{\Delta}}^{(j)}=\text{diag}[\widehat{\boldsymbol{\delta}}^{(j)}]$ with $\widehat{\boldsymbol{\delta}}^{(j)}=\langle\boldsymbol{\delta}\rangle_{q^{(j)}(\boldsymbol{\delta})}$. Hence, in the $(j+1)$th iteration, $\bar{\mathbf{h}}$ follows a complex Gaussian distribution with mean and covariance given by

[TABLE]

respectively. 2. 2.

Update of $q(\boldsymbol{\delta})$: From (18), (19) and (21), (23) becomes

[TABLE]

Hence, in the $(j+1)$th iteration, the factor $q(\delta_{l})$, the marginal distribution of the $l$th entry of $\boldsymbol{\delta}$, is a Gamma distribution with mean given by

[TABLE]

where $\langle|\bar{h}_{l}|^{2}\rangle_{q^{(j+1)}(\bar{\mathbf{h}})}=|\boldsymbol{\mu}^{(j+1)}_{\bar{\mathbf{h}}_{l}}|^{2}+\boldsymbol{\Sigma}^{(j+1)}_{\bar{\mathbf{h}}_{(l,l)}}$ for $\boldsymbol{\mu}^{(j+1)}_{\bar{\mathbf{h}}_{l}}$ and $\boldsymbol{\Sigma}^{(j+1)}_{\bar{\mathbf{h}}_{(l,l)}}$ being the $l$th entry of $\boldsymbol{\mu}^{(j+1)}_{\bar{\mathbf{h}}}$ and $(l,l)$th entry of $\mathbf{\Sigma}_{\bar{\mathbf{h}}}^{(j+1)}$, respectively. 3. 3.

Update of $q(\gamma)$: Similarly, (24) can be written as

[TABLE]

Hence, in the $(j+1)$th iteration, $\gamma$ follows a Gamma distribution with mean $\widehat{\gamma}^{(j+1)}=\langle\gamma\rangle_{q^{(j+1)}(\gamma)}$, where

[TABLE]

with $\langle||\mathbf{y}_{\mathrm{p}}-\bar{\mathbf{A}}_{\mathrm{p},\omega}^{(j)}\bar{\mathbf{h}}||^{2}\rangle_{q^{(j+1)}(\bar{\mathbf{h}})}=\|\mathbf{y}_{\mathrm{p}}-\bar{\mathbf{A}}_{\mathrm{p},\omega}^{(j)}\boldsymbol{\mu}^{(j+1)}_{\bar{\mathbf{h}}}\|^{2}\\ +\mathrm{Tr}(\bar{\mathbf{A}}_{\mathrm{p},\omega}^{(j)}\boldsymbol{\Sigma}^{(j+1)}_{\bar{\mathbf{h}}}\bar{\mathbf{A}}_{\mathrm{p},\omega}^{(j)^{H}})$.

Note that the updates to $q(\boldsymbol{\delta})$ and $q(\gamma)$ are independent of each other, so these steps can be executed in parallel.

To summarize, by iterating (22), (23) and (24), we obtain the posterior distribution of $\bar{\mathbf{h}}$ upon convergence. The mean of the posterior will have $\hat{P}$ nonzero entries and identifies the candidate grid locations near which the actual channel DS values lie. The next step is to refine the DS estimates by updating the dictionary, as detailed in the following subsection.

IV-A Dictionary Matrix Update

To update the dictionary matrix, we focus on the $\hat{P}$ paths identified by the VB-based CE algorithm. This involves the $\hat{P}$ columns of the dictionary matrix, ${\widetilde{\mathbf{A}}}_{\mathrm{p},\omega}\in\mathbb{C}^{M_{\mathrm{p}}\times\hat{P}}$, the corresponding $\hat{P}$ entries in the sparse channel vector $\widetilde{\mathbf{h}}\in\mathbb{C}^{\hat{P}}$, its mean $\boldsymbol{\mu}_{\widetilde{h}}\in\mathbb{C}^{\hat{P}}$, covariance matrix $\boldsymbol{\Sigma}_{\widetilde{\mathbf{h}}}\in\mathbb{C}^{\hat{P}\times\hat{P}}$, and the $\hat{P}$ entries in the precision vector, delay vector, and log-scale vector, $\widetilde{\boldsymbol{\delta}}\in\mathbb{R}_{+}^{\hat{P}}$, $\widetilde{\boldsymbol{\tau}}\in\mathbb{R}_{+}^{\hat{P}}$ and $\widetilde{\boldsymbol{\omega}}\in\mathbb{R}^{\hat{P}}$, respectively. Refinement of the grid-points that constitute the dictionary matrix is carried out using two approaches: the FOA and SOA methods, described in Sec.III-A and Sec.III-B, respectively. We now derive the corresponding update rules as follows:

IV-A1 FOA Method

We consider a uniform prior on the off-grid delay and log-scale variables, as follows

[TABLE]

We consider $\boldsymbol{\beta}_{\widetilde{\tau}}$ and $\boldsymbol{\beta}_{\widetilde{\omega}}$ to be latent variables and concatenate them with $\widetilde{\mathbf{z}}=\{\widetilde{\mathbf{h}},\widetilde{\boldsymbol{\delta}},\gamma\}$ to obtain the new set of latent variables $\mathbf{z}_{\text{new}}=\{\widetilde{\mathbf{z}},\boldsymbol{\beta}_{\widetilde{\tau}},\boldsymbol{\beta}_{\widetilde{\omega}}\}$. Assuming a factorized distribution $q(\mathbf{z}_{\text{new}})=q(\widetilde{\mathbf{z}})q(\boldsymbol{\beta}_{\widetilde{\tau}})q(\boldsymbol{\beta}_{\widetilde{\omega}})$, the update of the marginals $q(\boldsymbol{\beta}_{\widetilde{\tau}})$ and $q(\boldsymbol{\beta}_{\widetilde{\omega}})$ can be found from

[TABLE]

respectively. From (21) and (29), we can rewrite (30) as

[TABLE]

where $\mathbf{P}^{(j+1)}_{\widetilde{\tau}}$ is a positive semi-definite matrix given by

[TABLE]

and the vector $\mathbf{v}^{(j+1)}_{\widetilde{\tau}}$ is given by (35), at the top of the page, and other notation (such as $\mathbf{B}_{\mathrm{p}}$ and $\mathbf{C}_{\mathrm{p}}$) are defined after (13).

Note that, we obtain (33) from (32) by substituting the expression for $\widetilde{\mathbf{A}}_{\mathrm{p},\omega}^{(j)}$ from (13) and performing further simplifications following [27, 29]. Hence, if $\mathbf{P}^{(j+1)}_{\widetilde{\tau}}$ is invertible, $\boldsymbol{\beta}_{\widetilde{\tau}}$ follows a Gaussian distribution with mean

[TABLE]

Otherwise, the $\hat{p}$th element in $\boldsymbol{\beta}^{(j+1)}_{\widetilde{\tau}}$ follows a Gaussian distribution with mean

[TABLE]

where $\hat{p}\in\{1,\ldots,\hat{P}\}$, $[{\mathbf{v}}]_{\hat{p}}$ denotes the $\hat{p}$th entry of ${\mathbf{v}}$, and $[\mathbf{P}]_{\neq\hat{p}}$ is a matrix that is the same as $\mathbf{P}$, but with the $\hat{p}$th column discarded; $[\mathbf{P}]_{(:,\hat{p})}$ and $[\mathbf{P}]_{(\hat{p},:)}$ denote the $\hat{p}$th column and row of $\mathbf{P}$, respectively. Finally, if $[\mathbf{P}^{(j+1)}_{\widetilde{\tau}}]_{(\hat{p},\hat{p})}$ equals zero, we do not update the $\hat{p}$th entry of $\boldsymbol{\beta}_{\widetilde{\tau}}$ in the $(j+1)$th iteration.

Similarly, we can show that the off-grid log-scale variable follows a Gaussian distribution with mean

[TABLE]

if $\mathbf{P}^{(j+1)^{-1}}_{\widetilde{\omega}}$ is invertible, where $\mathbf{P}_{\widetilde{\omega}}$ is given by

[TABLE]

and $\mathbf{v}_{\widetilde{\omega}}$ is given by (40) on the next page.

Otherwise, $\hat{p}$th element in $\boldsymbol{\beta}^{(j+1)}_{\widetilde{\omega}}$ follows a Gaussian distribution with mean

[TABLE]

We now provide intuitive explanations for the update rule of the off-grid delay (36) and log-scale variables (38). Substitutng the FOA expression for $\widetilde{\mathbf{A}}_{\mathrm{p},\omega}^{(j)}$ from (13) into (12) yields

[TABLE]

which can be rewritten as

[TABLE]

where $\widetilde{\mathbf{h}}\sim\mathcal{CN}(\boldsymbol{\mu}_{\widetilde{\mathbf{h}}},\boldsymbol{\Sigma}_{\widetilde{\mathbf{h}}})$, and $\mathbf{w}_{\mathrm{p}}\sim\mathcal{CN}(\mathbf{0},\frac{1}{\gamma}\mathbf{I}_{M_{\mathrm{p}}})$, From (42), if we were to compute the Bayes’ optimal MMSE estimator for $\boldsymbol{\beta_{\widetilde{\tau}}}$, we get (36). In other words, (36) represents the Bayes’ optimal MMSE estimator of $\boldsymbol{\beta_{\widetilde{\tau}}}$ given $\boldsymbol{\beta_{\widetilde{\omega}}}$. Similarly, we obtain (38) as the Bayes’ optimal MMSE estimator of $\boldsymbol{\beta}_{\widetilde{\omega}}$ given the value of $\boldsymbol{\beta_{\widetilde{\tau}}}$.

IV-A2 SOA Method

As mentioned in Sec. III-B, here, we calculate element-wise updates of $\widetilde{\boldsymbol{\tau}}^{(j+1)}$ and $\widetilde{\boldsymbol{\omega}}^{(j+1)}$ using Newton’s method. The $\hat{p}$th element in $\widetilde{\boldsymbol{\tau}}^{(j+1)}$ and $\widetilde{\boldsymbol{\omega}}^{(j+1)}$ can be calculated by solving the optimization problem (43) and (44), respectively, on the next page.

We find the expressions for element-wise update using (16) and (17) as

[TABLE]

where $g_{\widetilde{\tau}}^{1}(\widetilde{\tau}_{\hat{p}}^{(j)}),g_{\widetilde{\tau}}^{2}(\widetilde{\tau}_{\hat{p}}^{(j)}),g_{\widetilde{\omega}}^{1}(\widetilde{\omega}_{\hat{p}}^{(j)}),$ and $g_{\widetilde{\omega}}^{2}(\widetilde{\omega}_{\hat{p}}^{(j)})$ are the first and second order derivates of the objective function with respect to the parameters, and are explicitly provided in (47), (49), (50), and (52), respectively, on the bottom of a page.

In the FOA method, the columns of the dictionary are approximated around the current grid points as a linear function of the off-grid parameters, $\boldsymbol{\beta_{\widetilde{\tau}}}$ and $\boldsymbol{\beta_{\widetilde{\omega}}}$, while ignoring higher-order effects. As a result, this method performs well only when the off-grid parameter values are small; however, if the true channel parameters (delay and scale) lie significantly off the grid, its accuracy degrades. In contrast, the SOA method uses a second-order approximation of the cost function with respect to the off-grid parameters and updates them using Newton iterations. By accounting for second-order terms, it can correct larger offsets and yield better updates.

Our algorithm for VB-based off-grid DS-spread CE is summarized in Algorithm 1. We iteratively update the marginal distributions of all the latent variables to obtain the approximate posterior distribution. The mean and covariance matrix of $\bar{\mathbf{h}}$ are updated using (25) and (26). The hyper-parameters are simultaneously updated according to (27) and (28). Then, the grid points of the delay and log-scale parameters are updated according to (14) or (45) and (15) or (46), respectively, for the FOA and SOA methods. This is followed by refining the dictionary matrix. The algorithm terminates when the normalized change in the mean of hyper-parameter $\boldsymbol{\delta}$ becomes smaller than a predefined convergence threshold $\epsilon$ or the iteration counter encounters the maximum number of iterations $J_{\max}$. We get the estimated sparse channel vector, delay vector, and log-scale vector as the algorithm’s outputs. Finally, the estimated channel parameters are used to compute the effective channel matrix using (9).

Note that, to solve the optimization problem (43) and (44), we use a single Newton update using (45) and (46), respectively, in Algorithm 1. This is justified because these optimization problems are intermediate steps in an overall iterative procedure. Furthermore, since $\widetilde{\boldsymbol{\tau}}$ and $\widetilde{\boldsymbol{\omega}}$ are updated in a block coordinate descent manner, it is not necessary to run the Newton updates to convergence at each step.

The VB-based iterative process follows the minorization-maximization principle of optimization, which ensures convergence to a stationary point [36]. The convergence of the first-order approximation-based algorithm in the VB framework (i.e., FVB) is discussed in [27, 29] for DD-spread CE; the approach readily extends to DS-spread CE also. The convergence of the second-order approximation-based Newtonized OMP (NOMP) algorithm, as established in [31], can also be extended to the VB framework (i.e., SVB) considered in this paper. Further, the per-iteration complexity of the FVB algorithm is $\mathcal{O}(M_{\mathrm{p}}^{3}+M_{\mathrm{p}}N_{\tau}M_{\alpha}+M_{\mathrm{p}}^{2}N_{\tau}M_{\alpha}+N_{\tau}^{2}M_{\alpha}^{2}M_{\mathrm{p}}+N_{\tau}M_{\alpha}+M_{\mathrm{p}}\hat{P}^{2}+\hat{P}^{3})$. The complexity of SVB is similar, except that the $\hat{P}^{3}$ term is absent due to coordinate descent updates on the individual entries of the grid points. To reduce the complexity of (26), the matrix inversion lemma is used. Since $M_{\tau}N_{\alpha}\gg M_{\mathrm{p}}$ and $M_{\tau}N_{\alpha}\gg\hat{P}$, the overall complexity of both FVB and SVB simplifies to $\mathcal{O}(N_{\tau}^{2}M_{\alpha}^{2}M_{\mathrm{p}})$.

V Performance Benchmark

We present a theoretical baseline to assess the performance of the proposed DS-spread channel estimators for on-grid channel parameters. To that end, we derive a CRLB for the estimation of $\mathbf{H}$ in (9). The bound is computed by first finding the Bayesian information matrix (BIM) [37, 38] corresponding to $\bar{\mathbf{h}}$ in (12). Since the entries of the effective channel matrix $\mathbf{H}\in\mathbb{C}^{M_{\mathrm{d}}\times M_{\mathrm{d}}}$ are a function of $\bar{\mathbf{h}}$ (see (9)), we subsequently use compound function theory [39] to find the required CRLB. Let $\widehat{\mathbf{H}}$ denote the estimated effective DS-spread channel matrix. The mean square error (MSE) matrix, ${\mathbf{E}}_{\text{MSE}}$, is defined as

[TABLE]

In order to derive the result, we impose a prior $\bar{\mathbf{h}}\sim\mathcal{CN}(0,\mathbf{P}_{\bar{\mathbf{h}}}^{-1})$ on $\bar{\mathbf{h}}$ in (12), where $\mathbf{P}_{\bar{\mathbf{h}}}\in\mathbb{R}_{+}^{N_{\tau}M_{\alpha}\times N_{\tau}M_{\alpha}}$ is a diagonal deterministic precision matrix containing the hyperparameters $\boldsymbol{\theta}=[\theta_{0},...,\theta_{N_{\tau}M_{\alpha}-1}]^{T}$. Hence, $p(\bar{\mathbf{h}};\mathbf{P}_{\bar{\mathbf{h}}})=\frac{|\mathbf{P}_{\bar{\mathbf{h}}}|}{\pi^{N}}\text{exp}(-\bar{\mathbf{h}}^{H}\mathbf{P}_{\bar{\mathbf{h}}}\bar{\mathbf{h}}).$ Now, to transform (12) from the complex field to the real field, we define

[TABLE]

[TABLE]

$\mathbf{w}_{\mathrm{R}}=\begin{bmatrix}\Re(\mathbf{w}_{\mathrm{p}})\\ \Im(\mathbf{w}_{\mathrm{p}})\end{bmatrix}\in\mathbb{R}^{2M_{\mathrm{p}}},$ so that the new system model is

[TABLE]

where $\mathbf{h}_{\mathrm{R}}\sim\mathcal{N}(0,\mathbf{P}_{\mathbf{h}_{\mathrm{R}}}^{-1})$ with real-valued precision matrix \mathbf{P}_{\mathbf{h}_{\mathrm{R}}}=\text{diag}\Bigl{[}[2\boldsymbol{\theta}^{T},2\boldsymbol{\theta}^{T}]^{T}\Bigr{]}\in\mathbb{R}_{+}^{2N_{\tau}M_{\alpha}\times 2N_{\tau}M_{\alpha}} and $\mathbf{w}_{\mathrm{R}}\sim\mathcal{N}(0,\sigma_{\mathrm{R}}^{2})$ with $\sigma_{\mathrm{R}}^{2}=\sigma_{\mathrm{p}}^{2}/2$. We then have the following result.

Theorem 1**.**

The CRLB on the MSE matrix ${\mathbf{E}}_{\mathrm{MSE}}$ is given by

[TABLE]

where $\mathbf{\Phi}_{\bar{\mathbf{h}}}=\frac{1}{4}([\mathbf{\Phi}_{\mathbf{h}_{\mathrm{R}}}]_{\Re\Re}+[\mathbf{\Phi}_{\mathbf{h}_{\mathrm{R}}}]_{\Im\Im})+\frac{j}{4}([\mathbf{\Phi}_{\mathbf{h}_{\mathrm{R}}}]_{\Re\Im}-[\mathbf{\Phi}_{\mathbf{h}_{\mathrm{R}}}]_{\Im\Re})\in\mathbb{C}^{{N_{\tau}M_{\alpha}}\times{N_{\tau}M_{\alpha}}}$ is the complex valued BIM corresponding to $\bar{\mathbf{h}}$, $\mathbf{\Phi}_{\mathbf{h}_{\mathrm{R}}}\!=\!\frac{\mathbf{A}_{\mathrm{R}}^{T}\mathbf{A}_{\mathrm{R}}}{\sigma_{\mathrm{R}}^{2}}+\mathbf{P}_{\mathbf{h}_{\mathrm{R}}}\in\mathbb{R}^{2N_{\tau}M_{\alpha}\times 2N_{\tau}M_{\alpha}}$, the matrices $[\mathbf{\Phi}_{\mathbf{h}_{\mathrm{R}}}]_{\Re\Re},[\mathbf{\Phi}_{\mathbf{h}_{\mathrm{R}}}]_{\Im\Im},[\mathbf{\Phi}_{\mathbf{h}_{\mathrm{R}}}]_{\Re\Im}$ and $[\mathbf{\Phi}_{\mathbf{h}_{\mathrm{R}}}]_{\Im\Re}$ are the top-left, bottom-right, bottom-left and top-right $N_{\tau}M_{\alpha}\times N_{\tau}M_{\alpha}$-sized submatrices of $\mathbf{\Phi}_{\mathbf{h}_{\mathrm{R}}}$, respectively. Also, $\mathbf{U}=[\mathbf{u}_{0},...,\mathbf{u}_{N_{\tau}M_{\alpha}-1}]\in\mathbb{C}^{M_{\mathrm{d}}^{2}\times N_{\tau}M_{\alpha}}$ with $\mathbf{u}_{l}=\mathrm{vec}(\mathbf{G}^{H}\sqrt{\bar{\alpha}_{l}}\mathbf{F}_{M_{\mathrm{d}}}^{H}\mathbf{\Gamma}_{l}\mathbf{F}_{M_{\mathrm{d}}}^{\cdot\frac{1}{\bar{\alpha}_{l}}}\mathbf{G})\in\mathbb{C}^{M_{\mathrm{d}}^{2}}$.

Proof.

See Appendix.

∎

Remark 1**.**

Since $M_{\mathrm{d}}$ data symbols are transmitted, the channel $\mathbf{H}$ is an $M_{\mathrm{d}}\times M_{\mathrm{d}}$ matrix, i.e., it consists of a total of $M_{\mathrm{d}}^{2}$ parameters. Hence, the CRLB corresponding to $\mathbf{H}$ is an $M_{\mathrm{d}}^{2}\times M_{\mathrm{d}}^{2}$ matrix, with its trace providing a lower bound on the MSE in the channel estimates for each channel coefficient. Also, note that this is the CRLB for all the waveforms due to the general framework under which it has been derived.

VI Data Detection

We use the above VB-based off-grid CE technique to estimate the DS-spread channel parameters from the preamble measurements given by (12). Then, we estimate the effective DS-spread channel matrix $\widehat{\mathbf{H}}$ according to (9) using the estimated parameters, and use it to detect the data symbols. In this section, we develop a novel VB-based technique for estimating the soft symbols and marginals of data bits from (8). Since $\widehat{\mathbf{H}}$ is not a diagonal matrix, exact computation of the posterior distribution of the data symbols involves the computationally expensive calculation of the so-called partition function $p(\mathbf{y})=\sum_{\mathbf{x}\in\mathbb{Q}^{M_{\mathrm{d}}}}p(\mathbf{y}|\mathbf{x})p(\mathbf{x})$, which requires summing $Q^{M_{\mathrm{d}}}$ exponential terms. This motivates the need to develop a low-complexity algorithm for soft-symbol detection. Let the data symbols be i.i.d. with a uniform prior $\mathrm{P}_{x_{i}}(x_{i})=\frac{1}{Q}$, so that their joint distribution is $\mathrm{P}_{\mathbf{x}}(\mathbf{x})=\prod_{i=1}^{M_{\mathrm{d}}}\mathrm{P}_{x_{i}}(x_{i})$. Instead of computing the exact posterior, we approximate it by minimizing the $\text{KL}(q(\mathbf{x})||p(\mathbf{x}|\mathbf{y}))$, where $q(\mathbf{x})$ is an arbitrary distribution. By imposing a factorized distribution on $q$, the optimal marginal distribution for the $i$th symbol $(i=1,2,\ldots,M_{\mathrm{d}})$ in $\mathbf{x}$ can be derived, similar to (20), as $\ln q^{*}(x_{i})\propto\mathbb{E}_{u\neq i}[\ln\mathrm{p}(\mathbf{x},\mathbf{y})].$ The approximate posterior distribution is obtained through an iterative update of the marginals using the minorization-maximization principle, ensuring convergence to a stationary point. Thus, the $(j+1)^{\text{th}}$ update for the marginal of $\mathbf{x}$ is given by

[TABLE]

where the exponential and division are element-wise, and $\mathbf{g}^{(j+1)}(z)\triangleq-\frac{1}{\sigma_{\mathrm{d}}^{2}}(\text{diag}[\mathbf{\widehat{H}}^{H}\mathbf{\widehat{H}}]|z|^{2}-2\Re\{[\mathbf{\widehat{H}}^{H}\mathbf{y}-(\mathbf{\widehat{H}}^{H}\mathbf{\widehat{H}}-\text{diag}[\text{diag}[\mathbf{\widehat{H}}^{H}\mathbf{\widehat{H}}]])\langle\mathbf{x}\rangle^{(j)}]z^{*}\}),$ with $\langle\mathbf{x}\rangle^{(j)}$ being the $j$th update of the mean of $\mathbf{x}$. The pseudocode for the proposed VSSD algorithm is presented in Algorithm 2.

Remark 2**.**

Algorithm 2 provides the approximate marginals $\mathbf{q}(z)$ and soft symbols $\langle\mathbf{x}\rangle$ as outputs. In uncoded communications, the transmitted data symbols can be decoded by slicing the soft symbols. For coded communications, the LLRs of the data bits are computed from the marginals and passed as the inputs to the channel decoder [40]. Also, Algorithm 2 avoids matrix inversions, which enhances its computational efficiency. Specifically, while the MMSE equalizer has a computational complexity of $\mathcal{O}(N^{3})$, the VSSD algorithm has a per-iteration complexity of $\mathcal{O}(N^{2})$.

VII Iterative Channel Estimation and Data Detection (ICED)

In this section, we introduce a data-aided ICED technique to improve channel estimation and data detection. Initially, the unknown data symbols are detected using the pilot-only based estimated off-grid DS-spread channel. Subsequently, these detected data symbols $\widehat{\mathbf{x}}$ are used as virtual pilots to formulate a new measurement model, similar to (10), as

[TABLE]

where the dictionary matrix $\mathbf{A}_{\mathrm{d}}(\boldsymbol{\tau},\boldsymbol{\alpha})\in\mathbb{C}^{M_{\mathrm{d}}\times P}$ is calculated using the virtual pilots as $\mathbf{A}_{\mathrm{d}}(\boldsymbol{\tau},\boldsymbol{\alpha})=[\mathbf{a}_{\mathrm{d}}(\tau_{1},\alpha_{1}),\ldots,\mathbf{a}_{\mathrm{d}}(\tau_{P},\alpha_{P})],$ with $\mathbf{a}_{\mathrm{d}}(\tau_{p},\alpha_{p})=\sqrt{\alpha_{p}}\mathbf{G}^{H}\mathbf{F}_{M_{\mathrm{d}}}^{H}\mathbf{\Gamma}_{p}\mathbf{F}_{M_{\mathrm{d}}}^{\cdot\frac{1}{\alpha_{p}}}\mathbf{G}\widehat{\mathbf{x}}\in\mathbb{C}^{M_{\mathrm{d}}}$. Next, we concatenate the two measurement models, (10) and (56), to form an extended measurement model as

[TABLE]

where $\mathbf{y}_{\mathrm{E}}=[\mathbf{y}_{\mathrm{p}}^{T},\mathbf{y}^{T}]^{T}\in\mathbb{C}^{(M_{\mathrm{p}}+M_{\mathrm{d}})}$, $\mathbf{A}_{\mathrm{E}}=[\mathbf{A}_{\mathrm{p}}^{T},\mathbf{A}_{\mathrm{d}}^{T}]^{T}\in\mathbb{C}^{(M_{\mathrm{p}}+M_{\mathrm{d}})\times P}$, and $\mathbf{w}_{\mathrm{E}}=[\mathbf{w}_{\mathrm{p}}^{T},\mathbf{w}^{T}]^{T}\in\mathbb{C}^{(M_{\mathrm{p}}+M_{\mathrm{d}})}$. Now, we use the extended data-aided measurement model (57) to reestimate the DS-spread off-grid channel parameters. We alternate between CE and data detection until convergence (i.e., the channel and data symbol estimates do not change in successive iterations) or a maximum number of iterations have elapsed. Within the ICED framework, we employ either FVB or SVB for CE and MMSE or VSSD equalizers for data detection, ensuring fair comparison across approaches.

VIII Simulation Results

In this section, we present simulation results to evaluate the performance of the proposed off-grid DS-spread CE methods.444Matlab code is available at https://shorturl.at/JJvQi. We also empirically compare the performance of different waveforms under PCSIR and ECSIR, in the context of UWA communications. We consider a DS-spread channel with parameters $\alpha_{\mathrm{max}}=1.001$, $\tau_{\max}=32$ ms, and $P=5$ paths. These parameters are typical for UWA channels such as rivers, estuaries, and harbors [41, 42, 43, 44]. We simulate the UWA DS-spread channel with the $p^{\text{th}}$ scattering path parameter: $h_{p}\overset{i.i.d.}{\sim}\mathcal{CN}(0,1),$ $\tau_{p}\overset{i.i.d.}{\sim}\mathcal{U}(0,\tau_{\max})$, and $\alpha_{p}\overset{i.i.d.}{\sim}\mathcal{U}(\alpha_{\max}^{-1},\alpha_{\max})$. Further, the preamble of the transmitted frame contains $M_{\mathrm{p}}=32$ pilot symbols of duration $T_{\mathrm{p}}=3.2$ ms. A guard interval of $T_{g}=10$ ms is used to avoid excessive interference from the preamble to the data symbols.

The virtual grid numbers along the delay axis and the scale axis are set to $N_{\tau}=50$ and $M_{\alpha}=5$, respectively, with $r_{\tau}=0.64$ ms and $q_{\alpha}=1.0005$. The delay grid points are uniformly spaced in $[0,\tau_{\max}]$, while the scale grid points are uniformly spaced in the log domain, i.e., over the interval $[-\log\alpha_{\max},\log\alpha_{\max}]$. This results in a preamble-based measurement matrix $\mathbf{A}_{\mathrm{p}}^{(0)}\in\mathbb{C}^{32\times 250}$. Moreover, to examine the effect of CE on the data payload, we consider $M_{\mathrm{d}}=128$ data symbols arranged in the 2D domain by choosing $N=2$ and $M=64$ for all the waveforms. The data symbols are chosen from the BPSK constellation. We consider UWA communication over the frequency band spanning $f_{L}=10$ kHz to $f_{H}=20$ kHz, i.e., the signal bandwidth is $B=10$ kHz. Hence, the subcarrier width for OFDM, OTFS, and OCDM becomes $\Delta f=156.25$ Hz resulting in $T_{s}=\frac{N}{\Delta f}=12.8$ ms. For ODSS, using the constant $q=1.001$, the base subcarrier width $W=146.61$ Hz satisfies $B=\sum_{m=0}^{M-1}q^{m}W$, and $T_{s}=\frac{N}{W}=13.6$ ms. Moreover, we consider a CP duration of $T_{\mathrm{CP}}=10$ ms, hence the total data duration becomes $T_{\mathrm{d}}=T_{s}+NT_{\mathrm{CP}}=32.8$ ms for OFDM, OTFS, and OCDM, and $33.6$ ms for ODSS.

We first present the NMSE of the estimated effective DS-spread channel matrix, defined as $\text{NMSE}=\frac{\|\mathbf{H}-\widehat{\mathbf{H}}\|_{\mathrm{F}}^{2}}{\|\mathbf{H}\|_{\mathrm{F}}^{2}}$, using ODSS modulated preamble under on-grid and off-grid channel parameters, respectively, in Fig. 3(a) and Fig. 3(b). The algorithms we consider are the state-of-the-art on-grid OMP and VB [1], off-grid NOMP [31], and the proposed off-grid FVB and SVB. The root parameters for the Algorithm 1 are set as $\epsilon_{1}=\epsilon_{2}=\epsilon_{3}=\epsilon_{4}=10^{-6}$, the convergence threshold $\epsilon=10^{-3}$, and the maximum number of iterations $J_{\max}=100$. To update the dictionary matrix, we threshold the sparse channel vector with $5\%$ level (i.e., we use $\widehat{P}=\lceil 0.05N_{\tau}M_{\alpha}\rceil=13>P=5$), with the thresholding being used only to identify the columns of the dictionary to be updated. In our experiments, we have seen that increasing the threshold level simply increases the time complexity without any appreciable improvement in the performance. The NMSE is averaged over $10^{4}$ Monte Carlo trials for each SNR value. We also include the CRLB for the estimation of the effective channel matrix under on-grid channel parameters using Theorem 1. Since the channel parameters take on-grid values, the off-grid and on-grid schemes exhibit similar performance in Fig. 3(a). Moreover, the VB-based CE schemes demonstrate monotonically decreasing NMSE with increasing SNR and eventually approaching the CRLB. In contrast, the OMP-based CE schemes floor at NMSE $=2.1\times 10^{-2}$. This result emphasizes the advantage of VB-based estimation, as it provides a close approximation to the posterior. However, in Fig. 3(b), where practical off-grid channel parameters are considered, the performance of on-grid CE algorithms degrade significantly, leading to a high NMSE floor. For instance, on-grid OMP floors at NMSE $=1.1\times 10^{-1}$, while on-grid VB floors at NMSE $=7.2\times 10^{-2}$. Further, off-grid NOMP also floors at NMSE $=9.3\times 10^{-2}$, indicating that on-grid VB outperforms off-grid NOMP and is a powerful approach for CE: even without considering the off-grid effects, it outperforms the OMP-based algorithms that estimate the off-grid parameter values. The proposed off-grid FVB further improves the performance, especially at high SNR and floors at NMSE $=4.1\times 10^{-2}$. Finally, the proposed off-grid SVB more effectively captures the impact of parameter perturbation in the dictionary, resulting in nearly an order of magnitude improvement in NMSE compared to the off-grid FVB. Unlike the on-grid channel parameters case, all the algorithms exhibit an error floor at high SNR; this is due to the residual error in approximating the off-grid components.

Next, we assess the NMSE of off-grid DS-spread CE for both pilot-only (using only the preamble) and data-aided ICED techniques, as shown in Fig. 4 and Fig. 5. The ICED process involves estimating the off-grid DS-spread channel using the preamble of $M_{\mathrm{p}}=32$ symbols, detecting $M_{\mathrm{d}}=128$ data symbols using VSSD equalizer, and iteratively refining the channel estimates through ICED. Both FVB or SVB can be used as channel estimators for the pilot-only and ICED techniques; resulting in four combinations: pilot-only FVB, ICED FVB, pilot-only SVB, and ICED SVB. In our experiments, we find that the ICED algorithm takes only $3$ iterations to converge, making its complexity comparable to the pilot-only CE. Fig. 4 highlights the comparison between pilot-only and ICED for the ODSS waveform. We see that both ICED FVB and ICED SVB outperform pilot-only CE significantly. For instance, ICED SVB achieves NMSE $=2\times 10^{-2}$ at SNR $=8$ dB and continues to decrease with increasing SNR, while pilot-only SVB floors, with NMSE $=2\times 10^{-2}$ at SNR $=20$ dB. In particular, the performance of the ICED FVB has dramatically improved and is comparable to ICED SVB. The enhanced performance of ICED FVB and ICED SVB is due to the effective data detection achieved by the VSSD equalizer, along with the use of the posterior channel distribution from the previous iteration as the prior for the next iteration.

For a fair comparison, in Fig. 5, we adopt the SVB-based channel estimator for all waveforms, since it delivers better NMSE for off-grid DS-spread channels than FVB while maintaining almost similar computational complexity. We see that, with a robust channel estimator like SVB, all four waveforms exhibit nearly the same NMSE regardless of whether the channel is estimated using pilot-only or ICED techniques.

We compare the BER of different waveforms for uncoded communications under PCSIR and ECSIR, as shown in Fig. 6 and Fig. 7, respectively. Fig. 6 presents the BER with a $1$-tap equalizer and a VSSD equalizer. With the $1$-tap equalizer, the BER of OTFS floors at $3\times 10^{-3}$, outperforming OFDM, which floors at a BER of $7\times 10^{-3}$. This is because OTFS mounts the symbols in the DD domain, and thereby gains a higher diversity. OCDM achieves a lower BER floor $(4\times 10^{-4})$, as its chirp-based subcarriers provide additional robustness. ODSS delivers the best performance with a BER floor of the order $10^{-4}$ beyond an SNR of $15$ dB, as it is designed to mitigate the time-scale effects in the DS-spread channels [8]. The VSSD equalizer, which iteratively operates on the entire channel matrix, achieves significantly lower BER than the $1$-tap equalizer. In particular, the VSSD equalizer reduces the SNR required to obtain a BER $=10^{-2}$ by at least $4$ dB. Under the VSSD equalizer, OTFS, OCDM, and ODSS exhibit similar performance; the gap between them vanishes. Moreover, they outperform OFDM by $11.5$ dB at a BER of $10^{-3}$.

Fig. 7 illustrates the impact of SVB CE with VSSD equalizer on the BER using pilot-only and ICED techniques. As observed in the NMSE results in Fig. 5, ICED improves the BER compared to the pilot-only technique. For instance, when using pilot-only estimated channels, an SNR of $=11.5$ dB is required to achieve a BER $=10^{-3}$ for ODSS. In contrast, ICED achieves the same BER at a lower SNR $=8.5$ dB, demonstrating the efficiency of the ICED scheme. Furthermore, comparing the results in Fig. 7 with those in Fig. 6, it is evident that CE only marginally degrades the performance. Specifically, at a BER $=10^{-3}$, ICED ODSS (under ECSIR) experiences a degradation of only $0.8$ dB compared to the VSSD ODSS (under PCSIR). Moreover, with the estimated channel, OTFS, OCDM, and ODSS maintain similar performance and consistently outperform OFDM.

We present the BER of different waveforms for coded communications under ECSIR with ICED SVB in Fig. 8(a) and Fig. 8(b). For channel coding, we adopt an LDPC code based on the 3GPP 5G NR standard [45]. Specifically, we utilize the parity-check matrix derived from LDPC base graph 1 with a lifting size of 8 and set index 0. This configuration produces code blocks containing 44 data bits and 136 encoded bits. Fig. 8(a) illustrates the BER for both uncoded and coded communications using MMSE equalizers, with MMSE soft symbols for channel decoding. The results demonstrate that channel coding enhances performance, providing at least $3$ dB improvement at BER $=10^{-3}$. The OTFS, OCDM, and ODSS waveforms show similar performance in both uncoded and coded cases. Although OFDM performs the worst among the waveforms in both uncoded and coded cases, the gap between OFDM and the other waveforms is considerably reduced in the coded case. For example, to achieve BER $=10^{-2}$, OTFS, OCDM, and ODSS require SNR $=8.9$ dB compared to $13.4$ dB for OFDM. In contrast, for BER $=10^{-4}$, OTFS, OCDM, and ODSS require an SNR of $12.5$ dB, which is only $1.4$ dB lower than OFDM. A similar behavior has been observed in [46], when comparing OFDM and OTFS in the case of DD-spread channels. In Fig. 8(b), we present results when the soft symbols from the VSSD algorithm (Algorithm 2) are used to detect the data symbols in the uncoded case, and the LLRs of the data bits, derived from the marginals of the VSSD algorithm, are passed to the channel decoder in the coded case. In the uncoded case, comparing results with Fig. 8(a), VSSD significantly outperforms the MMSE equalizer. For instance, in ODSS, VSSD achieves BER $=10^{-2}$ provides at $2.2$ dB lower SNR than the MMSE equalizer. Furthermore, the results confirm that channel coding significantly enhances the BER compared to the uncoded case, with VSSD LLRs offering a more pronounced improvement than MMSE soft symbols. Specifically, when using VSSD LLRs, channel coding requires $5.5$ dB lower SNR to achieve a BER of $10^{-2}$ than in the uncoded case, whereas the reduction is only $2$ dB with MMSE soft symbols. This highlights the importance of using LLRs obtained through the VSSD algorithm. Moreover, in the coded case, OFDM performs on par with the other waveforms, in sharp contrast to its relatively poor performance in the uncoded case. Also, the gap between OFDM and the other waveforms in the coded case with MMSE soft symbols completely disappears when VSSD LLRs are used. The overlapping of coded BER across all the waveforms can be attributed to the fact that the channel code effectively spreads each raw data bit in the TF domain, which gets fully exploited when VSSD LLRs are used as input to the channel decoder. Thus, when advanced channel estimation and data detection techniques such as VSSD and ICED are employed, all waveforms offer the same performance.

Finally, in Fig. 9, we present the BER of all the waveforms under ECSIR with SVB using a publicly available real-world measured channel, namely, the Watermark datasets [44]. We use NOF1 channel datasets measured in a shallow stretch of Oslofjorden between a stationary source and a stationary single-hydrophone receiver. The channel parameters correspond to the NOF1 channel described in [44], with a center frequency of $14$ kHz, a sounding duration of $32.9$ s, and $60$ soundings. The delay spread of the channel is $128$ ms, for which we insert a guard interval of $T_{g}=25$ ms between preamble and data as the power-delay profile decays by more than $20$ dB beyond $25$ ms. The data transmission schemes (modulation waveforms and data bits) and the channel coding scheme (for coded communication) are the same as those used in Fig. 8(b). For the uncoded case, we slice the soft symbols from the VSSD decoder, while for the coded case, we use LLRs from the VSSD output as input to the channel decoder. Similar to the results in Fig. 6 for the 1-tap equalizer, in the uncoded case, OFDM exhibits the worst performance, followed by OTFS, OCDM, and ODSS is the best. However, in coded case, all the waveforms exhibit similar performance, owing to the spreading of data bits over time and frequency, as also observed in Fig. 8(b).

IX Conclusions

We considered data transmission over UWA channels modeled using a DS-spread representation with off-grid channel parameters. We developed a unified, waveform-agnostic receiver processing framework for CE and data detection, enabling a fair performance comparison of different waveforms. Building on this framework, we proposed two VB-based off-grid CE algorithms, FVB and SVB, which achieve the CRLB for on-grid channel parameters while outperforming existing algorithms in off-grid CE, and with SVB outperforming FVB. We also introduced a data-aided ICED algorithm that leverages the detected data symbols as virtual pilots, enhancing CE and eventually improving data detection over the pilot-only technique. In addition, we designed a low-complexity VSSD algorithm for soft symbol detection and LLR computation, which consistently outperforms the MMSE equalizer in both uncoded and coded communications.

Regarding the relative behavior across different waveforms, all exhibit similar NMSE with SVB-based CE under both pilot-only and ICED techniques. In uncoded communications under PCSIR with a $1$-tap equalizer, ODSS achieves the lowest BER, followed by OCDM, OTFS, and OFDM. With the VSSD equalizer, ODSS, OCDM, and OTFS exhibit the same BER and outperform OFDM under both PCSIR and ECSIR. In coded communications, under ECSIR, the data detection (BER) performance gap between OFDM and the other waveforms reduces with MMSE soft symbols and disappears completely with VSSD LLRs. Thus, we conclude that if the receiver is complexity-constrained (or if latency constraints prevent channel coding over long time intervals), and the channels are DS-spread, waveform selection plays a critical role: ODSS offers the best performance, followed by OCDM and OTFS, with OFDM performing the worst. However, if more sophisticated receiver processing is feasible, or if the channels are benign, all waveforms perform equally well.

As future work, we can reduce the complexity of the off-grid UWA CE technique by further decoupling the delay and scale parameters [47, 48]. Further, we can study the interplay between the grid size and the efficacy of the off-grid refinement, extend the results to multiple antenna systems, and theoretically analyze the performance of different waveforms in coded communications.

Appendix A Proof of the Theorem 1

The BIM [37] of $\mathbf{h}_{\mathrm{R}}$ is calculated using (53) as

[TABLE]

By chain rule  [38], the complex-valued BIM can be found as $\mathbf{\Phi}_{\bar{\mathbf{h}}}=\frac{1}{4}([\mathbf{\Phi}_{\mathbf{h}_{\mathrm{R}}}]_{\Re\Re}+[\mathbf{\Phi}_{\mathbf{h}_{\mathrm{R}}}]_{\Im\Im})+\frac{j}{4}([\mathbf{\Phi}_{\mathbf{h}_{\mathrm{R}}}]_{\Re\Im}-[\mathbf{\Phi}_{\mathbf{h}_{\mathrm{R}}}]_{\Im\Re})$. We can write the vectorized version of $\mathbf{H}$ from (9) as $\text{vec}(\mathbf{H})=\mathbf{U}\bar{\mathbf{h}},$ where the $l$th column of $\mathbf{U}$ is $\mathbf{u}_{l}=\text{vec}(\mathbf{G}^{H}\sqrt{\bar{\alpha}_{l}}\mathbf{F}_{M_{\mathrm{d}}}^{H}\mathbf{\Gamma}_{l}\mathbf{F}_{M_{\mathrm{d}}}^{\cdot\frac{1}{\bar{\alpha}_{l}}}\mathbf{G})$. Finally, we use the transformation property [39] to obtain the CRLB of $\mathbf{H}$ as

[TABLE]

which completes the proof.