A Novel OTFS-based Massive Random Access Scheme in Cell-Free Massive MIMO Systems for High-Speed Mobility

We consider a cell-free massive MIMO system, as shown in Fig. 1, comprising $B$ APs and $U$ single-antenna UEs, which are randomly distributed over a large area. 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, denoted as $\mathcal{K_{A}}$, while the remaining UEs are silent. The channels of active UEs and APs experience doubly-selective fading. The maximum delay and Doppler shift are set to $\tau_{\max}$ and $\nu_{\max}$, respectively. The signal propagation from an active UE to an AP is characterized by a finite number of paths. The uplink transmission block consists of preamble sequences and data symbols. AUD and CE are performed based on received preamble sequences.

Consider a typical OTFS transceiver system occupies bandwidth ${M\Delta f}$ and time duration $NT$, where $M$ denotes the number of subcarriers with interval $\Delta f$ and $N$ denotes the number of time intervals $T$. 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:

On this basis, the transmitter applies the Heisenberg transform to convert it into time-domain:

where ${g_{tx}}(t)$ represents the rectangular window function in time domain with a duration of $T$. The delay-Doppler channel response model from UE $u$ to the $b$-th AP is defined as:

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}}/{P^{2}}}\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:

where ${n_{b}}\left(t\right)$ represents the additive Gaussian noise. Each AP is equipped with a UPA, where $N_{z}$ and $N_{y}$ represent the number of antennas in the $z$-direction and $y$-direction, respectively. As shown in Fig. 3, let the elevation angle and azimuth angle of the $i$-th incident signal from the $u$-th active UE to the $b$-th AP are denoted by ${\vartheta_{u,b,i}}$ and ${\varsigma_{u,b,i}}$, respectively. The antenna spacing is specified as half-wavelength, with $\rho_{u,b,i}^{y}=\cos{\vartheta_{u,b,i}}\cos{\varsigma_{u,b,i}}$ and $\rho_{u,b,i}^{z}=\sin{\vartheta_{u,b,i}}$. The response of UPA is given by:

where ${{\mathbf{a}}_{N}}\left(\rho\right)$ represents the spatial steering vector with dimension $N$, expressed as:

The combine matrix for each AP is denoted as ${\mathbf{W}}={{\mathbf{D}}_{{N_{y}}}}\otimes{{\mathbf{D}}_{{N_{z}}}}$, where ${{\mathbf{D}}_{{N}}}$ is the DFT matrix with dimension $N$. We define beam vector $\mathbf{a}_{u,b,i}=\mathbf{W}^{H}\mathbf{a}_{u,b,i}^{s}$, which has only one non-zero block. The local signal processing unit of AP performs Wigner transform on the time-domain received signal, yielding the received signal presented as

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

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:

Additionally, we introduce a function $\psi\left(x\right)$ defined as

Combining equations (1), (II-B), (8a), (8b) and (9), we obtain the received signal model in DD domain as $\mathbf{y}_{u,b}^{DD}\left[{k,l}\right]=\sum_{i=1}^{P}\mathbf{y}_{u,b,i}^{DD}\left[{k,l}\right]$, where $\mathbf{y}_{u,b,i}^{DD}\left[{k,l}\right]$ is shown in equation (11), $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.

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 for 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 for the maximum Doppler shift is ${k^{\prime}_{\max}}=\left\lceil{{\nu_{\max}}N^{\prime}T}\right\rceil$. Similar to equation (11), we obtain the reception model as:

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-{k^{\prime}_{u,b,i}}$ 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 small-value integer. 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}$, (a) holds approximately true.

For a concise expression, we define function:

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 (III-A) into matrix form:

where ${\mathbf{Y}}_{b}^{p1}=\left[\mathbf{y}_{b}^{DD}\left[{0,0}\right],...,\mathbf{y}_{b}^{DD}\left[{N^{\prime}-1,M^{\prime}-1}\right]\right]^{T}$, and

Here, ${{\mathbf{F}}_{M^{\prime},N^{\prime}}}$ is a matrix, where each element is given by $F_{M^{\prime},N^{\prime}}[k,l]=e^{j2\pi\frac{kl}{M^{\prime}N^{\prime}}}$. ${{\mathbf{p}}_{r1}}$ is the index vector that satisfies ${{\mathbf{p}}_{r1}}=[1:k^{\prime}_{\max}+\varepsilon^{\prime}+1]\cup[N^{\prime}-\varepsilon^{\prime}+1:N^{\prime}]$. ${\mathbf{H}}_{u,b}^{DD1}\in{\mathbb{C}^{({{k^{\prime}}_{\max}}+2\varepsilon^{\prime}+1)\times 1}}$ is expressed as

${{\mathbf{H}}_{u,b,i}^{DD1}}$ is presented in equation (17), where $-\varepsilon^{\prime}\leq t^{\prime}\leq\varepsilon^{\prime}$ and ${\hat{h}^{\prime}_{u,b,i}}={h_{u,b,i}}\psi_{M^{\prime}}({-\tilde{l}^{\prime}}_{u,b,i},0){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}}}}}$. Combining equations (14) and (17), ${\mathbf{A}}_{u}^{p1}={{\mathbf{X}}_{u}^{p1}\odot{\mathbf{\Phi^{\prime}}}}$ is regarded as the known measurement matrix at the AP, ${\mathbf{Y}}_{b}^{p1}$ as the observed matrix and ${\mathbf{H}}_{u,b}^{DD1}$ as an unknown 2-D block sparse matrix. In a multi-user scenario, the dimension of the sparse matrix expands, allowing us to utilize this model to detect the indices of non-zero entries in the sparse matrix 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 following subsections.

