Massive Random Access in Cell-Free Massive MIMO Systems for High-Speed Mobility with OTFS Modulation

Abstract

In the research of next-generation wireless communication technologies, orthogonal time frequency space (OTFS) modulation is emerging as a promising technique for high-speed mobile environments due to its superior efficiency and robustness in doubly selective channels. Additionally, the cell-free architecture, which eliminates the issues associated with cell boundaries, offers broader coverage for radio access networks. By combining cell-free network architecture with OTFS modulation, the system may meet the demands of massive random access required by machine-type communication devices in high-speed scenarios. This paper explores a massive random access scheme based on OTFS modulation within a cell-free architecture. A transceiver model for uplink OTFS signals involving multiple access points (APs) is developed, where channel estimation with fractional channel parameters is approximated as a block sparse matrix recovery problem. Building on existing superimposed and embedded preamble schemes, a hybrid preamble scheme is proposed. This scheme leverages superimposed and embedded preambles to respectively achieve rough and accurate active user equipment (UEs) detection (AUD), as well as precise channel estimation, under the condition of supporting a large number of access UEs. Moreover, this study introduces a generalized approximate message passing and pattern coupling sparse Bayesian learning with Laplacian prior (GAMP-PCSBL-La) algorithm, which effectively captures block sparse features after discrete cosine transform (DCT), delivering precise estimation results with reduced computational complexity. Simulation results demonstrate that the proposed scheme is effective and provides superior performance compared to other existing schemes.

I Introduction

The next-generation wireless communication will delve deeper into more ubiquitous Internet of Thing (IoT) scenarios in the coming decades, encompassing broader coverage areas and a significantly larger number of user equipment (UEs) [1]. Beyond human-type communication devices (HTCDs), there are numerous machine-type communication devices (MTCDs) that need to connect to wireless air interface to facilitate data transmission [2]. In high-speed massive machine-type communication (mMTC) scenarios, such as high-speed railways, Internet of Vehicles (IoV), unmanned aerial vehicle (UAV) communications, and high-speed integrated sensing and communication (ISAC) [3], the great number of MTCDs faces constraints on allocable resources and quality of service (QoS) [4]. Traditional coordinated access protocols, which require multiple handshake processes, not only cause delays but also generate substantial signaling overhead [5]. Moreover, coordinated orthogonal resources suffer severe orthogonality degradation in doubly-selective channels, thereby reducing system performance [6]. Unlike coordinated schemes, grant-free NOMA allows devices to transmit data without allocated resources. The receiver performs active UE detection (AUD) and channel estimation (CE) based on unique non-orthogonal preamble sequence assigned to each UE [7]. Therefore, grant-free NOMA in uncoordinated access schemes is considered one of the key technologies for mMTC [4].

Emerging machine-type wireless transmission services impose stringent demands on communication quality in high-mobility scenarios, such as intelligent driving, virtual reality, and UAV reconnaissance. Under this consideration, robust transceiver design becomes a critical challenge. Orthogonal frequency division multiplexing (OFDM), widely used in 4G and 5G, can eliminate inter-symbol interference caused by time dispersion using a cyclic prefix (CP), but struggles to mitigate frequency dispersion caused by Doppler shifts, leading to inter-carrier interference. Additionally, the length of the CP in OFDM is proportional to the number of subcarriers, impacting efficiency when the number of subcarriers is large [8]. To overcome these limitations, Hadani et al. proposed a novel two-dimensional modulation known as orthogonal time frequency space (OTFS) [9]. Compared to OFDM, OTFS has been shown to significantly improve transmission performance in doubly-selective channels with only a modest increase in system complexity [10]. Specifically, OTFS uses a two-dimensional inverse symplectic finite Fourier transform (ISFFT) to map signals from the Doppler-delay (DD) domain to the time-frequency (TF) domain. Unlike OFDM, each signal symbol in OTFS spans the entire TF domain channel, fully exploiting channel diversity and enhancing reliability [11]. Additionally, the number of reflectors is considerably smaller than the dimension of transmitted symbols, resulting in sparsity for channel parameters in the DD domain [11], which simplify the estimation of channel state information (CSI). Given these advantages, OTFS is considered a promising candidate for next-generation broadband communication modulation technology.

In addition, high-mobility communication inevitably requires wide coverage, as UEs may move considerable distances during communication intervals. Traditional cellular network models necessitate handovers for high-mobility UEs, which increases system processing complexity [12]. Moreover, boundary effects limit the QoS for UEs located at the cell edges [13]. Recently, a concept named cell-free massive MIMO has been proposed to support denser and wider network device coverage, significantly enhancing spectral efficiency and reliability [14]. Cell-free massive MIMO eliminates boundary effects by deploying numerous access points (APs) across the coverage area [15]. Each AP is equipped with an independent signal processing unit and connected to a central processing unit (CPU) via fronthaul links, providing a flexible networking [16]. Additionally, with UEs being closer to the receiving antennas, signal transmission and processing delays are significantly reduced. Mohammadi et al. theoretically demonstrated that OTFS modulation can achieve superior performance within a cell-free massive MIMO architecture [17]. However, there remain numerous challenges to be addressed for massive random access in high-mobility scenarios when integrating cell-free massive MIMO.

Current discussions on OTFS grant-free access schemes for high-mobility mainly focus on low Earth orbit (LEO) satellite communication [18]. Shen et al. approximated the OTFS channel as a sparse matrix and utilized the low-complexity pattern-coupled sparse Bayesian learning (PCSBL) for AUD and sparse CE [19]. Zhou et al. designed a novel training sequences aided OTFS (TS-OTFS) transmission protocol for LEO satellite IoT communication and proposed a two-stage AUD and CE method [20]. Besides, a high-speed railway IoT active detection method combining tandem spreading multiple access (TSMA) and OTFS was proposed in [21]. By pre-estimating propagation delays, a preamble transmission method was designed in [22], allowing UEs to perform pre-compensation. Considering the 3D channel sparsity characteristic of OTFS-massive MIMO, Shen et al. proposed a 3D-structured orthogonal matching pursuit (OMP) algorithm for CE [23]. However, existing researches lack the design of schemes for massive high-mobility MTCD access incorporating cell-free massive MIMO system. Moreover, the current CE methods, including the embedded [24] and superimposed [25] pilot schemes, have their limitations respectively: the former incurs high pilot overhead, while the latter has suboptimal estimation performance. There is a need for a balanced scheme that achieves good estimation performance while reducing overhead.

To address the aforementioned challenges, this paper investigates the AUD and CE schemes for massive random access in cell-free massive MIMO system with high-speed mobility scenarios. Firstly, we establish the OTFS uplink signal model in cell-free massive MIMO system and extend it to multi-UE scenarios. Secondly, we design a hybrid preamble scheme, where rough AUD is performed by superimposed preambles, and accurate AUD and CE are achieved based on embedded preambles. This scheme reduces the overall sparse signal dimension, allowing the system to accommodate more UEs. Finally, we propose a new block sparse matrix recovery algorithm for AUD and CE, named generalized approximate message passing and pattern coupling sparse Bayesian learning with Laplacian prior (GAMP-PCSBL-La). Numerical simulations demonstrate that this algorithm achieves more accurate estimation performance. In summary, our main contributions are as follows:

• We first analyze the OTFS SISO delay-Doppler domain transceiver signal model in this paper. Through mathematical approximation and selecting appropriate preamble sequence embedding positions, we model the preamble transceiver signals as a two-dimensional sparse compressed sensing model in the DD domains. We then extend this model to multi-UE and uniform planar antenna array, forming a four-dimensional sparse channel model in the delay-Doppler-UE-angle domains. By merging the delay, Doppler, and UE domains, we ultimately transform the AUD and CE problem into a two-dimensional block sparse matrix recovery problem.

• To address the scale constraints of high-dimensional sparse matrices in compressed sensing while minimizing the overhead of preambles, we propose a hybrid preamble scheme. Rough AUD is performed by superimposed preamble, followed by accurate AUD and CE based on embedded preamble. This approach reduces the sparse channel dimension for each estimation, enabling the system to support numerous UEs access.

• A novel GAMP-PCSBL-La algorithm is designed to recover the two-dimensional block sparse channel matrix. GAMP achieves good estimation performance while reducing computational complexity by avoiding matrix inversion [27]. PCSBL captures the block sparsity of two-dimensional matrix [28], and the Laplacian prior distribution has been proven to enhance reconstruction performance of sparse signals with discrete cosine transform (DCT) [29]. By combining these features, GAMP-PCSBL-La achieves excellent channel estimation accuracy with low computational complexity. Our simulation results further validate this conclusion.

The remainder of this paper is organized as follows. In Section II, we introduce the system model studied in this paper. Section III discusses the rough AUD and accurate AUD and CE strategies based on the hybrid preamble scheme. In Section IV, we present the novel block sparse matrix recovery algorithm, GAMP-PCSBL-La. Section V provides numerical simulations and analyzes the corresponding performance. Finally, we conclude the paper in the last section.

Notations: Bold lower letters and bold capital letters denote vectors and matrices, respectively. Normal lower letters and capital letters represent scalar variables and constants, respectively. $\mathbb{C}$ and $\mathbb{R}$ are complex number set and real number set, respectively. $\mathbf{A}_{a:b,c:d}$ represent a sliced matrix for $\mathbf{A}$ with $a$-th row to $b$-th row and $c$-th column to $d$-th column, while $\mathbf{a}_{a:b}$ is a sliced vector for $\mathbf{a}$ with $a$-th element to $b$-th element. Especially, $\mathbf{A}_{a:b,:}$ denotes the submatrix of $\mathbf{A}$ with $a$-th row to $b$-th row. $\mathbb{E}$ and $\mathbb{V}$ mean the expectation and variance, respectively. $\delta(\cdot)$ denotes a Dirac delta function, and $()^{H}$ is the conjugate transpose of a matrix or vector. $\mathbf{A}[a,b]$ means $(a,b)$-th element of matrix $\mathbf{A}$. $\otimes$ and $\odot$ represent Kronecker product and Hadamard product, respectively. Calligraphy letters are used to denote sets.

II System Model

II-A Cell-free Massive MIMO System

We consider a cell-free massive MIMO system, as shown in Fig. 1 below, comprising a total of $B$ APs and $U$ single-antenna UEs, randomly distributed over a large area. Each AP is equipped with a uniform planar antenna array with dimensions $N_{y}\times N_{z}$. Assume that each AP is connected to the CPU through fronthaul, allowing lossless data interaction. The UEs move within the area, and only a small portion of UEs transmit uplink data to APs in a specific transmission slot; these UEs are referred to as active UEs, while the remaining are silent. Active UEs form a set denoted as $\mathcal{K_{A}}$. The channels of active UEs and APs experience doubly selective fading. It is assumed that signals exceeding a maximum delay ${\tau_{\max}}$ at the receiver end are considered as noise. Furthermore, the range of UE velocity is $\left[{-{v_{\max}},{v_{\max}}}\right]$. For simplicity, we set the velocity zero point at $-{v_{\max}}$, which transforms the velocity range of UE to $\left[{0,2{v_{\max}}}\right]$. The signal propagation from an active UE to an AP is characterized by a finite number of paths. The uplink signal transmitted by UE consists of preamble sequences and data symbols. Each AP detects the activity of UE and performs CE based on received preamble sequences.

Figure 1: Massive random access in cell-free massive MIMO system.

II-B OTFS Modulation and Channel Model

Consider a typical OTFS transceiver system with $M$ subcarriers, subcarrier spacing of $\Delta f$, and symbol duration $T=1/\Delta f$ with $N$ symbols. Therefore, the system’s bandwidth is $M\Delta f$ and the duration of a block is $NT$. We assume the channel parameters remain constant within a transmission block. In DD domain, the resolutions of delay and Doppler parameters are $\frac{1}{{M\Delta f}}$ and $\frac{1}{{NT}}$, respectively. For a given active UE $u$, the modulated and power-allocated symbol $\left\{{X_{u}^{DD}\left[{k,l}\right],0\leq k\leq N-1,0\leq l\leq M-1}\right\}$ is assigned to the $(k,l)$-th grid point in $N\times M$ DD grid. Here, $k$ and $l$ represent the indices for Doppler domain and delay domain, respectively. By applying the ISFFT to ${\mathbf{X}}_{u}^{DD}\in{\mathbb{C}^{N\times M}}$ in DD domain, the $N\times M$ zero-mean symbols are transformed into TF domain:

$$\mathbf{X}_{u}^{TF}\left[{n,m}\right]=\frac{1}{{\sqrt{NM}}}\sum\limits_{k}{}\sum\limits_{l}{}\mathbf{X}_{u}^{DD}\left[{k,l}\right]{e^{-j2\pi\left({\frac{{ml}}{M}-\frac{{nk}}{N}}\right)}},$$ (1)

where $0\leq n\leq N-1,0\leq m\leq M-1$. On this basis, the transmitter applies the Heisenberg transform to the TF domain signal matrix ${\mathbf{X}}_{u}^{TF}\in{\mathbb{C}^{N\times M}}$, converting it into a time-domain signal, represented as:

$${s_{u}}\left(t\right)=\sum\limits_{n}{}\sum\limits_{m}{}\mathbf{X}_{u}^{TF}\left[{n,m}\right]{e^{j2\pi m\Delta f(t-nT)}}{g_{tx}}(t-nT).$$ (2)

Where ${g_{tx}}(t)$ represents the rectangular window function of the transmitted signal in time domain with a duration of $T$, defined as:

