Frequency Synchronization for Uplink Massive MIMO Systems

TL;DR

This paper introduces a frequency synchronization method for uplink massive MIMO OFDM systems that leverages angle information and FFT to efficiently estimate carrier frequency offsets for multiple users.

Contribution

It proposes a joint spatial-frequency alignment scheme utilizing angle information and FFT for efficient CFO estimation in massive MIMO uplink systems.

Findings

The scheme achieves accurate CFO estimation with reduced complexity.

A multi-branch beamforming design simplifies multi-user channel estimation.

Numerical results confirm the effectiveness of the proposed method.

Abstract

In this paper, we propose a frequency synchronization scheme for multiuser orthogonal frequency division multiplexing (OFDM) uplink with a large-scale uniform linear array (ULA) at base station (BS) by exploiting the angle information of users. Considering that the incident signal at BS from each user can be restricted within a certain angular spread, the proposed scheme could perform carrier frequency offset (CFO) estimation for each user individually through a \textit{joint spatial-frequency alignment} procedure and can be completed efficiently with the aided of fast Fourier transform (FFT). A multi-branch receive beamforming is further designed to yield an equivalent single user transmission model for which the conventional single-user channel estimation and data detection can be carried out. To make the study complete, the theoretical performance analysis of the CFO estimation is also conducted. We further develop a user grouping scheme to deal with the unexpected scenarios that some users may not be separated well from the spatial domain. Finally, various numerical results are provided to verify the proposed studies.

Equations

$${\mathbf{h}}_l^{(k)}=\sum _{p=1}^P{\alpha }_{l,p,k}\mathbf{a}({\theta }_{l,p,k}).$$

$$\displaystyle{\boldsymbol{a}}(\theta_{l,p,k})=\big{[}1,{\rm e}^{-{\bf j}\chi\cos\theta_{l,p,k}},\cdots,{\rm e}^{-{\bf j}\chi(M-1)\cos\theta_{l,p,k}}\big{]}^{T},$$

$$\mathbf{Y}=\sum _{k=1}^K\mathbf{E}({\varphi }_k){\mathbf{B}}_k{\mathbf{H}}_k+\mathbf{N},$$

$${\mathbf{Y}}_i^{\mathrm{d}}=\sqrt{N}\sum _{k=1}^K{\eta }_i({\varphi }_k)\mathbf{E}({\varphi }_k){\mathbf{F}}^H\text{diag}({\mathbf{s}}_i^{(k)}){\mathbf{F}}_L{\mathbf{H}}_k+{\mathbf{N}}_i^{\mathrm{d}},$$

$$\mathbf{Y}{\mathbf{a}}^{\ast }(\tilde{\theta })=\mathbf{E}({\varphi }_k){\mathbf{B}}_k{\mathbf{H}}_k{\mathbf{a}}^{\ast }(\tilde{\theta })+\underset{\text{MUI term}\simeq 0}{\underbrace{\sum _{q\ne k}\mathbf{E}({\varphi }_q){\mathbf{B}}_q{\mathbf{H}}_q{\mathbf{a}}^{\ast }(\tilde{\theta })}}+\mathbf{N}{\mathbf{a}}^{\ast }(\tilde{\theta }),$$

$$\mathbf{E}({\varphi }_k{)}^H\mathbf{Y}{\mathbf{a}}^{\ast }(\tilde{\theta })={\mathbf{B}}_k{\mathbf{H}}_k{\mathbf{a}}^{\ast }(\tilde{\theta })+\mathbf{E}({\varphi }_k{)}^H\mathbf{N}{\mathbf{a}}^{\ast }(\tilde{\theta }).$$

$$C_k({\tilde{\theta }}_k,{\tilde{\varphi }}_k)=\frac{{\int }_{{\tilde{\theta }}_k-{\theta }_{\mathrm{as}}}^{{\tilde{\theta }}_k+{\theta }_{\mathrm{as}}}\Vert {\mathbf{P}}_{{\mathbf{B}}_k}^{\perp }{\mathbf{E}}^H({\tilde{\varphi }}_k)\mathbf{Y}{\mathbf{a}}^{\ast }(\tilde{\theta }){\Vert }^{2}d\tilde{\theta }}{{\int }_{{\tilde{\theta }}_k-{\theta }_{\mathrm{as}}}^{{\tilde{\theta }}_k+{\theta }_{\mathrm{as}}}\Vert \mathbf{Y}{\mathbf{a}}^{\ast }(\tilde{\theta }){\Vert }^{2}d\tilde{\theta }},$$

$$\displaystyle\big{\{}\hat{\theta}_{k},\hat{\phi}_{k}\big{\}}=\arg\min_{\tilde{\theta}_{k},\tilde{\phi}_{k}}C_{k}(\tilde{\theta}_{k},\tilde{\phi}_{k}).$$

$$C_k({\theta }_k,{\varphi }_k)\simeq \frac{{\int }_{{\theta }_k-{\theta }_{\mathrm{as}}}^{{\theta }_k+{\theta }_{\mathrm{as}}}\Vert {\mathbf{P}}_{{\mathbf{B}}_k}^{\perp }\mathbf{E}({\varphi }_k{)}^H\mathbf{N}{\mathbf{a}}^{\ast }(\tilde{\theta }){\Vert }^{2}d\tilde{\theta }}{{\int }_{{\theta }_k-{\theta }_{\mathrm{as}}}^{{\theta }_k+{\theta }_{\mathrm{as}}}\Vert {\mathbf{B}}_k{\mathbf{H}}_k{\mathbf{a}}^{\ast }(\tilde{\theta }){\Vert }^{2}d\tilde{\theta }}\simeq \frac{2{\theta }_{\mathrm{as}}M(N-L){\sigma }_n^2}{{\int }_{{\theta }_k-{\theta }_{\mathrm{as}}}^{{\theta }_k+{\theta }_{\mathrm{as}}}\Vert {\mathbf{B}}_k{\mathbf{H}}_k{\mathbf{a}}^{\ast }(\tilde{\theta }){\Vert }^{2}d\tilde{\theta }},$$

$${\int }_{{\tilde{\theta }}_k-{\theta }_{\mathrm{as}}}^{{\tilde{\theta }}_k+{\theta }_{\mathrm{as}}}\Vert {\mathbf{P}}_{{\mathbf{B}}_k}^{\perp }\mathbf{E}(\tilde{\varphi })\mathbf{Y}{\mathbf{a}}^{\ast }(\tilde{\theta }){\Vert }^{2}d\tilde{\theta }$$

$$=$$

$$+{\int }_{{\tilde{\theta }}_k-{\theta }_{\mathrm{as}}}^{{\tilde{\theta }}_k+{\theta }_{\mathrm{as}}}\Vert {\mathbf{P}}_{{\mathbf{B}}_k}^{\perp }{\mathbf{Y}}_k^{(n)}{\mathbf{a}}^{\ast }(\tilde{\theta }){\Vert }^{2}d\tilde{\theta },$$

$$\Delta {\varphi }_k^{(n+1)}(\tilde{\theta })=\frac{\Re \left({\int }_{{\tilde{\theta }}_k-{\theta }_{\mathrm{as}}}^{{\tilde{\theta }}_k+{\theta }_{\mathrm{as}}}{\mathbf{a}}^T(\tilde{\theta })({\mathbf{Y}}_k^{(n)}{)}^H{\mathbf{P}}_{{\mathbf{B}}_k}^{\perp }{\mathbf{D}}^H{\mathbf{Y}}_k^{(n)}{\mathbf{a}}^{\ast }(\tilde{\theta })d\tilde{\theta }\right)}{{\int }_{{\tilde{\theta }}_k-{\theta }_{\mathrm{as}}}^{{\tilde{\theta }}_k+{\theta }_{\mathrm{as}}}\Vert {\mathbf{P}}_{{\mathbf{B}}_k}^{\perp }{\mathbf{D}}^H{\mathbf{Y}}_k^{(n)}{\mathbf{a}}^{\ast }(\tilde{\theta }){\Vert }^{2}d\tilde{\theta }}.$$

$${\hat{\theta }}_k^{(n+1)}=\arg \underset{{\tilde{\theta }}_k}{\min }\frac{1}{{\int }_{{\tilde{\theta }}_k-{\theta }_{\mathrm{as}}}^{{\tilde{\theta }}_k+{\theta }_{\mathrm{as}}}\Vert \mathbf{Y}{\mathbf{a}}^{\ast }(\tilde{\theta }){\Vert }^{2}d\tilde{\theta }}\times$$

$$\displaystyle\Bigg{(}\int_{\tilde{\theta}_{k}-\theta_{\rm as}}^{{\tilde{\theta}_{k}+\theta_{\rm as}}}\|{\bf P}_{{\bf B}_{k}}^{\bot}{\bf Y}^{(n)}_{k}{\boldsymbol{a}}^{*}(\tilde{\theta})\|^{2}d\tilde{\theta}-\frac{\Re\left(\int_{\tilde{\theta}_{k}-\theta_{\rm as}}^{{\tilde{\theta}_{k}+\theta_{\rm as}}}{\boldsymbol{a}}^{T}(\tilde{\theta})({\bf Y}^{(k)}_{n})^{H}{\bf P}_{{\bf B}_{k}}^{\bot}{\bf D}^{H}{\bf Y}^{(n)}_{k}{\boldsymbol{a}}^{*}(\tilde{\theta})d\tilde{\theta}\right)^{2}}{\int_{\tilde{\theta}_{k}-\theta_{\rm as}}^{{\tilde{\theta}_{k}+\theta_{\rm as}}}\|{\bf P}_{{\bf B}_{k}}^{\bot}{\bf D}^{H}{\bf Y}^{(n)}_{k}{\boldsymbol{a}}^{*}(\tilde{\theta})\|^{2}d\tilde{\theta}}\Bigg{)}.$$

$${\hat{\varphi }}_k^{(n+1)}={\hat{\varphi }}_k^{(n)}+\Delta {\varphi }_k^{(n+1)}({\hat{\theta }}_k^{(n+1)}).$$

$${\hat{\theta }}_k^{(n+1)}=$$

$$\displaystyle\kern 50.0pt+\frac{\Re\left(\int_{\tilde{\theta}_{k}-\theta_{\rm as}}^{{\tilde{\theta}_{k}+\theta_{\rm as}}}{\boldsymbol{y}}^{H}(\tilde{\theta}){\bf E}(\hat{\phi}^{(n)}_{k}){\bf P}_{{\bf B}_{k}}{\bf D}^{H}{\bf E}^{H}(\hat{\phi}^{(n)}_{k}){\boldsymbol{y}}(\tilde{\theta})d\tilde{\theta}\right)^{2}}{\int_{\tilde{\theta}_{k}-\theta_{\rm as}}^{{\tilde{\theta}_{k}+\theta_{\rm as}}}\|{\bf D}^{H}{\boldsymbol{y}}(\tilde{\theta})\|^{2}d\tilde{\theta}-\int_{\tilde{\theta}_{k}-\theta_{\rm as}}^{{\tilde{\theta}_{k}+\theta_{\rm as}}}\|{\bf P}_{{\bf B}_{k}}{\bf D}^{H}{\bf E}^{H}(\hat{\phi}^{(n)}_{k}){\boldsymbol{y}}(\tilde{\theta})\|^{2}d\tilde{\theta}}\Bigg{)}.$$

$${\hat{\theta }}_k^{(n+1)}=$$

$$\displaystyle\kern 50.0pt+\frac{\Re\left(\int_{\tilde{\theta}_{k}-\theta_{\rm as}}^{{\tilde{\theta}_{k}+\theta_{\rm as}}}{\boldsymbol{y}}^{H}(\tilde{\theta}){\bf E}(\hat{\phi}^{(n)}_{k}){\bf Q}_{k}{\bf Q}_{k}^{H}{\bf D}^{H}{\bf E}^{H}(\hat{\phi}^{(n)}_{k}){\boldsymbol{y}}(\tilde{\theta})d\tilde{\theta}\right)^{2}}{\int_{\tilde{\theta}_{k}-\theta_{\rm as}}^{{\tilde{\theta}_{k}+\theta_{\rm as}}}\|{\bf D}^{H}{\boldsymbol{y}}(\tilde{\theta})\|^{2}d\tilde{\theta}-\int_{\tilde{\theta}_{k}-\theta_{\rm as}}^{{\tilde{\theta}_{k}+\theta_{\rm as}}}\|{\bf Q}_{k}^{H}{\bf D}^{H}{\bf E}^{H}(\hat{\phi}^{(n)}_{k}){\boldsymbol{y}}(\tilde{\theta})\|^{2}d\tilde{\theta}}\Bigg{)}.$$

$$\displaystyle\Big{[}{\boldsymbol{y}}(\mathcal{D}_{1}),{\boldsymbol{y}}(\mathcal{D}_{2}),\cdots,{\boldsymbol{y}}(\mathcal{D}_{M_{\rm FFT}})\Big{]}$$

$${\mathcal{S}}_{{\tilde{\theta }}_k}=\left\{{\mathcal{D}}_i:|{\mathcal{D}}_i-{\tilde{\theta }}_k|<{\theta }_{\mathrm{as}}\right\}.$$

$${\hat{\theta }}_k^{(n+1)}=$$

$$\displaystyle\kern 50.0pt+\frac{\Re\left(\sum_{\tilde{\theta}\in\mathcal{S}_{\tilde{\theta}_{k}}}{\boldsymbol{y}}^{H}(\tilde{\theta}){\bf E}(\hat{\phi}^{(n)}_{k}){\bf Q}_{k}{\bf Q}_{k}^{H}{\bf D}^{H}{\bf E}^{H}(\hat{\phi}^{(n)}_{k}){\boldsymbol{y}}(\tilde{\theta})\right)^{2}}{\sum_{\tilde{\theta}\in\mathcal{S}_{\tilde{\theta}_{k}}}\|{\bf D}^{H}{\boldsymbol{y}}(\tilde{\theta})\|^{2}-\|{\bf Q}_{k}^{H}{\bf D}^{H}{\bf E}^{H}(\hat{\phi}^{(n)}_{k}){\boldsymbol{y}}(\tilde{\theta})\|^{2}}\Bigg{)}.$$

$${\mathbf{E}}^H({\hat{\varphi }}_k)\mathbf{Y}{\mathbf{a}}^{\ast }(\tilde{\theta })={\mathbf{B}}_k{\mathbf{H}}_k{\mathbf{a}}^{\ast }(\tilde{\theta }),$$

$${\eta }_i({\hat{\varphi }}_k{)}^{\ast }{\mathbf{E}}^H({\hat{\varphi }}_k){\mathbf{Y}}_i^{\mathrm{d}}{\mathbf{a}}^{\ast }(\tilde{\theta })=\sqrt{N}{\mathbf{F}}^H\text{diag}({\mathbf{s}}_i^{(k)}){\mathbf{F}}_L{\mathbf{H}}_k{\mathbf{a}}^{\ast }(\tilde{\theta }).$$

$${\mathcal{C}}_k(\tilde{\theta },{\hat{\varphi }}_k)=\frac{\Vert {\mathbf{P}}_{{\mathbf{B}}_k}{\mathbf{E}}^H({\hat{\varphi }}_k)\mathbf{Y}{\mathbf{a}}^{\ast }(\tilde{\theta }){\Vert }^2}{\Vert {\mathbf{P}}_{{\mathbf{B}}_k}^{\perp }{\mathbf{E}}^H({\hat{\varphi }}_k)\mathbf{Y}{\mathbf{a}}^{\ast }(\tilde{\theta }){\Vert }^2}$$

$${\vartheta }^{(k)}=\left\{{\mathcal{D}}_i:{\mathcal{C}}_k({\mathcal{D}}_i,{\hat{\varphi }}_k)\ge t_h\cdot L/(N-L)\right\},$$

$${\hat{\varphi }}_k=\arg \underset{\tilde{\varphi }}{\min }{G}_k(\tilde{\varphi }),$$

$$\mathrm{MSE}\{{\hat{\varphi }}_k\}={\frac{E[(\frac{\partial G_k(\tilde{\varphi })}{\partial \tilde{\varphi }}{)}^2]}{E[\frac{{\partial }^2{G}_k(\tilde{\varphi })}{\partial {\tilde{\varphi }}^2}{]}^2}|}_{\tilde{\varphi }={\varphi }_k}.$$

$${\frac{\partial G_k(\tilde{\varphi })}{\partial \tilde{\varphi }}|}_{\tilde{\varphi }={\varphi }_k}\simeq 2\Re \left({\int }_{{\theta }_k-{\theta }_{\mathrm{as}}}^{{\theta }_k+{\theta }_{\mathrm{as}}}{\mathbf{a}}^T(\tilde{\theta }){\mathbf{N}}^H{\mathbf{P}}_{{\mathbf{B}}_k}^{\perp }{\mathbf{D}}^H{\mathbf{B}}_k{\mathbf{H}}_k{\mathbf{a}}^{\ast }(\tilde{\theta })d\tilde{\theta }\right),$$

$$\displaystyle\left.\frac{\partial^{2}G_{k}(\tilde{\phi})}{\partial\tilde{\phi}^{2}}\right|_{\tilde{\phi}=\phi_{k}}\simeq 2\Big{(}\int_{\theta_{k}-\theta_{\rm as}}^{\theta_{k}+\theta_{\rm as}}{\boldsymbol{a}}^{T}(\tilde{\theta}){\bf H}_{k}^{H}{\bf B}_{k}^{H}{\bf D}{\bf P}_{{\bf B}_{k}}^{\bot}{\bf D}^{H}{\bf B}_{k}{\bf H}_{k}{\boldsymbol{a}}^{*}(\tilde{\theta})d\tilde{\theta}\Big{)}.$$

$$\displaystyle E\left[\left(\frac{\partial G_{k}(\tilde{\phi})}{\partial\tilde{\phi}}\right)^{2}\right]\Big{|}_{\tilde{\phi}=\phi_{k}}=2E\left[\left|\int_{\theta_{k}-\theta_{\rm as}}^{\theta_{k}+\theta_{\rm as}}{\boldsymbol{a}}^{T}(\tilde{\theta}){\bf N}^{H}{\bf P}_{{\bf B}_{k}}^{\bot}{\bf D}^{H}{\bf B}_{k}{\bf H}_{k}{\boldsymbol{a}}^{*}(\tilde{\theta})d\tilde{\theta}\right|^{2}\right]$$

$$=$$

Full text

I Introduction

Large-scale multiple-input multiple-output (MIMO) or “massive MIMO” systems have drawn considerable interests from both academia and industry [1, 2, 3, 4]. In a typical massive MIMO system, the base station (BS) is equipped with hundreds of antennas to simultaneously communicate with only a few number of single-antenna users at the same frequency band. Theoretically, such a massive MIMO system can almost perfectly relieve the inter-user interference in multiuser MIMO systems with simple linear transceivers. Other advantages, such as high spectral efficiency or robustness, also make massive MIMO a promising technique for the next generation wireless systems.

However, all these potential gains of massive MIMO systems rely heavily on perfect frequency synchronization among multiple users, which is rather challenging due to the coexistence of multiple CFOs at the BS. Both the multiuser CFO estimation and compensation problem are quite different from those in conventional single-user communication systems [8, 9, 10, 11, 12, 13]. In the past decade, there have been a number of works on the frequency synchronization for the conventional small-scale spatially multiplexed multiuser MIMO systems. For example, Besson and Stoica made the first attempt on multiple CFOs estimation [14] for flat fading channels. A joint CFO and channel estimation for multiuser cyclic-prefix (CP) MIMO orthogonal frequency division multiplexing (OFDM) systems was developed in [15] by exploiting the maximum likelihood (ML) criterion. Another suboptimal estimation algorithm was proposed in [16] based on long constant amplitude zero autocorrelation (CAZAC) training sequence. The optimal set of training Zadoff-Chu sequences for multiple CFOs estimation was studied in [17]. The work in [18] discretized the continuous-valued CFO parameters into a discrete set of bins and invoked the detection theory for CFO estimation.

The above works [14, 15, 16, 17, 18] for small-scale MIMO usually considered very few spatially multiplexed users. With the increase of the number of users, for example, in massive MIMO system, these methods may suffer from substantial performance degradation and may not be applied. Recently, the authors in [19] proposed an approximation to the joint maximum likelihood (ML) estimation for massive MIMO systems. However, the estimator [19] is designed for frequency flat channels and requires a multi-dimensional grid search. The work in [20] studied the impact of imperfect CFO estimation/compensation on the performance of zero-forcing (ZF) receiver for massive MIMO uplink transmissions. The authors in [21] addressed the CFO estimation for time division duplexed (TDD) massive MIMO system and proposed a simple uplink training scheme. The time-orthogonal training sequences from different users are guarded by zero interval in [21], such that the optimality criterion can be satisfied. However, the CFO estimates at BS must be fed back to multiple users over an error free control channel in [21], and thus the users can correct their oscillators accordingly before the data transmission.

The authors further proposed a spatially averaged periodogram based method for CFO estimation with a constant envelope pilot signal in [22]. The impact of frequency selectivity on the performance of CFO impaired single-carrier massive MU-MIMO uplink has been further studied in [23]. More recently, the authors in [24] proposed a blind frequency synchronization scheme for multiuser OFDM uplink with large number of antennas at the receiver. The scheme in [24] is blind and thus has the potential benefit of saving system resources. However, it may need a large number of blocks to obtain a satisfactory estimation performance, especially in the low signal-to-noise (SNR) regime.

An interesting array signal processing aided massive MIMO transmission scheme was proposed in [25, 34], where the angle information of the user is exploited to separate users and simplify the subsequent signal processing. It is shown that users can be grouped according to their angle information to release the inter-user interference, which can be named as angle division multiple access (ADMA). However, the work mainly focused on the channel estimation but does not provide the solution for frequency synchronization.

Motivated by [25, 34], in this paper, we design a new frequency synchronization scheme for multiuser OFDM uplink with a massive uniform linear array (ULA) at BS by exploiting the angle information of users. We consider that the incident signals at BS from each user can be restricted within a certain angular spread, since the massive BS is usually elevated at a very high altitude with few surrounding scatterers [5, 6, 7, 25, 34]. With sufficient spatial dimensions, the multiuser interference (MUI) effect can be substantially mitigated in angle domain via beamforming with the steering vectors pointing to the direction of the incoming signals. The proposed scheme could perform CFO estimation for each user individually through a joint spatial-frequency alignment procedure. A multi-branch beamforming is designed for each user and yield an equivalent single user transmission model for which the conventional single-user channel estimation and data detection can be further carried out. Moreover, the proposed scheme can be applied efficiently with the aided of fast Fourier transform (FFT) operation and therefore is of low complexity. To make the study complete, we also conduct the theoretical performance analysis for CFO estimation. Considering some users may not be separated well from the spatial domain, we further develop a user grouping based scheme, in which the users with similar DOAs are grouped together and a joint CFO estimation and data detection method is designed for each group. Finally, we provide numerical results to verify the proposed studies.

The rest of this paper is organized as follows. The system model is described in Section II. The proposed frequency synchronization scheme is presented in Section III. The improved scheme with user grouping is further developed in Section IV. Simulation results are provided in Section V and conclusions are drawn in Section VI.

Notations: Superscripts $(\cdot)^{*}$, $(\cdot)^{T}$, $(\cdot)^{H}$, $(\cdot)^{\dagger}$ and $E[\cdot]$ represent conjugate, transpose, Hermitian, pseudo-inverse and expectation, respectively; ${\bf j}=\sqrt{-1}$ is the imaginary unit; $\|\cdot\|$ denotes the Frobenius norm operator; for a vector $\boldsymbol{x}$, diag$(\boldsymbol{x})$ is a diagonal matrix with main diagonal of $\boldsymbol{x}$; $\otimes$ stands for the Kronecker product; $\textrm{vec}(\cdot)$ is the vectorization operator; ${\mathbb{C}}^{m\times n}$ defines the vector space of all $m\times n$ complex matrices; ${\bf I}_{N}$ is the $N\times N$ identity matrix; ${\bf 0}$ represents an all-zero matrix with appropriate dimension; $\textrm{Tr}(\cdot)$ denotes the trace operation. Parts of this paper have been presented in a conference paper [29].

II System Model

Consider a multiuser OFDM uplink system that consists of $K$ distributed single-antenna users and one BS with $M\gg 1$ antennas in the form of ULA, as shown in Fig. 1. The total number of subcarriers is $N$ and we assume prefect time synchronization among all users for the timing being. Denote the normalized CFO between the $k\rm{th}$ user and the BS by $\phi_{k}$, which is the ratio between the real CFO and the subcarrier spacing. Let $\Delta f_{k}$ and $T_{s}$ stand for real CFO of the $k\rm{th}$ user and the sampling interval of symbol-rate, respectively. We have $2\pi\phi_{k}=2\pi NT_{s}\Delta f_{k}$ [26], which represents the CFO-induced phase shift over one OFDM block. We consider the fractional CFO, i.e., $|\phi_{k}|<0.5$, which should be sufficient in the multiuser uplink [26, 27].

The classical one-ring channel propagation model is adopted [5, 6, 7], where each user is surrounded by a ring of $P\gg 1$ local scatterers. The multi-tap channel matrix between the $k\rm{th}$ user and the BS can be modeled as an $L\times M$ matrix {\bf H}_{k}=\big{[}{\bf h}^{(k)}_{1},{\bf h}^{(k)}_{2},\cdots,{\bf h}^{(k)}_{L}\big{]}^{T}, where the $l$-$\rm{th}$ tap channel vector ${\bf h}^{(k)}_{l}\in\mathbb{C}^{M\times 1}$ is composed of $P$ rays and can be expressed as:

[TABLE]

Here, $\theta_{l,p,k}$ represents the DOA of the $p\rm{th}$ ray in the $l$-$\rm{th}$ channel tap, while $\alpha_{l,p,k}\sim\mathcal{CN}(0,\sigma_{h,l}^{2}/P)$ represents the corresponding complex gain and is independent from each other [6]. Besides, $\sigma_{h,l}^{2}$, $l=1,2,\cdots,L$, models the power delay profile (PDP) of the propagation channel. We consider $\sum_{l=1}^{L}\sigma_{h,l}^{2}=1$ such that the total channel gain of each user at one receive antenna is normalized. Moreover, ${\boldsymbol{a}}(\theta_{l,p,k})\in\mathbb{C}^{M\times 1}$ is the steering vector and has the form

[TABLE]

where $\chi=\frac{2\pi d}{\lambda}$, $d$ is the antenna spacing, $\lambda$ denotes the signal wavelength. Similar to [5, 6, 7, 25, 34], we consider that the incident rays of each user can be constrained within a certain angular spread $\theta_{\rm as}$; Namely, the incident angles of the $k\rm{th}$ user is limited in the DOA region $(\theta_{k}-\theta_{\rm as},\theta_{k}+\theta_{\rm as})$ with uniform distribution, where $\theta_{k}$ represents the mean DOA of the $k\rm{th}$ user. Generally, this assumption has been supported in the following typical scenario [34]: BS equipped with a large number of antennas is always elevated at a very high altitude, say on the top of a high building, a dedicated tower, or an unmanned aerial vehicle platform, such that there are few surrounding scatterers. Hence, the angular spread seen by BS is quite small.

Define ${\bf F}$ as the $N\times N$ normalized DFT matrix with its ($i,j$)-$\rm{th}$ entry ${\bf F}_{i,j}=\frac{1}{\sqrt{N}}{\rm e}^{-{\bf j}\frac{2\pi(i-1)(j-1)}{N}}$. Let ${\bf F}_{L}$ stand for the submatrix of ${\bf F}$ containing the first $L$ column vectors and ${\bf E}(\phi_{k})$ represent the $N\times N$ diagonal phase rotation matrix: ${\bf E}(\phi_{k})=\textrm{diag}(1,{\rm e}^{{\bf j}\frac{2\pi\phi_{k}}{N}},\cdots,{\rm e}^{{\bf j}\frac{2\pi(N-1)\phi_{k}}{N}})$.

Each user transmits a whole training block at the beginning of the uplink frame, which is then followed by a number of data blocks. The uplink communication strategy is described in Fig. 2. Denote {\bf x}_{k}=\big{[}x_{k}(0),x_{k}(1),\cdots,x_{k}(N-1)\big{]}^{T} as the frequency-domain training block transmitted from the $k\rm{th}$ user. Then, the received time-domain training signal after CP removal can be written as the following $N\times M$ matrix:

[TABLE]

where ${\bf B}_{k}=\sqrt{N}{\bf F}^{H}\textrm{diag}({\bf x}_{k}){\bf F}_{L}\in\mathbb{C}^{N\times L}$, and ${\bf N}$ denotes the corresponding $N\times M$ additive white Gaussian noise (AWGN) matrix. We assume each element of ${\bf N}$ follows i.i.d. complex Gaussian distribution with variance $\sigma_{n}^{2}$, i.e., $E[{\bf N}{\bf N}^{H}]=M\sigma_{n}^{2}{\bf I}_{N}$.

Denote ${\bf s}_{i}^{(k)}\in\mathbb{C}^{N\times 1}$ as the frequency-domain data symbol vector transmitted from the $k\rm{th}$ user in the $i$-$\rm{th}$ data block. The corresponding received signal at BS can be expressed as the following $N\times M$ matrix:

[TABLE]

where $\eta_{i}(\phi_{k})={\rm e}^{{\bf j}\frac{2\pi i(N+N_{\rm cp})\phi_{k}}{N}}$ represents the accumulative phase rotation introduced by the CFO of the $k\rm{th}$ user with $N_{\rm cp}$ being the length of CP, and ${\bf N}_{i}^{d}\in\mathbb{C}^{N\times M}$ denotes the corresponding noise matrix with $E[{\bf N}_{i}^{\rm d}({\bf N}_{i}^{\rm d})^{H}]=M\sigma_{n}^{2}{\bf I}_{N}$.

III Joint Frequency Synchronization and Beamforming

III-A CFO Estimation

Denote ${\bf P}_{{\bf B}_{k}}={\bf B}_{k}({\bf B}_{k}^{H}{\bf B}_{k})^{-1}{\bf B}_{k}^{H}$ as the orthogonal projection operator onto the space spanned by the columns of ${\bf B}_{k}$, and define ${\bf P}_{{\bf B}_{k}}^{\bot}={\bf I}_{N}-{\bf P}_{{\bf B}_{k}}$. In order to guarantee the identifiability of the CFO estimation for the $k\rm{th}$ user, we consider the training sequence is designed such that ${\bf P}_{{\bf B}_{k}}^{\bot}{\bf E}(\tilde{\phi}){\bf B}_{k}\neq{\bf 0}$ holds for any $\tilde{\phi}\neq 0$. Meanwhile, assume ${\bf P}_{{\bf B}_{k}}^{\bot}{\bf E}(\tilde{\phi}){\bf B}_{q}\neq{\bf 0}$ holds for any $\tilde{\phi}$ and $k\neq q$.

The proposed multiuser frequency synchronization scheme is motivated by the following intuitive observation. With sufficient number of receive antennas, i.e., the sufficient spatial dimension, the steering vectors pointing to distinct DOAs would become nearly orthogonal to each other, which implies that we can simply employ the steering vector pointing to the DOA of the $k\rm{th}$ user as the beamforming vector to substantially suppress the interference from the other users in angle domain [25, 34]. This will yield an approximate single-user received signal model for the $k\rm{th}$ user, where many conventional single-user frequency synchronization algorithm can be carried out.

Let us first assume that the angular spread of all simultaneous users are disjoint through some proper user schedule strategy and consider $M\to\infty$ for the first glance. It can be readily checked that the steering vectors would become orthogonal to each other for any two distinct DOAs; namely ${\boldsymbol{a}}(\theta_{1})^{H}{\boldsymbol{a}}(\theta_{2})/M=0$ holds for any $\theta_{1}\neq\theta_{2}$ [6, 25]. Hence, for the $k\rm{th}$ user, when the angular spread of users are disjoint, the beamforming vector ${\boldsymbol{a}}(\tilde{\theta})$ with spatial alignment, i.e., $\tilde{\theta}\in(\theta_{k}-\theta_{\rm as},\theta_{k}+\theta_{\rm as})$, will eliminate the signals from the other users; namely, there are ${\bf H}_{k}{\boldsymbol{a}}^{*}(\tilde{\theta})\neq{\bf 0}$ and ${\bf H}_{q}{\boldsymbol{a}}^{*}(\tilde{\theta})={\bf 0}$ for any $q\neq k$. We then obtain

[TABLE]

which gives an equivalent frequency asynchronous single-user received signal model for the $k\rm{th}$ user.

The MUI term of (5) approaching zero arises from the following observation: The $l\rm{th}$ entry of ${\bf H}_{q}{\boldsymbol{a}}^{*}(\tilde{\theta})$, $l=1,2,\cdots,L$, equals $({\bf h}^{(q)}_{l})^{T}{\boldsymbol{a}}^{*}(\tilde{\theta})=\sum_{p=1}^{P}\alpha_{l,p,q}{\boldsymbol{a}}^{T}(\theta_{l,p,q}){\boldsymbol{a}}^{*}(\tilde{\theta})$. Once the angular spreads of the $k\rm{th}$ and the $q\rm{th}$ user are disjoint, ${\boldsymbol{a}}(\tilde{\theta})$ is nearly orthogonal to any ${\boldsymbol{a}}(\theta_{l,p,q})$ with sufficient large $M$ [6, 25].

We can further perform proper CFO compensation of (5) for purpose of frequency alignment as

[TABLE]

This implies that with the above spatial and frequency alignment, an equivalent frequency synchronous single-user received signal model will be obtained at BS corresponding to the $k\rm{th}$ user. The resultant signal in (6) would exactly lie in the column space of ${\bf B}_{k}$ in the absence of noise. It can be observed that, if $\tilde{\theta}\notin(\theta_{k}-\theta_{\rm as},\theta_{k}+\theta_{\rm as})$ or if $\tilde{\phi}\neq\phi_{k}$, then the multiuser interference or frequency asynchronism would deviate the column space of ${\bf E}(\tilde{\phi})^{H}{\bf Y}{\boldsymbol{a}}^{*}(\tilde{\theta})$ from that of $\textrm{Span}({\bf B}_{k})$.

Based on the above observations, the frequency synchronization for the $k\rm{th}$ user can be designed as the following joint spatial-frequency alignment procedure. Specifically, the joint CFO and DOA estimation for the $k\rm{th}$ user can be obtained by minimizing the cost function:

[TABLE]

where the trial DOA $\tilde{\theta}_{k}$ and the trial CFO $\tilde{\phi}_{k}$ target at the spatial and frequency alignment, respectively. Mathematically, the joint CFO and DOA estimation for the $k\rm{th}$ user can be expressed as

[TABLE]

We further obtain the following important Property:

Property: At high SNR region, the cost function in (8) achieves its minimum only in the case of both spatial and frequency alignment, i.e., $\tilde{\theta}_{k}=\theta_{k}$ and $\tilde{\phi}_{k}=\phi_{k}$.

Proof:

When $\tilde{\theta}_{k}=\theta_{k}$, all DOAs from the $k\rm{th}$ user will be included in the integration of (7) and then the beamforming vectors ${\boldsymbol{a}}(\tilde{\theta})$ in (7) will eliminate the signals from all users other than the $k\rm{th}$ user. In addition, when $\tilde{\phi}_{k}=\phi_{k}$, there is ${\bf E}(\tilde{\phi}_{k})^{H}{\bf Y}{\boldsymbol{a}}^{*}(\tilde{\theta})={\bf B}_{k}{\bf H}_{k}{\boldsymbol{a}}^{*}(\tilde{\theta})+{\bf E}(\tilde{\phi}_{k})^{H}{\bf N}{\boldsymbol{a}}^{*}(\tilde{\theta})$. As a result, the cost function is expressed as

[TABLE]

which tends to zero as SNR increases. On the other hand, the spatial and frequency misalignment includes the following four cases:

In the case of $\tilde{\theta}_{k}=\theta_{k}$ and $\tilde{\phi}_{k}\neq\phi_{k}$, the numerator of the cost function $C_{k}(\theta_{k},\tilde{\phi}_{k})$ approximately equals $\int_{\theta_{k}-\theta_{\rm as}}^{{\theta_{k}+\theta_{\rm as}}}\|{\bf P}_{{\bf B}_{k}}^{\bot}{\bf E}(\phi_{k}-\tilde{\phi}_{k}){\bf B}_{k}{\bf H}_{k}{\boldsymbol{a}}^{*}(\tilde{\theta})\|^{2}d\tilde{\theta}+2\theta_{\rm as}M(N-L)\sigma_{n}^{2}$, which indicates that additional signal power will be included in the numerator of (7). Then there is $C_{k}(\theta_{k},\tilde{\phi}_{k})>C_{k}(\theta_{k},\phi_{k})$. 2. 2.

When $0<|\tilde{\theta}_{k}-\theta_{k}|<2\theta_{\rm as}$, the integration interval $(\tilde{\theta}_{k}\!-\!\theta_{\rm as},\tilde{\theta}_{k}\!+\!\theta_{\rm as})$ intersects but not exactly coincides with the DOA region of the $k\rm{th}$ user. In this case, not all of the signal components of the $k\rm{th}$ user are included in the denominator of the cost function. Thus, there is $C_{k}(\tilde{\theta}_{k},\tilde{\phi}_{k})>C_{k}(\theta_{k},\phi_{k})$. 3. 3.

When $|\tilde{\theta}_{k}-\theta_{q}|>2\theta_{\rm as}$, $\forall q=1,2,\cdots,K$, the integration interval $(\tilde{\theta}_{k}\!-\!\theta_{\rm as},\tilde{\theta}_{k}\!+\!\theta_{\rm as})$ dose not intersect with the DOA region of any user. In this case, the denominator of the cost function approximately equals $\int_{\tilde{\theta}_{k}-\theta_{\rm as}}^{{\tilde{\theta}_{k}+\theta_{\rm as}}}\|{\bf N}{\boldsymbol{a}}^{*}(\tilde{\theta})\|^{2}d\tilde{\theta}\simeq 2\theta_{\rm as}MN\sigma_{n}^{2}$ and the cost function becomes $(N-L)/N$. We could obtain $C_{k}(\tilde{\theta}_{k},\tilde{\phi}_{k})>C_{k}(\theta_{k},\phi_{k})$ at high SNR condition. 4. 4.

When the integration interval intersects with the DOA region of one user (say the $q\rm{th}$ user) other than the $k\rm{th}$ user, the numerator of the cost function $C_{k}(\tilde{\theta}_{k},\tilde{\phi}_{k})$ approximately equals $\int_{\tilde{\theta}_{k}-\theta_{\rm as}}^{{\tilde{\theta}_{k}+\theta_{\rm as}}}\|{\bf P}_{{\bf B}_{k}}^{\bot}{\bf E}(\phi_{q}-\tilde{\phi}_{k}){\bf B}_{q}{\bf H}_{q}{\boldsymbol{a}}^{*}(\tilde{\theta})\|^{2}d\tilde{\theta}+2\theta_{\rm as}M(N-L)\sigma_{n}^{2}$. In this case, the cost function dose not tend to zero as SNR increases. Hence, $C_{k}(\tilde{\theta}_{k},\tilde{\phi}_{k})>C_{k}(\theta_{k},\phi_{k})$ holds at high SNR condition.

From the above discussions, we know that at high SNR region, the cost function achieves its minimum only in the case of both spatial and frequency alignment. This completes the proof. ∎

The above Property implies that, by finding the minimum point of (8), we can obtain a valid estimation of both CFO and DOA for the $k\rm{th}$ user. Note that directly solving the minimization problem (8) requires two dimensional search with respect to both the trial CFO and DOA. A low complexity iterative solution is further proposed by using the Taylor approximation. Specifically, denote the CFO estimation for the $k\rm{th}$ user in the $n\rm{th}$ iteration by $\hat{\phi}_{k}^{(n)}$. When $\tilde{\phi}$ is in the neighborhood of $\hat{\phi}_{k}^{(n)}$, we have the Taylor approximation of ${\bf E}(\tilde{\phi})\simeq{\bf E}(\hat{\phi}^{(n)}_{k})({\bf I}_{N}+{\bf D}\Delta\phi)$, where $\Delta\phi=\tilde{\phi}-\hat{\phi}^{(n)}_{k}$ and ${\bf D}=\frac{{\bf j}2\pi}{N}\textrm{diag}(0,1,\cdots,N-1)$. We further denote ${\bf Y}^{(n)}_{k}={\bf E}^{H}(\hat{\phi}^{(n)}_{k}){\bf Y}$, which corresponds to the signal matrix after the CFO compensation by $\hat{\phi}_{k}^{(n)}$. Then, there is

[TABLE]

which is a parabola function with respect to $\Delta\phi$. The above parabola function achieves its minimum value when $\Delta\phi$ equals

[TABLE]

By substituting $\Delta\phi_{k}^{(n+1)}(\tilde{\theta})$ for $\Delta\phi$ into (10), the DOA estimation in the ($n+1$)-$\rm{th}$ iteration for the $k\rm{th}$ user can be expressed as

[TABLE]

Accordingly, the CFO estimation can be updated as

[TABLE]

Note that optimization of (III-A) can be efficiently solved by exploiting FFT operation. Specifically, let ${\boldsymbol{y}}(\tilde{\theta})={\bf Y}{\boldsymbol{a}}^{*}(\tilde{\theta})$ represent the signal after the receive beamforming with steering vector ${\boldsymbol{a}}^{*}(\tilde{\theta})$. We can rewrite (III-A) into

[TABLE]

Note that ${\bf P}_{{\bf B}_{k}}$ has a dimension of $N\times N$, and thus directly computing (III-A) may suffer from a high computational burden. To further reduce the computational complexity, we can express ${\bf P}_{{\bf B}_{k}}$ as ${\bf P}_{{\bf B}_{k}}={\bf Q}_{k}{\bf Q}_{k}^{H}$ where ${\bf Q}_{k}\in\mathbb{C}^{N\times L}$ is a unitary matrix and has the same column space of ${\bf B}_{k}$, i.e., ${\rm Span}({\bf Q}_{k})={\rm Span}({\bf B}_{k})$. We further obtain

[TABLE]

Let both the trial DOA $\tilde{\theta}_{k}$ and the DOA samples $\tilde{\theta}$ in the numerical integration be drawn from $\mathcal{D}_{\rm FFT}=\{\mathcal{D}_{1},\mathcal{D}_{2},\cdots,\mathcal{D}_{M_{\rm FFT}}\}$ with $\mathcal{D}_{i}=\arccos\left(-\frac{2\pi(i-1)}{\chi M_{\rm FFT}}\right)$, where $M_{\rm FFT}$ represents the size of FFT. There is ${\boldsymbol{a}}^{*}(\mathcal{D}_{i})=[1,{\rm e}^{-{\bf j}2\pi(i-1)\frac{1}{M_{\rm FFT}}},\cdots,$ ${\rm e}^{-{\bf j}2\pi(i-1)\frac{(M_{\rm FFT}-1)}{M_{\rm FFT}}}]^{T}$ and thus,

[TABLE]

corresponds to the FFT conversion of ${\bf Y}$ along the spatial dimension. Denote the neighbor samples of $\tilde{\theta}_{k}$ by

[TABLE]

The numerical solution of (III-A) can be then expressed as

[TABLE]

Note that the FFT operation in (16) needs to be computed only once during the iterative estimation procedure for all users. Bearing in mind that the calculation of ${\boldsymbol{y}}(\mathcal{D}_{i})$ dominates the computational burden of the grid search in (III-A), the computational complexity can be greatly reduced by exploiting FFT in (16).

III-B Multi-Branch Beamforming

After the above joint CFO and DOA estimation, all steering vectors ${\boldsymbol{a}}^{*}(\tilde{\theta})$ with $\tilde{\theta}\in(\hat{\theta}_{k}-\theta_{\rm as},\hat{\theta}_{k}+\theta_{\rm as})$ can be utilized as the receive beamforming vectors for the $k\rm{th}$ user, and the CFO estimate $\hat{\phi}_{k}$ can be used for compensation. With perfect estimation and infinite receive antennas, i.e., $\hat{\phi}_{k}=\phi_{k}$, $\hat{\theta}_{k}=\theta_{k}$, $M\to\infty$ and $\forall\tilde{\theta}\in(\theta_{k}-\theta_{\rm as},\theta_{k}+\theta_{\rm as})$, we can obtain the signals after beamforming corresponding to the training and data blocks in noise-free environment as follows:

[TABLE]

The equivalent channel response ${\bf H}_{k}{\boldsymbol{a}}^{*}(\tilde{\theta})$ can be estimated from (19) by the conventional least square (LS) algorithm, which can be further utilized for subsequent data detection in (20).

However, in practice, the receiver has limited number of antennas and finite spatial resolution, and thus the MUI term in (5) may not be negligible. Hence, not all of the angles within the DOA region $(\theta_{k}-\theta_{\rm as},\theta_{k}+\theta_{\rm as})$ are suitable as the beamforming direction for the $k\rm{th}$ user. For example, once some users are close to each other in the spatial domain, beamforming towards the angles on the edge of DOA region may suffer from some non-negligible MUI effect. Therefore, instead of directly beamforming towards all angles within the DOA region $(\hat{\theta}_{k}-\theta_{\rm as},\hat{\theta}_{k}+\theta_{\rm as})$, we propose an alternative adaptive multi-branch beamforming scheme in this subsection. Note that after beamforming with ${\boldsymbol{a}}^{*}(\tilde{\theta})$ and CFO compensation with $\hat{\phi}_{k}$ on the received signal matrix ${\bf Y}$, the magnitude of

[TABLE]

actually stands for the ratio of the expected signal power to the resultant noise power. To some extent, we may consider the magnitude (21) as the SNR for the $k\rm{th}$ user with beamforming vector of ${\boldsymbol{a}}^{*}(\tilde{\theta})$. Let us draw the trial DOAs from set $\mathcal{D}_{\rm FFT}$. When $\tilde{\theta}$ does not correspond to the $k\rm{th}$ user, the magnitude of $\mathcal{C}_{k}(\tilde{\theta},\hat{\phi}_{k})$ can be approximately given by $L/(N\!-\!L)$. Then, we can define the set of qualified DOAs for the $k\rm{th}$ user as

[TABLE]

where $t_{h}$ is a pre-set threshold that can be larger than one. Here, the elements of $\vartheta^{(k)}$ are picked out from set $\mathcal{D}_{\rm FFT}$ in order to reuse the FFT operation on the received signal (16). The basic idea of the proposed multi-branch beamforming is to use beamforming vectors towards the qualified DOAs for each user. Mathematically, the corresponding receive beamforming matrix for the $k\rm{th}$ user can be obtained as ${\bf W}^{(k)}=\left[{\boldsymbol{a}}^{*}(\vartheta^{(k)}_{1}),{\boldsymbol{a}}^{*}(\vartheta^{(k)}_{2}),\cdots,{\boldsymbol{a}}^{*}(\vartheta^{(k)}_{M_{k}})\right]\in\mathbb{C}^{M\times M_{k}}$, where $M_{k}$ is the cardinality of $\vartheta^{(k)}$ and $\vartheta^{(k)}_{i}$ is the $i\rm{th}$ element of $\vartheta^{(k)}$.

III-C Performance Analysis of CFO Estimation

In this subsection, we provide the theoretical performance analysis for the proposed CFO estimation. For simplicity, we consider the case that BS has perfect knowledge of DOA information of each user, and utilize $\left(\theta_{k}-\theta_{\rm as},\theta_{k}+\theta_{\rm as}\right)$ as intergration interval in (7). Moreover, we consider the case that all the users have enough spatial separation between each other, and thus we can ignore the MUI effect during the analysis. In this case, the CFO estimation for the $k\rm{th}$ user is equivalent to

[TABLE]

where $G_{k}(\tilde{\phi})=\int_{\theta_{k}-\theta_{\rm as}}^{\theta_{k}+\theta_{\rm as}}{\boldsymbol{a}}^{T}(\tilde{\theta}){\bf Y}^{H}{\bf E}(\tilde{\phi}){\bf P}_{{\bf B}_{k}}^{\bot}{\bf E}^{H}(\tilde{\phi}){\bf Y}{\boldsymbol{a}}^{*}(\tilde{\theta})d\tilde{\theta}$. The mean square error (MSE) of CFO estimation under the high SNR can be expressed as [24]

[TABLE]

Bearing in mind of $\frac{\partial{\bf E}(\tilde{\phi})}{\partial\tilde{\phi}}={\bf D}{\bf E}(\tilde{\phi})$, we obtain

[TABLE]

Then, there is

[TABLE]

Define ${\bf R}_{k}=\int_{\theta_{k}-\theta_{\rm as}}^{\theta_{k}+\theta_{\rm as}}\frac{1}{2\theta_{\rm as}}{\boldsymbol{a}}^{*}(\tilde{\theta}){\boldsymbol{a}}^{T}(\tilde{\theta})d\tilde{\theta}$ with its ($i,j$)-$\rm{th}$ entry being ${\bf R}_{k}(i,j)=\frac{1}{2\theta_{\rm as}}\int_{\theta_{k}-\theta_{\rm as}}^{\theta_{k}+\theta_{\rm as}}{\rm e}^{{\bf j}\chi(i-j)\cos\tilde{\theta}}d\tilde{\theta}$. We can further simplify (III-C) as

[TABLE]

where ${\boldsymbol{\Sigma}}_{h}=\textrm{diag}(\sigma_{h,1},\sigma_{h,2},\cdots,\sigma_{h,L})$ corresponds to the channel PDP. The derivation of (28) arises from the fact that E\big{[}({\bf h}_{l}^{(k)})^{*}({\bf h}_{l}^{(k)})^{T}\big{]}=\sigma_{h,l}^{2}\cdot E\big{[}\sum_{p=1}^{P}\frac{1}{P}{\boldsymbol{a}}^{*}(\theta_{l,p,k}){\boldsymbol{a}}^{T}(\theta_{l,p,k})\big{]})=\sigma_{h,l}^{2}\cdot{\bf R}_{k}.

For the second order derivative in (26), we can obtain the following approximation:

[TABLE]

By substituting (28) and (III-C) into (24), the CFO estimation MSE for the $k\rm{th}$ user can be expressed as:

[TABLE]

Let $\sigma_{s}^{2}$ denote the transmit signal power from each user. With sufficiently large $M$ and $N$, the asymptotic CFO estimation MSE for the $k\rm{th}$ user can be expressed as

[TABLE]

indicating that the estimation MSE decreases as $1/M$. In other words, for a fixed desired MSE, the required transmit power of each user decreases as $1/M$. The detailed derivation of (31) can be found in Appendix A.

IV Improved Scheme with User Grouping

Due to user mobility, there exists the possibility that some users may have very close DOAs, or even overlapped DOAs. Meanwhile, the massive BS cannot have infinite number of antennas in practice, which indicates non-ideal spatial resolution. Thus, the proposed scheme in Section III may not work properly once the signal from multiple users cannot be separated well from the spatial domain, which can substantially deteriorate the CFO estimation and the subsequent data detection. Considering the random distribution of multiple users, we here propose an improved frequency synchronization scheme with the aid of user grouping, whose details are described as follows.

Here we should note that, the proposed scheme in Section III is originally designed for users with sufficient spatial separation; Nevertheless, for the users with a certain MUI effect, the previously obtained CFO estimate and qualified DOAs can still serve as valid coarse estimation, which can be employed to facilitate the following user grouping procedure. Specifically, according to the discussions on (21), the resultant SNR with beamforming matrix ${\bf W}^{(k)}$ can be expressed as $\rho_{k}=\sum_{\tilde{\theta}\in\vartheta^{(k)}}\mathcal{C}_{k}(\tilde{\theta},\hat{\phi}_{k})$. Moreover, the total received signal power at the BS can be approximately expressed as $\hat{\sigma}_{s}^{2}=\|{\bf Y}\|^{2}$. The noise power can be estimated as $\hat{\sigma}_{n}^{2}=\min\limits_{\tilde{\theta}\in\mathcal{D}_{\rm FFT}}\|{\bf Y}{\boldsymbol{a}}^{*}(\tilde{\theta})\|^{2}$, where the trial beamforming vector ${\boldsymbol{a}}^{*}(\tilde{\theta})$ aims at suppressing the signals from all users and thus the resultant ${\bf Y}{\boldsymbol{a}}^{*}(\tilde{\theta})$ mainly corresponds to the noise term. Then, the expected average SNR for each user can be evaluated as $\rho_{\rm ex}=\frac{\hat{\sigma}_{s}^{2}/K}{\hat{\sigma}_{n}^{2}}$. Once the estimated SNR of the $k\rm{th}$ user is below a threshold, say $\rho_{k}<\rho_{\rm th}\cdot\rho_{\rm ex}$ where $\rho_{\rm th}$ is a pre-set value, there may exist non-negligible MUI from other users. In this case, we denote the $k\rm{th}$ user as a critical user; otherwise, we denote it as a non-critical user.

The grouping strategy in angle domain includes the following two steps: First, we gather multiple critical users with similar DOAs into one group; namely, in one group the distance of the qualified DOAs between any user and the rest users should be less than a certain guard interval $\mathcal{G}$. Mathematically, denote the user index set of the $g\rm{th}$ group by $\mathcal{U}^{(g)}=\{c^{(g)}_{1},\cdots,c^{(g)}_{K_{g}}\}$, where $K_{g}$ is the number of the users in this group. Then, there should be $\min\limits_{j\neq i}{\rm dist}(\vartheta^{(c^{(g)}_{i})},\vartheta^{(c^{(g)}_{j})})<\mathcal{G}$ for any $i=1,2,\cdots,K_{g}$, where ${\rm dist}(\mathcal{U}_{1},\mathcal{U}_{2})\triangleq\max|c_{1}-c_{2}|$, $\forall c_{1}\in\mathcal{U}_{1}$ and $\forall c_{2}\in\mathcal{U}_{2}$. Second, for each group, we also take in the non-critical users whose minimal distance from all critical users within this group is smaller than the guard interval $\mathcal{G}$.

If one user is not grouped with all the other users, then this user is well separated from the others in the spatial domain, and the previously designed CFO estimation, beamforming, as well as data detection can be applied immediately for this user. On the contrary, for the user grouped with some other users, the non-negligible MUI would happen among the users in the same group. In this case, the CFO estimation and data detection should be jointly carried out for the multiple users in the same group.

Therefore, we here develop a joint CFO estimation and data detection scheme for the grouped users. Without loss of generality, we consider that the first $\kappa$ users stay in one group, i.e., $\mathcal{U}^{(1)}=\{1,2,\cdots,\kappa\}$, and illustrate the the CFO estimation and data detection for this group. Similar approach can be readily extended to other groups.

IV-A CFO Estimation

Suppose the DOA range for the considered group is [$\theta_{\min},\theta_{\max}$], where $\theta_{\min}$ and $\theta_{\max}$ are the minimum and maximum qualified DOAs from all $\kappa$ users, respectively, i.e., \theta_{\min}=\min\Big{(}\bigcup\nolimits_{k=1,\cdots,\kappa}\vartheta^{(k)}\Big{)} and \theta_{\max}=\max\Big{(}\bigcup\nolimits_{k=1,\cdots,\kappa}\vartheta^{(k)}\Big{)}. Let us then pick up the following $Q=\kappa L$ uniformly distributed DOAs within $[\theta_{\min},\theta_{\max}]$ to formulate the beamforming vectors for this group: $\psi_{i}=\theta_{\min}+(i-1)\frac{\theta_{\max}-\theta_{\min}}{Q-1}$, $i=1,2,\cdots,Q$. The corresponding beamforming matrix can be obtained as {\bf W}_{\rm Grp}=\big{[}{\boldsymbol{a}}(\psi_{1}),{\boldsymbol{a}}(\psi_{2}),\cdots,{\boldsymbol{a}}(\psi_{Q})\big{]}^{*}, and the beamforming is performed as

[TABLE]

Note that the user outside this group should have enough spatial separation from the users in this group, as defined by the grouping strategy, which says that the beamforming matrix ${\bf W}_{\rm Grp}$ could substantially mitigate the MUI from the users outside the group. Then, we can approximately consider that ${\bf Y}_{\rm Grp}$ is composed of only the $\kappa$ users in this group.

Define {\mathbb{B}}({\boldsymbol{\phi}})=\Big{[}{\bf E}(\phi_{1}){\bf B}_{1},{\bf E}(\phi_{2}){\bf B}_{2},\cdots,{\bf E}(\phi_{\kappa}){\bf B}_{\kappa}\Big{]}\in\mathbb{C}^{N\times\kappa L} where ${\boldsymbol{\phi}}=[\phi_{1},\phi_{2},\cdots,\phi_{\kappa}]^{T}$. Denote ${\bf H}_{{\rm eq}}^{(k)}={\bf H}_{k}{\bf W}_{\rm Grp}\in\mathbb{C}^{L\times Q}$ and {\bf H}_{\rm eq}=\Big{[}({\bf H}_{{\rm eq}}^{(1)})^{T},({\bf H}_{{\rm eq}}^{(2)})^{T},\cdots,({\bf H}_{{\rm eq}}^{\kappa})^{T}\Big{]}^{T}\in\mathbb{C}^{\kappa L\times Q}, and rewrite ${\bf Y}_{\rm Grp}$ as

[TABLE]

where ${\bf N}_{\rm Grp}={\bf N}{\bf W}_{\rm Grp}+\sum_{q=\kappa+1}^{K}{\bf E}(\phi_{q}){\bf B}_{q}{\bf H}_{q}{\bf W}_{\rm Grp}$ denotes the corresponding noise plus the mitigated MUI terms from the users outside the group.

Let $\tilde{\boldsymbol{\phi}}=[\tilde{\phi}_{1},\tilde{\phi}_{2},\cdots,\tilde{\phi}_{\kappa}]^{T}$ represent the trial value of ${\boldsymbol{\phi}}$. According to (33), the LS joint CFO estimation of the $\kappa$ users in this group can be given by

[TABLE]

where the cost function is C_{\rm Grp}(\tilde{\boldsymbol{\phi}})=\big{\|}{\bf P}^{\bot}_{{\mathbb{B}}}(\tilde{\boldsymbol{\phi}}){\bf R}_{\rm Grp}\big{\|}^{2}, ${\bf R}_{\rm Grp}={\bf Y}_{\rm Grp}{\bf Y}_{\rm Grp}^{H}$, ${\bf P}^{\bot}_{{\mathbb{B}}}(\tilde{\boldsymbol{\phi}})={\bf I}_{N}-{\bf P}_{{\mathbb{B}}}(\tilde{\boldsymbol{\phi}})$ and ${\bf P}_{{\mathbb{B}}}(\tilde{\boldsymbol{\phi}})=\mathbb{B}({\boldsymbol{\phi}})(\mathbb{B}^{H}({\boldsymbol{\phi}})\mathbb{B}({\boldsymbol{\phi}}))^{-1}\mathbb{B}^{H}({\boldsymbol{\phi}})$ is the orthogonal projection on the column space of $\mathbb{B}({\boldsymbol{\phi}})$.

Unfortunately, directly solving (34) requires multi-dimensional search, which is computationally inefficient. Note that although the MUI effect within the group cannot be easily mitigated from the spatial domain, previous proposed CFO estimate (8), i.e., $\hat{\phi}_{k}$, may still be able to serve as a valid coarse estimation. In the following, we propose to approximate the cost function $C_{\rm Grp}(\tilde{\boldsymbol{\phi}})$ as a high dimensional parabola in the neighborhood of coarse estimation $\hat{\boldsymbol{\phi}}=[\hat{\phi}_{1},\hat{\phi}_{2},\cdots,\hat{\phi}_{\kappa}]^{T}$, and the minimization of (34) can be then solved in a more efficient way. Specifically, there is

[TABLE]

We have the approximation ${\bf E}(\tilde{\phi}_{k}){\bf B}_{k}={\bf E}(\hat{\phi}_{k}){\bf B}_{k}+\Delta{\bf B}_{k}$, where $\Delta{\bf B}_{k}={\bf D}{\bf B}_{k}\Delta\phi_{k}+\frac{1}{2}{\bf D}^{2}{\bf B}_{k}\Delta\phi_{k}^{2}$ and $\Delta\phi_{k}=\tilde{\phi}_{k}-\hat{\phi}_{k}$. We can further obtain

[TABLE]

where

[TABLE]

Performing eigenvalue decomposition on ${\bf\Pi}(\hat{\boldsymbol{\phi}})$ gives

[TABLE]

where ${\bf U}_{\rm s}\in\mathbb{C}^{N\times\kappa L}$ and ${\bf U}_{\rm n}\in\mathbb{C}^{N\times(N-\kappa L)}$ correspond to the signal and noise subspace, respectively. Then, there is ${\bf P}^{\bot}_{{\mathbb{B}}}(\hat{\boldsymbol{\phi}})={\bf U}_{\rm n}{\bf U}_{\rm n}^{H}$. Denote the $i\rm{th}$ eigenvalue (in ascending order) and the corresponding eigenvector of ${\bf\Pi}(\hat{\boldsymbol{\phi}})$ by $\lambda_{i}$ and ${\boldsymbol{e}}_{i}$, $i=1,2,\cdots,L$, respectively. Then, the $i\rm{th}$ eigenvector of ${\bf\Pi}(\tilde{\boldsymbol{\phi}})$ can be expressed as ${\boldsymbol{e}}_{i}+\Delta{\boldsymbol{e}}_{i}$, where the perturbation term is given by [28]

[TABLE]

Bearing in mind that

[TABLE]

we further obtain the approximation

[TABLE]

By substituting (37) into (39), we obtain

[TABLE]

where

[TABLE]

With (41), we may rewrite the cost function as

[TABLE]

Letting the first derivative of $C_{\rm Grp}(\tilde{\boldsymbol{\phi}})$ with respect to $\Delta\phi_{k}$, $k=1,2,\cdots,\kappa,$ be zero yields

[TABLE]

Denote

[TABLE]

Then, (47) can be expressed as

[TABLE]

which would provide a linear estimation for $\Delta\phi_{k}$. Hence, the CFO estimation for the $\kappa$ users can be updated as $\hat{\boldsymbol{\phi}}_{\rm Grp}=[\hat{\phi}_{1,{\rm Grp}},\hat{\phi}_{2,{\rm Grp}},\cdots,\hat{\phi}_{\kappa,{\rm Grp}}]^{T}$ with $\hat{\phi}_{k,{\rm Grp}}=\hat{\phi}_{k}+\Delta\phi_{k}$, $k=1,2,\cdots,\kappa$.

IV-B Data Detection

After the CFO estimation, the equivalent channel response matrix can be obtained from (33) via the LS algorithm: $\hat{\bf H}_{\rm eq}=\mathbb{B}^{\dagger}(\hat{\boldsymbol{\phi}}_{\rm Grp}){\bf Y}_{\rm Grp}$. Similar to (32), the beamforming can also be performed for the data blocks as:

[TABLE]

where ${\bf N}_{i,{\rm Grp}}^{\rm d}\in\mathbb{C}^{N\times Q}$ denotes the corresponding noise plus the interference from the users outside the group. We rewrite the equivalent channel matrix for the $k\rm{th}$ user in (51) as {\bf H}_{{\rm eq}}^{(k)}=\big{[}{\bf h}_{1,{\rm eq}}^{(k)},{\bf h}_{2,{\rm eq}}^{(k)},\cdots,{\bf h}_{Q,{\rm eq}}^{(k)}\big{]}. Let $\mathcal{H}_{q,{\rm eq}}^{(k)}\in\mathbb{C}^{N\times N}$ stand for the circular convolution matrix corresponding to the equivalent channel response ${\bf h}^{(k)}_{q,{\rm eq}}$, i.e.,

[TABLE]

The vectorization of (51) is

[TABLE]

where

[TABLE]

Then, the ZF data detection for the $\kappa$ users can be performed as $\hat{\bf s}_{i}=({\boldsymbol{\eta}}_{i}^{*}({\boldsymbol{\phi}})\otimes{\bf F})\mathcal{H}_{\rm eq}^{{\dagger}}\textrm{vec}({\bf Y}^{\rm d}_{i,{\rm Grp}})$. In practice, ${\boldsymbol{\phi}}$ can be replaced by $\hat{\boldsymbol{\phi}}_{\rm Grp}$, while $\mathcal{H}_{\rm eq}$ can be constructed from both $\hat{\bf H}_{\rm eq}$ and $\hat{\boldsymbol{\phi}}_{\rm Grp}$ in the way as (52), (57) and (58).

V Simulations

In this section, we demonstrate the effectiveness of the proposed frequency synchronization scheme through numerical examples. 16-QAM constellation is adopted for data transmission. In the proposed scheme, the training symbols are randomly drawn from the QPSK constellations. The SNR is defined as $\sigma_{s}^{2}/\sigma_{n}^{2}$. The channel length is set as $L=10$ and the length of CP is taken as $N_{\rm cp}=L\!-\!1=9$. Unless otherwise stated, the total number of subcarriers is taken as $N=64$, while we consider $M=128$ and $d=\lambda/2$. The channel vectors of different users are formulated according to (1) with uniform PDP, i.e., $\sigma_{h,l}^{2}=1/L$. Both angular spreads of $\theta_{\rm as}=5^{\circ}$ and $\theta_{\rm as}=10^{\circ}$ are considered, while the guard interval is taken as $\mathcal{G}=4\theta_{\rm as}+5^{\circ}$. Unless otherwise stated, we assume perfect knowledge about angular spread is available at BS. The FFT size is set as $M_{\rm FFT}=2M$. The normalized CFO is randomly generated from $-\phi_{\max}$ to $\phi_{\max}$. The MSE of the normalized CFO is adopted as the figure of merit. The thresholds are taken as $t_{h}=10$ and $\rho_{\rm th}=2/3$. The proposed scheme and the improved scheme with user grouping are labelled as ‘FS-BEAM’ and ‘FS-BEAM-UG’. For comparison, we also include the results of the existing robust multi-CFO estimation scheme with iterative interference cancellation of [17], referred to as ‘RMCE-IIC’. Specifically, the ZC sequence in RMCE-IIC is adopted as $\exp({\bf j}\pi i^{2}/N)$, $i=0,1,\cdots,N-1$. The circular-shift set follows the constraint of [17, eq. (42)]. We consider that three iterations of iterative interference cancellation is adopted in RMCE-IIC.

First, we consider the case of $K=4$ and plot the iterative procedure of CFO estimation in the proposed FS-BEAM under different SNR conditions in Fig. 3. The angular spread is $\theta_{\rm as}=5^{\circ}$ and the mean DOAs of the four users are fixed as $\{30^{\circ},60^{\circ},120^{\circ},150^{\circ}\}$. We consider both the cases of $\phi_{\max}=0.1$ and $\phi_{\max}=0.2$ in this example. The CFOs are initialized as all zeros, i.e., $\hat{\phi}_{k}^{(0)}=0$, $k=1,2,\cdots,K$. The results clearly demonstrate the fast convergence rate of the proposed CFO estimation method. Especially, in the case with relatively small maximum CFO, i.e., $\phi_{\max}=0.1$, two iterations are sufficient for convergence under the moderate SNR condition, while only a few more iterations are required for the low SNR condition or a larger maximum CFO. In the following, we consider $\phi_{\max}=0.2$ and five iterations are employed in FS-BEAM.

We show the CFO estimation performance of FS-BEAM with $K=4$ users in Fig. 4(a) being compared with RMCE-IIC. The mean DOAs are set the same as Fig. 3. The corresponding analytical results from (30) are plotted as the curves with marker ‘$\times$’, while the asymptotic analysis results computed from (31) are represented by the dotted curve. Both two angular spreads $\theta_{\rm as}=5^{\circ}$ and $\theta_{\rm as}=10^{\circ}$ are included. We also include the following two time-division training based schemes for comparison, labelled as ‘FS-Time-Div-I’ and ‘FS-Time-Div-II’, as illustrated in Fig. 5.

• •

FS-Time-Div-I: Each user transmits a short training sequence of length-$N_{s}$ ($N_{s}\!<\!N$) in a time-division manner at the beginning of uplink frame. CP is also inserted in front of each sequence to eliminate the interference between the adjacent sequences in multipath propagation channel. After the $K$ time-division sequences, all users start data block transmission simultaneously. In FS-Time-Div-I, BS can perform single-user joint CFO and channel estimation for each user individually in the traditional way. The overall training length of FS-Time-Div-I can be expressed as $K(N_{\rm cp}+N_{s})$.

• •

FS-Time-Div-II: At the beginning of uplink frame, each user first transmits time-division short training sequences of length-$N_{s}$ in the same way as FS-Time-Div-I. Then, each user simultaneously transmits a length-$N$ training block as well as the data blocks. BS first perform DOA identification for each user individually by using the time-division training sequences. Denote ${\bf Y}_{k,{\rm TD}}$ as the received $N_{s}\times M$ signal matrix corresponding to the short training sequence from the $k\rm{th}$ user. The DOA estimation for the $k\rm{th}$ user can be performed as $\hat{\theta}_{k,{\rm TD}}=\arg\max\limits_{\tilde{\theta}_{k}}\int_{\tilde{\theta}_{k}-\theta_{\rm as}}^{\tilde{\theta}_{k}+\theta_{\rm as}}\|{\bf Y}_{k,{\rm TD}}{\boldsymbol{a}}^{*}(\tilde{\theta})\|^{2}d\tilde{\theta}$. Then, similarly to (5), to mitigate MUI for the $k\rm{th}$ user, BS can employ all steering vectors ${\boldsymbol{a}}(\tilde{\theta})$ with $\tilde{\theta}\in(\hat{\theta}_{k,{\rm TD}}-\theta_{\rm as},\hat{\theta}_{k,{\rm TD}}+\theta_{\rm as})$ as the receive beamforming vectors on the subsequent length-$N$ training block as well as data blocks. Afterwards, the conventional single-user synchronization and channel estimation can be performed for each user based on the length-$N$ training block. The overall training length of FS-Time-Div-II can be expressed as $K(N_{\rm cp}+N_{s})+N_{\rm cp}+N$.

We consider $N_{s}=16$ in this example. The following observations can be made from Fig. 4(a): 1) We see that the proposed method almost touches the corresponding analytical results above SNR $=0$ dB. This indicates that our estimation method can effectively suppress the MUI effect by exploiting the high spatial resolution from the large-scale antenna array. In this example, the performance of our method is basically unchanged with different magnitude of angular spreads, unless under the low SNR condition. 2) It is also seen that our method could behave much better than RMCE-IIC and a performance gap as large as 10 dB is observed. Moreover, the proposed method also substantially outperforms FS-Time-Div-I scheme. This is not unexpected and because the training block in FS-BEAM could span a much wider time interval as compared to the time-division sequence for each user in FS-Time-Div-I, which provides the potential of higher resolution for CFO estimation [35]. In the meanwhile, the proposed scheme could work similarly as FS-Time-Div-II. Since FS-Time-Div-II decouples the DOA identification and CFO estimation problems, it should have lower computational complexity than FS-BEAM. Nevertheless, we should note that, the training overhead of FS-Time-Div-II is about $1+\frac{K(N_{s}+N_{\rm cp})}{N+N_{\rm cp}}=237\%$ of the proposed scheme in this example. In other words, we may conclude that as compared to FS-Time-Div-II, the proposed scheme could remarkably reduce the training overhead without obviously sacrificing the estimation performance.

The corresponding symbol error rate (SER) performance with $K=4$ users by using our proposed FS-BEAM scheme is depicted in Fig. 4(b). The result of FS-Time-Div-I is not included in this example due to its poor synchronization performance. Moreover, for comparison, we also include the ideal SER results without the multiuser CFOs as the benchmark, referred to as ‘Sync-Perfect-ZF’ and ‘Sync-Ch-Est-ZF’. Specifically, for the ideal benchmark schemes, we consider that all users are perfectly synchronized with BS in carrier frequency and conventional ZF data detection is employed. Moreover, BS has perfect knowledge of channel information in Sync-Perfect-ZF, while obtains channel information from conventional LS algorithm in Sync-Ch-Est-ZF.

Note that after the initial frequency synchronization, the residual frequency offset would lead to accumulative phase shift over the data blocks and may result in large detection error rate. Hence, phase tracking algorithms can be further employed to deal with this issue [30]. For example, the authors in [30] have analytically evaluated the effect of phase impairment due to phase noise plus a residual frequency offset and studied the corresponding compensation scheme. However, in this work we mainly focus on the initial frequency synchronization performance, and thus, for simplicity we adopt only the first data block to evaluate the SER performance in simulations. The following observations can be made:

1. The results clearly demonstrate the validity of our proposed frequency synchronization scheme. We also see that the SER curves of both Sync-Perfect-ZF and the proposed scheme drop more quickly with a larger angular spread under high SNR condition. This phenomenon can be explained as follows. A larger angular spread may increase the rank of the channel covariance matrix from each user [5], and thus can potentially provide a higher diversity gain of data detection performance [32].
2. It is seen that the performance gap between our scheme and Sync-Perfect-ZF is less than 2 dB with $K=4$ users. Our scheme can substantially outperform Sync-Ch-Est-ZF. This is not unexpected and because Sync-Ch-Est-ZF does not first perform user separation from spatial domain and needs to estimate a total of $KL=40$ channel parameters based on the $N=64$ received training samples. The insufficiency of training samples of Sync-Ch-Est-ZF would lead to both poor channel estimation and detection performance.
3. It is observed that SER curve of the proposed scheme basically touches the corresponding results of FS-Time-Div-II. This implies that as compared to FS-Time-Div-II, the proposed scheme could remarkably reduce the training overhead without obvious sacrifice on detection performance.

Next, we evaluate the CFO estimation and SER performance in Fig. 6 with $K=10$ users and $\theta_{\rm as}=5^{\circ}$. The mean DOAs of ten users are fixed as $\{30^{\circ},45^{\circ},120^{\circ},150^{\circ},80^{\circ},95^{\circ},20^{\circ},60^{\circ},110^{\circ},130^{\circ}\}$. We include both results of $M=128$ and $M=256$ in this example. The user separation in spatial domain becomes smaller with $K=10$ and thus lead to the increased MUI effect. We can observe some discrepancy between FS-BEAM and the corresponding analytical results under high SNR from Fig. 6(a). As expected, more receive antennas at BS bring higher spatial resolution and thus improve the system performance in both CFO estimation and data detection. Especially, FS-BEAM shows some observable SER performance floor with $M=128$, while it would work far more better when $M$ increases to $256$. Moreover, we see that FS-BEAM outperform FS-Time-Div-II under moderate and high SNR condition. This should be attributed to the adaptive beamforming scheme in FS-BEAM, which employs only the qualified DOAs for beamforming. In comparison, FS-Time-Div-II employs all DOAs within the angular spread of each user for beamforming, which suffers from a certain magnitude of MUI effect with $K=10$ users.