Given a sufficiently large dimensions of transmission block, the embedded preamble can be adopted for joint accurate AUD and CE. Let ${k_{\max}}=\left\lfloor{{\nu_{\max}}NT}\right\rfloor$ and ${l_{\max}}={\tau_{\max}}M\Delta f$. It can be observed from (11) that the $(k,l)$-th DD domain received symbol is affected by the transmitted symbols 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. 4.

Assuming the starting coordinates of the preamble are $(k_{p},l_{p})$, with the 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}+L_{p}<M$, for $l\in\left[{{l_{p}}-{l_{\max}},{l_{p}}+{L_{p}}}\right]$, only the case when $l\geq{l_{u,b,i}}$ (since $l\geq{l_{\max}}\geq{l_{u,b,i}}$) in (11) are considerable. We denote $M_{p}=L_{p}+l_{\max}$, $N_{p}=K_{p}+k_{\max}$, and let ${{\mathbf{X}}_{u,p}}=({{\mathbf{X}}_{u}})_{{k_{p}}:{k_{p}}+{N_{p}}+\varepsilon-1,{l_{p}}:{l_{p}}+M_{p}-1}$, ${\mathbf{Y}}_{u,b}^{p2}=({{\mathbf{Y}}_{u,b}})_{\text{vec}({k_{p}}:{k_{p}}+{N_{p}}-1,{l_{p}}:{l_{p}}+M_{p}-1),:}$, where $\mathbf{X}_{u}\in\mathbb{C}^{N\times M}$ and $\mathbf{Y}_{u,b}\in\mathbb{C}^{NM\times N_{z}N_{y}}$ are DD domain transmitted symbols of $u$-th UE and received symbols of $b$-th AP, respectively. $\text{vec}({k_{p}}:{k_{p}}+{N_{p}}-1,{l_{p}}:{l_{p}}+M_{p}-1)$ denotes the vectorized matrix indices $[{k_{p}}:{k_{p}}+{N_{p}}-1,{l_{p}}:{l_{p}}+M_{p}-1]$. Similar to (14), the received embedded preamble signal is expressed as:

where ${\mathbf{X}}_{u}^{p2}$ is described in equation (19), and

${{\mathbf{p}}_{r2}}$ is the index vector that satisfies ${{\mathbf{p}}_{r2}}=[1:k_{\max}+\varepsilon+1]\cup[N_{p}-\varepsilon+1:N_{p}]$. ${\mathbf{H}}_{u,b}^{DD2}\in{\mathbb{C}^{({k_{\max}}+2\varepsilon+1)({l_{\max}}+1)\times N_{z}N_{y}}}$ is expressed as

${{\mathbf{H}}_{u,b,i}^{DD2}}$ is presented as in equation (22), where $-\varepsilon\leq t^{\prime}\leq\varepsilon$ and ${\hat{h}_{u,b,i}}={h_{u,b,i}}{e^{-j2\pi\frac{{{{l}_{u,b,i}}\left({{{k}_{u,b,i}}+{{\tilde{k}}_{u,b,i}}}\right)}}{{NM}}}}$. According to equation (18), ${\mathbf{A}}_{u}^{p2}={{\mathbf{X}}_{u}^{p2}\odot{\mathbf{\Phi}}}$ can be considered as the known measurement matrix at the AP, ${\mathbf{Y}}_{u,b}^{p2}$ as the observed matrix and ${\mathbf{H}}_{u,b}^{DD2}$ as an unknown 2-D block sparse matrix. Since in accurate estimation, $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 estimation, the dimension of the sparse vector for accurate estimation is larger, making it more difficult to recover the sparse vector in multi-UE scenarios. Therefore, a hybrid preamble scheme is designed to achieve precise detection and estimation with lower overhead and complexity.

In the $N\times M$ DD domain, we superimpose the superimposed preamble, denoted as preamble1 ${\mathbf{X}}_{u,p1}^{DD}$, and the block symbols ${\mathbf{X}}_{u,2}^{DD}$, which includes the embedded preamble denoted as preamble2 ${\mathbf{X}}_{u,p2}$, and data symbols ${\mathbf{X}}_{u,d}$. Different power levels are allocated to ${\mathbf{X}}_{u,p1}^{DD}$ and ${\mathbf{X}}_{u,2}^{DD}$, ensuring a significant difference in energy domain between these two types of signals. The superimposed result forms transmission block, structured as shown in the Fig. 5.

Then we have

where $\mathcal{PG}$ area represents the grids designated for placing preamble2 and the guard intervals. We set the preamble1 is placed at intervals of $\frac{1}{\beta}$ along the delay axis while being placed continuously along the Doppler axis. 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 guard interval. This arrangement ensures that the received signal of preamble1 does not interfere with preamble2.

Building on this,the received signal in the TF domain can be expressed as:

where ${\mathbf{Y}^{TF}_{b1}}$ and ${\mathbf{Y}^{TF}_{b2}}$ represent $b$-th AP’s received ${\mathbf{X}}_{1}^{DD}$ and ${\mathbf{X}}_{2}^{DD}$ signals in TF domain, respectively. $\mathbf{\tilde{Z}}^{TF}_{b}={\mathbf{Y}^{TF}_{b2}}+\mathbf{N}^{TF}_{b}$ 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}\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}M^{\prime}\times N_{z}N_{y}}}$ and ${\mathbf{Y}}_{u,b,1}^{TF}\in{\mathbb{C}^{NM\times N_{z}N_{y}}}$ respectively. Based on equations (1), (9), and (22), we can derive: