{'fraction_non_alphanumeric': 0.043174700960008215, 'fraction_numerical': 0.009677088146208738, 'mean_word_length': 4.299102285092492, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Frequency Synchronization for Uplink Massive MIMO Systems

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.

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]- [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]- [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, 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]- [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 March 14, 2017 DRAFT The mean DOA and angular spread of the kth user are θ k and θas, respectively.
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 (·) * , (·) T , (·) H , (·) † and E[·] represent conjugate, transpose, Hermitian, pseudo-inverse and expectation, respectively; j = √ −1 is the imaginary unit; · denotes the Frobenius norm operator; for a vector x, diag(x) is a diagonal matrix with main diagonal of x; ⊗ stands for the Kronecker product; vec(·) is the vectorization operator; C m×n defines the vector space of all m × n complex matrices; I N is the N × N identity matrix; 0 represents an all-zero matrix with appropriate dimension; Tr(·) 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 ≫ 1 antennas in the form of ULA, as shown in Fig. 1 which represents the CFO-induced phase shift over one OFDM block. We consider the fractional CFO, i.e., |φ k | < 0.5, which should be sufficient in the multiuser uplink [26], [27].
The classical one-ring channel propagation model is adopted [5]- [7], where each user is surrounded by a ring of P ≫ 1 local scatterers. The multi-tap channel matrix between the kth user and the BS can be modeled as an is composed of P rays and can be expressed as: Here, θ l,p,k represents the DOA of the pth ray in the l-th channel tap, while α l,p,k ∼ CN (0, σ 2 h,l /P ) represents the corresponding complex gain and is independent from each other [6]. Besides, σ 2 h,l , l = 1, 2, · · · , L, models the power delay profile (PDP) of the propagation channel. We consider L l=1 σ 2 h,l = 1 such that the total channel gain of each user at one receive antenna is normalized. Moreover, a(θ l,p,k ) ∈ C M ×1 is the steering vector and has the form a(θ l,p,k ) = 1, e −jχ cos θ l,p,k , · · · , e −jχ(M −1) cos θ l,p,k T , where χ = 2πd λ , d is the antenna spacing, λ denotes the signal wavelength. Similar to [5]- [7], [25], [34], we consider that the incident rays of each user can be constrained within a certain angular spread θ as ; Namely, the incident angles of the kth user is limited in the DOA region (θ k − θ as , θ k + θ as ) with uniform distribution, where θ k represents the mean DOA of the kth 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.
Let F L stand for the submatrix of F containing the first L column vectors and E(φ k ) represent the .
Each user transmits a whole training block at the beginning of the uplink frame, which is then March 14, 2017 DRAFT followed by a number of data blocks. The uplink communication strategy is described in Fig. 2.
T as the frequency-domain training block transmitted from the kth user. Then, the received time-domain training signal after CP removal can be written as the following N × M matrix: where where η i (φ k ) = e j 2πi(N+Ncp )φ k N represents the accumulative phase rotation introduced by the CFO of the kth user with N cp being the length of CP, and N d i ∈ C N ×M denotes the corresponding noise

A. CFO Estimation
H k as the orthogonal projection operator onto the space spanned by the columns of B k , and define P ⊥ B k = I N − P B k . In order to guarantee the identifiability of the CFO estimation for the kth user, we consider the training sequence is designed such that 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 kth 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 kth 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 → ∞ 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 a(θ 1 ) H a(θ 2 )/M = 0 holds for any θ 1 = θ 2 [6], [25]. Hence, for the kth user, when the angular spread of users are disjoint, the beamforming vector a(θ) with spatial alignment, i.e.,θ ∈ (θ k − θ as , θ k + θ as ), will eliminate the signals from the other users; namely, there are H k a * (θ) = 0 and H q a * (θ) = 0 for any q = k. We then obtain which gives an equivalent frequency asynchronous single-user received signal model for the kth user. The MUI term of (5) approaching zero arises from the following observation: The lth entry of H q a * (θ), l = 1, 2, · · · , L, equals (h (q) l ) T a * (θ) = P p=1 α l,p,q a T (θ l,p,q )a * (θ). Once the angular spreads of the kth and the qth user are disjoint, a(θ) is nearly orthogonal to any a(θ l,p,q ) with sufficient large M [6], [25].
We can further perform proper CFO compensation of (5) for purpose of frequency alignment as 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 kth user.
The resultant signal in (6) would exactly lie in the column space of B k in the absence of noise.
It can be observed that, ifθ / ∈ (θ k − θ as , θ k + θ as ) or ifφ = φ k , then the multiuser interference or frequency asynchronism would deviate the column space of E(φ) H Ya * (θ) from that of Span(B k ).
Based on the above observations, the frequency synchronization for the kth user can be designed as the following joint spatial-frequency alignment procedure. Specifically, the joint CFO and DOA estimation for the kth user can be obtained by minimizing the cost function: where the trial DOAθ k and the trial CFOφ k target at the spatial and frequency alignment, respectively. Mathematically, the joint CFO and DOA estimation for the kth user can be expressed as θ k ,φ k = arg miñ 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.,θ k = θ k andφ k = φ k .
Proof: Whenθ k = θ k , all DOAs from the kth user will be included in the integration of (7) and then the beamforming vectors a(θ) in (7) will eliminate the signals from all users other than the kth user. In addition, As a result, the cost function is expressed as which tends to zero as SNR increases. On the other hand, the spatial and frequency misalignment includes the following four cases: 1) In the case ofθ k = θ k andφ k = φ k , the numerator of the cost function C k (θ k ,φ k ) approximately equals n , which indicates that additional signal power will be included in the numerator of (7). Then there is 2) When 0 < |θ k − θ k | < 2θ as , the integration interval (θ k − θ as ,θ k + θ as ) intersects but not exactly coincides with the DOA region of the kth user. In this case, not all of the signal components of the kth user are included in the denominator of the cost function. Thus, there 3) When |θ k − θ q | > 2θ as , ∀q = 1, 2, · · · , K, the integration interval (θ k −θ as ,θ k +θ as ) dose not intersect with the DOA region of any user. In this case, the denominator of the cost function approximately equals θ k +θas θ k −θas Na * (θ) 2 dθ ≃ 2θ as MNσ 2 n and the cost function becomes