Considering the user locations are randomly distributed, we next assume the mean DOA of each user is randomly generated in each simulation trial. Specifically, we assume the mean DOAs of half of the users follow uniform distribution from $30^{\circ}$ to $60^{\circ}$ where the rest half follow uniform distribution from $120^{\circ}$ to $150^{\circ}$. We plot the CFO estimation performance with $K=4$ of FS-BEAM and FS-BEAM-UG in Fig. 7(a). Both angular spreads of $\theta_{\rm as}=5^{\circ}$ and $\theta_{\rm as}=10^{\circ}$ are considered in this example. From the results, we can observe a performance error floor of the proposed FS-BEAM, especially with larger angular spread. This is not unexpected and mainly due to the cases where some users have very similar mean DOAs and the spatial separation may not be enough. Nevertheless, it is observed that our proposed FS-BEAM-UG could work with randomly distributed users and can achieve much better performance as compared to FS-BEAM, which indicates the robustness of FS-BEAM-UG. Moreover, we see that the proposed FS-BEAM-UG also substantially outperforms the competitor of RMCE-IIC and a performance gap more than 5 dB can be observed. The corresponding SER performance with $K=4$ is depicted in Fig. 7(b). Both the ideal benchmarks of Sync-Perfect-ZF and Sync-Ch-Est-ZF are included for comparison. As expected, the proposed FS-BEAM is sensitive to the random distribution of user DOAs and cannot work properly in this example. This suggests that our FS-BEAM may rely on some proper user schedule strategy to avoid the case where some users do not have enough spatial separation. Nevertheless, we see that our proposed FS-BEAM-UG still would work properly and substantially outperform both the Sync-Ch-Est-ZF and FS-BEAM scheme.

Next, we evaluate the CFO estimation and SER performance of $K=10$ users with randomly distributed DOAs in Fig. 8. We include additional results of $M=256$ and $N=128$ in this example. The angular spread is assumed as $\theta_{\rm as}=5^{\circ}$. As expected, FS-BEAM shows little change when $M$ or $N$ increases due to the severe MUI effect. In comparison, we can observe performance improvement of the proposed FS-BEAM-UG by increasing $M$ or $N$. Specifically, we can still observe error performance floor of FS-BEAM-UG in terms of both CFO estimation and data detection with $N=64$ even when $M$ increases to 256, whereas it can behave much better when $N$ increases to 128. This is because FS-BEAM-UG should perform joint estimation for the grouped users and the training sequence length becomes a bottleneck in this situation. Moreover, Fig. 8(a) shows that, with a larger training block length of $N=128$, the CFO estimation MSE of RMCE-IIC could also decrease as SNR increases. Nevertheless, the performance gap between RMCE-IIC and the proposed FS-BEAM-UG is still more than 10 dB.

In the next example, we investigate the sensitivity of the proposed scheme with respect to the imperfect knowledge of angular spread. Suppose BS takes the angular spread as $\hat{\theta}_{\rm as}$. In Fig. 9, we evaluate the impact of the assumed angular spread $\hat{\theta}_{\rm as}$ on the CFO estimation performance under SNR $=5$ dB. We consider $\theta_{\rm as}=5^{\circ}$ in this example. Both the results of fixed DOAs and random DOAs are considered, represented by the solid and dashed curves, respectively. Specifically, the fixed DOAs are generated in the same way as Fig. 4 and Fig. 6 for $K=4$ and $K=10$, respectively. The random DOAs are generated in the same way as Fig. 7 and Fig. 8. We can observe that with $K=4$ users, only a little performance degradation when the assumed angular spread $\hat{\theta}_{\rm as}$ deviates from the true value, which says that our schemes are robust to the imperfect knowledge of angular spread in this situation. However, from the results with $K=10$ users, we see that the proposed FS-BEAM scheme shows much more sensitivity to the imperfect knowledge of angular spread. This is not unexpected and because more users would lead to fewer ‘white space’ in spatial domain and decrease the robustness to imperfect knowledge of angular spread.

In Fig. 10, we plot the probability density function (PDF) of the CFO estimation error of the proposed FS-BEAM with different number of receive antennas. In this example, we include both the cases for $K=4$ and $K=10$, where the DOAs are generated in the same way as Fig. 4 and Fig. 6, respectively. As expected, the error distribution tends to be more concentrated to zero with more receive antennas.

It is also interesting to evaluate the performance of the proposed method as $M$ increases. In Fig. 11, we plot the SER performance of the proposed FS-BEAM scheme with $K=4$ users. We consider the following two different DOA configurations: $\{30^{\circ},60^{\circ},120^{\circ},150^{\circ}\}$ and $\{30^{\circ},32^{\circ},120^{\circ},150^{\circ}\}$, represented by the solid and dashed curves, respectively. The overall SNR of the multiple receive antennas can be expressed as $M\sigma_{s}^{2}/\sigma_{n}^{2}$, which is the ratio between the overall average received signal power at multiple receive antennas and the noise power. The aforementioned SNR defined as $\sigma_{s}^{2}/\sigma_{n}^{2}$ can be interpreted as the SNR per receive antenna, since the average channel gain from each user at one receive antenna is normalized, i.e., $\sum_{l=1}^{L}\sigma_{h,l}^{2}=1$. In this example, we vary the number of antennas with fixed overall SNR $M\sigma_{s}^{2}/\sigma_{n}^{2}=19$ dB. The performance of maximum ratio combining (MRC) detection with perfect knowledge of channel information and perfect synchronization is also plotted for comparison, labelled as ‘Sync-Perfect-MRC’. The following observations can be made:

1. Regarding the solid curves, since each user has enough spatial separation from the others, we see that Sync-Perfect-MRC closely approaches the corresponding performance of Sync-Perfect-ZF, especially as $M$ increases. Moreover, note that SNR $\sigma_{s}^{2}/\sigma_{n}^{2}$ in fact decreases when $M$ increases with fixed overall SNR $M\sigma_{s}^{2}/\sigma_{n}^{2}$. The proposed FS-BEAM exhibits no much SER performance improvement as $M$ increases with a smaller training sequence length $N$. Nevertheless, with enough large $N$, the proposed FS-BEAM also approaches Sync-Perfect-ZF.
2. In comparison, from the dashed curves, we see that due to the insufficient spatial separation, the performance of both Sync-Perfect-MRC and FS-BEAM suffers from substantial deterioration and becomes far away from that of Sync-Perfect-ZF.
3. We also plot the single-user analytical SER [31, eq. (45)] under SNR= 19 dB and AWGN channel as a benchmark, which can be considered as the asymptotic performance when $M$ tends to be infinity. Owing to the MUI interference as well as the channel fading effect, a performance gap between Sync-Perfect-ZF and the asymptotic benchmark is still observable.

Let us now evaluate the computational complexity of the proposed CFO estimation in terms of complex multiplications. We denote $\alpha$ as the number of required iterations in FS-BEAM. Then, the complexity of the proposed FS-BEAM is approximately in the order of \mathcal{O}\big{(}NM_{\rm FFT}\log_{2}M_{\rm FFT}+2\alpha KLNM_{\rm FFT}\big{)}\simeq\mathcal{O}\big{(}2\alpha KLNM_{\rm FFT}\big{)}. In addition, considering the multiuser estimation in one group with user number of $\kappa$ in FS-BEAM-UG, SVD on $\mathbb{B}(\hat{\boldsymbol{\phi}})$ requires around $\mathcal{O}(\kappa^{2}L^{2}N)$, and computing ${\bf R}_{\rm Grp}$ requires $\mathcal{O}(\kappa LN^{2})$. Calculating (IV-A) requires around \mathcal{O}\big{(}\kappa LN^{2}+2\kappa^{2}LN^{2}+2\kappa N^{3}\big{)}. Then, the total complexity for this group can be expressed as \mathcal{O}\big{(}2\kappa N^{3}+2(\kappa^{2}+\kappa)LN^{2}+\kappa^{2}L^{2}N\big{)}\simeq\mathcal{O}(4\kappa N^{3}). Considering the extreme case where all users are grouped, as compared to FS-BEAM, the FS-BEAM-UG further requires the complexity of $\mathcal{O}(4KN^{3})$. On the other side, the complexity of RMCE-IIC can be expressed as $\mathcal{O}(\eta\beta KLNM)$, where $\eta$ is the required iterations and $\beta$ denotes the number of trial CFOs in [17, eq. (21)]. Let us consider the example of above $K=6$ scenario. The required complexities of the proposed FS-BEAM, FS-BEAM-UG and RMCE-IIC can be approximately expressed as $\mathcal{O}(5.9\times 10^{6})$, $\mathcal{O}(1.2\times 10^{7})$ and $\mathcal{O}(1.8\times 10^{8})$, respectively. This indicates that the proposed schemes also have the advantage of lower complexity as compared to RMCE-IIC.

VI Conclusions

In this paper, we developed a new frequency synchronization scheme for massive multiuser uplink transmissions by exploiting the angle domain information of different users. We designed both the CFO estimation and beamforming method, and convert the original complicated multi-user problem into equivalent single user transmission model. We further design an improved user grouping scheme to deal with the unexpected scenarios that some users may not be separated well from the spatial domain. Both the simulation and theoretical results are provided to corroborate the proposed studies.

Appendix A Asymptotic CFO Estimation Performance

We have ${\bf R}_{k}^{2}=\frac{1}{4\theta_{\rm as}^{2}}\int_{\theta_{k}-\theta_{\rm as}}^{\theta_{k}+\theta_{\rm as}}\int_{\theta_{k}-\theta_{\rm as}}^{\theta_{k}+\theta_{\rm as}}{\boldsymbol{a}}^{*}(\tilde{\theta}_{1}){\boldsymbol{a}}^{T}(\tilde{\theta}_{1}){\boldsymbol{a}}^{*}(\tilde{\theta}_{2}){\boldsymbol{a}}^{T}(\tilde{\theta}_{2})d\tilde{\theta}_{1}d\tilde{\theta}_{2}$. With small angular spread $\theta_{\rm as}$, we further obtain

[TABLE]

The derivation of (59) arises from the following approximation: For any function $f(\cdot)$ and $a>0$, $b>0$, there holds $\lim\limits_{M\to\infty}\int_{x_{0}-a}^{x_{0}+b}\frac{\sin\frac{\chi M(x-x_{0})}{2}}{\chi(x-x_{0})/2}f(x)dx\simeq\lim\limits_{M\to\infty}\int_{-\infty}^{\infty}\frac{\sin\frac{\chi M(x-x_{0})}{2}}{\chi(x-x_{0})/2}f(x)dx=\frac{2\pi}{\chi}f(x_{0})$. Then, based on (59), we can rewrite ${\bf R}_{k}^{2}$ as

[TABLE]

Hence, with sufficiently large $M$, we have the approximation: $\textrm{Tr}({\bf R}_{k}^{2})\simeq\frac{\pi}{\theta_{\rm as}\chi\sin\theta_{k}}\textrm{Tr}({\bf R}_{k})\frac{M\pi}{\theta_{\rm as}\chi\sin\theta_{k}}$ and $\textrm{Tr}({\bf R}_{k}^{3})=\frac{\pi^{2}}{\theta_{\rm as}^{2}\chi^{2}\sin\theta_{k}^{2}}\textrm{Tr}({\bf R}_{k})=\frac{M\pi^{2}}{\theta_{\rm as}^{2}\chi^{2}\sin\theta_{k}^{2}}$. Then, the asymptotic CFO estimation MSE with sufficiently large $M$ for the $k\rm{th}$ user can be expressed as

[TABLE]

Next, with sufficiently large $N$, we have the approximation: ${\bf B}^{H}{\bf B}\simeq N\sigma_{s}^{2}{\bf I}_{L}$. Then, there holds

[TABLE]

After substituting (A) into (61), we arrive at (31).