Frequency Synchronization Algorithms for MIMO-OFDM Systems with Periodic Preambles

Abstract

This paper addresses the problem of frequency synchronization in multiple-input multiple-output (MIMO) orthogonal frequency-division multiplexing (OFDM) systems with periodic preambles. Two Maximum Likelihood Estimators (MLEs) and four low complexity Best Linear Unbiased Estimators (BLUEs) are proposed. When cyclic shift delay is applied to periodic preambles, MLEs can adopt a unique structure which employs diversity and correlation calculation in the estimator. We first propose a multiple-antenna estimator named as MLE-MA, which has an estimation range of the entire signal bandwidth. MLE-MAS is a simplified version of MLE-MA, which reduces the computational load at the cost of a narrower estimation range. Moreover, phases or differential phases obtained from the autocorrelation of periodic sequences with different delays can be used to estimate carrier frequency offset (CFO). The weighting coefficients of all BLUEs are derived in closed form and shown to be independent from the signal-to-noise ratio (SNR). BLUEs have asymptotically consistent variances. Theoretical and simulation results demonstrate that all estimators can asymptotically approach the theoretical Cramer Rao Lower bound (CRLB) at different SNRs.

1. Introduction

Multiple-input multiple-output (MIMO) systems have greater capacities than single-input and single-output (SISO) systems. This capacity increases linearly with the minimum of transmitter and receiver antenna [1, 2] if uncorrelated Tx-Rx antenna pairs are assumed. Diversity and multiplexing techniques can be employed to improve the reliability and spectral efficiency of MIMO systems [3], respectively. For more practical case, one can refer to some useful MIMO channel models [4, 5] for capacity analysis and simulations. Orthogonal frequency-division multiplexing (OFDM) splits a wide frequency band into a bank of narrow subbands. Although the whole bandwidth may be subject to frequency selective fading, each subband normally undergoes flat fading. This flat fading effect can be easily compensated by a single-tap frequency domain equalizer. Therefore, OFDM can overcome inter-symbol interference (ISI) very well [6] compared with single carrier systems. MIMO combined with OFDM [3] is seen as an essential technique for advanced communication systems.

Frequency synchronization is crucial to any communication system. In particular, OFDM is sensitive to carrier frequency offset (CFO), which is caused by oscillator mismatch between transmitter (Tx) and receiver (Rx) and/or Doppler shift. If CFO is not compensated at the Rx, accumulated phase rotation, amplitude degradation, and intercarrier interference (ICI) will severely degrade the overall system performance [7]. Thus, accurate CFO estimation is indispensable for OFDM systems. In MIMO scenarios, signals from different transmit antennas are superimposed at the Rx, so that the accuracy and reliability of CFO estimation can be improved using diversity techniques.

Synchronization techniques for SISO systems have been well documented. Due to its low complexity, periodic preambles are widely utilized in wireless standards, such as IEEE 802.11a/n [8]. Schmidl and Cox [9] employed two training symbols to accomplish the synchronization task. Timing and fractional CFO estimations were obtained by the first symbol, which is composed of two identical short sequences. CFO was calculated as the phase of timing metric at the correct timing position. The second symbol was exploited to estimate the integer CFO. However, the use of two data symbols caused extra overhead. Two Maximum Likelihood Estimators (MLEs) and an ad hoc CFO estimator for signal transmission over frequency selective channels were considered in [10]. The MLEs were based on Fast Fourier Transform (FFT) implementation and the ad hoc method was based on best linear unbiased estimator (BLUE). The sensitivity of the BLUE algorithms to quantization errors of the arg function (the argument of a complex number) was analyzed and quantified in [11]. Further discussions on BLUEs were presented by Minn et al. in [12], where three BLUEs were proposed. However, the weighting coefficients were calculated based on simulations. Therefore, the weighting coefficients were not closed form and depend on the signal-to-noise ratio (SNR). The above papers all concentrated on single antenna systems. For synchronization in MIMO-OFDM systems, some contributions have been reported [13, 14]. Through solving the roots of a real polynomial, fractional CFO was estimated in [13]. A practical and low complexity algorithm without solving polynomial was discussed in [14]. After a coarse compensation, fractional CFO estimation was found by trigonometric approximation and integer CFO was also obtained by a search method.

This paper proposes novel frequency synchronization methods for MIMO-OFDM systems with periodic preambles. Unlike the above-mentioned schemes [9, 10, 12–14], the presence of cyclic delay is further considered in this paper. After joint estimation of channel and CFO, we first propose two MLEs which employ spatial diversity and weighted correlation calculation in the estimator. Weighting coefficient for correlation calculation only relies on the number of pilots, which is totally independent of modulation signals. Although different cyclic delays may be introduced in the preambles, the proposed estimators have a definite and unique structure. Compared with nonperiodic preambles, correlation based estimator for periodic preambles has lower complexity and simpler structure. Similar to the estimators in [10], FFT implementation and numerical interpolation are exploited in the proposed MLEs. The first MLE has a very large estimation range up to the signal bandwidth of OFDM but requires a large point FFT implementation. The second MLE can reduce complexity but has a smaller estimation range. Furthermore, four different low complexity BLUEs are proposed. The accumulated phase increases linearly as time elapses in OFDM, so the phase of time delay autocorrelation for periodic sequences is mainly overwhelmed by CFO. This characteristic is utilized to estimate CFO. BLUE-PA and BLUE-PA-N take the phase of correlation metrics directly and a coarse estimation is needed to enlarge the estimation range, while BLUE-DPA and BLUE-DPA-N utilize differential phase of correlation metrics to avoid coarse estimation steps. In particular, BLUE-DPA-N is a straightforward extension of the estimator in [15]. BLUE-PA and BLUE-DPA perform a little better than BLUE-PA-N and BLUE-DPA-N by considering noise-noise product items when deriving covariance. We note that weighting coefficients in [12] are based on simulations and related to the SNRs. In contrast, all weighting coefficients in this paper have closed-form expressions and are independent from SNRs.

The imperfect CFO estimation and compensation will degrade the performance of SISO-OFDM and MIMO-OFDM systems. There are some literatures to discuss the impact of CFO on the channel estimation or system BER performance [16, 17]. We do not discuss the question in this paper due to page limits. The interested reader can refer to these literatures.

The rest of this paper is organized as follows. The system model is described in Section 2. ML CFO estimation algorithms for periodic preambles are proposed in Section 3. Four different BLUEs are presented in Section 4. In Section 5, computational complexity of MLEs and BLUEs are discussed. Section 6 shows the simulation results of frequency synchronization algorithms. At last, conclusions are drawn in Section 7.

Notations. Superscripts ${(·)}^{*}$ , ${(·)}^T$ , ${(·)}^H$ , and ${(·)}^{+}$ denote conjugate, transpose, Hermitian transpose, and Moore-penrose inversion, respectively. $x\hspace{0.17em}|\hspace{0.17em}P$ denotes x module P. Matrices ${\mathbf{F}}_N$ , ${\mathbf{I}}_N$ , and ${\mathbf{1}}_{N\times N}$ are N-by-N Fourier matrix, identity matrix, and all-ones matrix, respectively. ${\mathbf{F}}_{N,L}$ takes the first L columns of ${\mathbf{F}}_N$ , and ${\mathbf{1}}_N$ is a N-by-1 all-ones column vector. Moreover, $\lfloor ·\rfloor$ and ${\parallel ·\parallel }^2$ denote floor operation and square of Euclidean norm, respectively, ⊙ and ⊗ denote Hadamard and Kronecker product operators, respectively, and $\Re \{·\}$ , $\Im \{·\}$ , and $\arg \{·\}$ denote the real part, imaginary part, and principal argument of a complex number, respectively.

2. System Model

Considering a MIMO-OFDM system with $N_t$ transmit antennas and $N_r$ receive antennas. The number of subcarriers is N. The frequency domain pilot sequence of the ith Tx branch is ${\mathbf{c}}_i={[c_i(0),c_i(1),\dots ,c_i(N-1)]}^T$ , where $1\le i\le N_t$ . At first, we do not apply any constraints, for example, periodic property, on the pilot sequences. This means that we first consider the general case. Usually, some virtual subcarriers are inserted at high frequency components. Let $N_u$ be the number of subcarriers that can be utilized, that is, except virtual subcarriers and DC component.

The channel impulse response (CIR) of the frequency selective fading channel between the ith Tx and the jth Rx antenna is modeled as a L order finite impulse response (FIR) filter, ${\mathbf{h}}_{i,j}={[h_{i,j}(0),h_{i,j}(1),\dots ,h_{i,j}(L-1)]}^T$ , where $1\le j\le N_r$ . Channels between different antenna pairs are assumed to be statistically independent. Each tap coefficient is a complex Gaussian random variable with zero mean. Cyclic prefix (CP) is inserted before each OFDM symbol with a length of $N_g$ and we assume that $L\le N_g$ . At the Rx end, after discarding CP, the received signal vector at the jth Rx branch is ${\mathbf{r}}_j={[r_j(0),r_j(1),\dots ,r_j(N-1)]}^T$ . CFO is denoted by ν, which is normalized to the subcarrier spacing; then ${\mathbf{r}}_j$ is given by

$$\begin{array}{cc}& {\mathbf{r}}_j=\mathbf{\Gamma }{\displaystyle \sum }_{i=1}^{N_t}‍{\mathbf{F}}_N^H\mathrm{diag}\left\{{\mathbf{c}}_i\right\}{\mathbf{F}}_{N,L}{\mathbf{h}}_{i,j}+{\mathbf{v}}_j,\end{array}$$

(1)

where $\mathbf{\Gamma }=\mathrm{diag}\{1,e^{j2\pi \nu /N},\dots ,e^{j2\pi \nu (N-1)/N}\}$ . The entries of Fourier matrix ${\mathbf{F}}_N$ are ${[{\mathbf{F}}_N]}_{n,k}=1/\sqrt{N}{e}^{-j2\pi nk/N}$ , where $0\le n$ , $k\le N-1$ . AWGN noise vector at the jth Rx branch is ${\mathbf{v}}_j={[v_j(0),v_j(1),\dots ,v_j(N-1)]}^T$ . The noise samples are independent circularly symmetric zero-mean complex Gaussian random variables with variance $1/2{\sigma }_v^2$ per dimension.

Due to the presence of CP, channel linear convolution is converted into circular convolution. The noise-free received signal vector is ${\mathbf{x}}_j={\sum }_{i=1}^{N_T}‍{\mathbf{S}}_i{\mathbf{h}}_{i,j}$ , where ${\mathbf{S}}_i$ is transmit signal matrix at the ith antenna, ${\mathbf{S}}_i={\mathbf{F}}_N^H\mathrm{diag}\{{\mathbf{c}}_i\}{\mathbf{F}}_{N,L}$ . The received signal power is defined as ${\sigma }_x^2=1/N{\sum }_{n=0}^{N-1}‍{\left|x_j(k)\right|}^2$ .

We denote $\mathbf{S}=[{\mathbf{S}}_1,{\mathbf{S}}_2,\dots ,{\mathbf{S}}_{N_t}]$ , $\mathbf{h}={[{\mathbf{h}}_1^T,{\mathbf{h}}_2^T,\dots ,{\mathbf{h}}_m^T,\dots ,{\mathbf{h}}_{N_r}^T]}^T$ , ${\mathbf{h}}_m={[{\mathbf{h}}_{1,m}^T,{\mathbf{h}}_{2,m}^T,\dots ,{\mathbf{h}}_{N_t,m}^T]}^T$ , and $\mathbf{v}={[{\mathbf{v}}_1^T,{\mathbf{v}}_2^T,\dots ,{\mathbf{v}}_{N_r}^T]}^T$ . After stacking all the received vectors from different Rx branches ${\mathbf{r}}_j$ , $1\le j\le N_r$ , the new received vector, $\mathbf{r}={[{\mathbf{r}}_1^T,{\mathbf{r}}_2^T,\dots ,{\mathbf{r}}_{N_r}^T]}^T$ , is given by

$$\begin{array}{cc}& \mathbf{r}=\left({\mathbf{I}}_{N_r}\otimes \mathbf{\Gamma }\mathbf{S}\right)\mathbf{h}+\mathbf{v}=\mathbf{A}\mathbf{h}+\mathbf{v}\mathbf{,}\end{array}$$

(2)

where $\mathbf{A}={\mathbf{I}}_{N_r}\otimes \mathbf{\Gamma }\mathbf{S}$ . The mean of v is zero and the covariance matrix of v is ${\sigma }_v^2{\mathbf{I}}_{N{N}_r}$ .

3. ML Carrier Frequency Offset Estimation Algorithms

For known A and h, the received vector r is Gaussian with mean $\mathbf{A}\mathbf{h}$ and covariance matrix ${\sigma }_w^2{\mathbf{I}}_{N{N}_r}$ . The likelihood function of r is a $N{N}_r$ dimension complex Gaussian function:

$$\begin{array}{cc}& \Lambda (\mathbf{r};\mathbf{h},\nu )=\frac{1}{{\left(\pi {\sigma }_v^2\right)}^{N{N}_r}}\exp \left\{-\frac{{\parallel \mathbf{r}-\mathbf{A}\mathbf{h}\parallel }^2}{{\sigma }_v^2}\right\}.\end{array}$$

(3)

ML algorithm searches for the optimum value of h and ν, which can maximize the likelihood function in (3). The maximization is equivalent to the minimization of the cost function below

$$\begin{array}{cc}& C\left(\mathbf{r};\mathbf{h},\nu \right)={\parallel \mathbf{r}-\mathbf{A}\mathbf{h}\parallel }^2={\displaystyle \sum }_{j=1}^{N_r}‍{\parallel {\mathbf{r}}_j-\mathbf{\Gamma }\mathbf{S}{\mathbf{h}}_j\parallel }^2.\end{array}$$

(4)

Assuming ν is fixed, we take the first order partial derivative of $C(\mathbf{r};\mathbf{h},\nu )$ with respect to h and set it to zero. The least square (LS) estimation of CIR becomes [18]

$$\begin{array}{cc}& {\widehat{\mathbf{h}}}_j={\left(\mathbf{\Gamma }\mathbf{S}\right)}^{+}{\mathbf{r}}_j,\hspace{1em}1\le j\le N_r,\end{array}$$

(5)

where ${\mathbf{X}}^{+}$ is the Moore-Penrose generalized inverse of X. Substituting (5) into (4) and discarding the irrelevant terms, (4) becomes

$$\begin{array}{cc}& C\left(\mathbf{r};\nu \right)={\displaystyle \sum }_{j=1}^{N_r}‍{\mathbf{r}}_j^H\left(\mathbf{\Gamma }\mathbf{S}\right){\left(\mathbf{\Gamma }\mathbf{S}\right)}^{+}{\mathbf{r}}_j={\displaystyle \sum }_{j=1}^{N_r}‍{\mathbf{r}}_j^H\mathbf{\Gamma }\mathbf{B}{\mathbf{\Gamma }}^H{\mathbf{r}}_j,\end{array}$$

(6)

where $\mathbf{B}=\mathbf{S}{\mathbf{S}}^{+}$ is called project matrix. Notice that the projection matrix B is Hermitian symmetrical and $C(\mathbf{r};\nu )$ can be further derived as

$$\begin{array}{cc}& C\left(\mathbf{r};\nu \right)=\rho \left(0\right)-2\Re \left\{{\displaystyle \sum }_{l=0}^{N-1}‍\rho \left(l\right)e^{-j2\pi \nu l/N}\right\},\end{array}$$

(7)

where $\rho (l)$ is given by

$$\begin{array}{cc}& \rho \left(l\right)={\displaystyle \sum }_{j=1}^{N_r}‍\hspace{0.17em}{\displaystyle \sum }_{k=l}^{N-1}‍{\mathbf{B}}_{k-l,k}{r}_j\left(k\right)r_j^{*}(k-l).\end{array}$$

(8)

Compared with SISO systems, there are diversity effect in $\rho (l)$ . The ML CFO estimator for nonperiodic preambles is

$$\begin{array}{cc}& \widehat{\nu }=\arg \underset{\nu }{\max }\left\{\Re \left\{{\displaystyle \sum }_{l=0}^{N-1}‍\rho \left(l\right)e^{-j2\pi \nu l/N}\right\}\right\}.\end{array}$$

(9)

The last term of right hand of (7) has the same structure as Discrete Fourier Transform (DFT). So this term can be implemented by FFT effectively.

3.1. MLE-MA

As we know, for system with periodic preambles, its packet detectors and CFO estimator have simpler form and lower complexity [9, 10]. A periodic preamble has multiple identical short slots and each slot has a length of P. For simplicity, the length of the periodic sequence is first set to N, the same as the FFT size. The periodic sequence includes $M=N/P$ identical short slots. Actually, more repeated slots may be used in practice, such as in 802.11a/n [8], the short preamble contains ten identical short slots, and each slot has a length of 16. To generate periodic preambles, nonzero pilots are inserted with a fixed spacing M; namely, the entries of ${\mathbf{c}}_i$ are

$$\begin{array}{cc}& {\mathit{c}}_i\left(n\right)=\{\begin{array}{cc}{d}_i\left(n^{\prime }\right),\hfill & n=n^{\prime }M,\hfill \\ 0,\hfill & \text{otherwise},\hfill \end{array}\end{array}$$

(10)

where pilot $d_i(n^{\prime })$ is usually constant power modulation symbols with $\left|d_i(n^{\prime })\right|=\sqrt{1/(N_D{N}_t)}$ , and $n^{\prime }\in \{\varphi :[1,N_D/2]\cup [P-N_D/2,\hspace{0.17em}P-1]\}$ , where the number of pilots $N_D$ is an even number and $N_D\le \lfloor N_u/M\rfloor$ . In 802.11n, the frequency domain pilot sequences used in each Tx antenna are the same but with different cyclic delay. Now the pilot sequences ${\mathbf{c}}_i$ is represented as ${\mathbf{c}}_i={\mathbf{c}}_1\odot \mathbf{u}$ , where $\mathbf{u}={[1,e^{j2\pi {\tau }_i/N},\dots ,e^{j2\pi (N-1){\tau }_i/N}]}^T$ . The time delay of the ith Tx branch is ${\tau }_i$ . We have ${\tau }_1=0$ . In the case of multiple (>P) repeated slots, assuming that timing uncertainty μ existed, but there are enough samples to assure the receiving vector not containing other parts of the signal, thus the received signal vector at the jth Rx branch is given by

$$\begin{array}{cc}& {\mathbf{r}}_j=\mathbf{\Gamma }{\displaystyle \sum }_{i=1}^{N_t}‍{\mathbf{F}}_N^H\mathbf{\Gamma }\left({\tau }_i+\mu \right)\mathrm{diag}\left\{{\mathbf{c}}_1\right\}{\mathbf{F}}_{N,L}{\mathbf{h}}_{i,j}+{\mathbf{v}}_j,\end{array}$$

(11)

where $\mathbf{\Gamma }({\tau }_i+\mu )=\mathrm{diag}\{1,e^{j2\pi ({\tau }_i+\mu )/N},\dots ,e^{j2\pi ({\tau }_i+\mu )(N-1)/N}\}$ . Now, the transmit signal matrix at the ith antenna becomes ${\mathbf{S}}_i={\mathbf{F}}_N^H\mathbf{\Gamma }({\tau }_i+\mu )\mathrm{diag}\{{\mathbf{c}}_1\}{\mathbf{F}}_{N,L}$ .

We observe that $\mathrm{diag}\{{\mathbf{c}}_1\}$ is rank deficiency and ${\mathbf{F}}_N^H\mathrm{diag}\{{\mathbf{c}}_1\}$ has $N_D$ nonzero columns. We define a new N-by- $N_D$ column selection matrix E, with one at $(n^{\prime },i_{n^{\prime }})$ th entries and zero at other entries, to extract the valid columns of ${\mathbf{F}}_N^H\mathrm{diag}\{{\mathbf{c}}_1\}$ , where $i_{n^{\prime }}$ denotes the index of $n^{\prime }$ in ϕ. Correspondingly, ${\mathbf{E}}^T$ is a row selection matrix. Then, the total transmit signal matrix S can be expressed as $\mathbf{S}={\mathbf{F}}_N^H\mathrm{diag}\{{\mathbf{c}}_1\}\mathbf{E}{\mathbf{E}}^T[\mathbf{\Gamma }({\tau }_1+\mu ){\mathbf{F}}_{N,L},\dots ,\mathbf{\Gamma }({\tau }_{N_t}+\mu ){\mathbf{F}}_{N,L}]$ . The rank of matrix ${\mathbf{F}}_N^H\mathrm{diag}\{{\mathbf{c}}_1\}\mathbf{E}$ (denoted by U) and matrix ${\mathbf{E}}^T[\mathbf{\Gamma }({\tau }_1+\mu ){\mathbf{F}}_{N,L},\dots ,\mathbf{\Gamma }({\tau }_{N_t}+\mu ){\mathbf{F}}_{N,L}]$ (denoted by V) is $N_D$ . Both matrices constitute a full-rank decomposition of $\mathbf{S},\text{that is}\mathrm{\hspace{0.17em}\hspace{0.17em}}\mathbf{S}=\mathbf{U}\mathbf{V}$ , based on which we can use the following expression to calculate the Moore-Penrose generalized inverse ${\mathbf{S}}^{+}$ [19],

$$\begin{array}{cc}& {\mathbf{S}}^{+}={\mathbf{V}}^{+}{\mathbf{U}}^{+}={\mathbf{V}}^H{\left(\mathbf{V}{\mathbf{V}}^H\right)}^{-1}{\left({\mathbf{U}}^H\mathbf{U}\right)}^{-1}{\mathbf{U}}^H,\end{array}$$

(12)

and the projection matrix B ,

$$\begin{array}{cc}& \mathbf{B}=\mathbf{S}{\mathbf{S}}^{+}=\mathbf{U}{\left({\mathbf{U}}^H\mathbf{U}\right)}^{-1}{\mathbf{U}}^H=\mathbf{U}{\mathbf{U}}^{+},\end{array}$$

(13)

where we utilize two matrix properties, $\mathbf{V}{\mathbf{V}}^H{(\mathbf{V}{\mathbf{V}}^H)}^{-1}=\mathbf{I}$ , and the Moore-Penrose generalized inverse of N-by- $N_D$ matrix U with rank $N_D$ is ${\mathbf{U}}^{+}={({\mathbf{U}}^H\mathbf{U})}^{-1}{\mathbf{U}}^H$ . It is interesting to observe that B is independent of cyclic delay $\{{\tau }_i,\hspace{0.17em}\hspace{0.17em}1\le i\le N_t\}$ and timing uncertainty μ as well. If we proceed to another full-rank decomposition of  U as $\mathbf{U}=({\mathbf{F}}_N^H)(\mathrm{diag}\{{\mathbf{c}}_1\}\mathbf{E})$ and use the first equation in (12), the projection matrix B is further given by

$$\begin{array}{l}\mathbf{B}={\mathbf{F}}_N^H\left(\mathrm{diag}\left\{{\mathbf{c}}_1\right\}\mathbf{E}\right){\left(\mathrm{diag}\left\{{\mathbf{c}}_1\right\}\mathbf{E}\right)}^{+}{\mathbf{F}}_N\\ ={\mathbf{F}}_N^H\mathrm{diag}\left\{{\mathbf{c}}_1\right\}\mathrm{diag}{\left\{{\mathbf{c}}_1\right\}}^{+}{\mathbf{F}}_N\\ ={\mathbf{F}}_N^H\mathrm{diag}\left\{\mathbf{\Theta }\right\}{\mathbf{F}}_N\\ =\frac{1}{M}{\mathbf{1}}_{M\times M}\otimes {\mathbf{Q}}_P,\end{array}$$

(14)

where $\mathbf{\Theta }={[0,\dots ,\mathrm{0,1},0,\dots ,\mathrm{0,1},\dots ,0]}^T$ and Θ takes 1 at nonzero pilot position in frequency domain of the preamble, which is only related to the position and number of pilots and independent of the pilot modulation symbols. Let q be the time-domain transform of Θ, $\mathbf{q}={\mathbf{F}}_N^H\mathbf{\Theta }$ . q is periodic with period of P. ${\mathbf{Q}}_P$ is a circulant matrix and its entries are ${[{\mathbf{Q}}_P]}_{m,n}=q(m-n\hspace{0.17em}|\hspace{0.17em}P)$ , where $0\le m,n\le P-1$ . After some straightforward manipulations, $\rho (l)$ can be written as

$$\begin{array}{cc}& \rho \left(l\right)=\frac{1}{M}{q}_{l\hspace{0.17em}|\hspace{0.17em}P}{\displaystyle \sum }_{j=1}^{N_r}\hspace{0.17em}{\displaystyle \sum }_{k=l}^{N-1}‍r_j\left(k\right)r_j^{*}(k-l).\end{array}$$

(15)

It is obvious that (15) is the weighted correlation of the received signal. For $N=64$ , $P=16$ , and $N_D\ge 8$ , the weighting vector q is illustrated in Table 1. In 802.11a/n systems, $N_D=12$ . Spatial diversity effect is seen in (15) by accumulating correlation of the signals from multiple receiving antennas. The weighted correlation structure has a reduced complexity compared with (8). Moreover, B is correlated with $N_D$ and independent from the specific preamble. Such a special structure makes the estimator simpler. Thus, the first MLE for periodic preambles takes the following form:

$$\begin{array}{cc}& {\widehat{\nu }}_{\text{MA}}=\arg \underset{\nu }{\max }\left\{\Re \left\{{\displaystyle \sum }_{l=1}^{N-1}‍\rho \left(l\right)e^{-j2\pi \nu l/N}\right\}\right\}.\end{array}$$

(16)

Table 1

Pilots number and q vector.

$N_D$	q vector
8	[0.5000, 0.2517, −0.1250, −0.1560, 0, −0.0207, −0.1250, −0.0749, 0, −0.0749, −0.1250, −0.0207, 0, −0.1560, −0.1250, 0.2517 ${]}^T$
10	[0.6250, 0.2039, −0.2134, −0.0406, 0, −0.1362, −0.0366, −0.0271, −0.1250, −0.0271, −0.0366, −0.1362, 0, −0.0406, −0.2134, 0.2039 ${]}^T$
12	[0.7500, 0.1155, −0.2134, 0.0478, −0.1250, −0.0478, −0.0366, −0.1155, 0, −0.1155, −0.0366, −0.0478, −0.1250, 0.0478, −0.2134, 0.1155 ${]}^T$
14	[0.8750, 0, −0.1250, 0, −0.1250, 0, −0.1250, 0, −0.1250, 0, −0.1250, 0, −0.1250, 0, −0.1250, 0 ${]}^T$
16	[1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 ${]}^T$

Note that the Direct-Current (DC) component, $\rho (0)$ , can be removed from the estimation equation above because it contributes nothing to estimate the frequency offset. This estimator is a multiple-antenna CFO estimator, which is called MLE-MA. Although different cyclic delay may be applied at the Tx, this estimator has a definite and unique expression in (16). Similar to (9), this estimator can be implemented by FFT. After zero padding to $\rho (l)$ , the actual FFT length $K_{P}N$ can improve the accuracy and resolution [10], where $K_P$ is a designed factor.

If the index of maximum FFT value is λ, then the estimated CFO is $\widehat{\nu }=\lambda /K_P$ . The maximum value of FFT may not locate at the maximum of spectrum. Numerical method, such as Lagrangian function interpolation [20], can be used to further improve the accuracy.

A simple three-point quadratic interpolation is employed here. The three frequency values near index λ are $y(\lambda -1)=(\lambda -1)/K_P$ , $y(\lambda )=\lambda /K_P$ , and $y(\lambda +1)=(\lambda +1)/K_P$ . The corresponding FFT values are denoted as $f(\lambda -1)$ , $f(\lambda )$ , and $f(\lambda +1)$ . The Lagrangian interpolation function is

$$\begin{array}{l}{L}_2\left(y\right)=\{\frac{1}{2}\left(y-y\left(\lambda \right)\right)\left(y-y\left(\lambda +1\right)\right)f\left(\lambda -1\right)\\ \phantom{\text{K}}-(y-y(\lambda -1))(y-y(\lambda +1))f\left(\lambda \right)\\ \phantom{\text{K}}+\frac{1}{2}(y-y(\lambda -1))(y-y\left(\lambda \right))f(m+1)\}.\end{array}$$

(17)

Taking the derivation of $L_2(y)$ , $L_2^{\prime }(y)$ , to be 0, the refined CFO value is obtained as

$$\begin{array}{cc}& {\widehat{\nu }}_{\text{MA}}=\frac{(2\lambda +1)f(\lambda -1)-4\lambda f\left(\lambda \right)+(2\lambda -1)f(\lambda +1)}{2K_P(f(\lambda -1)-2f\left(\lambda \right)+f(\lambda +1))}.\end{array}$$

(18)

MLE-MA has an estimation range of $\left|{\nu }_{\text{MA}}\right|\le N/2$ , which is up to the entire signal bandwidth of OFDM.

3.2. MLE-MAS

For preamble which has a pilot number of $N_D=P$ , the project matrix B equals $(1/M){\mathbf{1}}_{M\times M}\otimes {\mathbf{I}}_P$ , which is shown in the last row of Table 1. $\rho (l)$ has nonzero values only when l is the multiple of P. Denote $\xi (l)$ as

$$\begin{array}{cc}& \xi \left(l\right)=\frac{1}{M}{\displaystyle \sum }_{j=1}^{N_r}\hspace{0.17em}{\displaystyle \sum }_{k=lP}^{N-1}‍r_j\left(k\right)r_j^{*}\left(k-lP\right),\end{array}$$

(19)

where $l=\mathrm{0,1},\dots ,M-1$ . Only M correlation values are needed, so the computation load is further reduced. A new CFO estimator becomes

$$\begin{array}{cc}& {\widehat{\nu }}_{\text{MAS}}=\arg \underset{\nu }{\max }\left\{\Re \left\{{\displaystyle \sum }_{l=0}^{M-1}‍\xi \left(l\right)e^{-j\left(2\pi \nu l/M\right)}\right\}\right\}.\end{array}$$

(20)

Similarly, the DC component $\xi (0)$ is removed because of zero contribution to estimate the frequency offset. This estimator above is named as MLE-MAS. Notice that when the number of receive antenna $N_r=1$ , MLE-MAS coincides with MLE #2 of [10] in SISO system. If the FFT maximum value index is ${\lambda }_{\text{MAS}}$ , estimated CFO is $\widehat{\nu }={\lambda }_{\text{MAS}}/(K_{P}P)$ . The same interpolation method used for MLE-MA can be used here to improve the accuracy of MLE-MAS. Finally, we note that MLE-MAS has an estimation range of $\left|{\nu }_{\text{MAS}}\right|\le M/2$ .

4. Best Linear Unbiased Estimators (BLUEs)

Usually, periodic sequences passing through multipath channels still have a periodic structure except for the phase rotation induced by CFO. Let $\xi (l)$ be the time delay autocorrelation at integer periods. It can be expanded as

$$\begin{array}{l}\xi \left(l\right)=e^{j\left(2\pi l\nu /M\right)}\\ \times \frac{1}{M}{\displaystyle \sum }_{j=1}^{N_r}‍\hspace{0.17em}{\displaystyle \sum }_{k=lP}^{N-1}\{{\left|x_j\left(k\right)\right|}^2+x_j\left(k\right){\tilde{v}}_j^{*}(k-lP)‍\\ \phantom{kkkkkkkkkkkkk}+\hspace{0.17em}{x}_j^{*}\left(k\right){\tilde{v}}_j\left(k\right)+{\tilde{v}}_j\left(k\right){\tilde{v}}_j^{*}(k-lP)\phantom{{\left|x_j\left(k\right)\right|}^2}\}\\ =e^{j\left(2\pi l\nu /M\right)}\frac{1}{M}{\sigma }_x^2\left\{N_r(N-lP)(1+{\gamma }_l)\right\},\end{array}$$

(21)

where $l=\mathrm{0,1},\dots ,M-1$ and the equivalent noise term is ${\tilde{v}}_j(k)=v_j(k)e^{j2\pi k\nu /N}$ . The components excluding signal in (21) are viewed as interference ${\gamma }_l$ , which is given by

$$\begin{array}{l}{\gamma }_l=\frac{1}{N_r(N-lP){\sigma }_x^2}\\ \times {\displaystyle \sum }_{j=1}^{N_r}‍\hspace{0.17em}{\displaystyle \sum }_{k=lP}^{N-1}\{x_j\left(k\right){\tilde{v}}_j^{*}(k-lP)‍\\ \phantom{kkkkkkkkkk}+\hspace{0.17em}{x}_j^{*}\left(k\right){\tilde{v}}_j\left(k\right)+{\tilde{v}}_j\left(k\right){\tilde{v}}_j^{*}(k-lP)\}.\end{array}$$

(22)

Due to the fact that $\xi (l)$ is the autocorrelation of repeated sequences, the argument of $\xi (l)$ is mainly affected by $lP$ times of CFO with mean $E\{\arg (\xi (l))\}=2\pi l\nu /M$ , so this characteristic can be exploited to estimate CFO. A typical estimator is BLUE, which is firstly proposed by Kay [18].

BLUE is linear, unbiased and the CFO estimation variance is minimized. For CFO estimation, BLUE has the form of

$$\begin{array}{cc}& \widehat{\nu }=\frac{M}{2\pi }{\displaystyle \sum }_{l=1}^H‍\varphi \left(l\right)w\left(l\right),\end{array}$$

(23)

where $\varphi (l)$ is the lth CFO estimation derived from corresponding estimation metric, $w(l)$ is the entry of weighting coefficient vector w, and $\mathbf{w}={[w(1),w(2),\dots ,w(H)]}^T$ . w is calculated as

$$\begin{array}{cc}& \mathbf{w}=\frac{{\mathbf{C}}_{\varphi }^{-1}\mathbf{1}}{{\mathbf{1}}^T{\mathbf{C}}_{\varphi }^{-1}\mathbf{1}},\end{array}$$

(24)

where $\varphi ={[\varphi (1),\varphi (2),\dots ,\varphi (H)]}^T$ and ${\mathbf{C}}_{\varphi }$ is the covariance matrix of ϕ. $\mathbf{1}={[\mathrm{1,1},\dots ,1]}^T$ is a H-by-1 all-ones vector, and H is a designed parameter to estimate CFO. The variance of BLUE is

$$\begin{array}{cc}& \text{var}\left(\nu \right)=\frac{M^2}{4{\pi }^2}\frac{1}{{\mathbf{1}}^T{\mathbf{C}}_{\varphi }^{-1}\mathbf{1}}.\end{array}$$

(25)

4.1. BLUE-PA

The simplest method is taking the phase of $\xi (l)$ directly to estimate CFO. Because of the limitation of $\left|\arg (\xi (l))\right|\le \pi$ , the estimation range of this algorithm decreases with the increases of  l. The largest estimation range is $\left|\nu \right|\le M/2$ when $l=1$ and the total range depends on the maximum of l, which is $\left|\nu \right|\le M/(2(M-1))$ . In order to enlarge the range to $\left|\nu \right|\le M/2$ , a coarse estimation is first made as $l=1$ . This coarse estimation value is denoted as ${\widehat{\nu }}_1=(M/2\pi )\arg (\xi (1))$ . Multiplying a counter-rotation factor $e^{-j2\pi l{\widehat{\nu }}_1/M}$ to each $\xi (l)$ , we obtain ${\xi }^{\prime }(l)=\xi (l)e^{-j2\pi l{\widehat{\nu }}_1/M}$ . After this coarse compensation, the estimation metric is given by

$$\begin{array}{cc}& {\varphi }_1\left(l\right)=\frac{1}{l}\arg \left({\xi }^{\prime }\left(l\right)\right)=\frac{2\pi \tilde{\nu }}{M}+\frac{1}{l}\arctan \left(\frac{{\gamma }_{I,l}}{1+{\gamma }_{R,l}}\right),\end{array}$$

(26)

where $\tilde{\nu }=\nu -{\widehat{\nu }}_1$ and ${\gamma }_{R,l}$ and ${\gamma }_{I,l}$ are the real and imaginary part of ${\gamma }_l$ . At moderate and high SNR values, ${\gamma }_{R,l}\ll 1$ ; then ${\varphi }_1(l)$ is approximated as

$$\begin{array}{cc}& {\varphi }_1\left(l\right)\approx \frac{2\pi \tilde{\nu }}{M}+\frac{1}{l}\arctan \left({\gamma }_{I,l}\right)\approx \frac{2\pi \tilde{\nu }}{M}+\frac{1}{l}{\gamma }_{I,l}.\end{array}$$

(27)

The mean of ${\varphi }_1(l)$ is $2\pi \tilde{\nu }/M$ , so the estimator is unbiased. The covariance matrix ${\mathbf{C}}_{{\varphi }_1}$ is ${\mathbf{C}}_{{\varphi }_1}(m,l)=(1/ml)c(m,l)$ , where $c(m,l)=E\{{\gamma }_{I,m}{\gamma }_{I,l}\}$ is the covariance of ${\gamma }_{I,m}$ and ${\gamma }_{I,l}$ , which is expanded as in

$$\begin{array}{l}c\left(m,l\right)=\frac{1}{N_r^2{\sigma }_x^4{P}^2(M-m)(M-l)}\\ \times {\displaystyle \sum }_{j=1}^{N_r}‍E\{\Im \{\{\phantom{{\displaystyle \sum }_{p=lP}^{N-1}‍}{\displaystyle \sum }_{k=mP}^{N-1}‍x_j\left(k\right){\tilde{v}}_j^{*}\left(k-mP\right)\\ \phantom{kkkkkkkkkkkkkkkk}+x_j^{*}\left(k\right){\tilde{v}}_j\left(k\right)\\ \phantom{kkkkkkkkkkkkkkkk}+{\tilde{v}}_j\left(k\right){\tilde{v}}_j^{*}\left(k-mP\right)\phantom{{\displaystyle \sum }_{p=lP}^{N-1}‍}\}\\ \phantom{kkkkkkkkkkkkk}\times \{{\displaystyle \sum }_{p=lP}^{N-1}‍x_j\left(p\right){\tilde{v}}_j^{*}\left(p-lP\right)\\ \phantom{\text{kkkkkkkkkkkkkkkkkkkk}}+x_j^{*}\left(p\right){\tilde{v}}_i\left(p\right)\\ \phantom{\text{kkkkkkkkkkkkkkkkkkk}}+{\tilde{v}}_j\left(p\right){\tilde{v}}_j^{*}\left(p-lP\right)\phantom{{\displaystyle \sum }_{p=lP}^{N-1}‍}\}\}\}\\ =\frac{1}{\text{SNR}{N}_{r}P\left(M-m\right)\left(M-l\right)}\\ \times \{\begin{array}{ll}m+(M-2m)u(\frac{M}{2}-m)\hfill & \hfill \\ \hspace{1em}+\frac{M-m}{\text{SNR}},\hfill & \text{if}\hspace{0.17em}\hspace{0.17em}m=l,\hfill \\ M-\max (m,l)\hfill & \hfill \\ \hspace{1em}-\left(M-m-l\right)u\left(M-m-l\right),\hfill & \text{if}\hspace{0.17em}\hspace{0.17em}m\ne l.\hfill \end{array}\end{array}$$

(28)

Matrix ${\mathbf{C}}_{{\varphi }_1}$ is full rank and $H=M-1$ . The inversion of ${\mathbf{C}}_{{\varphi }_1}$ cannot be represented in a closed form. However, we observe that the numerator and denominator in (24) are the summation of each row and the summation of all elements of entire matrix, respectively. We denote ${\mathbf{a}}_1={\mathbf{C}}_{{\varphi }_1}^{-1}\mathbf{1}$ and $b_1={\mathbf{1}}^T{\mathbf{C}}_{{\varphi }_1}^{-1}\mathbf{1}$ , respectively. After some lengthy computation, ${\mathbf{a}}_1$ and $b_1$ can be precisely solved as

$$\begin{array}{cc}& a_1\left(l\right)=\frac{2N_r\text{SN}{\text{R}}^{2}P{l}^2\left(M-l\right)}{M\text{SNR}+1},\hspace{1em}l=\mathrm{1,2},\dots ,M-1,\\ b_1=\frac{N_r\text{SN}{\text{R}}^{2}P{M}^2(M^2-1)}{6(M\text{SNR}+1)}.\end{array}$$

(29)

The entries of weighting coefficient vector w can be represented in a closed-form expression as

$$\begin{array}{cc}& w_1\left(l\right)=\frac{12l^2(M-l)}{M^2(M^2-1)}.\end{array}$$

(30)

Both ${\mathbf{a}}_1$ and $b_1$ are correlated with SNR. However, they have common components in the denominators, which makes the weighting coefficients uncorrelated with SNR. This estimator is named as BLUE-PA for utilizing the phase of autocorrelation. Its variance is

$$\begin{array}{cc}& \text{va}{\text{r}}_1\left(\nu \right)=\frac{3(M\text{SNR}+1)}{2{\pi }^2{N}_r\text{SN}{\text{R}}^{2}P(M^2-1)}.\end{array}$$

(31)

$\text{va}{\text{r}}_1(\nu )$ is a function of SNR. At high SNR values, it can be simplified as

$$\begin{array}{cc}& \overline{\text{var}}\left(\nu \right)=\frac{3M}{2{\pi }^2{N}_r\text{SNR}P(M^2-1)}.\end{array}$$

(32)

The approximate variance $\overline{\text{var}}$ coincides with the theoretical Cramer-Rao lower bound (CRLB) given in [14]. Considering there is a coarse estimation before BLUE-PA, the actual CFO estimation value takes ${\nu }_{\text{actual}}={\widehat{\nu }}_1+\widehat{\nu }$ .

4.2. BLUE-PA-N

In this part, we make a further approximation in (21). By discarding the noise-noise product ${\tilde{v}}_j(k){\tilde{v}}_j^{*}(k-lP)$ , ${\gamma }_{I,l}$ can be approximated as ${\tilde{\gamma }}_{I,l}$

$$\begin{array}{l}{\tilde{\gamma }}_{I,l}\\ \hspace{0.5em}=\frac{1}{N_r{\sigma }_x^{2}P(M-m)}\\ \hspace{1em}\hspace{0.5em}\times {\displaystyle \sum }_{j=1}^{N_r}\{{\displaystyle \sum }_{k=m}^{N-1}‍‍\left\{\Re \left\{x_j\left(k\right)\right\}\left\{\Im \left\{{\tilde{v}}_j\left(k\right)\right\}-\Im \left\{{\tilde{v}}_j^{*}\left(k-mP\right)\right\}\right\}\right\}\\ \phantom{kkkkkkk}\phantom{\text{k}}\hspace{1em}\hspace{0.17em}\hspace{0.17em}+\{\phantom{-\Re \left\{{\tilde{v}}_j^{*}\left(k\right)\right\}\}}\Im \left\{x_i\left(k\right)\right\}\{\phantom{\left\{{\tilde{v}}_j^{*}\left(k\right)\right\}}\Re \left\{{\tilde{v}}_j\left(k-mP\right)\right\}\\ \phantom{kkkkkkkkkkkkkkkkkkkk}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}-\hspace{0.17em}\Re \left\{{\tilde{v}}_j^{*}\left(k\right)\right\}\}\}\phantom{{\displaystyle \sum }_{k=m}^{N-1}‍}\}.\end{array}$$