4) When the integration interval intersects with the DOA region of one user (say the qth user)
other than the kth user, the numerator of the cost function C k (θ k ,φ k ) approximately equals In this case, the cost function dose not tend to zero as SNR increases. Hence, C k (θ k ,φ k ) > C k (θ k , φ 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 March 14, 2017 DRAFT estimation of both CFO and DOA for the kth 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 kth user in the nth iteration byφ k )Y, which corresponds to the signal matrix after the CFO compensation byφ which is a parabola function with respect to ∆φ. The above parabola function achieves its minimum value when ∆φ equals By substituting ∆φ (10), the DOA estimation in the (n + 1)-th iteration for the kth user can be expressed aŝ Accordingly, the CFO estimation can be updated aŝ Note that optimization of (12) can be efficiently solved by exploiting FFT operation. Specifically, let y(θ) = Ya * (θ) represent the signal after the receive beamforming with steering vector a * (θ).
We can rewrite (12) intô Note that P B k has a dimension of N × N, and thus directly computing (14) may suffer from a high computational burden. To further reduce the computational complexity, we can express P B k is a unitary matrix and has the same column space of B k , i.e., Span(Q k ) = Span(B k ). We further obtain Let both the trial DOAθ k and the DOA samplesθ in the numerical integration be drawn from M FFT ] T and thus, corresponds to the FFT conversion of Y along the spatial dimension. Denote the neighbor samples ofθ k by The numerical solution of (15) can be then expressed aŝ 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 y(D i ) dominates the computational burden of the grid search in (18), the computational complexity can be greatly reduced by exploiting FFT in (16).