$${g_{tx}}\left(t\right)=\left\{{\begin{array}[]{*{20}{c}}{\frac{1}{{\sqrt{T}}}}&{0\leq t\leq T}\\ 0&{{\text{otherwise}}}\end{array}}\right..$$ (3)

To facilitate the analysis of the received signal, we assume that each AP has a single antenna and there is only one active UE. In fact, the multi-antenna and multi-UE scenarios can be easily extended from this results, as shown in the subsequent derivations. The delay-Doppler channel response model from UE $u$ to the $b$-th AP is defined as:

$${h_{u,b}}\left({\tau,\nu}\right)=\sum\limits_{i=1}^{P}{{h_{u,b,i}}\delta\left({\tau-{\tau_{u,b,i}}}\right)\delta\left({\nu-{\nu_{u,b,i}}}\right)}.$$ (4)

Here, ${h_{u,b,i}}$, ${\tau_{u,b,i}}$ and ${\nu_{u,b,i}}$ represent the gain, delay, and Doppler shift, respectively, of the $i$-th path from UE $u$ to the $b$-th AP. $P$ is the number of path. Consider the path loss and shadow fading, we have ${h_{u,b,i}}\sim{\mathcal{C}\mathcal{N}}\left({0,{\lambda_{u,b}}}\right)$, with ${\lambda_{u,b}}$ representing the large-scale fading coefficient of the channels from UE $u$ to $b$-th AP. The corresponding received time-domain signals for $b$-th AP is presented as:

$${r_{b}}\left(t\right)=\iint{}{h_{u,b}}\left({\tau,\nu}\right){s_{u}}\left({t-\tau}\right){e^{j2\pi\left({t-\tau}\right)\nu}}d\tau d\nu+{n_{b}}\left(t\right)$$ (5)

Where ${n_{b}}\left(t\right)$ represents the additive Gaussian noise signal in the time domain, with a mean of 0. The local signal processing unit of AP performs Wigner transform on the time-domain received signal, yielding the received signal in TF domain presented as follows:

$$\begin{gathered}{\mathbf{Y}_{b}}\left[{n,m}\right]=\int{}{r_{b}}\left(t\right)g_{rx}^{*}\left({t-nT}\right){e^{-j2\pi m\Delta f(t-nT)}}dt\hfill\\ =\int{}\sum\limits_{i}{}{h_{u,b,i}}{s_{u}}\left({t-{\tau_{u,b,i}}}\right){e^{j2\pi\left({t-{\tau_{u,b,i}}}\right){\nu_{u,b,i}}}}g_{rx}^{*}\left({t-nT}\right)\hfill\\ {e^{-j2\pi m\Delta f(t-nT)}}dt+{\mathbf{N}_{b}}\left[{n,m}\right]\hfill\\ =\int{}\sum\limits_{i}{}{h_{u,b,i}}\sum\limits_{n^{\prime}}{}\sum\limits_{m^{\prime}}{}\mathbf{X}_{u}^{TF}\left[{n^{\prime},m^{\prime}}\right]{e^{j2\pi m^{\prime}\Delta f(t-{\tau_{u,b,i}}-n^{\prime}T)}}\hfill\\ {g_{tx}}(t-{\tau_{u,b,i}}-n^{\prime}T){e^{j2\pi\left({t-{\tau_{u,b,i}}}\right){\nu_{u,b,i}}}}g_{rx}^{*}\left({t-nT}\right)\hfill\\ {e^{-j2\pi m\Delta f(t-nT)}}dt+{\mathbf{N}_{b}}\left[{n,m}\right].\hfill\\ \end{gathered}$$ (6)

Where ${g_{rx}}(t)$ denotes the rectangular window function in receiver, which is defined identically to equation (3). ${\mathbf{N}_{b}}$ is the noise in TF domain. According to the properties of the window function, in equation (6), ${g_{tx}}(t-{\tau_{u,b,i}}-n^{\prime}T)$ and $g_{rx}^{*}\left({t-nT}\right)$ can take nonzero values when $0\leq t-{\tau_{u,b,i}}-n^{\prime}T\leq T$ and $0\leq t-nT\leq T$. When $n=n^{\prime}$, the integration range in equation (6) is $\left[{nT+{\tau_{u,b,i}},(n+1)T}\right]$, and when $n=n^{\prime}+1$, it becomes $\left[{nT,nT+{\tau_{u,b,i}}}\right]$. When $n$ takes other values, the integral result in equation (6) is 0. Therefore, equation (6) can be rewritten as:

$$\begin{gathered}{\mathbf{Y}_{b}}\left[{n,m}\right]=\frac{1}{T}\sum\limits_{i}{}{h_{u,b,i}}\sum\limits_{m^{\prime}}{}{\mathbf{X}^{TF}}\left[{n,m^{\prime}}\right]{e^{-j2\pi m^{\prime}\Delta f{\tau_{u,b,i}}}}\hfill\\ {e^{-j2\pi{\nu_{u,b,i}}{\tau_{u,b,i}}}}{e^{j2\pi{\nu_{u,b,i}}nT}}\int_{{\tau_{u,b,i}}}^{T}{}{e^{-j2\pi\Delta ft(m-m^{\prime}-\frac{{{\nu_{u,b,i}}}}{{\Delta f}})}}dt\hfill\\ +\frac{1}{T}\sum\limits_{i}{}{h_{u,b,i}}\sum\limits_{m^{\prime}}{}{\mathbf{X}^{TF}}\left[{n-1,m^{\prime}}\right]{e^{-j2\pi m^{\prime}\Delta f{\tau_{u,b,i}}}}{e^{j2\pi m^{\prime}\Delta fT}}\hfill\\ {e^{-j2\pi{\nu_{u,b,i}}{\tau_{u,b,i}}}}{e^{j2\pi{\nu_{u,b,i}}nT}}\int_{0}^{{\tau_{u,b,i}}}{}{e^{-j2\pi\Delta ft(m-m^{\prime}-\frac{{{\nu_{u,b,i}}}}{{\Delta f}})}}dt\hfill\\ +{\mathbf{N}_{b}}\left[{n,m}\right].\hfill\\ \end{gathered}$$ (7)

Usually, $M$ is greater than $N$. We assume that each delay parameter is an integer multiple of the resolution, i.e.,

$${\tau_{u,b,i}}=\frac{{{l_{u,b,i}}}}{{M\Delta f}},$$ (8)

$${\nu_{u,b,i}}=\frac{{{k_{u,b,i}}+{{\tilde{k}}_{u,b,i}}}}{{NT}}.$$ (9)

Where both ${l_{u,b,i}}$ and ${k_{u,b,i}}$ are integers, and ${\tilde{k}_{u,b,i}}$ is a fraction value between -0.5 and 0.5. Using symplectic finite Fourier transform (SFFT), the received signal in the TF domain is transformed into the DD domain:

$$\mathbf{Y}_{b}^{DD}[k,l]=\frac{1}{{\sqrt{NM}}}\sum\limits_{n}{}\sum\limits_{m}{}\mathbf{Y}_{b}^{TF}\left[{n,m}\right]{e^{j2\pi\left({\frac{{ml}}{M}-\frac{{nk}}{N}}\right)}}.$$ (10)

Combining equations (1), (7), (8), (9) and (10), we obtain the received signal model in DD domain as equation (11), where $k^{\prime\prime}\in\left[{{k_{u,b,i}}-\varepsilon,{k_{u,b,i}}+\varepsilon}\right]$ is defined as the neighborhood of integer Doppler parameters, and $\varepsilon$ is a very small integer. The derivation of equation (11) can be found in Appendix A.

$$Y_{b}^{DD}\left[{k,l}\right]\approx\left\{{\begin{array}[]{*{20}{c}}{\sum\limits_{i}{}{h_{u,b,i}}\sum\limits_{k^{\prime\prime}}{}X_{u}^{DD}\left[{k-k^{\prime\prime},l-{l_{u,b,i}}}\right]{e^{-j2\pi\frac{{k-k^{\prime\prime}}}{N}}}{e^{j2\pi\frac{{\left({l-{l_{u,b,i}}}\right)\left({{k_{u,b,i}}+{{\tilde{k}}_{u,b,i}}}\right)}}{{NM}}}}\frac{1}{N}\frac{{1-{e^{j2\pi{{\tilde{k}}_{u,b,i}}}}}}{{1-{e^{-j2\pi\frac{{k^{\prime\prime}-{k_{u,b,i}}-{{\tilde{k}}_{u,b,i}}}}{N}}}}}}&{l<{l_{u,b,i}}}\\ {\sum\limits_{i}{}{h_{u,b,i}}\sum\limits_{k^{\prime\prime}}{}X_{u}^{DD}\left[{k-k^{\prime\prime},l-{l_{u,b,i}}}\right]{e^{j2\pi\frac{{\left({l-{l_{u,b,i}}}\right)\left({{k_{u,b,i}}+{{\tilde{k}}_{u,b,i}}}\right)}}{{NM}}}}\frac{1}{N}\frac{{1-{e^{j2\pi{{\tilde{k}}_{u,b,i}}}}}}{{1-{e^{-j2\pi\frac{{k^{\prime\prime}-{k_{u,b,i}}-{{\tilde{k}}_{u,b,i}}}}{N}}}}}}&{l\geq{l_{u,b,i}}}\end{array}}\right.$$ (11)

III Hybrid Preamble-based AUD and CE Scheme

Inspired by the existing OTFS superimposed pilots scheme and embedded pilot-aided scheme, we propose a hybrid preamble strategy that achieves precise AUD and CE with reduced preamble overhead. Firstly, we perform sparse channel recovery based on the superimposed preamble (preamble1) received signals to obtain a rough set of active UEs. We then perform another sparse channel recovery on the embedded preamble (preamble2) received signals to achieve a more accurate result based on the rough set. Given the smaller dimension of preamble1, this scheme can support a larger number of user devices. Preamble2, with its higher resolution, enables more accurate estimation for the detected rough active UE set. Furthermore, by leveraging the characteristics of planar antenna arrays, we obtain a sparser block-sparse channel matrix model, yielding precise detection and estimation results. A more detailed description is provided in the subsequence of this section.

III-A Rough AUD

We are going to consider the signal model with both delay domain and Doppler domain dimensions are relatively small. Assuming $N^{\prime}=\alpha N$ and $M^{\prime}=\beta M$ are both integers, where $0<\alpha,\beta<1$. The quantization value corresponding to the maximum delay ${\tau_{\max}}$ is ${\tilde{l}^{\prime}_{\max}}={\tau_{\max}}M^{\prime}\Delta f=\beta{\tau_{\max}}M\Delta f\ll 1$ when $\beta$ is particularly small, which implies that any delay parameter $0<{\tilde{l}^{\prime}_{u,b,i}}\leq{\tilde{l}^{\prime}_{\max}}\ll 1$ is a fractional value. The quantization value corresponding to the maximum Doppler shift is defined as ${k^{\prime}_{\max}}=\left\lfloor{{\nu_{\max}}N^{\prime}T}\right\rfloor$. Similar to equation (11), we obtain the reception model as:

$$\begin{gathered}\mathbf{Y}_{b}^{DD}\left[{k^{\prime},l^{\prime}}\right]\mathop{\approx}\limits^{(a)}\frac{1}{{N^{\prime}M^{\prime}}}\sum\limits_{i}{}{h_{u,b,i}}\sum\limits_{k^{\prime\prime}}{}\mathbf{X}_{u}^{DD}\left[{k^{\prime}-k^{\prime\prime},l^{\prime}}\right]{e^{j2\pi\frac{{l^{\prime}{{k^{\prime}}_{u,b,i}}}}{{N^{\prime}M^{\prime}}}}}\hfill\\ {e^{-j2\pi\frac{{{{\tilde{l}^{\prime}}_{u,b,i}}\left({{{k^{\prime}}_{u,b,i}}+{{\tilde{k}^{\prime}}_{u,b,i}}}\right)}}{{N^{\prime}M^{\prime}}}}}\frac{{1-{e^{j2\pi{{\tilde{k}^{\prime}}_{u,b,i}}}}}}{{1-{e^{-j2\pi\frac{{k^{\prime\prime}-{{k^{\prime}}_{u,b,i}}-{{\tilde{k}^{\prime}}_{u,b,i}}}}{{N^{\prime}}}}}}}\frac{{1-{e^{-j2\pi{{\tilde{l}^{\prime}}_{u,b,i}}}}}}{{1-{e^{-j2\pi\frac{{{{\tilde{l}^{\prime}}_{u,b,i}}}}{{M^{\prime}}}}}}},\hfill\\ \end{gathered}$$ (12)

where ${\tilde{l}^{\prime}_{u,b,i}}={\tau_{u,b,i}}M^{\prime}\Delta f$, ${k^{\prime}_{u,b,i}}=\left[{{\nu_{u,b,i}}N^{\prime}T}\right]_{\text{R}}$, ${\tilde{k}^{\prime}_{u,b,i}}={\nu_{u,b,i}}N^{\prime}T-\left[{{\nu_{u,b,i}}N^{\prime}T}\right]$ and $k^{\prime\prime}\in\left[{{{k^{\prime}}_{u,b,i}}-\varepsilon^{\prime},{{k^{\prime}}_{u,b,i}}+\varepsilon^{\prime}}\right]$ is defined as the neighborhood of ${{k^{\prime}}_{u,b,i}}$ with $\varepsilon^{\prime}$ is a integer with very small value. Due to $l^{\prime}<M^{\prime}$, $-0.5\leq{\tilde{k}^{\prime}_{u,b,i}}<0.5$, and especially when $M^{\prime}$ is very small and $N^{\prime}$ is larger compared to $M^{\prime}$, equation (a) holds approximately true. The derivation process of the above equation can be found in Appendix B. Unlike equation (11), this reception model is not represented in segments. We define functions:

$$\phi^{\prime}\left({k,l}\right)={e^{j2\pi\frac{{kl}}{{N^{\prime}M^{\prime}}}}},$$ (13)

$$\varphi( {k,{k_{i}},{{\tilde{k}}_{i}},{{\tilde{l}}_{i}},N,M})=\frac{1}{{NM}}\frac{{1-{e^{j2\pi{{\tilde{k}}_{i}}}}}}{{1-{e^{-j2\pi\frac{{k-{k_{i}}-{{\tilde{k}}_{i}}}}{N}}}}}$$ (14)

$$\frac{{1-{e^{-j2\pi{{\tilde{l}}_{i}}}}}}{{1-{e^{-j2\pi\frac{{{{\tilde{l}}_{i}}}}{M}}}}}{e^{-j2\pi\frac{{{{\tilde{l}}_{i}}\left({{k_{i}}+{{\tilde{k}}_{i}}}\right)}}{{NM}}}},$$

$$\begin{gathered}C\left({{\mathbf{x}},k,\varepsilon}\right)=[{\text{circ(}}{\mathbf{x}},0{\text{),circ(}}{\mathbf{x}},1{\text{),}}\ldots{\text{,circ(}}{\mathbf{x}},k+\varepsilon{\text{),}}\hfill\\ {\text{circ(}}{\mathbf{x}},-\varepsilon{\text{),circ(}}{\mathbf{x}},-\varepsilon+1{\text{),}}\ldots{\text{,circ(}}{\mathbf{x}},-1{\text{)}}].\hfill\\ \end{gathered}$$ (15)

Where ${\text{circ(}}{\mathbf{x}},i{\text{)}}$ represents the vector obtained by circularly shifting vector $\mathbf{x}$ by $i$ positions. Based on the above definitions, we transform equation (12) into vector form:

$$\begin{gathered}{\mathbf{y}}_{b}^{p1}\approx\left({{\mathbf{X}}_{u}^{p1}\odot{\mathbf{\Phi^{\prime}}}}\right){\mathbf{h}}_{u,b}^{DD1}+{\mathbf{n}}_{b}^{DD1}\hfill\\ ={\mathbf{A}}_{u}^{p1}{\mathbf{h}}_{u,b}^{DD1}+{\mathbf{n}}_{b}^{DD1}\in{\mathbb{C}^{N^{\prime}M^{\prime}\times 1}},\hfill\\ \end{gathered}$$ (16)

where ${\mathbf{y}}_{b}^{p1}$ is the vectorization of $\mathbf{Y}_{b}^{DD}$, and

$${\mathbf{X}}_{u}^{p1}=\left[{\begin{array}[]{*{20}{c}}{C\left({{\mathbf{X^{\prime}}}_{u}^{DD}\left({:,1}\right),{{k^{\prime}}_{\max}},\varepsilon^{\prime}}\right)}\\ {C\left({{\mathbf{X^{\prime}}}_{u}^{DD}\left({:,2}\right),{{k^{\prime}}_{\max}},\varepsilon^{\prime}}\right)}\\ \vdots\\ {C\left({{\mathbf{X^{\prime}}}_{u}^{DD}\left({:,M^{\prime}}\right),{{k^{\prime}}_{\max}},\varepsilon^{\prime}}\right)}\end{array}}\right],$$ (17)

$${\mathbf{\Phi^{\prime}}}={\left[{{\mathbf{\Phi^{\prime}}}_{1}^{H},\ldots,{\mathbf{\Phi^{\prime}}}_{l}^{H},\ldots,{\mathbf{\Phi^{\prime}}}_{M^{\prime}}^{H}}\right]^{H}},$$ (18)

$${{\mathbf{\Phi^{\prime}}}_{l}}={\bm{\phi}^{\prime}_{l}}\otimes{{\mathbf{1}}^{N^{\prime}\times 1}},$$ (19)

$$\begin{gathered}{\bm{\phi}^{\prime}_{l}}=[\phi^{\prime}(l,0),\ldots,\phi^{\prime}(l,{k_{\max}}+\varepsilon),\hfill\\ \phi^{\prime}(l,N-\varepsilon-1),\ldots,\phi^{\prime}(l,N-1){]^{T}}.\hfill\\ \end{gathered}$$ (20)

${\mathbf{h}}_{u,b}^{DD1}\in{\mathbb{C}^{({{k^{\prime}}_{\max}}+2\varepsilon^{\prime}+1)\times 1}}$, it expressed as

$${\mathbf{h}}_{u,b}^{DD1}=\sum\limits_{i}{{\mathbf{h}}_{u,b,i}^{DD1}}.$$ (21)

$$h_{u,b,i}^{DD1}\left(t\right)=\left\{{\begin{array}[]{*{20}{c}}{{h_{u,b,i}}\varphi\left({{{k^{\prime}}_{i}}+t^{\prime},{{k^{\prime}}_{i}},{{\tilde{k}^{\prime}}_{i}},{{\tilde{l}^{\prime}}_{i}},N^{\prime},M^{\prime}}\right)}&{\left\{{\begin{array}[]{*{20}{c}}{{\text{if }}{{k^{\prime}}_{i}}+t^{\prime}<0{\text{, }}t=({{k^{\prime}}_{\max}}+2\varepsilon^{\prime}+1)+{{k^{\prime}}_{i}}+t^{\prime}+1}\\ {{\text{if }}{{k^{\prime}}_{i}}+t^{\prime}\geq 0{\text{, }}t={{k^{\prime}}_{i}}+t^{\prime}+1}\end{array}}\right.}\\ 0&{{\text{otherwise}}}\end{array}}\right.$$ (22)

${{\mathbf{h}}_{u,b,i}^{DD1}}$ is presented as in equation (22), where $t^{\prime}$ is an integer and satisfies $-\varepsilon^{\prime}\leq t^{\prime}\leq\varepsilon^{\prime}$. Combining equations (16) and (22), ${\mathbf{A}}_{u}^{p1}={{\mathbf{X}}_{u}^{p1}\odot{\mathbf{\Phi^{\prime}}}}$ can be considered as the known measurement matrix at the AP, ${\mathbf{y}}_{b}^{p1}$ as the observed vector and ${\mathbf{h}}_{u,b}^{DD1}$ as an unknown sparse vector that can be recovered using compressive sensing methods. In a multi-user scenario, the dimension of the sparse vector expands, allowing us to utilize this model to detect the indices of non-zero entries in the sparse vector and thereby identify potential active UEs. Given the coarse approximations made during rough AUD, especially under conditions where $M^{\prime}$ is notably small, accurate estimation of channel parameters becomes challenging. Therefore, it necessitates further refinement based on initial rough detection for accurate AUD and CE. Detailed elaboration on this can be found in subsequent subsections.

III-B Accurate AUD

Inspired by embedded pilot scheme, within a transmission block, a portion of the Doppler-delay grid resources will be utilized as preamble sequences for accurate AUD and CE. The dimensions $N$ and $M$ corresponding to this transmission block can be sufficiently large, allowing the received signal model to be formulated as shown in equation (11). Let ${k_{\max}}=\left\lfloor{{\nu_{\max}}NT}\right\rfloor$ and ${l_{\max}}={\tau_{\max}}M\Delta f$. Observing equation (11), we can see that the $(k,l)$-th DD domain received signal is affected by the transmitted signals with range of $[k-{k_{\max}}-\varepsilon:k+\varepsilon,l-{l_{\max}}:l]$. Therefore, to avoid interference between preamble and data, a guard interval needs to be estabilshed, where symbols within this interval are set to zero, as illustrated in Fig. 2.

Figure 2: Symbols arrangement for active UE.

Assuming the starting coordinates of the preamble symbol are $(k_{p},l_{p})$, with the preamble having a dimension of $L_{p}$ on the delay axis and $K_{p}$ on the Doppler axis. If we set $l_{p}-l_{max}\geq l_{max}$ and $l_{p}+Lp<M$, for the received model in equation (11) and $l\in\left[{{l_{p}}-{l_{\max}},{l_{p}}+{L_{p}}}\right]$, we only need to consider the case when $l\geq{l_{u,b,i}}$ (since $l\geq{l_{\max}}\geq{l_{u,b,i}}$). Let ${{\mathbf{X}}_{u,p}}={{\mathbf{X}}_{u}}\left[{{k_{p}}:{k_{p}}+{K_{p}}+{k_{\max}}+\varepsilon,{l_{p}}:{l_{p}}+{L_{p}}+{l_{\max}}}\right]$, ${\mathbf{Y}}_{b}^{DDp}={{\mathbf{Y}}_{b}}\left[{{k_{p}}:{k_{p}}+{K_{p}}+{k_{\max}},{l_{p}}:{l_{p}}+{L_{p}}+{l_{\max}}}\right]$, $\mathbf{X}_{u}\in\mathbb{C}^{N\times M}$ and $\mathbf{Y}_{b}\in\mathbb{C}^{N\times M}$ are DD domain transmitted symbols of $u$-th UE and received symbols of $b$-th AP, respectively. ${\mathbf{y}}_{b}^{p2}=\text{vec}({{\mathbf{Y}}_{b}^{DDp}})$, $N_{p}=K_{p}+k_{max}$, $M_{p}=L_{p}+l_{max}$ and we define:

$$\psi\left({k,{k_{i}},{{\tilde{k}}_{i}},{l_{i}},N,M}\right) =\frac{1}{N}\frac{{1-{e^{j2\pi{{\tilde{k}}_{i}}}}}}{{1-{e^{-j2\pi\frac{{k-{k_{i}}-{{\tilde{k}}_{i}}}}{N}}}}}{e^{-j2\pi\frac{{{l_{i}}\left({{k_{i}}+{{\tilde{k}}_{i}}}\right)}}{{NM}}}},$$ (23a)

$$\phi\left({k,l}\right) ={e^{j2\pi\frac{{kl}}{{NM}}}}.$$ (23b)

With equation (15), we can similarly obtain:

$$\begin{gathered}{\mathbf{y}}_{b}^{p2}\approx\left({{\mathbf{X}}_{u}^{p2}\odot{\mathbf{\Phi}}}\right){\mathbf{h}}_{u,b}^{DD2}+{\mathbf{n}}_{b}^{DD2}\hfill\\ ={\mathbf{A}}_{u}^{p2}{\mathbf{h}}_{u,b}^{DD2}+{\mathbf{n}}_{b}^{DD2}\in{\mathbb{C}^{{N_{p}}{M_{p}}\times 1}},\hfill\\ \end{gathered}$$ (24)

where ${\mathbf{X}}_{u}^{p2}$ is expressed as in equation (25), and

$$\begin{gathered}{\mathbf{X}}_{u}^{p2}=\left[{\begin{array}[]{*{20}{c}}{{C_{c}}\left({{{\mathbf{X}}_{u,p}}\left({:,1}\right),{k_{\max}},\varepsilon}\right)}\\ {{C_{c}}\left({{{\mathbf{X}}_{u,p}}\left({:,2}\right),{k_{\max}},\varepsilon}\right)}\\ \vdots\\ {{C_{c}}\left({{{\mathbf{X}}_{u,p}}\left({:,{L_{p}}+{l_{\max}}}\right),{k_{\max}},\varepsilon}\right)}\end{array}\begin{array}[]{*{20}{c}}{{C_{c}}\left({{{\mathbf{X}}_{u,p}}\left({:,{L_{p}}+{l_{\max}}}\right),{k_{\max}},\varepsilon}\right)}\\ {{C_{c}}\left({{{\mathbf{X}}_{u,p}}\left({:,1}\right),{k_{\max}},\varepsilon}\right)}\\ \vdots\\ {{C_{c}}\left({{{\mathbf{X}}_{u,p}}\left({:,{L_{p}}+{l_{\max}}-1}\right),{k_{\max}},\varepsilon}\right)}\end{array}\begin{array}[]{*{20}{c}}{}\hfil\\ {}\hfil\\ \cdots\\ {}\hfil\end{array}\begin{array}[]{*{20}{c}}{{C_{c}}\left({{{\mathbf{X}}_{u,p}}\left({:,{L_{p}}+1}\right),{k_{\max}},\varepsilon}\right)}\\ {{C_{c}}\left({{{\mathbf{X}}_{u,p}}\left({:,{L_{p}}+2}\right),{k_{\max}},\varepsilon}\right)}\\ \vdots\\ {{C_{c}}\left({{{\mathbf{X}}_{u,p}}\left({:,{L_{p}}}\right),{k_{\max}},\varepsilon}\right)}\end{array}}\right]\hfill\\ \end{gathered}$$ (25)

$${C_{c}}\left({{\mathbf{x}},k,\varepsilon}\right)={\left[{C\left({{\mathbf{x}},k,\varepsilon}\right)}\right]_{1:\dim({\mathbf{x}})-\varepsilon,:}},$$ (26)

$${\mathbf{\Phi}}={\left[{{\mathbf{\Phi}}_{1}^{H},\ldots,{\mathbf{\Phi}}_{l}^{H},\ldots,{\mathbf{\Phi}}_{{M_{p}}}^{H}}\right]^{H}},$$ (27)

$${{\mathbf{\Phi}}_{l}}={{\mathbf{1}}^{1\times({l_{\max}}+1)}}\otimes\bm{\phi}_{l}^{T}\otimes{{\mathbf{1}}^{{N_{p}}\times 1}},$$ (28)

$$\begin{gathered}{\bm{\phi}_{l}}=[\phi(l+{l_{p}}-1,0),\ldots,\phi(l+{l_{p}}-1,{k_{\max}}+\varepsilon),\hfill\\ \phi(l+{l_{p}}-1,{N_{p}}-\varepsilon-1),\ldots,\phi(l+{l_{p}}-1,{N_{p}}-1){]^{T}}.\hfill\\ \end{gathered}$$ (29)

${\mathbf{h}}_{u,b}^{DD2}\in{\mathbb{C}^{({k_{\max}}+2\varepsilon+1)({l_{\max}}+1)\times 1}}$ expressed as

$${\mathbf{h}}_{u,b}^{DD2}=\sum\limits_{i}{{\mathbf{h}}_{u,b,i}^{DD2}}.$$ (30)

$$h_{u,b,i}^{DD2}\left(t\right)=\left\{{\begin{array}[]{*{20}{c}}{{h_{u,b,i}}\psi\left({{k_{i}}+t^{\prime},{k_{i}},{{\tilde{k}}_{i}},{l_{i}},N,M}\right)}&{\left\{{\begin{array}[]{*{20}{c}}{{\text{if }}{k_{i}}+t^{\prime}<0{\text{, }}t=({l_{i}}+1)({k_{\max}}+2\varepsilon+1)+{k_{i}}+t+1}\\ {{\text{if }}{k_{i}}+t^{\prime}\geq 0{\text{, }}t={l_{i}}({k_{\max}}+2\varepsilon+1)+{k_{i}}+t^{\prime}+1}\end{array}}\right.}\\ 0&{{\text{otherwise}}}\end{array}}\right.$$ (31)

${{\mathbf{h}}_{u,b,i}^{DD2}}$ is presented as in equation (31), where $t^{\prime}$ is an integer which satisfies $-\varepsilon\leq t^{\prime}\leq\varepsilon$. According to equation (24), ${\mathbf{A}}_{u}^{p2}={{\mathbf{X}}_{u}^{p2}\odot{\mathbf{\Phi}}}$ can be considered as the known measurement matrix at the AP, ${\mathbf{y}}_{b}^{p2}$ as the observed vector and ${\mathbf{h}}_{u,b}^{DD2}$ as an unknown sparse vector that can be recovered using compressive sensing methods. Since in accurate AUD, $M>M^{\prime}$ and $N>N^{\prime}$, which make a higher resolution in delay and Doppler, resulting in more precise quantization. However, compared to rough AUD, because ${k_{\max}}>{k^{\prime}_{\max}}$ and ${l_{\max}}\gg{\tilde{l}^{\prime}_{\max}}$, the dimension of the sparse vector in accurate AUD is larger, making it more difficult to recover the sparse vector in multi-UE scenarios. Therefore, we need to utilize both rough and accurate AUD to achieve precise detection and estimation with lower computational complexity.

III-C Hybrid Preamble Transmission Block Structure and Multi-UE Joint Active Detection and Channel Estimation

In the $N\times M$ DD domain, we superimpose the preamble1, denoted as ${\mathbf{X}}_{u,p1}^{DD}$, for rough AUD, and the block symbols ${\mathbf{X}}_{u,2}^{DD}$, which includes the preamble2 ${\mathbf{X}}_{u,p2}$ and data symbols ${\mathbf{X}}_{u,d}$, for accurate AUD and CE. The superimposed result forms transmission block, structured as shown in the Fig. 3.

Figure 3: The hybrid preamble transmission block structure.

Then we have

$$\begin{gathered}{\mathbf{X}}_{u}^{DD}={\mathbf{X}}_{u,1}^{DD}+{\mathbf{X}}_{u,2}^{DD},\hfill\\ {\mathbf{X}}_{u,1}^{DD}(k,l)=\left\{{\begin{array}[]{*{20}{c}}{{\mathbf{X}}_{u,p1}^{DD}(k^{\prime},l^{\prime})}&{k=k^{\prime},l=\frac{{l^{\prime}}}{\beta}}\\ 0&{{\text{otherwise}}}\end{array}}\right.,\hfill\\ {\mathbf{X}}_{u,2}^{DD}(k,l)=\left\{{\begin{array}[]{*{20}{c}}{{\mathbf{X}}_{u,p2}^{DD}(k^{\prime},l^{\prime})}&{k=k^{\prime}+{k_{p}},l=l^{\prime}+{l_{p}}}\\ {{\mathbf{X}}_{u,d}^{DD}(k,l)}&{\left({k,l}\right){\text{ not in }}{\mathcal{P}\mathcal{G}}{\text{ area}}}\\ 0&{{\text{otherwise}}}\end{array}}\right.,\hfill\\ \end{gathered}$$ (32)

The $\mathcal{PG}$ area represents the grids designated for placing preamble2 and the guard intervals, as illustrated by the light green area in Fig. 3. As can be seen, the preamble1 is placed at intervals of $\frac{1}{\beta}$ along the delay dimension while being placed continuously along the Doppler dimension. Since the delay dimension $M^{\prime}=\beta M$ of preamble1 is assumed to be very small, and the Doppler dimension satisfies $N^{\prime}=\alpha N\leq\frac{N}{2}$, there is sufficient space within DD dimension to place preamble2 and its corresponding guard interval. This arrangement ensures that the received signal of preamble1 does not interfere with the received signal of preamble2.

Building on this, we allocate different power levels to ${\mathbf{X}}_{u,1}^{DD}$ and ${\mathbf{X}}_{u,2}^{DD}$, ensuring a significant difference in energy domain between these two types of signals. The received signal in the time-frequency domain can be expressed as:

$$\begin{gathered}{\mathbf{Y}_{b}}\left[{n,m}\right]={\mathbf{Y}_{b1}}\left[{n,m}\right]+{\mathbf{Y}_{b2}}\left[{n,m}\right]+\mathbf{N}\left[{n,m}\right]\hfill\\ \mathop{=}\limits^{a}{\mathbf{Y}_{b1}}\left[{n,m}\right]+\mathbf{\tilde{Z}}\left[{n,m}\right],\hfill\\ \end{gathered}$$ (33)

where $\mathbf{\tilde{Z}}\left[{n,m}\right]={\mathbf{Y}_{b2}}\left[{n,m}\right]+\mathbf{N}\left[{n,m}\right]$ is treated as noise. Assuming that applying ISFFT to ${\mathbf{X}}_{u,p1}^{DD}\in{\mathbb{C}^{N^{\prime}\times M^{\prime}}}$ results in a TF domain signal ${\mathbf{X}}_{u,p1}^{TF}\in{\mathbb{C}^{N^{\prime}\times M^{\prime}}}$, and applying ISFFT to ${\mathbf{X}}_{u,1}^{DD}(k,l)\in{\mathbb{C}^{N\times M}}$ to obtain a TF domain signal ${\mathbf{X}}_{u,1}^{TF}\in{\mathbb{C}^{N\times M}}$, both signals pass through the same channel to arrive at the $b$-th AP. After performing the Wigner transform, the TF domain received signals are ${\mathbf{Y}}_{u,b,p1}^{TF}\in{\mathbb{C}^{N^{\prime}\times M^{\prime}}}$ and ${\mathbf{Y}}_{u,b,1}^{TF}\in{\mathbb{C}^{N\times M}}$ respectively. Based on equations (1), (10), and (31), we can derive:

$${\mathbf{X}}_{u,p1}^{TF}[n^{\prime},m^{\prime}]=\frac{1}{{\sqrt{\alpha\beta}}}{\mathbf{X}}_{u,1}^{TF}[\frac{{n^{\prime}}}{\alpha},m^{\prime}],$$ (34)

$$\begin{gathered}{\mathbf{Y}}_{u,b,1}^{TF}[\frac{{n^{\prime}}}{\alpha},m^{\prime}]=\sqrt{\alpha\beta}\frac{1}{T}\sum\limits_{i}{}{h_{u,b,i}}\sum\limits_{m^{\prime\prime}}{}{\mathbf{X}}_{u,p1}^{TF}[n^{\prime},m^{\prime\prime}]\hfill\\ {e^{-j2\pi m^{\prime\prime}\Delta f{\tau_{u,b,i}}}}{e^{-j2\pi{\nu_{u,b,i}}{\tau_{u,b,i}}}}{e^{j2\pi{\nu_{u,b,i}}\frac{{n^{\prime}}}{\alpha}T}}\hfill\\ \int_{{\tau_{u,b,i}}}^{T}{}{e^{-j2\pi\Delta ft(m^{\prime}-m^{\prime\prime}-\frac{{{\nu_{u,b,i}}}}{{\Delta f}})}}dt+\mathbf{N}_{u,b,1}^{TF}[\frac{{n^{\prime}}}{\alpha},m^{\prime}],\hfill\\ \end{gathered}$$ (35)

$$\begin{gathered}{\mathbf{Y}}_{u,b,p1}^{TF}[n^{\prime},m^{\prime}]=\frac{1}{T}\sum\limits_{i}{}{h_{u,b,i}}\sum\limits_{m^{\prime\prime}}{\mathbf{X}}_{u,p1}^{TF}[n^{\prime},m^{\prime\prime}]\hfill\\ {e^{-j2\pi m^{\prime\prime}\Delta f{\tau_{u,b,i}}}}{e^{-j2\pi{\nu_{u,b,i}}{\tau_{u,b,i}}}}{e^{j2\pi{\nu_{u,b,i}}n^{\prime}T}}\hfill\\ \int_{{\tau_{u,b,i}}}^{T}{}{e^{-j2\pi\Delta ft(m^{\prime}-m^{\prime\prime}-\frac{{{\nu_{u,b,i}}}}{{\Delta f}})}}dt+\mathbf{N}_{u,b,p1}^{TF}[n^{\prime},m^{\prime}].\hfill\\ \end{gathered}$$ (36)

The detailed derivation can be referred to Appendix C. By comparing equations (35) and (36), it is apparent that for $0\leq n^{\prime}\leq N^{\prime}$ and $0\leq m^{\prime}\leq M^{\prime}$, ${\mathbf{Y^{\prime}}}_{u,b,1}^{TF}=\left[{\sqrt{\alpha\beta}{\mathbf{Y}}_{u,b,1}^{TF}(n^{\prime\prime},m^{\prime})|n^{\prime\prime}=\frac{{n^{\prime}}}{\alpha}}\right]\in{\mathbb{C}^{N^{\prime}\times M^{\prime}}}$ can be approximated as the received signal of ${\mathbf{X}}_{u,p1}^{TF}(n^{\prime},m^{\prime})$ through a channel with the same parameters, except that the Doppler parameter is $\frac{1}{\alpha}$ times the original parameter. Therefore, the maximum Doppler quantization parameter ${k^{\prime}_{\max}}$ will also become $\frac{1}{\alpha}$ times the original value. Based on this inference, we apply an $N^{\prime}\times M^{\prime}$ SFFT to ${\mathbf{Y^{\prime}}}_{u,b,1}^{TF}$ and perform rough AUD on this signal with the maximum Doppler quantization parameter ${k^{\prime}_{\max}}=[\frac{1}{\alpha}N^{\prime}\Delta f{\nu_{\max}}]_{\text{R}}$. In the multi-UE scenario, the reception model for preamble1, as described in equation (16), can be written as:

$${\mathbf{y}}_{b}^{p1}\approx{{\mathbf{A}}^{p1}}{\mathbf{h}}_{b}^{DD1}+{\mathbf{n}}_{b}^{DD}\in{\mathbb{C}^{N^{\prime}M^{\prime}\times 1}},$$ (37)

where ${\mathbf{y}}_{b}^{p1}$ is the vectorization of SFFT result of ${\mathbf{Y^{\prime}}}_{u,b,1}^{TF}$, and

$${{\mathbf{A}}^{p1}}=[{\mathbf{A}}_{1}^{p1},{\mathbf{A}}_{2}^{p1},\ldots,{\mathbf{A}}_{U}^{p1}],$$ (38)

$$\begin{gathered}{\mathbf{h}}_{b}^{DD1}={\left[{{{\left({{\mathbf{h}}_{b,1}^{DD1}}\right)}^{H}},{{\left({{\mathbf{h}}_{b,2}^{DD1}}\right)}^{H}},\ldots,{{\left({{\mathbf{h}}_{b,U}^{DD1}}\right)}^{H}}}\right]^{H}}.\hfill\\ \end{gathered}$$ (39)

After completing the rough AUD, the receiver obtains information about active UEs. Each AP transmits the detection results, representing the set of active UEs, to CPU. The CPU merges these results to form the system-wide rough active UEs set, as ${\mathcal{U}_{a}}=\bigcup\nolimits_{b}{{\mathcal{U}_{b,a}}}$, where ${\mathcal{U}_{b,a}}$ denotes the set of active UEs detected by the $b$-th AP. Assuming that for $1\leq i\leq\left|{{\mathcal{U}_{a}}}\right|$, we have $u_{i}\in\mathcal{U}_{a}$, and with the presence of guard intervals, preamble1 does not interfere with the reception of preamble2. Thus, equation (24) is rewritten as:

$${\mathbf{y}}_{b}^{p2}\approx{{\mathbf{A}}^{p2}}{\mathbf{h}}_{b}^{DD2}+{\mathbf{n}}_{b}^{DD2}\in{\mathbb{C}^{{N_{p}}{M_{p}}\times 1}},$$ (40)

where

$$\begin{gathered}{{\mathbf{A}}^{p2}}=[{\mathbf{A}}_{{u_{1}}}^{p2},{\mathbf{A}}_{{u_{2}}}^{p2},\ldots,{\mathbf{A}}_{{u_{\left|{{\mathcal{U}_{a}}}\right|}}}^{p2}],\hfill\\ \end{gathered}$$ (41)

$$\begin{gathered}{\mathbf{h}}_{b}^{DD2}={\left[{{{\left({{\mathbf{h}}_{b,{u_{1}}}^{DD2}}\right)}^{H}},{{\left({{\mathbf{h}}_{b,{u_{2}}}^{DD2}}\right)}^{H}},\ldots,{{\left({{\mathbf{h}}_{b,{u_{\left|{{\mathcal{U}_{a}}}\right|}}}^{DD2}}\right)}^{H}}}\right]^{H}}.\hfill\\ \end{gathered}$$ (42)

Since $\left|{{\mathcal{U}_{a}}}\right|\ll U$, the dimension of the sparse vector to be recovered is smaller than that of the estimated vector in a scheme that solely performs accurate AUD with the same sparsity (i.e., the same number of non-zero elements) and observation vector. Therefore, a more accurate estimation can be achieved by the hybrid preamble scheme. Moreover, it is necessary to appropriately choose the values for parameters such as $N$, $M$, $\alpha$, $\beta$.

After obtaining active UEs and their corresponding channels, the influence of preamble1 on the received signal can be removed by using successive interference cancellation (SIC). Based on the residual received signal and estimated channel parameters, the data signal can be recovered using algorithms such as belief propagation. The system’s signal processing flow can be seen in the diagram below. We won’t explore this part of data recovery in this paper.

Figure 4: The system’s signal processing flow.

III-D Multi-Antenna Reception Model

We extend the single-antenna reception model to a multi-antenna reception model. Suppose each AP is equipped with a uniform planar antenna array of dimension $N_{y}\times N_{z}$, with an antenna spacing of half a wavelength. The received signal matrix for preamble1 is represented as:

$$\begin{gathered}{\mathbf{Y}}_{b}^{DDS1}=\sum\limits_{u}{}\sum\limits_{i}{}{\mathbf{A}}_{u}^{p1}{\mathbf{h}}_{u,b,i}^{DD1}{{\mathbf{a}}_{s}}({\theta_{u,b,i,e}},{\theta_{u,b,i,a}})+{\mathbf{N}}_{b}^{DDS1}\hfill\\ ={{\mathbf{A}}^{p1}}{\mathbf{H}}_{b}^{DDS1}+{\mathbf{N}}_{b}^{DDS1}\in{\mathbb{C}^{NM\times{N_{y}}{N_{z}}}},\hfill\\ \end{gathered}$$ (43)

where

$$\begin{gathered}{{\mathbf{a}}_{s}}({\theta_{u,b,i,e}},{\theta_{u,b,i,a}})=[{e^{j\pi\sin{\theta_{u,b,i,e}}\sin{\theta_{u,b,i,a}}}}{e^{j\pi\cos{\theta_{u,b,i,e}}}},\hfill\\ \ldots,{e^{j\pi{n_{y}}\sin{\theta_{u,b,i,e}}\sin{\theta_{u,b,i,a}}}}{e^{j\pi{n_{z}}\cos{\theta_{u,b,i,e}}}},\ldots,\hfill\\ {e^{j\pi{N_{y}}\sin{\theta_{u,b,i,e}}\sin{\theta_{u,b,i,a}}}}{e^{j\pi{N_{z}}\cos{\theta_{u,b,i,e}}}}]\in{\mathbb{C}^{1\times{N_{y}}{N_{z}}}},\hfill\\ \end{gathered}$$ (44)

$$\begin{gathered}{\mathbf{H}}_{b}^{DDS}={\left[{{{\left({{\mathbf{H}}_{1,b}^{DDS}}\right)}^{H}},{{\left({{\mathbf{H}}_{2,b}^{DDS}}\right)}^{H}},\ldots,{{\left({{\mathbf{H}}_{U,b}^{DDS}}\right)}^{H}}}\right]^{H}},\hfill\\ \end{gathered}$$ (45)

$${\mathbf{H}}_{u,b}^{DDS}=\sum\limits_{i}{}{\mathbf{h}}_{u,b,i}^{DD}{{\mathbf{a}}_{s}}({\theta_{u,b,i,e}},{\theta_{u,b,i,a}}),$$ (46)

with ${\theta_{u,b,i,e}}$ and ${\theta_{u,b,i,a}}$ denote the elevation angle and azimuth angle, respectively, of the $i$-th path from $u$-th UE to $b$-th AP. The definition of ${\mathbf{h}}_{u,b,i}^{DD}$ can be found in equation (22). Applying the discrete fourier transform along the rows of ${\mathbf{Y}}_{b}^{DDS}$ with dimensions $N_{y}$ and $N_{z}$ results in a new matrix ${\mathbf{Y}}_{b}^{DDA}$. This allows transforming equation (43) into the angular domain. Specifically, for $0\leq{n_{y}}\leq{N_{y}}-1$ and $0\leq{n_{z}}\leq{N_{z}}-1$, we have

$${\mathbf{Y}}_{b}^{DDA}(t ,{n_{y}}+{N_{y}}{n_{z}}+1)=\frac{1}{{\sqrt{{N_{y}}{N_{z}}}}}\sum\limits_{{{n^{\prime}}_{y}}=0}^{{N_{y}}-1}{}\sum\limits_{{{n^{\prime}}_{z}}=0}^{{N_{z}}-1}{}$$ (47)

$${\mathbf{Y}}_{b}^{DDS}(t,{n^{\prime}_{y}}+{N_{y}}{n^{\prime}_{z}}+1){e^{-j\frac{{2\pi{{n^{\prime}}_{y}}{n_{y}}}}{{{N_{y}}}}}}{e^{-j\frac{{2\pi{{n^{\prime}}_{z}}{n_{z}}}}{{{N_{z}}}}}}.$$

Where $t$ is the row index. Consequently, we obtain:

$$\begin{gathered}{\mathbf{Y}}_{b}^{DDA1}=\sum\limits_{u}{}\sum\limits_{i}{}{\mathbf{A}}_{u}^{p1}{\mathbf{h}}_{u,b,i}^{DD1}{{\mathbf{a}}_{a}}({\theta_{u,b,i,e}},{\theta_{u,b,i,a}})+{\mathbf{N}}_{b}^{DDA1}\hfill\\ ={{\mathbf{A}}^{p1}}{\mathbf{H}}_{b}^{DDA1}+{\mathbf{N}}_{b}^{DDA1}\in{\mathbb{C}^{NM\times{N_{y}}{N_{z}}}},\hfill\\ \end{gathered}$$ (48)

where

$$\begin{gathered}{{\mathbf{a}}_{a}}({\theta_{u,b,i,e}},{\theta_{u,b,i,a}})=[{a_{a}}(0,0,{\theta_{u,b,i,e}},{\theta_{u,b,i,a}}),\ldots,{a_{a}}({n_{y}},\hfill\\ {n_{z}},{\theta_{u,b,i,e}},{\theta_{u,b,i,a}}),\ldots,{a_{a}}({N_{y}}-1,{N_{z}}-1,{\theta_{u,b,i,e}},{\theta_{u,b,i,a}})],\hfill\\ \end{gathered}$$ (49)

$${a_{a}}({n_{y}},{n_{z}}, {\theta_{u,b,i,e}},{\theta_{u,b,i,a}})=\frac{1}{{\sqrt{{N_{y}}{N_{z}}}}}\sum\limits_{{{n^{\prime}}_{y}}=0}^{{N_{y}}-1}{}\sum\limits_{{{n^{\prime}}_{z}}=0}^{{N_{z}}-1}{}$$ (50)

$${e^{-j\frac{{2\pi{{n^{\prime}}_{y}}}}{{{N_{y}}}}\left({{n_{y}}-{{\rm O}_{u,b,i}}}\right)}}{e^{-j\frac{{2\pi{{n^{\prime}}_{z}}}}{{{N_{z}}}}\left({{n_{z}}-{\Omega_{u,b,i}}}\right)}},$$

$${{\rm O}_{u,b,i}}=\frac{{\sin{\theta_{u,b,i,e}}\sin{\theta_{u,b,i,a}}{N_{y}}}}{2},$$ (51)

$${\Omega_{u,b,i}}=\frac{{\cos{\theta_{u,b,i,e}}{N_{z}}}}{2},$$ (52)

$$\begin{gathered}{\mathbf{H}}_{b}^{DDA}={\left[{{{\left({{\mathbf{H}}_{1,b}^{DDA}}\right)}^{H}},{{\left({{\mathbf{H}}_{2,b}^{DDA}}\right)}^{H}},\ldots,{{\left({{\mathbf{H}}_{U,b}^{DDA}}\right)}^{H}}}\right]^{H}},\hfill\\ \end{gathered}$$ (53)

$${\mathbf{H}}_{u,b}^{DDA}=\sum\limits_{i}{}{\mathbf{h}}_{u,b,i}^{DD}{{\mathbf{a}}_{a}}({\theta_{u,b,i,e}},{\theta_{u,b,i,a}}).$$ (54)

From equations (49) and (50), we can see that the row vector ${{\mathbf{a}}_{a}}({\theta_{u,b,i,e}},{\theta_{u,b,i,a}})$ is a block sparse vector, and ${\mathbf{h}}_{u,b,i}^{DD1}$ is also a block sparse vector in the DD domain. Additionally, given the inherent sparsity of active UEs, ${\mathbf{H}}_{b}^{DDA1}$ exhibits block sparse characteristics, as illustrated in the Fig. 5.

Figure 5: (a) The row sparsity of channel matrix ${\mathbf{H}}_{u,b}^{DDS}$ for an active UE $u$ in the DD-spatial domain and (b) The block sparsity of channel matrix ${\mathbf{H}}_{u,b}^{DDA}$ for an active UE $u$ in the DD-angular domain.

Similarly, the received signal matrix for preamble2 in a massive MIMO system can also be expressed in the angular domain:

$${\mathbf{Y}}_{b}^{DDA2}={{\mathbf{A}}^{p2}}{\mathbf{H}}_{b}^{DDA2}+{\mathbf{N}}_{b}^{DDA2}\in{\mathbb{C}^{{N_{p}}{M_{p}}\times{N_{y}}{N_{z}}}},$$ (55)

$${\mathbf{H}}_{b}^{DDA2}={\left[{{{\left({{\mathbf{H}}_{{u_{1}},b}^{DDA2}}\right)}^{H}},{{\left({{\mathbf{H}}_{{u_{2}},b}^{DDA2}}\right)}^{H}},\ldots,{{\left({{\mathbf{H}}_{{u_{\left|{{U_{a}}}\right|}},b}^{DDA2}}\right)}^{H}}}\right]^{H}},$$ (56)

$${\mathbf{H}}_{u,b}^{DDA2}=\sum\limits_{i}{}{\mathbf{h}}_{u,b,i}^{DD2}{{\mathbf{a}}_{a}}({\theta_{u,b,i,e}},{\theta_{u,b,i,a}}).$$ (57)

Also, ${\mathbf{H}}_{b}^{DDA2}$ exhibits block sparsity. Therefore, to achieve both rough AUD and accurate AUD along with CE, it is necessary to implement block sparse matrix recovery for the models in equations (48) and (55), respectively. Based on the estimated indices and values of the non-zero elements, we will detect the active UEs and estimate the channel matrix. The hybrid preamble based AUD and CE scheme is summarized as in Algorithm 1.

Algorithm 1 Hybrid Preamble Based AUD and CE Scheme

0:  $\left\{{{\mathbf{Y}}_{b}^{DDA1}}\right\}$, ${{\mathbf{A}}^{p1}}$, $\left\{{{\mathbf{Y}}_{b}^{DDA2}}\right\}$

0:  ${\bar{\mathcal{U}_{a}}}=\bigcup\nolimits_{b}{{{\bar{\mathcal{U}}}_{b,a}}}$, $\left\{{{\mathbf{\bar{H}}}_{b}^{DDA2}|1\leq b\leq B}\right\}$

1:  %Rough AUD

2:  for $b=1$ to $B$ do

3:     Recover ${\mathbf{H}}_{b}^{DDA1}$ based on ${\mathbf{Y}}_{b}^{DDA1}$ and ${\mathbf{A}}^{p1}$ by block sparse matrix recovery algorithm (such as GAMP-PCSBL-La proposed in Section IV);

4:     Obtain ${\mathcal{U}_{b,a}}$ based on non-zero entries of estimated ${\mathbf{H}}_{b}^{DDA1}$;

5:  end for

6:  Form ${{\mathbf{A}}^{p2}}$ based on ${\mathcal{U}_{a}}=\bigcup\nolimits_{b}{{\mathcal{U}_{b,a}}}$;

7:  %Accurate AUD and CE

8:  for $b=1$ to $B$ do

9:     Recover ${\mathbf{H}}_{b}^{DDA2}$ based on ${\mathbf{A}}^{p2}$ and ${\mathbf{Y}}_{b}^{DDA2}$ by block sparse matrix recovery algorithm (such as GAMP-PCSBL-La proposed in Section IV);

10:     Obtain accurate detected active UEs’ set ${\bar{\mathcal{U}}_{a}}$ and channel matrix $\left\{{{\mathbf{\bar{H}}}_{b,i}^{DDA2}|i\in{{\bar{\mathcal{U}}}_{a}}}\right\}$ based on non-zero entries of estimated ${\mathbf{H}}_{b}^{DDA2}$.

11:  end for

IV GAMP-PCSBL-La Algorithm

IV-A GAMP-PCSBL and Laplacian Prior

To achieve block sparse vector recovery, $\mathbf{y}=\mathbf{Ax}+\mathbf{n}$, Fang et.al proposed the PCSBL algorithm [28]. In this algorithm, akin to conventional sparse Baysian learning (SBL), the atoms of the block sparse vector are modeled using a two-layer hierarchical probability model. Specifically, it is assumed that there exists hyperparameter ${\bm{\alpha}}$, and the prior distribution of the sparse vector is given by:

$$p({\mathbf{x}})=\int{}p({\mathbf{x}}|{\bm{\alpha}})p({\bm{\alpha}})d{\bm{\alpha}}.$$ (58)

In the conventional SBL framework, the likelihood distribution $p({\mathbf{x}}|{\bm{\alpha}})$ is defined as a conditional Gaussian distribution, given by:

$$p({\mathbf{x}}|{\mathbf{\alpha}})=\prod\limits_{i}{\mathcal{C}\mathcal{N}\left({{x_{i}};0,\alpha_{i}^{-1}}\right)}.$$ (59)

At this point, the vector ${\bm{\alpha}}$ is defined as the precision of $\mathbf{x}$, where $\alpha_{i}^{-1}$ represents the variance of ${x_{i}}$. By assuming ${\alpha_{i}}$ follows a Gamma distribution, ${\bm{\alpha}}$ is ensured to be non-negative, which in turn can control the sparsity of vector ${\mathbf{x}}$. Specifically, when ${\alpha_{i}}$ is a large value, the corresponding $x_{i}$ tends to be zero. By maximizing the posterior probability $p({\bm{\alpha}}|{\mathbf{x}},{\mathbf{y}})$, the hyperparameters can be learned based on the estimated values of the sparse vector.

Unlike conventional SBL, PCSBL leverages the fact that the sparse patterns of adjacent coefficients exhibit a certain degree of correlation. Therefore, the likelihood probability is defined not as a Gaussian distribution dependent on a single hyperparameter, but as a conditional Gaussian distribution related to the hyperparameters of ”neighboring regions” [28]. Specifically, the likelihood distribution of ${\mathbf{x}}$ given ${\bm{\alpha}}$ is expressed as:

$$p({\mathbf{x}}|{\mathbf{\alpha}})=\prod\limits_{i}{\mathcal{C}\mathcal{N}\left({{x_{i}};0,\tau_{i}^{-1}}\right)},$$ (60)

$${\tau_{i}}={\alpha_{i}}+\beta{\alpha_{i-1}}+\beta{\alpha_{i+1}}.$$ (61)

$\beta$ represents the coupling factor, which characterizes the degree of correlation between the current element ${x_{i}}$ and the elements in its neighboring region $\left\{{{x_{i-1}},{x_{i+1}}}\right\}$. At $t$-th iteration, we first fix the current estimate of the hyperparameter ${{\mathbf{\alpha}}{(t)}}$. Then using Bayes’ theorem, the posterior distribution of ${\mathbf{x}}$ is given by:

$$\begin{gathered}p({\mathbf{x}}|{\mathbf{y}},{\mathbf{\alpha}})\propto p({\mathbf{y}}|{\mathbf{x}})p({\mathbf{x}}|{{\mathbf{\alpha}}^{(t)}})\hfill\\ =\mathcal{C}\mathcal{N}\left({{\mathbf{x}};{\mathbf{\mu}},{\mathbf{\Sigma}}}\right),\hfill\\ \end{gathered}$$ (62)

with

$$\mathbf{\Sigma}={\left({{\gamma^{-1}}{{\mathbf{A}}^{H}}{\mathbf{A}}+{\mathbf{D}}}\right)^{-1}},$$ (63)

$${\mathbf{D}}={\text{diag}}\left({{\tau_{1}},\ldots,{\tau_{N}}}\right),$$ (64)

$${\mathbf{\mu}}={\gamma^{-1}}{\mathbf{\Sigma}}{{\mathbf{A}}^{H}}{\mathbf{y}}.$$ (65)

Here, $\gamma$ and $\mathbf{A}$ represent the noise variance and measurement matrix, respectively. Subsequently, according to the expectation-maximization (EM) algorithm, hyperparameters ${{\mathbf{\alpha}}{(t+1)}}$ is updated. Specifically, we have:

$${{\mathbf{\alpha}}{(t+1)}}=\mathop{\arg\max}\limits_{\mathbf{\alpha}}{{\mathbb{E}}_{p({\mathbf{x}}|{\mathbf{y}},{\mathbf{\alpha}})}}\left[{\log\left({p\left({{\mathbf{x}}|{\mathbf{\alpha}}}\right)p\left({\mathbf{\alpha}}\right)}\right)}\right],$$ (66)

where ${{\mathbb{E}}_{p}}\left[\cdot\right]$ denotes expectation according to the distribution $p$. It should be noted that equation (63) involves a matrix inversion process, which makes the algorithm computationally intensive. To achieve lower complexity in computing the posterior distribution, Fang et.al. proposed the pattern-coupled sparse Bayesian learning via generalized approximate message passing (GAMP-PCSBL) algorithm [30]. Its core idea is to use the GAMP algorithm with $p({\mathbf{x}}|{{\bm{\alpha}}{(t)}})$ as the prior distribution, replacing the computational process of equation (63). The GAMP algorithm was introduced by Rangan et.al [31]. Based on the BP algorithm, it leverages the central limit theorem and quadratic approximation to avoid the computationally intensive marginal posterior probability calculations, especially effective in high-dimensional settings. Literature [30] concludes that GAMP-PCSBL performs similarly to conventional PCSBL in terms of performance but significantly reduces the average execution time.

Furthermore, inspired by reference [29], the Laplacian distribution can be utilized to capture the sparsity of signals after undergoing DCT. The probability density function of the Laplacian distribution with mean 0 and scale parameter $b$ is represented as follows:

$$\mathcal{L}\mathcal{A}(x;b)=\frac{1}{{2b}}\exp\left({-\frac{{\left|x\right|}}{{2b}}}\right).$$ (67)

Compared to the Gaussian distribution, the Laplacian distribution enhances sparsity constraints by placing more posterior probability mass along the axes, as described in references [32] and [33]. In the field of image processing, the Laplacian distribution is widely used as a prior distribution for approximating sparsity generated by natural images after DCT, achieving more accurate estimation results than Gaussian mixture model priors [34]. Considering the presence of fractional Doppler and arrival angles in our system model, DFT in respective domains also produces sparsity, similar to that in image DCT. Combining these superiorities, we propose a novel GAMP-PCSBL-La algorithm for two-dimensional block sparse matrix recovery. Compared to existing algorithms, our proposed algorithm achieves superior estimation performance with lower complexity.

IV-B Probability Model

Given an AWGN channel model:

$${\mathbf{\tilde{Y}}}={\mathbf{\tilde{A}\tilde{X}}}+{\mathbf{\tilde{N}}},$$ (68)

where ${\mathbf{\tilde{Y}}}\in{\mathbb{C}^{L\times M}}$ is the observed matrix, ${\mathbf{\tilde{A}}}\in{\mathbb{C}^{L\times K}}$ is the measurement matrix (its values are known at the receiver), ${\mathbf{\tilde{X}}}\in{\mathbb{C}^{K\times M}}$ is the block sparse matrix to be estimated, and ${\mathbf{\tilde{N}}}\in{\mathbb{C}^{L\times M}}$ is the additive noise matrix. Since the Laplacian distribution is defined only for real-valued random variables, we need to convert the complex form model of equation (68) into the following real equivalent model:

$$\begin{gathered}{\mathbf{Y}}={\mathbf{AX}}+{\mathbf{N}}\\ {\mathbf{Y}}\buildrel\Delta\over{=}\left[{\begin{array}[]{*{20}{c}}{\mathcal{R}\left\{{{\mathbf{\tilde{Y}}}}\right\}}\\ {\mathcal{I}\left\{{{\mathbf{\tilde{Y}}}}\right\}}\end{array}}\right]\in{\mathbb{R}^{2L\times M}}\\ {\mathbf{A}}\buildrel\Delta\over{=}\left[{\begin{array}[]{*{20}{c}}{\mathcal{R}\left\{{{\mathbf{\tilde{A}}}}\right\}}&{-\mathcal{I}\left\{{{\mathbf{\tilde{A}}}}\right\}}\\ {\mathcal{I}\left\{{{\mathbf{\tilde{A}}}}\right\}}&{\mathcal{R}\left\{{{\mathbf{\tilde{A}}}}\right\}}\end{array}}\right]\in{\mathbb{R}^{2L\times 2K}}\\ {\mathbf{X}}\buildrel\Delta\over{=}\left[{\begin{array}[]{*{20}{c}}{\mathcal{R}\left\{{{\mathbf{\tilde{X}}}}\right\}}\\ {\mathcal{I}\left\{{{\mathbf{\tilde{X}}}}\right\}}\end{array}}\right]\in{\mathbb{R}^{2K\times M}}\\ {\mathbf{N}}\buildrel\Delta\over{=}\left[{\begin{array}[]{*{20}{c}}{\mathcal{R}\left\{{{\mathbf{\tilde{N}}}}\right\}}\\ {\mathcal{I}\left\{{{\mathbf{\tilde{N}}}}\right\}}\end{array}}\right]\in{\mathbb{R}^{2L\times M}}\\ \end{gathered}$$ (69)

Here, $\mathcal{R}\left\{\cdot\right\}$ and $\mathcal{I}\left\{\cdot\right\}$ represent the operations of taking the real and imaginary parts of a complex matrix, respectively. In practical systems, the noise variance is often unpredictable. We assume that the communication between the transmitter and receiver occurs over an AWGN channel, i.e.:

$$p({\mathbf{Y}}|{\mathbf{Z}})=\prod\nolimits_{l,j}{\mathcal{N}\left({{y_{l,j}};{z_{l,j}},\gamma}\right)},$$ (70)

where ${z_{l,j}}$ is $(i,j)$-th element of matirx ${\mathbf{Z}}$, ${\mathbf{Z}}={\mathbf{AX}}$ and $\gamma$ denotes the noise variance. Referencing the two-layer hierarchical probabilistic model of SBL, we assume the hyperparameters $\left\{{{\alpha_{i,j}}}\right\}$ and establish the probability distribution of ${\mathbf{X}}$ as:

$$p({\mathbf{X}}|{\mathbf{\alpha}}) ={\prod\limits_{0<i<N+1,j}{\mathcal{L}\mathcal{A}\left({{x_{i,j}};\tau_{i,j}^{-1}}\right)}}{\prod\limits_{N<i<2N,j}{\mathcal{L}\mathcal{A}\left({{x_{i,j}};\tau_{i-N,j}^{-1}}\right)}}$$ (71)

$$={\prod\limits_{0<i<N+1,j}{\frac{1}{2}}{\tau_{i,j}}\exp\left({-\frac{1}{2}{\tau_{i,j}}\left|{{x_{i,j}}}\right|}\right)}$$

$${\prod\limits_{N<i<2N,j}{\frac{1}{2}}{\tau_{i-N,j}}\exp\left({-\frac{1}{2}{\tau_{i-N,j}}\left|{{x_{i,j}}}\right|}\right)}\hfill$$

$${\tau_{i,j}}={\alpha_{i,j}}+\beta{\alpha_{i-1,j}}+\beta{\alpha_{i+1,j}}+\beta{\alpha_{i,j-1}}+\beta{\alpha_{i,j+1}},$$ (72)

$$p({\alpha_{i,j}})=\mathcal{G}\mathcal{A}\left({{\alpha_{i,j}};a,b}\right).$$ (73)

Among them, $\mathcal{G}\mathcal{A}\left({{\alpha_{i,j}};a,b}\right)=\Gamma{\left(a\right)^{-1}}{b^{a}}\alpha_{i,j}^{a}{e^{-b{\alpha_{i,j}}}}$ denotes the Gamma distribution with parameters $a$ and $b$. $\Gamma\left(a\right)=\int_{0}^{\infty}{{t^{a-1}}{e^{-t}}dt}$ is the Gamma funcation. Equation (71) shows that the real and imaginary parts of ${\tilde{x}_{i,j}}$ share the same hyperparameter ${\tau_{i,j}}$, as defined in equation (72).

IV-C GAMP Algorithm for Sparse Matrix Recovery

As previously explained, given a prior distribution, the GAMP algorithm can achieve sparse signal recovery with relatively low computational complexity. Following the sum-product and max-sum forms of the BP algorithm, the GAMP algorithm uses Gaussian and quadratic approximations to provide the minimum mean square error (MMSE) estimation and maximum a posteriori (MAP) estimation of the sparse matrix, respectively. By defining scalar estimation functions, ${g_{in}}\left(\cdot\right)$ and ${g_{out}}\left(\cdot\right)$, the GAMP algorithm iteratively performs scalar operations at the input and output nodes to decompose the vector-valued estimation problem. Assuming that in $t$-th iteration, the prior distribution of the sparse matrix is expressed as $p({\mathbf{X}}|{{\bm{\alpha}}^{(t)}})$, with ${{\bm{\alpha}}^{(t)}}$ is the hyperparameter obtained in the $t$-th iteration. For AWGN channel, the GAMP algorithm process is shown from line 3 to line 13 in Algorithm 2.

Algorithm 2 GAMP-PCSBL-La

0:  ${\mathbf{Y}}$, ${\mathbf{A}}$, $p({\mathbf{X}}|{\bm{\alpha}})$, $p({\bm{\alpha}})$, $\beta$

0:  ${\mathbf{\hat{X}}}\left({t+1}\right)$, ${\bm{\alpha}}(t+1)$, and $\gamma(t+1)$

1:  Initialize: ${\bm{\alpha}}(1)$, ${\mathbf{\hat{X}}}\left(1\right)=\mathbf{0}$, ${\mathbf{S}}\left(0\right)=\mathbf{0}$, $\gamma(1)$, ${u_{i,j}^{x}(1)}$;

2:  for $t=1$ to $T$ do

3:     $\forall l,j$, $u_{l,j}^{p}(t)=\sum\nolimits_{i}{{{\left|{{a_{l,i}}}\right|}^{2}}}u_{i,j}^{x}(t)$

4:     $\forall l,j$, ${\hat{p}_{l,j}}(t)=\sum\nolimits_{i}{{a_{l,i}}{{\hat{x}}_{i,j}}(t)}-u_{l,j}^{p}(t){\hat{s}_{l,j}}(t-1)$

5:     $\forall l,j$, $u_{l,j}^{z}(t)=\frac{{u_{l,j}^{p}(t)\gamma(t)}}{{u_{l,j}^{p}(t)+\gamma(t)}}$

6:     $\forall l,j$, ${\hat{z}_{l,j}}(t)=\frac{{u_{l,j}^{p}(t){y_{l,j}}+\gamma(t){{\hat{p}}_{l,j}}(t)}}{{u_{l,j}^{p}(t)+\gamma(t)}}$

7:     $\forall l,j$, ${\hat{s}_{l,j}}(t)={g_{out}}\left({t,{{\hat{p}}_{l,j}}(t),{y_{l,j}},u_{l,j}^{p}(t)}\right)$

8:     $\forall l,j$, $u_{l,j}^{s}(t)=-\frac{{\partial{g_{out}}\left({t,{{\hat{p}}_{l,j}}(t),{y_{l,j}},u_{l,j}^{p}(t)}\right)}}{{\partial{{\hat{p}}_{l,j}}(t)}}$

9:     $\forall i,j$, $u_{i,j}^{r}(t)={\left[{\sum\nolimits_{l}{{{\left|{{a_{l,i}}}\right|}^{2}}u_{l,j}^{s}(t)}}\right]^{-1}}$

10:     $\forall i,j$, ${\hat{r}_{i,j}}(t)={\hat{x}_{i,j}}(t)+u_{i,j}^{r}(t)\sum\nolimits_{l}{{a_{l,i}}{{\hat{s}}_{l,j}}(t)}$

11:     $\forall i,j$, ${\tau_{i,j}}(t)={\alpha_{i,j}}(t)+\beta{\alpha_{i-1,j}}(t)+\beta{\alpha_{i+1,j}}(t)+\beta{\alpha_{i,j-1}}(t)+\beta{\alpha_{i,j+1}}(t)$

12:     $\forall i,j$, ${\hat{x}_{i,j}}(t+1)={g_{in}}\left({t,{{\hat{r}}_{i,j}}(t),{\tau_{i,j}}(t),u_{i,j}^{r}(t)}\right)$

13:     $\forall i,j$, $u_{i,j}^{x}(t+1)=u_{i,j}^{r}(t)\frac{{\partial{g_{in}}\left({t,{{\hat{r}}_{i,j}}(t),{\tau_{i,j}}(t),u_{i,j}^{r}(t)}\right)}}{{\partial{{\hat{r}}_{i,j}}(t)}}$

14:     $\forall i,j$, ${\alpha_{i,j}}(t+1)=\frac{a}{{b+{\omega_{i,j}}(t+1)+{\omega_{N+i,j}}(t+1)}}$

15:     Update $\gamma(t+1)=\frac{{\sum\nolimits_{l,j}{}{{\left\|{{y_{l,j}}-{{\hat{z}}_{l,j}}(t)}\right\|}^{2}}+u_{l,j}^{z}(t)}}{{2MN}}$

16:     If $\frac{{\left\|{{\mathbf{\hat{X}}}\left({t+1}\right)-{\mathbf{\hat{X}}}\left(t\right)}\right\|_{F}^{2}}}{{\left\|{{\mathbf{\hat{X}}}\left({t+1}\right)}\right\|_{F}^{2}}}<\varepsilon$, break

17:  end for

For sum-product GAMP, ${g_{out}}\left({t,{{\hat{p}}_{l,j}}(t),{y_{l,j}},u_{l,j}^{p}(t)}\right)$ is defined as

$${g_{out}}\left({t,{{\hat{p}}_{l,j}}(t),{y_{l,j}},u_{l,j}^{p}(t)}\right)=\frac{{{{\hat{z}}_{l,j}}(t)-{{\hat{p}}_{l,j}}(t)}}{{u_{l,j}^{p}(t)}}=\frac{{{y_{l,j}}-{{\hat{p}}_{l,j}}(t)}}{{u_{l,j}^{p}(t)+\gamma(t)}}.$$ (74)

Then $u_{l,j}^{s}(t)=-\frac{{\partial{g_{out}}\left({t,{{\hat{p}}_{l,j}}(t),{y_{l,j}},u_{l,j}^{p}(t)}\right)}}{{\partial{{\hat{p}}_{l,j}}(t)}}=\frac{1}{{u_{l,j}^{p}(t)+\gamma(t)}}$. Based on MMSE estimation, in input node, we have

$${g_{in}}\left({t,{{\hat{r}}_{i,j}}(t),{\tau_{i,j}}(t),u_{i,j}^{r}(t)}\right)={\mathbb{E}}_{p(x|r,\tau,u^{r})}\{{x_{i,j}}\},$$ (75)

$$\frac{{\partial{g_{in}}\left({t,{{\hat{r}}_{i,j}}(t),{\tau_{i,j}}(t),u_{i,j}^{r}(t)}\right)}}{{\partial{{\hat{r}}_{i,j}}(t)}}{\tau_{i,j}}(t)={\mathbb{V}}_{p(x|r,\tau,u^{r})}\{{x_{i,j}}\}.$$ (76)

${\hat{x}_{i,j}}(t+1)$ and $u_{i,j}^{x}(t+1)$ represent the mean and variance of the approximate posterior distribution of $(i,j)$-th element of the matrix to be estimated, respectively. In the sum-product derivation, the messages from the factor node $p\left({y|x}\right)$ to the variable node ${x_{i,j}}$ are approximated like:

$${\vec{m}_{{x_{i,j}}}}(t)\approx\mathcal{N}\left({{x_{i,j}};{{\hat{r}}_{i,j}}(t),u_{i,j}^{r}(t)}\right).$$ (77)

As previously mentioned, the prior of ${x_{i,j}}$ is defined as:

$$p({x_{i,j}}|{\tau_{\tilde{i},j}}(t))=\mathcal{L}\mathcal{A}\left({{x_{i,j}};{{\left({\tau_{\tilde{i},j}(t)}\right)}^{-1}}}\right),$$

with

$$\tilde{i}=\left\{{\begin{array}[]{*{20}{c}}i&{1\leq i\leq N}\\ {i-N}&{N+1\leq i\leq 2N}\end{array}}\right..$$

Therefore, the approximate posterior distribution of ${x_{i,j}}$ can be expressed as:

$$\begin{gathered}p\left({{x_{i,j}}|{{\hat{r}}_{i,j}}(t),{\tau_{i,j}}(t),u_{i,j}^{r}(t)}\right)\propto{{\vec{m}}_{{x_{i,j}}}}(t)p({x_{i,j}}|{\tau_{\tilde{i},j}}(t))\hfill\\ =\mathcal{N}\left({{x_{i,j}};{{\hat{r}}_{i,j}}(t),u_{i,j}^{r}(t)}\right)\mathcal{L}\mathcal{A}\left({{x_{i,j}};{{\left({\tau_{\tilde{i},j}(t)}\right)}^{-1}}}\right)\hfill\\ =\frac{1}{{{\psi_{i,j}}(t)}}\exp\left\{{-{\xi_{i,j}}\left(t\right)-\frac{1}{{2u_{i,j}^{r}(t)}}{{\left({{x_{i,j}}-{\varphi_{i,j}}\left(t\right)}\right)}^{2}}}\right\},\hfill\\ \end{gathered}$$ (78)

where

$${\psi_{i,j}}(t)=\int{\exp\left\{{-{\xi_{i,j}}\left(t\right)-\frac{1}{{2u_{i,j}^{r}(t)}}{{\left({{x_{i,j}}-{\varphi_{i,j}}\left(t\right)}\right)}^{2}}}\right\}d{x_{i,j}}}$$

$$=\int_{-\infty}^{0}{}\exp\left\{{{\tau_{\tilde{i},j}}(t){{\hat{r}}_{i,j}}(t)+\frac{1}{2}u_{i,j}^{r}(t){{\left({{\tau_{\tilde{i},j}}(t)}\right)}^{2}}}\right\}$$

$$\exp\left\{{-\frac{1}{{2u_{i,j}^{r}(t)}}{{\left({{x_{i,j}}-\left({{{\hat{r}}_{i,j}}(t)+u_{i,j}^{r}(t){\tau_{\tilde{i},j}}(t)}\right)}\right)}^{2}}}\right\}d{x_{i,j}}+$$

$$\int_{0}^{\infty}{}\exp\left\{{\frac{1}{2}u_{i,j}^{r}(t){{\left({{\tau_{\tilde{i},j}}(t)}\right)}^{2}}-{\tau_{\tilde{i},j}}(t){{\hat{r}}_{i,j}}(t)}\right\}$$

$$\exp\left\{{-\frac{1}{{2u_{i,j}^{r}(t)}}{{\left({{x_{i,j}}-\left({{{\hat{r}}_{i,j}}(t)-u_{i,j}^{r}(t){\tau_{\tilde{i},j}}(t)}\right)}\right)}^{2}}}\right\}d{x_{i,j}}$$

$$=\sqrt{2\pi u_{i,j}^{r}(t)}\left[\exp\left\{{-\xi_{i,j}^{-}\left(t\right)}\right\}Q\left({{\varphi_{i,j}^{-}(t)}}/{{\sqrt{u_{i,j}^{r}(t)}}}\right)+\right.$$

$$\left.\exp\left\{{-\xi_{i,j}^{+}\left(t\right)}\right\}Q\left({{{\varphi_{i,j}^{+}(t)}}/{{\sqrt{u_{i,j}^{r}(t)}}}}\right)\right],$$ (79)

$${\xi_{i,j}}\left(t\right)={\tau_{\tilde{i},j}}(t){\hat{r}_{i,j}}(t){\text{sign}}\left({{x_{i,j}}}\right)-\frac{1}{2}u_{i,j}^{r}(t){\left({{\tau_{\tilde{i},j}}(t)}\right)^{2}},$$ (80)

$${\varphi_{i,j}}\left(t\right)={\hat{r}_{i,j}}(t)-u_{i,j}^{r}(t){\tau_{\tilde{i},j}}(t){\text{sign}}\left({{x_{i,j}}}\right),$$ (81)

$$\xi_{i,j}^{-}\left(t\right)=-{\tau_{\tilde{i},j}}(t){\hat{r}_{i,j}}(t)-\frac{1}{2}u_{i,j}^{r}(t){\left({{\tau_{\tilde{i},j}}(t)}\right)^{2}},$$ (82)

$$\xi_{i,j}^{+}\left(t\right)={\tau_{\tilde{i},j}}(t){\hat{r}_{i,j}}(t)-\frac{1}{2}u_{i,j}^{r}(t){\left({{\tau_{\tilde{i},j}}(t)}\right)^{2}},$$ (83)

$$\varphi_{i,j}^{-}\left(t\right)={\hat{r}_{i,j}}(t)+u_{i,j}^{r}(t){\tau_{\tilde{i},j}}(t),$$ (84)

$$\varphi_{i,j}^{+}\left(t\right)={\hat{r}_{i,j}}(t)-u_{i,j}^{r}(t){\tau_{\tilde{i},j}}(t),$$ (85)

$${\text{sign}}({x_{i,j}})=\left\{{\begin{array}[]{*{20}{c}}1&{{x_{i,j}}>0}\\ 0&{{x_{i,j}}=0}\\ {-1}&{{x_{i,j}}<0}\end{array}}\right..$$ (86)

According to equation (78), the posterior mean and variance of ${x_{i,j}}$ can be calculated as in equations (87) and (88), where $Q\left(\cdot\right)$ is the standard Q-function, representing the tail probability of the normal distribution, defined as:

$${\hat{x}_{i,j}}(t+1)=\frac{{\sqrt{2\pi u_{i,j}^{r}(t)}}}{{{\psi_{i,j}}(t)}}\left[{{e^{-\xi_{i,j}^{-}(t)}}\varphi_{i,j}^{-}(t)Q\left({{{\varphi_{i,j}^{-}(t)}}/{{\sqrt{u_{i,j}^{r}(t)}}}}\right)+{e^{-\xi_{i,j}^{+}(t)}}\varphi_{i,j}^{+}(t)Q\left({{{\varphi_{i,j}^{+}(t)}}/{{\sqrt{u_{i,j}^{r}(t)}}}}\right)}\right]$$ (87)

$$\begin{gathered}u_{i,j}^{x}(t+1)=\frac{{\sqrt{2\pi u_{i,j}^{r}(t)}}}{{{\psi_{i,j}}(t)}}\left(\left({{{\left({\varphi_{i,j}^{+}(t)}\right)}^{2}}+u_{i,j}^{r}(t)}\right){e^{-\xi_{i,j}^{+}(t)}}Q\left({-{{\varphi_{i,j}^{+}(t)}}/{{\sqrt{u_{i,j}^{r}(t)}}}}\right)+\left({{{\left({\varphi_{i,j}^{-}(t)}\right)}^{2}}+u_{i,j}^{r}(t)}\right){e^{-\xi_{i,j}^{-}(t)}}\right.\hfill\\ Q\left({{{\varphi_{i,j}^{-}(t)}}/{{\sqrt{u_{i,j}^{r}(t)}}}}\right)\left.-\frac{{2{\tau_{i,j}}(t)u_{i,j}^{r}{{(t)}^{2}}}}{{\sqrt{2\pi u_{i,j}^{r}(t)}}}{e^{-{{{{\hat{r}}_{i,j}}{{(t)}^{2}}}}/{{2u_{i,j}^{r}(t)}}}}\right)-{\left({{{\hat{x}}_{i,j}}(t)}\right)^{2}}\hfill\\ \end{gathered}$$ (88)

$$Q\left(x\right)=\frac{1}{{\sqrt{2\pi}}}\int_{x}^{\infty}{{e^{-\frac{{{u^{2}}}}{2}}}du}.$$ (89)

The derivation process of equations (87) and (88) can be found in Appendix D. This completes the GAMP portion of Algorithm 2.

IV-D Learning Hyperparameters via EM Algorithm

After obtaining the posterior distribution of ${\mathbf{X}}$, our objective shifts to finding appropriate hyperparameters ${\bm{\alpha}}$ and $\gamma$ that maximize the posterior probability of them. A direct strategy is to use the EM algorithm, where ${\mathbf{X}}$ is treated as a hidden variable. In the E-step, the log-posterior mean is computed, and in the M-step, the log-posterior is maximized. The iterative process of these two steps is summarized as follows.

E Step: Given the posterior distribution of ${\mathbf{X}}$ and the observed matrix ${\mathbf{Y}}$, we compute the mean of the log-posterior of the hyperparameters ${\bm{\alpha}}$ with respect to the hidden variable ${\mathbf{X}}$. Let ${\bm{\theta}}=\left\{{{\bm{\alpha}},\gamma}\right\}$ and we define R function as:

$$\begin{gathered}R\left({{\bm{\theta}}|{\bm{\theta}}(t)}\right)={\mathbb{E}_{p({\mathbf{X}}|{\mathbf{Y}},{{\bm{\theta}}{\left(t\right)}})}}\left\{{\log p\left({{\bm{\theta}}|{\mathbf{X}},{\mathbf{Y}},{\bm{\theta}}(t)}\right)}\right\}\hfill\\ =R\left({{\bm{\alpha}}|{\bm{\theta}}(t)}\right)+R\left({\gamma|{\bm{\theta}}(t)}\right)+c,\hfill\\ \end{gathered}$$ (90)

where $c$ represents a constant that is independent of ${\bm{\theta}}$. Next, we calculate $R\left({{\bm{\alpha}}|{\bm{\theta}}(t)}\right)$ and $R\left({\gamma|{\bm{\theta}}(t)}\right)$ as follow.

$$\begin{gathered}R\left({{\bm{\alpha}}|{\bm{\theta}}(t)}\right)={\mathbb{E}_{p({\mathbf{X}}|{\mathbf{Y}},{\bm{\theta}}(t))}}\left\{{\log p\left({{\mathbf{X}}|{\bm{\alpha}}}\right)+\log p\left({\bm{\alpha}}\right)}\right\}\hfill\\ =\sum\nolimits_{i,j}{}2\ln({\alpha_{i,j}}+\beta{\alpha_{i-1,j}}+\beta{\alpha_{i+1,j}}+\beta{\alpha_{i,j-1}}+\beta{\alpha_{i,j+1}})\hfill\\ -({\alpha_{i,j}}+\beta{\alpha_{i-1,j}}+\beta{\alpha_{i+1,j}}+\beta{\alpha_{i,j-1}}+\beta{\alpha_{i,j+1}})\hfill\\ \left\langle{\left|{{x_{i,j}}(t)}\right|+\left|{{x_{i+N,j}}(t)}\right|}\right\rangle+a\ln{\alpha_{i,j}}-b{\alpha_{i,j}},\hfill\\ \end{gathered}$$ (91)

$$\begin{gathered}R\left({\gamma|{\mathbf{\theta}}(t)}\right)={E_{p({\mathbf{X}}|{\mathbf{Y}},{\mathbf{\theta}}(t))}}\left\{{\log p\left({{\mathbf{Y}}|{\mathbf{Z}},\gamma}\right)}\right\}\hfill\\ =-MN\ln\gamma-\frac{{\left\|{{\mathbf{Y}}-{\mathbf{\hat{Z}}}\left(t\right)}\right\|_{F}^{2}+\sum\nolimits_{i,j}{u_{i,j}^{z}(t)}}}{{2\gamma}}.\hfill\\ \end{gathered}$$ (92)

$\left\langle{\left|{{x_{i,j}}(t)}\right|}\right\rangle$ represents the mean of the absolute value of ${x_{i,j}}(t)$, that is

$$\begin{gathered}\left\langle{\left|{{x_{i,j}}(t)}\right|}\right\rangle=\int{\left|{{x_{i,j}}(t)}\right|p({x_{i,j}}(t)|{\mathbf{Y}},{\mathbf{\theta}}(t))d}{x_{i,j}}(t)\hfill\\ =\frac{1}{{{\psi_{i,j}}(t)}}\left[{e^{-\xi_{i,j}^{+}(t)}}\sqrt{2\pi u_{i,j}^{r}(t)}\varphi_{i,j}^{+}(t)Q\left({\frac{{\varphi_{i,j}^{+}(t)}}{{\sqrt{u_{i,j}^{r}(t)}}}}\right)-\right.\hfill\\ \left.{e^{-\xi_{i,j}^{-}(t)}}\sqrt{2\pi u_{i,j}^{r}(t)}\varphi_{i,j}^{-}(t)Q\left({\frac{{\varphi_{i,j}^{-}(t)}}{{\sqrt{u_{i,j}^{r}(t)}}}}\right)+2u_{i,j}^{r}(t){e^{-\frac{{{{\left({{{\hat{r}}_{i,j}}(t)}\right)}^{2}}}}{{2u_{i,j}^{r}(t)}}}}\right].\hfill\\ \end{gathered}$$ (93)

M Step: We update the hyperparameters ${\bm{\alpha}}$ and $\gamma$ by maximizing R funcation:

$${\bm{\alpha}}(t+1)=\mathop{\arg\max}\limits_{\bm{\alpha}}R\left({{\bm{\alpha}}|{\bm{\theta}}(t)}\right),$$ (94)

$$\gamma(t+1)=\mathop{\arg\max}\limits_{\gamma}R\left({\gamma|{\bm{\theta}}(t)}\right).$$ (95)

First, we consider $\bm{\alpha}$. Unlike conventional SBL, in PCSBL, the hyperparameters are interdependent, meaning that the element-wise estimation of parameters cannot be performed independently. Directly solving the result of equation (94) is challenging. To address this, we refer to the derivation process in [29] and consider an alternative suboptimal solution that achieves good estimation accuracy while simplifying the computation process. Assuming ${{\bm{\alpha}}^{*}}$ is the optimal solution to equation (94), the first-order derivative of R funcation with respect to ${\bm{\alpha}}$ equals zero at ${{\bm{\alpha}}^{*}}$. That is, for any $i$, $j$, the following condition holds:

$$\begin{gathered}\frac{{\partial Q\left({{\mathbf{\alpha}}|{{\mathbf{\theta}}^{(t)}}}\right)}}{{\partial{\alpha_{i,j}}}}{|_{{\mathbf{\alpha}}={{\mathbf{\alpha}}^{*}}}}=\frac{a}{{\alpha_{i,j}^{*}}}+2\left({\upsilon_{i,j}}+\beta{\upsilon_{i-1,j}}+\beta{\upsilon_{i+1,j}}+\right.\hfill\\ \left.\beta{\upsilon_{i,j-1}}+\beta{\upsilon_{i,j+1}}\right)-b-{\omega_{i,j}}(t)-{\omega_{N+i,j}}(t)=0,\hfill\\ \end{gathered}$$ (96)

where

$$\begin{gathered}{\omega_{i,j}}(t)\buildrel\Delta\over{=}\left\langle{\left|{{x_{i,j}}(t)}\right|}\right\rangle+\beta\left\langle{\left|{{x_{i-1,j}}(t)}\right|}\right\rangle\hfill\\ +\beta\left\langle{\left|{{x_{i+1,j}}(t)}\right|}\right\rangle+\beta\left\langle{\left|{{x_{i,j-1}}(t)}\right|}\right\rangle+\beta\left\langle{\left|{{x_{i,j+1}}(t)}\right|}\right\rangle,\hfill\\ \end{gathered}$$ (97)

$${\upsilon_{i,j}}\buildrel\Delta\over{=}\frac{1}{{\alpha_{i,j}^{*}+\beta\alpha_{i-1,j}^{*}+\beta\alpha_{i+1,j}^{*}+\beta\alpha_{i,j-1}^{*}+\beta\alpha_{i,j+1}^{*}}}.$$ (98)

In our model, the parameters $\beta\geq 0$ and ${\alpha_{i,j}}\geq 0$ hold true for any $i,j$. Building on this, based on equation (98), ${\upsilon_{i,j}}$ satisfies the following inequality constraint:

$$\begin{gathered}0\leq{\upsilon_{i,j}}\leq\frac{1}{{\alpha_{i,j}^{*}}},\\ 0\leq{\upsilon_{i,j}}\leq\frac{1}{{\beta\alpha_{i-1,j}^{*}}},\\ 0\leq{\upsilon_{i,j}}\leq\frac{1}{{\beta\alpha_{i+1,j}^{*}}},\\ 0\leq{\upsilon_{i,j}}\leq\frac{1}{{\beta\alpha_{i,j-1}^{*}}},\\ 0\leq{\upsilon_{i,j}}\leq\frac{1}{{\beta\alpha_{i,j+1}^{*}}}.\\ \end{gathered}$$ (99)

Substituting the above results into equation (96), we obtain:

$$\frac{a}{{\alpha_{i,j}^{*}}}\leq b+{\omega_{i,j}}(t)+{\omega_{N+i,j}}(t)\leq\frac{{a+10}}{{\alpha_{i,j}^{*}}}.$$ (100)

Then $\alpha_{i,j}^{*}\in\left[{\frac{a}{{b+{\omega_{i,j}}(t)+{\omega_{N+i,j}}(t)}},\frac{{a+10}}{{b+{\omega_{i,j}}(t)+{\omega_{N+i,j}}(t)}}}\right]$ is held. Therefore, a simple suboptimal solution for equation (94) can be given by:

$${\alpha_{i,j}}(t+1)=\frac{a}{{b+{\omega_{i,j}}(t)+{\omega_{N+i,j}}(t)}},$$ (101)

where $a>0$ and $b={10^{-4}}$ are shape parameter and scale parameter of Gamma distribution, respectively. As analyzed in [29], typically $a$ takes much larger values than $b$. In our scenario, considering the estimated matrix as a two-dimensional block sparse matrix with a Laplacian distribution as the prior, $a$ is empirically set to 2. This value has been found through simulation validation to achieve good convergence of the algorithm.

Then we focus on noise variance $\gamma$. Suppose ${\gamma^{*}}$ is the optimal solution of equation (95), according to equation (92), ${\gamma^{*}}$ satisfies:

$$\frac{{\partial R\left({\gamma|{\mathbf{\theta}}(t)}\right)}}{{\partial\gamma}}{|_{\gamma={\gamma^{*}}}}=-\frac{{MN}}{{{\gamma^{*}}}}+\frac{{\left\|{{\mathbf{Y}}-{\mathbf{\hat{Z}}}\left(t\right)}\right\|_{F}^{2}+\sum\nolimits_{i,j}{u_{i,j}^{z}(t)}}}{{2{{\left({{\gamma^{*}}}\right)}^{2}}}}=0.$$ (102)

It is easy to obtain the expression of $\gamma(t+1)$ is:

$$\gamma(t+1)={\gamma^{*}}=\frac{{\left\|{{\mathbf{Y}}-{\mathbf{\hat{Z}}}\left(t\right)}\right\|_{F}^{2}+\sum\nolimits_{i,j}{u_{i,j}^{z}(t)}}}{{2MN}}.$$ (103)

Thus completes the update process for ${\bm{\theta}}$. The conclusions from equations (101) and (103) serve as the output of the EM algorithm, reflected in lines 14 and 15 of Algorithm 2. With this, we have completed the entire derivation process of the GAMP-PCSBL-La algorithm. In Section V, we validate that the proposed GAMP-PCSBL-La algorithm can more accurately estimate block sparse matrix with DCT sparse properties. Furthermore, compared to other algorithms, it exhibits lower computational complexity. We will employ this algorithm for rough AUD, accurate AUD, and CE.

V Simulations

Figure 6: Performance comparisons for block sparse matrix recovery algorithms.

To validate the effectiveness and superiority of the proposed scheme, we conducted numerical simulations. The specific simulation parameters are detailed in the table below. We consider the 3GPP vehicular models, namely extended vehicular A (EVA) with number of path is 9 and ${\tau_{\max}}=9\mu s$ [35]. Assuming that the channel gain for each path follows a Gaussian distribution with zero mean and a variance equal to the corresponding LSFC. The delay and Doppler parameters are randomly generated within the range of 0 to their respective maximum values. For the uniform planar antenna array at the receiver, the elevation angle ${\theta_{u,b,i,e}}$ is uniformly distributed within the range $\left[{0,\pi}\right]$, while the azimuth angle ${\theta_{u,b,i,a}}$ is uniformly distributed within $\left[{-\frac{\pi}{2},\frac{\pi}{2}}\right]$.

TABLE I: Simulation Parameters

Parameters	Value
Doppler dimension for a block $N$	128
Delay dimension for a block $M$	512
Number of paths $P$	9
Coupling factor $\beta$	0.3
Shape parameter for Gamma distribution $a$	1.5
Scale parameter for Gamma distribution $b$	1e-3
UEs’ maximum velocity $v_{\max}$	300Km/h
Maximum path delay $\tau_{\max}$	2.5us
Number of total UEs $U$	1000
Carrier frequency $f_{c}$	4GHz
Subcarrier spacing $\Delta f$	15KHz
Doppler dimension for preamble 1 $N^{\prime}$	64
Delay dimension for preamble 1 $M^{\prime}$	4
Doppler dimension for preamble 2 $K_{p}$	20
Delay dimension for preamble 2 $L_{p}$	20
LSFC $\lambda_{u,b}$ at distance $d_{u,b}$ in km	-128.1-37.6log${}_{10}(d_{u,b})$dB
Background noise power	-174dBm/Hz
Transmission power	10dBm

To evaluate the performance of the massive random access scheme, we use the detection error rate (DER) and the normalized mean squared error (NMSE) as performance metrics for AUD and CE, respectively. They are defined as follows:

$${\text{DER}}=\frac{{\left|{{{\mathcal{K}}_{a}}\backslash{{\bar{\mathcal{U}}}_{a}}}\right|+\left|{{{\bar{\mathcal{U}}}_{a}}\backslash{{\mathcal{K}}_{a}}}\right|}}{K},$$ (104)

$${\text{NMSE}}=10{\log_{10}}\frac{{\left\|{{\mathbf{\bar{H}}}^{DDA2}-{\mathbf{H}}^{DDA2}}\right\|_{F}^{2}}}{{\left\|{{\mathbf{H}}^{DDA2}}\right\|_{F}^{2}}}.$$ (105)

${\mathcal{A}}\backslash{\mathcal{B}}$ represents a set whose elements are in ${\mathcal{A}}$ but not in ${\mathcal{B}}$. $\left|{\mathcal{A}}\right|$ denotes the cardinality of set $\mathcal{A}$. $\left\|\cdot\right\|_{F}$ denotes Frobenius norm. A smaller DER or NMSE indicates more accurate detection and estimation results, corresponding to better AUD and CE performance.

The performance of the proposed GAMP-PCSBL-La algorithm are compared with other existing algorithms in block sparse matrix recovery at first. we set the dimensions of the block sparse matrix to $256\times 64$, the observation matrix to $64\times 64$, and the sensing matrix to $64\times 256$. A block sparse matrix was generated by randomly creating non-zero values and applying a DCT. The elements of the sensing matrix and noise matrix followed a Gaussian distribution. Compared algorithms including generalized approximate message passing and pattern coupling sparse Bayesian learning with Gaussian prior (GAMP-PCSBL-Gs) [30], PCSBL [28], Orthogonal Matching Pursuit (OMP) [36], Block OMP (BOMP), Turbo-Variational Bayesian Inference-Markov Random Field (Turbo-VBI-MRF) [37], and GAMP-SBL [38]. The GAMP-based algorithms avoid matrix inversion or optimization processes, resulting in lower computational complexity.

In Fig. 6 (a), with the number of non-zero blocks fixed at 5, we compared the performance of various algorithms in recovering block sparse matrix under different signal-to-noise ratios (SNRs). The simulation curves show that as the SNR increases, the NMSE performance of all algorithms improves. In Fig. 6 (b), with the SNR fixed at 12.5 dB, we analyzed the impact of varying the column dimensions of the block sparse matrix on the performance of each algorithm. It is evident that as the dimensions of the sparse matrix increase, the estimation accuracy of all algorithms gradually declines. Additionally, in Fig. 6 (c), with the SNR fixed at 12.5 dB and the sparse matrix dimensions set to $256\times 64$, we compared the performance trends of each algorithm under different numbers of non-zero blocks. This figure implies that as the number of non-zero blocks increases, the estimated accuracy decreases across all algorithms. The simulation results in Fig. 6 (b) and (c) are consistent with the relevant conclusions of compressed sensing theory. These simulation curves also demonstrate that algorithms utilizing PCSBL outperform other algorithms in block sparse matrix recovery. Moreover, the proposed GAMP-PCSBL-La algorithm outperforms the other algorithms, showcasing its unique performance advantages in recovering block sparse matrices formed through DCT.

Figure 7: Convergence trends of iterative algorithms.

In Fig. 7, we compare the convergence trends of several iterative algorithms, with the simulation settings being consistent with those in Fig. 6. The ”estimated error” in the figure is defined as the non-logarithmic form of NMSE, i.e., $10^{\text{NMSE}/10}$. The figure shows that the PCSBL and turbo-VBI-MRF algorithms, which are based on direct matrix inversion, converge faster than the GAMP-based algorithms, reaching convergence in approximately 5 iterations. In contrast, the GAMP-based algorithms converge after about 20 iterations. This simulation result demonstrates that the proposed GAMP-PCSBL-La algorithm exhibits good convergence performance and reliability.

To compare the proposed hybrid preamble scheme, we set up control groups using only the superimposed preamble and only the embedded preamble. Additionally, to avoid the complex computation of large matrix inversions, we limited our comparison to low-complexity GAMP-based algorithms to evaluate their performance in massive random access, thereby verifying the superiority of the proposed scheme. For simplicity, we set the magnitude of non-zero entries of ${\mathbf{X}}_{u,1}^{DD}$ to be ten times that of ${\mathbf{X}}_{u,2}^{DD}$. Especially, we use ‘HP’, ‘SP’ and ‘EP’ to denote hybrid preamble, superimposed preamble and embedded preamble schemes, respectively. In Fig. 8, with the number of active UEs fixed at 30, we simulated the impact of antenna array dimensions on the performance of massive random access schemes. The results indicate that as the number of antennas increases, both the hybrid preamble scheme and the superimposed preamble scheme exhibit improved AUD and CE performance. Due to the absence of a rough activity detection step to reduce the dimensionality of the matrix to be estimated, the embedded preamble scheme alone, with its excessively large block-sparse channel matrix, fails to achieve effective AUD and CE results. It is evident that in massive random access, the proposed hybrid preamble scheme significantly outperforms the schemes that utilize either the superimposed or embedded preamble alone. Additionally, the simulation curves demonstrate that, compared to the GAMP-PCSBL-Gs and GAMP-SBL algorithms, the proposed GAMP-PCSBL-La algorithm more effectively captures the block-sparsity caused by fractional channel parameters, resulting in superior AUD and CE performance.

Figure 8: Performance comparisons for massive random access schemes versus dimensions of antenna array: (a) DER; (b) NMSE.

In Fig. 9, with the number of active UEs fixed at 30 and the antenna array dimension set at $8\times 8$, we experimented with different dimensions of preamble 2 sequence, ranging from 12 12 to 24 24. The simulation curves indicate that, due to the increased dimensionality of the received signals, a larger preamble sequence dimension in the hybrid preamble scheme leads to more accurate AUD and CE performance. The performance of the superimposed preamble scheme remains unaffected by the preamble 2 sequence dimension. Similarly, although the embedded preamble scheme shows limited improvement as the preamble2 dimension increases, it still lacks reliability. Consistent with previous simulation results, among the three algorithms tested, the GAMP-PCSBL-La algorithm exhibits superior performance compared to GAMP-PCSBL-Gs and GAMP-SBL.

Figure 9: Performance comparisons for massive random access schemes versus dimensions of embedded preamble sequence: (a) DER; (b) NMSE.

Figure 10: Performance comparisons for massive random access schemes versus number of active UEs: (a) DER; (b) NMSE.

Similar conclusions can be drawn from the simulation curves in Fig. 10. With the antenna array dimension fixed at $8\times 8$, we varied the number of active UEs from 10 to 40. The simulation results show that as the number of active UEs increases, the performance of all the massive random access schemes declines. This decline is attributed to the fact that more active UEs correspond to more non-zero elements, making the matrix less sparse. The proposed hybrid preamble scheme achieves significantly better AUD and CE performance compared to the other two preamble schemes when addressing the demands of massive random access. The larger dimension of the channel matrix to be estimated causes the embedded preamble scheme to fail when used alone. Moreover, the simulation curves in Figure 3 demonstrate that the proposed GAMP-PCSBL-La algorithm outperforms other iterative algorithms in block-sparse matrix recovery.

VI Conclusion

This paper proposes a hybrid preamble scheme for massive machine-type random access in high-mobility scenarios within cell-free massive MIMO systems using OTFS modulation. This scheme employs a superimposed preamble for rough AUD and then performs accurate AUD and CE based on the rough detected UE set and embedded preamble. By leveraging the advantages of both preamble schemes, the proposed hybrid preamble scheme achieves more precise detection and estimation with reduced preamble overhead. Additionally, a GAMP-PCSBL-La algorithm is introduced to estimate the channel matrix, effectively capturing the block-sparse characteristics of the channel caused by fractional channel parameters, while maintaining low computational complexity. Simulation results demonstrate that the proposed hybrid preamble scheme better meets the requirements for massive random access in cell-free massive MIMO systems, and that the GAMP-PCSBL-La algorithm is particularly well-suited for this scheme.

Appendix A

By substituting equations (1) and (8), (9) into equation (7), we obtain:

$$\begin{gathered}Y\left[{n,m}\right]=\frac{1}{T}\sum\limits_{i}{}{h_{i}}\sum\limits_{m^{\prime}}{}\frac{1}{{\sqrt{NM}}}\sum\limits_{k^{\prime}}{}\sum\limits_{l^{\prime}}{}{X^{DD}}\left[{k^{\prime},l^{\prime}}\right]\hfill\\ {e^{-j2\pi\left({\frac{{m^{\prime}l^{\prime}}}{M}-\frac{{nk^{\prime}}}{N}}\right)}}{e^{-j2\pi\frac{{m^{\prime}{l_{i}}}}{M}}}{e^{-j2\pi\frac{{{l_{i}}\left({{k_{i}}+{{\tilde{k}}_{i}}}\right)}}{{NM}}}}{e^{j2\pi\frac{{n({k_{i}}+{{\tilde{k}}_{i}})}}{N}}}\hfill\\ \sum\limits_{p={l_{i}}}^{M}{}\frac{1}{{M\Delta f}}{e^{-j2\pi\frac{p}{M}(m-m^{\prime}-\frac{{{k_{i}}+{{\tilde{k}}_{i}}}}{N})}}+\frac{1}{T}\sum\limits_{i}{}{h_{i}}\sum\limits_{m^{\prime}}{}\frac{1}{{\sqrt{NM}}}\sum\limits_{k^{\prime}}{}\hfill\\ \sum\limits_{l^{\prime}}{}{X^{DD}}\left[{k^{\prime},l^{\prime}}\right]{e^{-j2\pi\left({\frac{{m^{\prime}l^{\prime}}}{M}-\frac{{nk^{\prime}}}{N}}\right)}}{e^{-j2\pi\frac{{m^{\prime}{l_{i}}}}{M}}}{e^{-j2\pi\frac{{{l_{i}}\left({{k_{i}}+{{\tilde{k}}_{i}}}\right)}}{{NM}}}}\hfill\\ {e^{j2\pi\frac{{\left({n-1}\right)({k_{i}}+{{\tilde{k}}_{i}})}}{N}}}\sum\limits_{p=0}^{{l_{i}}}{}\frac{1}{{M\Delta f}}{e^{-j2\pi\frac{p}{M}(m-m^{\prime}-\frac{{{k_{i}}+{{\tilde{k}}_{i}}}}{N})}}+N[n,m].\hfill\\ \end{gathered}$$ (106)

Furthermore, substituting equation (10) into the above equation then we get:

$${Y^{DD}}\left[{k,l}\right]=\sum\limits_{n}{}\sum\limits_{m}{}\sum\limits_{i}{}{h_{i}}\sum\limits_{m^{\prime}}{}\frac{1}{{NM}}\sum\limits_{k^{\prime}}{}{\sum\limits_{l^{\prime}}X^{DD}}\left[{k^{\prime},l^{\prime}}\right]$$

$${e^{-j2\pi\left({\frac{{m^{\prime}l^{\prime}}}{M}-\frac{{nk^{\prime}}}{N}}\right)}}{e^{-j2\pi\frac{{m^{\prime}{l_{i}}}}{M}}}{e^{-j2\pi\frac{{{l_{i}}\left({{k_{i}}+{{\tilde{k}}_{i}}}\right)}}{{NM}}}}{e^{j2\pi\frac{{n({k_{i}}+{{\tilde{k}}_{i}})}}{N}}}\sum\limits_{p={l_{i}}}^{M}{}$$

$$\frac{1}{M}{e^{-j2\pi\frac{p}{M}(m-m^{\prime}-\frac{{{k_{i}}+{{\tilde{k}}_{i}}}}{N})}}{e^{j2\pi\left({\frac{{ml}}{M}-\frac{{nk}}{N}}\right)}}+\sum\limits_{n}{}\sum\limits_{m}{}\sum\limits_{i}{}{h_{i}}\sum\limits_{m^{\prime}}{}$$

$$\frac{1}{{NM}}\sum\limits_{k^{\prime}}{}{\sum\limits_{l^{\prime}}X^{DD}}\left[{k^{\prime},l^{\prime}}\right]{e^{-j2\pi\left({\frac{{m^{\prime}l^{\prime}}}{M}-\frac{{\left({n-1}\right)k^{\prime}}}{N}}\right)}}{e^{-j2\pi\frac{{m^{\prime}{l_{i}}}}{M}}}$$

$${e^{-j2\pi\frac{{{l_{i}}\left({{k_{i}}+{{\tilde{k}}_{i}}}\right)}}{{NM}}}}{e^{j2\pi\frac{{n({k_{i}}+{{\tilde{k}}_{i}})}}{N}}}\sum\limits_{p=0}^{{l_{i}}}{}\frac{1}{M}{e^{-j2\pi\frac{p}{M}(m-m^{\prime}-\frac{{{k_{i}}+{{\tilde{k}}_{i}}}}{N})}}$$

$${e^{j2\pi\left({\frac{{ml}}{M}-\frac{{nk}}{N}}\right)}}$$

$$=\sum\limits_{i}{}{h_{i}}\sum\limits_{k^{\prime}}{}{\sum\limits_{l^{\prime}}X^{DD}}\left[{k^{\prime},l^{\prime}}\right]{e^{-j2\pi\frac{{{l_{i}}\left({{k_{i}}+{{\tilde{k}}_{i}}}\right)}}{{NM}}}}\sum\limits_{p={l_{i}}}^{M}{}{e^{j2\pi\frac{p}{M}(\frac{{{k_{i}}+{{\tilde{k}}_{i}}}}{N})}}$$

$$\sum\limits_{m^{\prime}}{}\frac{1}{M}{e^{-j2\pi\frac{{m^{\prime}}}{M}(l^{\prime}+{l_{i}}-p)}}\sum\limits_{m}{}\frac{1}{M}{e^{-j2\pi\frac{m}{M}(p-l)}}\sum\limits_{n}{}$$

$$\frac{1}{N}{e^{-j2\pi\frac{n}{N}(k-{k_{i}}-{{\tilde{k}}_{i}}-k^{\prime})}}+\sum\limits_{i}{}{h_{i}}\sum\limits_{k^{\prime}}{}{\sum\limits_{l^{\prime}}X^{DD}}\left[{k^{\prime},l^{\prime}}\right]{e^{-j2\pi\frac{{k^{\prime}}}{N}}}$$

$$\sum\limits_{p=0}^{{l_{i}}}{}{e^{j2\pi\frac{{p-{l_{i}}}}{M}(\frac{{{k_{i}}+{{\tilde{k}}_{i}}}}{N})}}\sum\limits_{m^{\prime}}{}\frac{1}{M}{e^{-j2\pi\frac{{m^{\prime}}}{M}(l^{\prime}+{l_{i}}-p)}}\sum\limits_{m}{}\frac{1}{M}{e^{-j2\pi\frac{m}{M}(p-l)}}$$

$$\sum\limits_{n}{}\frac{1}{N}{e^{-j2\pi\frac{n}{N}(k-{k_{i}}-{{\tilde{k}}_{i}}-k^{\prime})}}.$$ (107)

We define ${\prod_{N}}\left({x-a}\right)=\frac{1}{N}\sum\limits_{n=0}^{N-1}{}{e^{-j2\pi\frac{n}{N}\left({x-a}\right)}}$, thus, the above equation can be expressed as:

$$\begin{gathered}{Y^{DD}}\left[{k,l}\right]=\sum\limits_{i}{}{h_{i}}\sum\limits_{k^{\prime}}{}{\sum\limits_{l^{\prime}}X^{DD}}\left[{k^{\prime},l^{\prime}}\right]\hfill\\ \sum\limits_{p={l_{i}}}^{M}{}{e^{j2\pi\frac{{\left({p-{l_{i}}}\right)\left({{k_{i}}+{{\tilde{k}}_{i}}}\right)}}{{MN}}}}\delta\left({l^{\prime}+{l_{i}}-p}\right)\delta\left({p-l}\right)\hfill\\ {{\prod}_{N}}\left({k-{k_{i}}-{{\tilde{k}}_{i}}-k^{\prime}}\right)+\sum\limits_{i}{}{h_{i}}\sum\limits_{k^{\prime}}{}{\sum\limits_{l^{\prime}}X^{DD}}\left[{k^{\prime},l^{\prime}}\right]\hfill\\ {e^{-j2\pi\frac{{k^{\prime}}}{N}}}\sum\limits_{p=0}^{{l_{i}}-1}{}{e^{j2\pi\frac{{\left({p-{l_{i}}}\right)\left({{k_{i}}+{{\tilde{k}}_{i}}}\right)}}{{MN}}}}\delta\left({l^{\prime}+{l_{i}}-p}\right)\delta\left({p-l}\right)\hfill\\ {{\prod}_{N}}\left({k-{k_{i}}-{{\tilde{k}}_{i}}-k^{\prime}}\right).\hfill\\ \end{gathered}$$ (108)

Analyzing the equation (108), and using the properties of the Dirac delta function, it can be written as a segment function. When ${l_{i}}\leq l$, we have:

$$\begin{gathered}{Y^{DD}}\left[{k,l}\right]=\sum\limits_{i}{}{h_{i}}\sum\limits_{k^{\prime}}{}{X^{DD}}\left[{k^{\prime},l-{l_{i}}}\right]{e^{j2\pi\frac{{\left({l-{l_{i}}}\right)\left({{k_{i}}+{{\tilde{k}}_{i}}}\right)}}{{NM}}}}\hfill\\ {{\prod}_{N}}\left({k-{k_{i}}-{{\tilde{k}}_{i}}-k^{\prime}}\right)\hfill\\ \mathop{\approx}\limits^{a}\sum\limits_{i}{}{h_{i}}\sum\limits_{k^{\prime\prime}}{}{X^{DD}}\left[{k-k^{\prime\prime},l-{l_{i}}}\right]{e^{j2\pi\frac{{\left({l-{l_{i}}}\right)\left({{k_{i}}+{{\tilde{k}}_{i}}}\right)}}{{NM}}}}\hfill\\ \frac{1}{N}\frac{{1-{e^{j2\pi{{\tilde{k}}_{i}}}}}}{{1-{e^{-j2\pi\frac{{k^{\prime\prime}-{k_{i}}-{{\tilde{k}}_{i}}}}{N}}}}}\hfill\\ \end{gathered}$$ (109)

Here, the approximate equality $a$ retains only the $2\varepsilon+1$ integer points near the extremum to approximate ${\prod_{N}}\left({k-{k_{i}}-{{\tilde{k}}_{i}}-k^{\prime}}\right)$, with $k^{\prime\prime}\in\left[{{k_{i}}-\varepsilon,{k_{i}}+\varepsilon}\right]$. Similarly, when ${l_{i}}>l$, we have

$${Y^{DD}}\left[{k,l}\right]=\sum\limits_{i}{}{h_{i}}\sum\limits_{k^{\prime}}{}{X^{DD}}\left[{k^{\prime},l-{l_{i}}}\right]{e^{-j2\pi\frac{{k^{\prime}}}{N}}}$$

$${e^{j2\pi\frac{{\left({l-{l_{i}}}\right)\left({{k_{i}}+{{\tilde{k}}_{i}}}\right)}}{{NM}}}}{{\prod}_{N}}\left({k-{k_{i}}-{{\tilde{k}}_{i}}-k^{\prime}}\right)$$

$$\approx\sum\limits_{i}{}{h_{i}}\sum\limits_{k^{\prime\prime}}{}{X^{DD}}\left[{k-k^{\prime\prime},l-{l_{i}}}\right]{e^{-j2\pi\frac{{k-k^{\prime\prime}}}{N}}}$$

$${e^{j2\pi\frac{{\left({l-{l_{i}}}\right)\left({{k_{i}}+{{\tilde{k}}_{i}}}\right)}}{{NM}}}}\frac{1}{N}\frac{{1-{e^{j2\pi{{\tilde{k}}_{i}}}}}}{{1-{e^{-j2\pi\frac{{k^{\prime\prime}-{k_{i}}-{{\tilde{k}}_{i}}}}{N}}}}}.$$ (110)

Combining the above derivations, we obtain equation (11).

Appendix B

In the case where $M$ is very small, the received signal in the time-frequency domain can be expressed as:

$$\begin{gathered}Y\left[{n,m}\right]=\frac{1}{T}\sum\limits_{i}{}{h_{i}}\sum\limits_{m^{\prime}}{}\frac{1}{{\sqrt{NM}}}\sum\limits_{k^{\prime}}{}\sum\limits_{l^{\prime}}{}{X^{DD}}\left[{k^{\prime},l^{\prime}}\right]\hfill\\ {e^{-j2\pi\left({\frac{{m^{\prime}l^{\prime}}}{M}-\frac{{nk^{\prime}}}{N}}\right)}}{e^{-j2\pi\frac{{m^{\prime}{{\tilde{l}}_{i}}}}{M}}}{e^{-j2\pi\frac{{{{\tilde{l}}_{i}}\left({{k_{i}}+{{\tilde{k}}_{i}}}\right)}}{{NM}}}}{e^{j2\pi\frac{{n({k_{i}}+{{\tilde{k}}_{i}})}}{N}}}\hfill\\ \left({\sum\limits_{p=0}^{M}{}\frac{1}{{M\Delta f}}{e^{-j2\pi\frac{p}{M}(m-m^{\prime}-\frac{{{k_{i}}+{{\tilde{k}}_{i}}}}{N})}}-\int_{0}^{{\tau_{i}}}{}{e^{-j2\pi\Delta ft(m-m^{\prime}-\frac{{{\nu_{i}}}}{{\Delta f}})}}dt}\right)\hfill\\ \mathop{\approx}\limits^{a}\frac{1}{T}\sum\limits_{i}{}{h_{i}}\sum\limits_{m^{\prime}}{}\frac{1}{{\sqrt{NM}}}\sum\limits_{k^{\prime}}{}\sum\limits_{l^{\prime}}{}{X^{DD}}\left[{k^{\prime},l^{\prime}}\right]{e^{-j2\pi\left({\frac{{m^{\prime}l^{\prime}}}{M}-\frac{{nk^{\prime}}}{N}}\right)}}\hfill\\ {e^{-j2\pi\frac{{m^{\prime}{{\tilde{l}}_{i}}}}{M}}}{e^{-j2\pi\frac{{{{\tilde{l}}_{i}}\left({{k_{i}}+{{\tilde{k}}_{i}}}\right)}}{{NM}}}}{e^{j2\pi\frac{{n({k_{i}}+{{\tilde{k}}_{i}})}}{N}}}\sum\limits_{p=0}^{M}{}\frac{1}{{M\Delta f}}\hfill\\ {e^{-j2\pi\frac{p}{M}(m-m^{\prime}-\frac{{{k_{i}}+{{\tilde{k}}_{i}}}}{N})}}.\hfill\\ \end{gathered}$$ (111)

In our system, assuming that the delay parameter is much smaller than the duration of one symbol, we can roughly establish the approximation $a$ in the above equation. Similar to the derivation in Appendix A, we can substitute equations (1), (8), (9), and (10) to obtain:

$${Y^{DD}}\left[{k,l}\right]\approx\sum\limits_{n}{}\sum\limits_{m}{}\sum\limits_{i}{}{h_{i}}\sum\limits_{m^{\prime}}{}\frac{1}{{NM}}\sum\limits_{k^{\prime}}{}{\sum\limits_{l^{\prime}}X^{DD}}\left[{k^{\prime},l^{\prime}}\right]$$

$${e^{-j2\pi\left({\frac{{m^{\prime}l^{\prime}}}{M}-\frac{{nk^{\prime}}}{N}}\right)}}{e^{-j2\pi\frac{{m^{\prime}{{\tilde{l}}_{i}}}}{M}}}{e^{-j2\pi\frac{{{{\tilde{l}}_{i}}\left({{k_{i}}+{{\tilde{k}}_{i}}}\right)}}{{NM}}}}{e^{j2\pi\frac{{n({k_{i}}+{{\tilde{k}}_{i}})}}{N}}}$$

$$\sum\limits_{p={l_{i}}}^{M}{}\frac{1}{M}{e^{j2\pi\left({\frac{{ml}}{M}-\frac{{nk}}{N}}\right)}}{e^{-j2\pi\frac{p}{M}(m-m^{\prime}-\frac{{{k_{i}}+{{\tilde{k}}_{i}}}}{N})}}$$

$$=\sum\limits_{i}{}{h_{i}}\sum\limits_{k^{\prime}}{}{\sum\limits_{l^{\prime}}X^{DD}}\left[{k^{\prime},l^{\prime}}\right]{e^{-j2\pi\frac{{{{\tilde{l}}_{i}}\left({{k_{i}}+{{\tilde{k}}_{i}}}\right)}}{{NM}}}}\sum\limits_{p={l_{i}}}^{M}{}{e^{j2\pi\frac{p}{M}(\frac{{{k_{i}}+{{\tilde{k}}_{i}}}}{N})}}$$

$$\sum\limits_{m^{\prime}}{}\frac{1}{M}{e^{-j2\pi\frac{{m^{\prime}}}{M}(l^{\prime}+{{\tilde{l}}_{i}}-p)}}\sum\limits_{m}{}\frac{1}{M}{e^{-j2\pi\frac{m}{M}(p-l)}}$$

$$\sum\limits_{n}{}\frac{1}{N}{e^{-j2\pi\frac{n}{N}(k-{k_{i}}-{{\tilde{k}}_{i}}-k^{\prime})}}$$

$$=\sum\limits_{i}{}{h_{i}}\sum\limits_{k^{\prime}}{}{\sum\limits_{l^{\prime}}X^{DD}}\left[{k^{\prime},l^{\prime}}\right]{e^{-j2\pi\frac{{{{\tilde{l}}_{i}}\left({{k_{i}}+{{\tilde{k}}_{i}}}\right)}}{{NM}}}}\sum\limits_{p={l_{i}}}^{M}{}{e^{j2\pi\frac{p}{M}(\frac{{{k_{i}}+{{\tilde{k}}_{i}}}}{N})}}$$

$${{\prod}_{M}}\left({l^{\prime}+{{\tilde{l}}_{i}}-p}\right)\delta\left({p-l}\right){{\prod}_{N}}\left({k-{k_{i}}-{{\tilde{k}}_{i}}-k^{\prime}}\right)$$

$$=\sum\limits_{i}{}{h_{i}}\sum\limits_{k^{\prime}}{}{\sum\limits_{l^{\prime}}X^{DD}}\left[{k^{\prime},l^{\prime}}\right]{e^{j2\pi\frac{{(l-{{\tilde{l}}_{i}})\left({{k_{i}}+{{\tilde{k}}_{i}}}\right)}}{{NM}}}}{{\prod}_{M}}\left({l^{\prime}+{{\tilde{l}}_{i}}-l}\right)$$

$${{\prod}_{N}}\left({k-{k_{i}}-{{\tilde{k}}_{i}}-k^{\prime}}\right)$$

$$\approx\frac{1}{{NM}}\sum\limits_{i}{}{h_{i}}\sum\limits_{k^{\prime\prime}}{}{X^{DD}}\left[{k-k^{\prime\prime},l}\right]{e^{j2\pi\frac{{\left({l-{{\tilde{l}}_{i}}}\right)\left({{k_{i}}+{{\tilde{k}}_{i}}}\right)}}{{NM}}}}$$

$$\frac{{1-{e^{j2\pi{{\tilde{k}}_{i}}}}}}{{1-{e^{-j2\pi\frac{{k^{\prime\prime}-{k_{i}}-{{\tilde{k}}_{i}}}}{N}}}}}\frac{{1-{e^{-j2\pi{{\tilde{l}}_{i}}}}}}{{1-{e^{-j2\pi\frac{{{{\tilde{l}}_{i}}}}{M}}}}}.$$

Then we can derive equation (12) based on above results.

Appendix C

For the case where $M$ is small, the derivation of equation (36) can be found in Appendix B. For the case where $M$ is large, according to equation (7), we have:

$${Y_{b}}\left[{\frac{{n^{\prime}}}{\alpha},m^{\prime}}\right]=\frac{1}{T}\sum\limits_{i}{}{h_{u,b,i}}\sum\limits_{m^{\prime\prime}}{}{\mathbf{X}}_{u,1}^{TF}\left[\frac{{n^{\prime}}}{\alpha},m^{\prime\prime}\right]$$

$${e^{-j2\pi m^{\prime\prime}\Delta f{\tau_{u,b,i}}}}{e^{-j2\pi{\nu_{u,b,i}}{\tau_{u,b,i}}}}{e^{j2\pi{\nu_{u,b,i}}\frac{{n^{\prime}}}{\alpha}T}}$$

$$\int_{{\tau_{u,b,i}}}^{T}{}{e^{-j2\pi\Delta ft(m^{\prime}-m^{\prime\prime}-\frac{{{\nu_{u,b,i}}}}{{\Delta f}})}}dt$$

$$+\frac{1}{T}\sum\limits_{i}{}{h_{u,b,i}}\sum\limits_{m^{\prime\prime}}{}{\mathbf{X}}_{u,1}^{TF}\left[\frac{{n^{\prime}}}{\alpha}-1,m^{\prime\prime}\right]{e^{-j2\pi m^{\prime\prime}\Delta f{\tau_{u,b,i}}}}$$

$${e^{-j2\pi m^{\prime\prime}\Delta fT}}{e^{-j2\pi{\nu_{u,b,i}}{\tau_{u,b,i}}}}{e^{j2\pi{\nu_{u,b,i}}\frac{{n^{\prime}}}{\alpha}T}}$$

$$\int_{0}^{{\tau_{u,b,i}}}{}{e^{-j2\pi\Delta ft(m^{\prime}-m^{\prime\prime}-\frac{{{\nu_{u,b,i}}}}{{\Delta f}})}}dt$$

$$\mathop{=}\limits^{a}\frac{1}{T}\sum\limits_{i}{}{h_{u,b,i}}\sum\limits_{m^{\prime\prime}}{}{\mathbf{X}}_{u,1}^{TF}\left[\frac{{n^{\prime}}}{\alpha},m^{\prime\prime}\right]{e^{-j2\pi m^{\prime\prime}\Delta f{\tau_{u,b,i}}}}$$

$${e^{-j2\pi{\nu_{u,b,i}}{\tau_{u,b,i}}}}{e^{j2\pi{\nu_{u,b,i}}\frac{{n^{\prime}}}{\alpha}T}}\sum\limits_{p=0}^{M-1}{\frac{1}{{M\Delta f}}{e^{-j2\pi\frac{p}{M}(m^{\prime}-m^{\prime\prime}-\frac{{{\nu_{u,b,i}}}}{{\Delta f}})}}}$$

$$+\frac{1}{T}\sum\limits_{i}{}{h_{u,b,i}}\sum\limits_{m^{\prime\prime}}{}\sum\limits_{p=0}^{{l_{i}}}{\frac{1}{{M\Delta f}}{e^{-j2\pi\frac{p}{M}(m^{\prime}-m^{\prime\prime}-\frac{{{\nu_{u,b,i}}}}{{\Delta f}})}}}$$

$${e^{-j2\pi m^{\prime\prime}\Delta f{\tau_{u,b,i}}}}{e^{-j2\pi{\nu_{u,b,i}}{\tau_{u,b,i}}}}{e^{j2\pi{\nu_{u,b,i}}\frac{{n^{\prime}}}{\alpha}T}}$$

$$\left({{\mathbf{X}}_{u,1}^{TF}\left[\frac{{n^{\prime}}}{\alpha}-1,m^{\prime\prime}\right]{e^{-j2\pi m^{\prime\prime}\Delta fT}}-{\mathbf{X}}_{u,1}^{TF}\left[\frac{{n^{\prime}}}{\alpha},m^{\prime\prime}\right]}\right).$$

Assuming the time-frequency domain symbols ${\mathbf{X}}_{u,1}^{TF}[n^{\prime},m^{\prime}]$ follow a zero-mean Gaussian distribution, according to the central limit theorem, the ratio of variances between the first and second terms on the right-hand side of the equation $a$ is $\frac{M}{{2\bar{l}}}$, where $\bar{l}$ is the expected value of the delay quantization value ${l_{u,b,i}}$. Typically, delays are assumed to be uniformly randomly distributed, so $2\bar{l}={l_{\max}}$ and ${l_{\max}}\ll M$ It can be considered that the first term on the right-hand side of equation $a$ dominates the numerical value. By placing the second term of equation $a$ into the noise, we obtain the equation (36).

Appendix D

According to equations (78)-(86), the posterior mean of ${x_{i,j}}$ can be expressed as:

$$\begin{gathered}{{\hat{x}}_{i,j}}(t+1)=\int{{x_{i,j}}p\left({{x_{i,j}}|{\mathbf{Y}},{\mathbf{\tau}},\gamma}\right)d}{x_{i,j}}\hfill\\ =\frac{{\sqrt{2\pi u_{i,j}^{r}(t)}}}{{{\psi_{i,j}}(t)}}\left[{e^{-\xi_{i,j}^{+}\left(t\right)}}{\psi_{1}}\left({\varphi_{i,j}^{+}\left(t\right),u_{i,j}^{r}(t)}\right)-{e^{-\xi_{i,j}^{-}\left(t\right)}}\right.\hfill\\ \left.{\psi_{1}}\left({-\varphi_{i,j}^{-}\left(t\right),u_{i,j}^{r}(t)}\right)\right],\hfill\\ \end{gathered}$$ (112)

where

$$\begin{gathered}{\psi_{1}}\left({\varphi,u}\right)=\frac{1}{{\sqrt{2\pi u}}}\int_{0}^{+\infty}{}t\exp\left\{{-\frac{{{{\left({t-\varphi}\right)}^{2}}}}{{2u}}}\right\}dt\hfill\\ =\varphi Q\left({-\frac{\varphi}{{\sqrt{u}}}}\right)+\frac{u}{{\sqrt{2\pi u}}}\exp\left\{{-\frac{{{\varphi^{2}}}}{{2u}}}\right\}.\hfill\\ \end{gathered}$$ (113)

Then we have

$${\hat{x}_{i,j}}(t+1)=\frac{{\sqrt{2\pi u_{i,j}^{r}(t)}}}{{{\psi_{i,j}}(t)}}\left[{e^{-\xi_{i,j}^{+}\left(t\right)}}\varphi_{i,j}^{+}\left(t\right)Q\left({{{-\varphi_{i,j}^{+}\left(t\right)}}/{{\sqrt{u_{i,j}^{r}(t)}}}}\right)\right.$$

$$+{e^{-\xi_{i,j}^{-}\left(t\right)}}\varphi_{i,j}^{-}\left(t\right)Q\left({{{\varphi_{i,j}^{-}\left(t\right)}}/{{\sqrt{u_{i,j}^{r}(t)}}}}\right)$$

$$\left.+\frac{{u_{i,j}^{r}(t)}}{{\sqrt{2\pi u_{i,j}^{r}(t)}}}\left({{e^{-\xi_{i,j}^{+}\left(t\right)-\frac{{{{\left({\varphi_{i,j}^{+}\left(t\right)}\right)}^{2}}}}{{2u_{i,j}^{r}(t)}}}}-{e^{-\xi_{i,j}^{-}\left(t\right)-\frac{{{{\left({\varphi_{i,j}^{-}\left(t\right)}\right)}^{2}}}}{{2u_{i,j}^{r}(t)}}}}}\right)\right].$$ (114)

From equations (82) to (85), we can get that

$$\xi_{i,j}^{+}\left(t\right)+\frac{{{{\left({\varphi_{i,j}^{+}\left(t\right)}\right)}^{2}}}}{{2u_{i,j}^{r}(t)}}=\xi_{i,j}^{-}\left(t\right)+\frac{{{{\left({\varphi_{i,j}^{-}\left(t\right)}\right)}^{2}}}}{{2u_{i,j}^{r}(t)}}=\frac{{{{\left({{{\hat{r}}_{i,j}}(t)}\right)}^{2}}}}{{2u_{i,j}^{r}(t)}}.$$ (115)

The last two terms of equation (D) can be eliminated, resulting in equation (87). We define

$$\begin{gathered}{\chi_{i,j}}(t+1)=\int{}x_{i,j}^{2}p\left({{x_{i,j}}|{\mathbf{Y}},{\mathbf{\tau}},\gamma}\right)d{x_{i,j}}\hfill\\ =\frac{{\sqrt{2\pi u_{i,j}^{r}(t)}}}{{{\psi_{i,j}}(t)}}\left[{e^{-\xi_{i,j}^{+}\left(t\right)}}{\psi_{2}}\left({\varphi_{i,j}^{+}\left(t\right),u_{i,j}^{r}(t)}\right)+{e^{-\xi_{i,j}^{-}\left(t\right)}}\right.\hfill\\ \left.{\psi_{2}}\left({-\varphi_{i,j}^{-}\left(t\right),u_{i,j}^{r}(t)}\right)\right],\hfill\\ \end{gathered}$$ (116)

where

$${\psi_{2}}\left({\varphi,u}\right)=\frac{1}{{\sqrt{2\pi u}}}\int_{0}^{+\infty}{}{t^{2}}\exp\left\{{-\frac{{{{\left({t-\varphi}\right)}^{2}}}}{{2u}}}\right\}dt.$$ (117)

First we have

$$\begin{array}[]{*{20}{c}}{g\left(t\right)=\exp\left\{{-\frac{{{{\left({t-\varphi}\right)}^{2}}}}{{2u}}}\right\}}&\to&{g^{\prime}\left(t\right)=-\frac{{t-\varphi}}{u}g\left(t\right)}\\ {f\left(t\right)=t}&\to&{f^{\prime}\left(t\right)=1}\end{array},$$ (118)

using the fact that

$$\int_{0}^{+\infty}{f\left(t\right)g^{\prime}\left(t\right)dt}=\left.{f\left(t\right)g\left(t\right)}\right|_{0}^{+\infty}-\int_{0}^{+\infty}{f^{\prime}\left(t\right)g\left(t\right)dt},$$ (119)

and $\left.{f\left(t\right)g\left(t\right)}\right|_{0}^{+\infty}=0$ to get

$$\int_{0}^{+\infty}{}\frac{{t\left({t-\varphi}\right)}}{u}\exp\left\{{-\frac{{{{\left({t-\varphi}\right)}^{2}}}}{{2u}}}\right\}dt=\int_{0}^{+\infty}{}\exp\left\{{-\frac{{{{\left({t-\varphi}\right)}^{2}}}}{{2u}}}\right\}dt.$$ (120)

In the right-hand side of equation (120), we set $x=\left({t-\varphi}\right)/\sqrt{u}$ and substitute the definitions of ${\psi_{1}}\left({\varphi,u}\right)$ and ${\psi_{2}}\left({\varphi,u}\right)$ into the left-hand side, yielding:

$$\frac{{\sqrt{2\pi u}}}{u}{\psi_{2}}\left({\varphi,u}\right)-\frac{{\varphi\sqrt{2\pi u}}}{u}{\psi_{1}}\left({\varphi,u}\right)=\sqrt{2\pi u}Q\left({\frac{{-\gamma}}{{\sqrt{u}}}}\right).$$ (121)

Then we get

$${\psi_{2}}\left({\varphi,u}\right)=\varphi{\psi_{1}}\left({\varphi,u}\right)+uQ\left({\frac{{-\gamma}}{{\sqrt{u}}}}\right).$$ (122)

With the definition of ${\psi_{1}}\left({\varphi,u}\right)$, it can be obtained that

$$\begin{gathered}{e^{-\xi_{i,j}^{+}\left(t\right)}}{\psi_{2}}\left({\varphi_{i,j}^{+}\left(t\right),u_{i,j}^{r}(t)}\right)=\hfill\\ \left({{{\left({\varphi_{i,j}^{+}\left(t\right)}\right)}^{2}}+u_{i,j}^{r}(t)}\right){e^{-\xi_{i,j}^{+}\left(t\right)}}Q\left({\frac{{-\varphi_{i,j}^{+}\left(t\right)}}{{\sqrt{u_{i,j}^{r}(t)}}}}\right)\hfill\\ +\frac{{u_{i,j}^{r}(t)\varphi_{i,j}^{+}\left(t\right)}}{{\sqrt{2\pi u_{i,j}^{r}(t)}}}\exp\left\{{-\frac{{{{\left({{{\hat{r}}_{i,j}}(t)}\right)}^{2}}}}{{2u_{i,j}^{r}(t)}}}\right\},\hfill\\ \end{gathered}$$ (123)

$$\begin{gathered}{e^{-\xi_{i,j}^{-}\left(t\right)}}{\psi_{2}}\left({-\varphi_{i,j}^{-}\left(t\right),u_{i,j}^{r}(t)}\right)=\hfill\\ \left({{{\left({\varphi_{i,j}^{-}\left(t\right)}\right)}^{2}}+u_{i,j}^{r}(t)}\right){e^{-\xi_{i,j}^{-}\left(t\right)}}Q\left({\frac{{-\varphi_{i,j}^{-}\left(t\right)}}{{\sqrt{u_{i,j}^{r}(t)}}}}\right)\hfill\\ +\frac{{u_{i,j}^{r}(t)\varphi_{i,j}^{-}\left(t\right)}}{{\sqrt{2\pi u_{i,j}^{r}(t)}}}\exp\left\{{-\frac{{{{\left({{{\hat{r}}_{i,j}}(t)}\right)}^{2}}}}{{2u_{i,j}^{r}(t)}}}\right\}.\hfill\\ \end{gathered}$$ (124)

Combining equations (84) and (85), and substituting (123) and (124) into (116), finally using the variance definition $u_{i,j}^{x}(t+1)={\chi_{i,j}}(t+1)-{\left({{{\hat{x}}_{i,j}}(t+1)}\right)^{2}}$, we obtain the result of equation (88).