(33)

Taking the same step as in BLUE-PA, a coarse estimation is carried out for $l=1$ . After coarse CFO compensation, the estimation metric is

$$\begin{array}{cc}& {\varphi }_2\left(l\right)=\frac{2\pi \tilde{\nu }}{M}+\frac{1}{l}{\tilde{\gamma }}_{I,l},\end{array}$$

(34)

The mean of ${\varphi }_2$ is $2\pi \tilde{\nu }/M$ , so the estimator is also unbiased. The covariance between ${\tilde{\gamma }}_{I,m}$ and ${\tilde{\gamma }}_{I,l}$ is

$$\begin{array}{cc}& \tilde{c}(m,l)=\frac{M-\max (m,l)-(M-m-l)u(M-m-l)}{\text{SNR}{N}_r^2{P}^2(M-m)(M-l)}.\end{array}$$

(35)

Notice that the mth row/column of ${\mathbf{C}}_{{\varphi }_2}$ is the same as its $(M-m)$ th row/column, so ${\mathbf{C}}_{{\varphi }_2}$ does not have full rank. For $m,l>M/2$ , ${\mathbf{C}}_{{\varphi }_2}^{-1}$ is singular, then a designed parameter H is selected so that the first $H\times H$ element of ${\mathbf{C}}_{{\varphi }_2}$ is full rank, where $K\le M/2$ . We have

$$\begin{array}{cc}& \tilde{c}(m,l)=\frac{\min (m,l)}{\text{SNR}{N}_r^2{P}^2(M-m)(M-l)}.\end{array}$$

(36)

The covariance matrix ${\mathbf{C}}_{{\varphi }_2}$ is Hermitian. The inversion lemma of Hermitian matrix [21] is used to solve ${\mathbf{C}}_{{\varphi }_2}^{-1}$ . We can readily check that the entries of ${\mathbf{C}}_{{\varphi }_2}^{-1}$ are given by

$$\begin{array}{l}{\left[{\mathbf{C}}_{{\varphi }_2}^{-1}\right]}_{m,l}\\ \hspace{1em}=\{\begin{array}{ll}-N_{r}mlP\left(M-m\right)\left(M-l\right)\text{SNR},\hfill & m=l+1,\hfill \\ 2N_r{m}^{2}P{\left(M-m\right)}^2\text{SNR},\hfill & m=l\ne H,\hfill \\ N_r{m}^{2}P{\left(M-m\right)}^2\text{SNR},\hfill & m=l=H,\hfill \\ 0,\hfill & \text{otherwise}.\hfill \end{array}\end{array}$$

(37)

Then the weighting coefficient is obtained to be

$$\begin{array}{cc}& w_2\left(l\right)=\{\begin{array}{cc}\frac{6\left(M-l\right)l}{H\left(4H^2-6MH+3M^2-1\right)},\hfill & 1\le l<H,\hfill \\ \frac{3(M-H)(M-2H+1)}{(4H^2-6MH+3M^2-1)},\hfill & l=H.\hfill \end{array}\end{array}$$

(38)

The second estimator is named as BLUE-PA-N for utilizing the phase of autocorrelation without considering the noise-noise product. Its variance is

$$\begin{array}{cc}& \text{va}{\text{r}}_2\left(\nu \right)=\frac{3M^2}{4{\pi }^2{N}_{r}PH(4H^2-6MH+3M^2-1)\text{SNR}}.\end{array}$$

(39)

Taking the derivative of $\text{va}{\text{r}}_2(\nu )$ and setting it to 0, the minimization is checked at $H=M/2$ . When $H=M/2$ , $\text{va}{\text{r}}_2(\nu )$ is identical with (32). Similar to BLUE-PA, the actual CFO estimation value of BLUE-PA-N should also be ${\nu }_{\text{actual}}={\widehat{\nu }}_1+\widehat{\nu }$ .

4.3. BLUE-DPA

In order to avoid the coarse estimation steps in BLUE-PA and BLUE-PA-N while keeping the same estimation range, we take the differentiation between adjacent $\xi (l)$ . Define $R(l)$ as

$$\begin{array}{l}R\left(l\right)=\xi \left(l\right)\xi {\left(l-1\right)}^{*}\\ =e^{j\left(2\pi \nu /M\right)}\frac{1}{M^2}{\sigma }_x^4{N}_r^2(N-lP)(N-(l-1)P)\\ \times \left(1+{\gamma }_l\right){\left(1+{\gamma }_{l-1}\right)}^{*},\end{array}$$

(40)

where $l=\mathrm{1,2},\dots ,M-1$ , and the differential phase is represented as the third estimation metric

$$\begin{array}{cc}& {\varphi }_3\left(l\right)=\arg \left\{R\left(l\right)\right\}\approx \frac{2\pi \nu }{M}+{\gamma }_{I,l}-{\gamma }_{I,l-1}.\end{array}$$

(41)

The mean value of ${\varphi }_3(l)$ is $2\pi \nu /M$ , so this estimator is unbiased, too. The covariance matrix ${\mathbf{C}}_{{\varphi }_3}$ is given by

$$\begin{array}{l}{\mathbf{C}}_{{\varphi }_3}\left(m,l\right)=E\left\{\left[{\gamma }_{I,m}-{\gamma }_{I,m-1}\right]\left[{\gamma }_{I,l}-{\gamma }_{I,l-1}\right]\right\}\\ =c(m,l)+c(m-1,l-1)\\ -c\left(m-1,l\right)-c\left(m,l-1\right).\end{array}$$

(42)

The covariance matrix ${\mathbf{C}}_{{\varphi }_3}$ is a full rank matrix, and $H=M-1$ . Similar to BLUE-PA, the inversion of ${\mathbf{C}}_{{\varphi }_3}$ is difficult to be solved. But the numerator and denominator can be represented as simple expressions through computation. Denoting ${\mathbf{a}}_3={\mathbf{C}}_{{\varphi }_3}^{-1}\mathbf{1}$ and $b_3={\mathbf{1}}^T{\mathbf{C}}_{{\varphi }_3}^{-1}\mathbf{1}$ , we have

$$\begin{array}{l}{a}_3\left(l\right)=\frac{2N_{r}P\text{SN}{\text{R}}^{2}l\left(M-l\right)\left(M-l+1\right)}{\left(M+2\right)\text{SNR}+1},\\ \phantom{kkkkk}\phantom{kkkkkk}l=\mathrm{1,2},\dots ,M-1,\\ b_3=\frac{N_{r}P\text{SN}{\text{R}}^{2}M(M^2-1)(M+2)}{6((M+2)\text{SNR}+1)}.\end{array}$$

(43)

The third estimator is named as BLUE-DPA for utilizing differential phase of autocorrelation. The weighting coefficients are

$$\begin{array}{cc}& w_3\left(l\right)=\frac{12l(M-l)(M-l+1)}{M(M^2-1)(M+2)},\end{array}$$

(44)

which is also closed form and not related to SNR. BLUE-DPA has a variance of

$$\begin{array}{cc}& \text{va}{\text{r}}_3\left(\nu \right)=\frac{3M((M+2)\text{SNR}+1)}{2{\pi }^{2}P\text{SN}{\text{R}}^2(M^2-1)(M+2)}.\end{array}$$

(45)

At large SNR values, $\text{va}{\text{r}}_3(\nu )$ also approaches the CRLB asymptotically.

4.4. BLUE-DPA-N

Discarding the noise-noise product terms in ${\gamma }_l$ , the estimation metric of the fourth estimator is represented as

$$\begin{array}{cc}& {\varphi }_4\left(l\right)=\arg \left\{R\left(l\right)\right\}\approx \frac{2\pi \nu }{M}+{\tilde{\gamma }}_{I,l}-{\tilde{\gamma }}_{I,l-1}.\end{array}$$

(46)

The mean of ${\varphi }_4(l)$ is $2\pi \nu /M$ and estimator is still unbiased. The covariance matrix ${\mathbf{C}}_{{\varphi }_4}$ is given by

$$\begin{array}{l}{\mathbf{C}}_{{\varphi }_4}\left(m,l\right)=\tilde{c}(m,l)+\tilde{c}(m-1,l-1)-\tilde{c}(m-1,l)\\ -\tilde{c}(m,l-1).\end{array}$$

(47)

Similar to BLUE-PA-N, ${\mathbf{C}}_{{\varphi }_4}$ is not full rank for $m,l>M/2$ , so a design parameter H is also chosen here. Using the matrix inversion lemma [21], the inversion of ${\mathbf{C}}_{{\varphi }_4}$ is found as

$$\begin{array}{cc}& {\left[{\mathbf{C}}_{{\varphi }_4}^{-1}\right]}_{m,n}=\{\begin{array}{cc}{N}_{r}P\left(M-H\right)\text{SNR},\hfill & \text{if}\mathrm{\hspace{0.17em}\hspace{0.17em}}m\ne n,\hfill \\ N_{r}P\{M-H+(H-m)\hfill & \hfill \\ \phantom{\text{kkkk}}\times (H-m+1)\}\text{SNR},\hfill & m=n.\hfill \end{array}\end{array}$$

(48)

The fourth estimator is named as BLUE-DPA-N for utilizing differential phase of autocorrelation without considering the noise-noise product. The weighting coefficient of BLUE-DPA-N is

$$\begin{array}{cc}& w_4\left(l\right)=\frac{3(M-l)(M-l+1)-H(M-H)}{H(4H^2-6HM+3M^2-1)}.\end{array}$$

(49)

This estimator, namely BLUE-DPA-N, can be seen as the multiple-antenna extension of its SISO form in [15]. It can be easily validated that BLUE-DPA-N has the same variance as BLUE-PA-N, which is given in (39). This indicates that at high SNR, four estimators have asymptotically consistent variances, which asymptotically approaches the CRLB in (32).

Compared with MLEs, the four BLUEs have lower complexity with certain performance loss, but the estimation range is the same as MLE-MAS. BLUE-PA takes the phase of autocorrelation between periodic sequences with different delays directly, at the same time, noise-noise product components are taken into account in covariance derivation. When the noise-noise product effects are discarded, the estimator becomes BLUE-PA-N. Coarse estimation and compensation is introduced in the first two BLUEs to enlarge the estimation range. An alternative to attain this range is to consider differential phase. BLUE-DPA is based on differential phase and the noise-noise products terms are considered. However, BLUE-DPA-N discards the noise-noise product components.

It is worth noting that the length of the periodic sequence used for estimation is not limited to the OFDM symbol length. More repeated slots may be utilized to obtain better MSE performance or accuracy as desired.

5. Complexity Comparison

In this section, we evaluate and compare the computational complexities of the proposed frequency estimators with those of the existing estimators.

In MLE-MA, correlation and diversity are utilized, and the weighting coefficients for correlation calculation are only correlated with pilots number, which can be calculated offline. The two MLEs are based on FFT implementation and numerical interpolation. The main computational complexity of MLEs is caused by correlation calculation and FFT implementation, while numerical interpolation can be neglected compared with them. For $\rho (l)$ , $l=\mathrm{1,2},\dots ,N-1$ , there are $N_r(N-1)N/2$ complex products, $2N$ real products, and $N_r(N-1)(N-2)/2+(N_r-1)(N-1)$ complex additions, while $\xi (l)$ , $l=\mathrm{1,2},\dots ,M-1$ , need $N_r(M-1)N/2$ complex products and $N_r(M-1)N/2-M+1$ complex additions. When computing FFT, $\eta K_{P}N/2{\log }_2(K_{P}N)$ complex products and $\eta K_{P}N/2{\mathrm{\hspace{0.17em}\hspace{0.17em}}\log }_2(K_{P}N)$ complex additions are introduced, where factor η represents the computational saving gained from skipping the operations on the zeros in the FFT [10].

Though the coarse frequency estimation is required for BLUE-PA, we can take differentiating $\xi (l),l=\mathrm{2,3},\dots ,M-1$ , with $\xi (1)$ and its power series followed by phase angle calculations to replace the operations of coarse estimation and compensation. Thus the steps are similar to those in the BLUE-DPA. Keep in mind that BLUE-PA-N and BLUE-DPA-N only need to calculate first $M/2\hspace{0.17em}\hspace{0.17em}\xi (l)$ . The differentiations take $2(M-2)$ , $2(M/2-1)$ , $M-2$ , and $M/2-1$ complex products for BLUE-PA, BLUE-DPA, BLUE-PA-N, and BLUE-DPA-N, respectively. In hardware implementation, the phase angle calculation can be conducted by CORDIC algorithm [22]. BLUE-PA and BLUE-DPA require $M-1$ CORDIC computations, while BLUE-PA-N and BLUE-DPA-N require $M/2$ CORDIC computations. Extra computation consumed on the weighted adding of individual estimation, namely (23), is easier to obtain.

Meanwhile, by applying simple average operations on the phases ${\varphi }_2(l)$ in (34) and their differentiation ${\varphi }_4(l)$ in (46), we can obtain two Averaged estimators, named as Averaged-PA-N and Averaged-DPA-N, respectively:

$$\begin{array}{cc}& {\widehat{\nu }}_{\text{Averaged-PA-N}}={\widehat{\nu }}_1+\frac{M}{2\pi H}{\displaystyle \sum }_{l=1}^H{\varphi }_2\left(l\right),\end{array}$$

(50a)

$$\begin{array}{cc}& {\widehat{\nu }}_{\text{Averaged-DPA-N}}=\frac{M}{2\pi H}{\displaystyle \sum }_{l=1}^H{\varphi }_4\left(l\right).\end{array}$$

(50b)

Besides the computation of $\xi (l)$ , the Averaged-PA-N takes extra $2(M/2-1)$ complex products, $M/2$ real products, and $M/2-1$ real additions. While Averaged-DPA-N takes $M/2-1$ complex products, one real products, and $M/2-1$ real additions.

The CFO estimator based on the trigonometric approximation of likelihood function, namely, Morelli's algorithm, is given as [14]

$$\begin{array}{cc}& {\widehat{\nu }}_{\text{Morelli}}={\widehat{\nu }}_1+\frac{{\sum }_{j=1}^{N_r}\hspace{0.17em}{\sum }_{l=1}^{M-1}l\left|{{\rm Y}}_j\left(l\right)\right|{\psi }_j\left(l\right)}{2\pi {\sum }_{j=1}^{N_r}\hspace{0.17em}{\sum }_{l=1}^{M-1}{l}^2\left|{{\rm Y}}_i\left(l\right)\right|},\end{array}$$

(51)

where ${\widehat{\nu }}_1$ is the coarse estimation the same as in BLUE-PA, ${{\rm Y}}_j(l)={\sum }_{k=lP}^{N-1}‍r_j(k)r_j^{*}(k-lP)$ , and ${\psi }_j(l)=\arg \{{{\rm Y}}_j(l)e^{-j2\pi l{\widehat{\nu }}_1/M}\}$ .

The complexities of the proposed MLEs, BLUEs, Averaged methods and Morelli's algorithm [14] are outlined in Table 2. As shown in this Table, the two MLEs have the highest computational load, while MLE-MA has higher complexity than MLE-MAS. The proposed BLUEs have lower complexity than Morelli's algorithm because they require less CORDIC computations. BLUE-PA-N and BLUE-DPA-N algorithm have slightly lower complexity than BLUE-PA and BLUE-DPA. Averaged-PA-N has the same complexity as BLUE-PA-N, while Averaged-DPA-N is the simplest one.

Table 2

Computational complexity comparison for different algorithms.

CFO estimator	Real products	Real additions	Phase derivations	Absolute values
MLE-MA	$2N_{r}N(N-1)+2N+2C$	$2\left(N_{r}N-1\right)\left(N-1\right)+3C$	0	0
MLE-MAS	$M_{\xi ,1}+2C$	$A_{\xi ,1}+3C$	0	0
BLUE-PA	$M_{\xi ,1}+9M-16$	$A_{\xi ,1}+5M-9$	$M-1$	0
BLUE-PA-N	$M_{\xi ,2}+9M/2-8$	$A_{\xi ,2}+5M/2-5$	$M-1$	0
BLUE-DPA	$M_{\xi ,1}+5M-8$	$A_{\xi ,1}+3M-5$	$M/2$	0
BLUE-DPA-N	$M_{\xi ,2}+5M/2-4$	$A_{\xi ,2}+3M/2-3$	$M/2$	0
Averaged-PA-N	$M_{\xi ,2}+9M/2-8$	$A_{\xi ,2}+5M/2-5$	$M/2$	0
Averaged-DPA-N	$M_{\xi ,2}+2M-3$	$A_{\xi ,2}+3M/2-3$	$M/2$	0
Morelli	$N_r(2N(M-1)+8(M-1))$	$2N_r(N(M-1)+M-2)$	$N_r(M-1)+1$	$N_r(M-1)$

Consider $C=\eta K_{P}N{\log }_{\text{2}}(K_{P}N)$ , $M_{\xi ,1}=2N_{r}N(M-1)$ , and $A_{\xi ,1}=2(M-1)(N_{r}N-1)$ ;

$M_{\xi ,2}=N_{r}N(3M/2-1/2)$ , $A_{\xi ,2}=N_{r}N(3M/4-1/2)-M$ .

6. Simulation Results

In this paper, the legacy short training field (L-STF) preamble in 802.11n system is taken as an example in our simulations. It has ten identical short slots in a 20 M bandwidth. There are 16 samples in each slot. Cyclic delay is employed for different transmit antennas to avoid unintentional beamforming.

The multipath channel between each pair of Tx and Rx branches is L paths Rayleigh channel with quasistatic fading. Tap coefficients are kept constant in a frame but may vary from frames to frames. The channel power delay profile decays exponentially. That is, $E\{{\left|h_{i,j}(l_c)\right|}^2\}=\beta \exp (-l_c/L)$ , where $l_c=\mathrm{0,1},\dots ,L-1$ and β is a factor to keep the total channel power normalized to 1. In this paper, the channel has 6 independent paths; that is, $L=6$ .

The theoretical CRLB is calculated based on [23]

$$\begin{array}{cc}& \text{CRLB}=\frac{{\sigma }_v^2{N}^2}{8{\pi }^2{\sum }_{j=1}^{N_r}{\mathbf{h}}_j{\mathbf{S}}^H\mathbf{\Lambda }{\mathbf{\Pi }}_S\mathbf{\Lambda }\mathbf{S}{\mathbf{h}}_j^H},\end{array}$$

(52)

where $\mathbf{\Lambda }=\mathrm{diag}\{\mathrm{0,1},\dots ,N-1\}$ , ${\mathbf{\Pi }}_S$ is the orthogonal projection matrix of S and ${\mathbf{\Pi }}_S={\mathbf{I}}_N-\mathbf{B}$ . As CRLB depends on h, we obtain CRLB by averaging (52) on $10^5$ realizations of the channel. Zeros are padded to channel realizations to assure $L=N_D$ . Usually, virtual subcarriers are inserted to high frequency components of the training sequence and CRLB in (52) is suitable for these training sequences, but CRLB in [14] is for the special case without virtual subcarriers.

The mean square error (MSE) of MLE-MA is shown in Figure 1. The numbers in the bracket are marked as $(N_t,N_r)$ . CFO is randomly generated within interval $[-\mathrm{30,30}]$ with uniform distribution. In this figure, $K_P=8$ , $M=4$ , and the actual FFT length is 512. It is observed that MLE-MA approaches CRLB at SNRs of 22 dB, 12 dB, and 4 dB. Configurations with more antennas can reach the CRLB quicker. We find that MIMO systems can obtain diversity gains compared to SISO systems, which have been clearly shown in (15). This estimator needs a large $K_P$ to improve accuracy. For very large $K_P$ , the estimator becomes impractical. However, this estimator can cope with very large CFO up to the whole signal bandwidth of OFDM. From Figure 1, we observe that the MSE curves fall sharply and reach the CRLBs at certain SNRs. The reason is that there exists integer frequency ambiguity at some channel realizations, which is of rare occurrence, to make the MSEs departing from the CRLBs. Moreover the occurrence times incline to zero as SNR increases due to the limited realization of channel in the simulation. Frequency ambiguity problem is the topic of integer frequency offset estimation and it is beyond the topic of this paper.

Figure 1

MSE performance for MLE-MA in MIMO-OFDM systems with different M, $N_t$ , and $N_r$ ( $N_t,N_r=1,2,4$ , $N=64$ , $M=4$ , $K_P=8$ , and the FFT length is 512).

The performance of MLE-MAS is illustrated in Figure 2. In this figure, $M=4$ , $K_P=0.5$ , and the actual FFT length is 32. CFO is randomly generated within interval $[-1.9,1.9]$ . MLE-MA has an estimation range of $\left|\nu \right|\le N/2$ , while MLE-MAS has a range of $\left|\nu \right|\le M/2$ . Parameter ${\rho }_B(l)$ is decimated from ${\rho }_A(l)$ with a fixed length P, then FFT resolution of MLE-MAS is $1/P$ of MLE-MA. Small $K_P$ may behave excellently. Due to the smaller CFO estimation range, this algorithm has a smaller fluctuation and MSE. After linear interpolation, all three antenna configurations can approach CRLB at 14 dB, 8 dB, and 6 dB, respectively.

Figure 2

MSE performance for MLE-MAS in MIMO-OFDM systems with different M, $N_t$ , and $N_r$ ( $N_t,\hspace{0.17em}\hspace{0.17em}{N}_r=1,2,4$ , $N=64$ , $M=4$ , $K_P=0.5$ , and the FFT length is 32).

Performance comparison of MLE-MAS and Morelli's method when $N_t=2$ and $N_r=2$ is shown in Figure 3. In MLE-MAS, the K factor is set to 1 with a FFT length 64. It can be seen that MLE-MAS approaches CRLB at 4 dB and 2 dB, and Morelli's method approaches CRLB at the same SNR of 4 dB for both $M=4$ and $M=6$ , while MLE-MAS performs better than Morelli's algorithm below 4 dB.

Figure 3

MSE performance comparison between MLE-MAS and Morelli's estimator in MIMO-OFDM systems ( $N_t=2$ , $N_r=2$ , $N=64$ , and $M=4,6$ ).

The mean square error (MSE) of BLUE-PA is shown in Figure 4. BLUE-PA approaches CRLB at 8 dB, 6 dB, and 4 dB for different antennas configurations. Due to the presence of the diversity effect, larger antenna configuration can get lower MSE and approach CRLB more quickly.

Figure 4

MSE performance of BLUE-PA in MIMO-OFDM systems with three different antenna configurations ( $N_t,N_r=\mathrm{1,2},4$ , $N=64$ , and $M=4$ ).

MSE comparison between two MLEs, four BLUEs, two Averaged methods, and Morelli's algorithm [14] for a 2-by-2 MIMO system is shown in Figure 5. Here, H is set to $M/2=2$ for BLUE-PA-N, BLUE-DPA-N, and the two Averaged methods. All the estimators except Averaged-DPA-N can asymptotically approach CRLB at different SNRs. The MSEs at SNRs ranging from $-6$  dB to 2 dB are magnified to show details. It can be observed that the MLE-MA algorithm achieves the best performance but the performance gap between MLE-MA and MLE-MAS is small. Surprisingly, at SNRs below 10 dB, BLUE-DPA and BLUE-DPA-N have worse performance than BLUE-PA, BLUE-PA-N, and even Averaged-PA-N. This is because the differentiation items between correlations have amplified noise, leading to larger variances on the estimation of frequency offset. More differentiation items involved make BLUE-PA-N get poorer performance than BLUE-DPA-N. The operations of coarse frequency estimation and compensation make four estimations perform similarly. They are BLUE-PA, BLUE-PA-N, Averaged-PA-N, and Morelli's method. Notice that BLUE-PA has a lower computation complexity than Morelli's method, but they perform almost the same. Averaged-PA-N has lower complexity and suitable performance for all SNRs. If the complexity is an important factor for implementation, it can be considered. In the simulations, all the estimation results are processed to remove the possible frequency ambiguity for clearance.

Figure 5

MSE performance comparison between different CFO estimators for 2-by-2 MIMO-OFDM systems ( $N_t=2$ , $N_r=2$ , $N=64$ , $K=2$ , and $M=4$ ).

MLE-MAS and BLUEs have the same estimation range of $\left|\nu \right|\le M/2$ , but, under some circumstances, CFO is very large; then MLE-MA with small $K_P$ combined with BLUEs is a good choice to deal with these worse cases. Performance of MLE-MA combined with BLUE-PA for 2-by-2 system is given by Figure 6. CFO is randomly generated within interval $[-30,30]$ . We can see that combined algorithm can still approach CRLB at SNR of 12 dB for $M=4,6$ and $K_P=2,4$ , respectively. These estimators can be flexibly combined to satisfy different performance requirements, such as MSE, complexity, and SNR requirements.

Figure 6

MSE performance of MLE-MA combined with BLUE-PA for 2-by-2 MIMO-OFDM systems ( $N_t=2$ , $N_r=2$ , $N=64$ , $M=\mathrm{4,6}$ , $K_P=2$ , and the FFT length is 128).

7. Conclusions

In this paper, two MLEs and four BLUEs have been proposed as frequency offset estimators for MIMO-OFDM systems with periodic preambles and cyclic delays. The two proposed MLEs are MLE-MA and MLE-MAS. The MLE-MA can cope with a very large CFO up to the entire signal bandwidth but has a very high computational complexity. In contrast, MLE-MAS has shown to achieve a good MSE performance with low complexity, at the cost of narrow estimation range. In addition, four low complexity BLUEs with high accuracy have been designed and analyzed in detail. Due to the inherent spatial diversity, synchronization in MIMO-OFDM systems can achieve better performance than that of SISO-OFDM systems. The proposed CFO estimators can approach the theoretical bound given enough diversity order and adequate SNR. Simulation results have shown that the proposed schemes are effective in making synchronization in MIMO-OFDM systems more robust than in SISO-OFDM systems.