B. Multi-Branch Beamforming
After the above joint CFO and DOA estimation, all steering vectors a * (θ) withθ ∈ (θ k −θ as ,θ k + θ as ) can be utilized as the receive beamforming vectors for the kth user, and the CFO estimateφ k March 14, 2017 DRAFT can be used for compensation. With perfect estimation and infinite receive antennas, i.e.,φ k = φ k , θ k = θ k , M → ∞ and ∀θ ∈ (θ k − θ as , θ k + θ as ), we can obtain the signals after beamforming corresponding to the training and data blocks in noise-free environment as follows: The equivalent channel response H k a * (θ) 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 (θ k − θ as , θ k + θ as ) are suitable as the beamforming direction for the kth 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 (θ k − θ as ,θ k + θ as ), we propose an alternative adaptive multi-branch beamforming scheme in this subsection. Note that after beamforming with a * (θ) and CFO compensation withφ k on the received signal matrix Y, the magnitude of 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 kth user with beamforming vector of a * (θ). Let us draw the trial DOAs from set D FFT . Whenθ does not correspond to the kth user, the magnitude of C k (θ,φ k ) can be approximately given by L/(N −L). Then, we can define the set of qualified DOAs for the kth user as where t h is a pre-set threshold that can be larger than one. Here, the elements of ϑ (k) are picked out from set D 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 kth user can be obtained as

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 (θ k − θ as , θ k + θ as ) 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 kth user is equivalent where The mean square error (MSE) of CFO estimation under the high SNR can be expressed as [24] MSE{φ Then, there is Define R k = θ k +θas θ k −θas 1 2θas a * (θ)a T (θ)dθ with its (i, j)-th entry being R k (i, j) = 1 2θas θ k +θas θ k −θas e jχ(i−j) cosθ dθ. We can further simplify (27) as where Σ h = diag(σ h,1 , σ h,2 , · · · , σ h,L ) corresponds to the channel PDP. The derivation of (28) arises from the fact that E (h 1 P a * (θ l,p,k )a T (θ l,p,k ) ) = σ 2 h,l · R k . For the second order derivative in (26), we can obtain the following approximation: By substituting (28) and (29) into (24), the CFO estimation MSE for the kth user can be expressed as: Let σ 2 s denote the transmit signal power from each user. With sufficiently large M and N, the asymptotic CFO estimation MSE for the kth user can be expressed as 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 W (k) can be expressed as . Once the estimated SNR of the kth user is below a threshold, say ρ k < ρ th · ρ ex where ρ th is a pre-set value, there may exist non-negligible MUI from other users. In this case, we denote the kth 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 G.
Mathematically, denote the user index set of the gth group by U (g) = {c Kg }, where K g is the number of the users in this group. Then, there should be min max |c 1 − c 2 |, ∀c 1 ∈ U 1 and ∀c 2 ∈ 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 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 κ users stay in one group, i.e., U (1) = {1, 2, · · · , κ}, and illustrate the the CFO estimation and data detection for this group. Similar approach can be readily extended to other groups.

A. CFO Estimation
Suppose the DOA range for the considered group is [θ min , θ max ], where θ min and θ max are the minimum and maximum qualified DOAs from all κ users, respectively, i.e., θ min = min k=1,··· ,κ ϑ (k) and θ max = max k=1,··· ,κ ϑ (k) . Let us then pick up the following Q = κL uniformly distributed DOAs within [θ min , θ max ] to formulate the beamforming vectors for this group: ψ i = θ min + (i − 1) θmax−θ min Q−1 , i = 1, 2, · · · , Q. The corresponding beamforming matrix can be obtained as March 14, 2017 DRAFT W Grp = a(ψ 1 ), a(ψ 2 ), · · · , a(ψ Q ) * , and the beamforming is performed as 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 W Grp could substantially mitigate the MUI from the users outside the group. Then, we can approximately consider that Y Grp is composed of only the κ users in this group.
where N Grp = NW Grp + K q=κ+1 E(φ q )B q H q W Grp denotes the corresponding noise plus the mitigated MUI terms from the users outside the group. Letφ = [φ 1 ,φ 2 , · · · ,φ κ ] T represent the trial value of φ. According to (33), the LS joint CFO estimation of the κ users in this group can be given bŷ where the cost function is is the orthogonal projection on the column space of B(φ).
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.,φ k , may still be able to serve as a valid coarse estimation. In the following, we propose to approximate the cost function C Grp (φ) as a high dimensional parabola in the neighborhood of coarse estimationφ = [φ 1 ,φ 2 , · · · ,φ κ ] T , and the minimization of (34) can be then solved in a more efficient way. Specifically, there is We have the approximation and ∆φ k =φ k −φ k . We can further obtain where Performing eigenvalue decomposition on Π(φ) gives where U s ∈ C N ×κL and U n ∈ C N ×(N −κL) correspond to the signal and noise subspace, respectively.

B. Data Detection
After the CFO estimation, the equivalent channel response matrix can be obtained from (33) via the LS algorithm:Ĥ eq = B † (φ Grp )Y Grp . Similar to (32), the beamforming can also be performed for the data blocks as: where N d i,Grp ∈ C N ×Q denotes the corresponding noise plus the interference from the users outside the group. We rewrite the equivalent channel matrix for the kth user in (51) as H q,eq ∈ C N ×N stand for the circular convolution matrix corresponding to the equivalent channel response h (k) q,eq , i.e., The vectorization of (51) is where Then, the ZF data detection for the κ users can be performed asŝ i = (η * i (φ)⊗F)H † eq vec(Y d i,Grp ). In practice, φ can be replaced byφ Grp , while H eq can be constructed from bothĤ eq andφ Grp in the way as (52), (57)  The normalized CFO is randomly generated from −φ max to φ max . The MSE of the normalized CFO is adopted as the figure of merit. The thresholds are taken as t h = 10 and ρ 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(jπi 2 /N), i = 0, 1, · · · , 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.   Denote Y k,TD as the received N s × M signal matrix corresponding to the short training sequence from the kth user. The DOA estimation for the kth user can be performed aŝ Then, similarly to (5), to mitigate MUI for the kth user, BS can employ all steering vectors a(θ) withθ ∈ (θ k,TD − θ as ,θ k,TD + θ 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 cp + N s ) + N cp + N.
We consider N s = 16 in this example. The following observations can be made from Fig. 4(a  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 + K(Ns+Ncp) N +Ncp = 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  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 • to 60 • where the rest half follow uniform distribution from 120 • to 150 • . We plot the CFO estimation performance with K = 4 of FS-BEAM and FS-BEAM-UG in Fig. 7(a). Both angular spreads of θ as = 5 • and θ as = 10 • 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 θ as = 5 • . As expected, FS-BEAM shows little change when M or   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θ as . In Fig. 9, we evaluate the impact of the assumed angular spreadθ as on the CFO estimation performance under SNR = 5 dB. We consider θ as = 5 • 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θ 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 ig. 11. SER performance of the proposed schemes as a function of M . 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 represented by the solid and dashed curves, respectively. The overall SNR of the multiple receive antennas can be expressed as Mσ 2 s /σ 2 n , which is the ratio between the overall average received signal power at multiple receive antennas and the noise power. The aforementioned SNR defined as σ 2 s /σ 2 n 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., L l=1 σ 2 h,l = 1. In this example, we vary the number of antennas with fixed overall SNR Mσ 2 s /σ 2 n = 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 σ 2 s /σ 2 n in fact decreases when M increases with fixed overall SNR Mσ 2 s /σ 2 n . 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.

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 multiuser 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.