Maximum likelihood-based adaptive iteration algorithm design for joint CFO and channel estimation in MIMO-OFDM systems

In this paper, we present a joint time-variant carrier frequency offset (CFO) and frequency-selective channel response estimation scheme for multiple input-multiple output-orthogonal frequency-division multiplexing (MIMO-OFDM) systems for mobile users. The signal model of the MIMO-OFDM system is introduced, and the joint estimator is derived according to the maximum likelihood criterion. The proposed algorithm can be separated into three major parts. In the first part of the proposed algorithm, an initial CFO is estimated using derotation, and the result is used to apply a frequency-domain equalizer. In the second part, an iterative method is employed to locate the fine frequency peak for better CFO estimation. An adaptive process is used in the third part of the proposed algorithm to obtain the updated CFO estimation and track parameter time variations, including the time-varying CFOs and time-varying channels. The computational complexity of the proposed algorithm is considerably lower than that of the maximum likelihood-based grid search method. In a simulation, the mean squared error performance of the proposed algorithm was close to the Cramer-Rao lower bound. The simulation results indicate that the proposed novel joint estimation algorithm provides a bit error rate performance close to that in the perfect channel estimation condition. The results also suggest that the proposed method has reliable tracking performance in Jakes’ channel models.

channel into a set of frequency-flat fading channels for easier equalization [3]. From another perspective, OFDM can resist inter-symbol interference (ISI) by appending cyclic prefix (CP) samples in front of OFDM samples.
In the case of single-antenna OFDM systems or other multicarrier systems, three aspects should generally be managed: frequency synchronization, timing synchronization, and channel estimation. Multiple input-multiple output-orthogonal frequency-division multiplexing (MIMO-OFDM) is extremely sensitive to frequency synchronization and channel estimation errors [4]. Carrier frequency offsets (CFOs) are induced by local oscillator mismatches between transmitters and receivers as well as Doppler shifts [5]. CFOs break subcarrier orthogonality, which results in intercarrier interference and the possible performance degradation of the system [6]. Compared with frequency synchronization and channel estimation, timing synchronization is less critical because CP insertion has some tolerance to timing errors [7]. In this paper, we assume the timing offsets are equal to zero and focus on frequency synchronization and channel estimation.
To achieve high transmission quality, a frequency offset must be estimated and compensated at the receiver [8]. The combined estimation of the channel and the frequency offset causes complex problems in MIMO-OFDM systems due to the number of unknown parameters. In some estimation algorithms, CFO and channel estimations are treated separately by using different training sequences [9][10][11][12]. In these schemes, channel estimation is performed assuming zero CFO or frequency synchronization is achieved. However, such assumptions are rarely valid in practice with the presence of noise [9,10].
In other schemes, to save bandwidth, joint estimation of the channel and frequency offset is attempted. A prohibitively high computational complexity is required to obtain maximum likelihood (ML) solutions for both the frequency offset and channel impulse response (CIR) [13]. These estimators in [3,[5][6][7][9][10][11][12][13] are designed for single-input single-output OFDM systems and not MIMO-OFDM systems. In [14,15], the channel estimation of MIMO-OFDM systems was performed assuming no CFO imbalance or perfect frequency synchronization with the training sequences. In [16,17], the proposed algorithms for CFO estimation were executed assuming an estimated channel and negligible channel effect. In [18], a method based on the alternating projection algorithm was proposed for ML synchronization and channel estimation in orthogonal frequency-division multiple access uplink transmission. In [19][20][21][22][23], joint estimation algorithms have been proposed for time-frequency synchronization with channel identification in MIMO-OFDM systems. However, all these algorithms have rather high complexity. In [24][25][26], an iterative method was used to reduce the complexity. However, in these schemes in [24][25][26], the CFO estimator must use point search, which may assume a search range with a small interval to achieve low complexity. Moreover, the methods presented in [18][19][20][21][22][23][24][25][26] do not provide the ability to track time-varying CFOs and time-varying channels simultaneously.
This paper discusses the problem of joint CFO and channel estimation in MIMO-OFDM systems for mobile users, in which all subcarriers can be utilized simultaneously. The complexity of joint ML estimation through a grid search procedure motivated us to develop a new iteration algorithm [27] that can provide ML solutions for joint estimation problems in an iterative manner. The proposed algorithm can find the initial estimated CFO by employing a derotation algorithm and use this result as an initial value. It then uses the initial value to apply the frequency-domain equalizer and uses the derotation operation again for better estimation. Subsequently, an iterative method is adopted to improve estimation accuracy. To benchmark the performance of the proposed scheme, Cramer-Rao lower bounds (CRBs) were derived for the CFO and the CIR. Finally, the CFO was estimated through a proposed parameterized adaptive process by tracking channels and CFOs simultaneously. The simulation results indicated the suitable performance of the proposed adaptive iteration estimator. The computational complexity of the proposed algorithm is lower than that of the grid search-based method. Our contribution and new ideas are as follows: (a) a joint CFO and channel estimation in MIMO-OFDM systems for mobile users, (b) a new iteration algorithm with lower complexity than the grid search-based method, and (c) a newly designed mechanism with an adaptive mode to simultaneously track the time-varying CFOs and time-varying channels. Up to the present time, the proposed adaptive iteration algorithm is currently a state-of-the art approach for achieving near-optimal performance for this joint CFO and channel estimation problem. To the best of our knowledge, in other studies such as [18][19][20][21][22][23][24][25][26], no report has investigated the joint CFO and channel estimation problem with simultaneous tracking of the time-varying CFOs and time-varying channels.
The remainder of this paper is organized as follows. Section 2 presents the system model. Section 3 describes the development of the adaptive iterative CFO estimation algorithm. Section 4 presents the simulation results, which indicate that the performance of the proposed algorithm is close to the CRBs. Finally, Section 5 reports the conclusions.

System model
Consider a MIMO-OFDM communication system with K subcarriers as well as N transmitting (TX) and M receiving (RX) antennas (Fig. 1).
The model of each transmitter structure of the MIMO-OFDM system is illustrated in Fig. 2 [28].
At the jth transmitter, the preamble OFDM signal in the baseband time domain obtained after performing inverse discrete Fourier transform and CP insertion can be denoted as follows [29]: where d(k) is the frequency-domain data symbol for each subcarrier (k = 0, 1, …, K) and N g represents the CP samples of the OFDM symbol, which is used for ISI resistance. The discrete-time composite CIR between the jth transmitter and the ith received antenna with L multipaths is represented as follows: where h i, j (l), l = 0, 1, …, L − 1 is the complex Gaussian gain of the lth multipath. We assume the complex CIR is time-invariant over one OFDM symbol. Thus, the superposition of signals from all TX antennas plus noise is received at each RX antenna. Assume that all the received signals are down-converted to baseband with the same local oscillator centered at f c . In the presence of CFO ε, the samples at the ith RX antenna can be represented as follows: where v i (m) is an additive white complex Gaussian distributed noise with a mean of 0 and variance of σ 2 v i . Moreover, ε denotes the CFO normalized to the subcarrier spacing [5]. To avoid ISI, we assume that the timing offset is equal to 0 and that the number of multipaths (L) is smaller than the CP length.
According to (3), the received signal vector r i = [y i (0), y i (1), …, y i (K − 1)] Τ at the i th RX antenna can be used to rewrite the input-output relationship as follows: where D(ε) = diag (1, e j2πε/K , …, e j2π(K − 1)ε/K ) K × K is a diagonal CFO matrix, X j is a K × L  circular matrix formed by the j th transmitted signal x j (k), and v i is a K × 1 noise vector that can be expressed using the covariance matrix σ 2 v i I K ÂK . Then, we define two matrices X and h i as follows: Let D(ε)Xh i = s i and s i = [s i (0), …, s i (K − 1)] T so that the signal-to-noise ratio (SNR) is defined [10] as follows: 3 Proposed joint CFO and channel iterative estimation algorithm

Receiver design
The processing diagram at the receiver for the proposed algorithm is displayed in Fig. 3.
In the common RX structure used in OFDM transmission, a frequency-domain equalizer is adopted to eliminate multipath interference. The coefficients of the equalizer are the results obtained through channel estimation. These results are fed back to the time-domain signal, and the steps are run again to obtain more accurate CFO estimates. The main idea is that the proposed algorithm combines the CFOs estimated in the previous processes as a rough estimated value for the followed iteration. Through small-step iterative searching, we can find the maximum peak near the actual CFO as the fine-tuned estimate. Finally, the CFO adaptation process is performed using the final fine-tuned CFO estimate to improve the performance of the proposed estimator for tracking the time-varying parameters. The aforementioned are described in detail as the following text.

Initial CFO estimation
According to (4) and the Gaussian probability density function, the ML function can be written as follows: Then, the log-likelihood function for two unknown parameters can be expressed as follows: To evaluate the log-likelihood function in (9), the constant is eliminated to minimize the newly defined objective function, which is the second term of the log-likelihood function.
We can expand (10) and take the partial differentiation operator for h i ; then, let the first-order derivative be equal to 0 to obtainh i as follows: By substituting (11) into (10), for the purpose of minimizing the ML objective function (10), it becomes to maximize The estimation of ε for the i th RX antenna can be obtained as follows: Equation (13) can be used to perform a grid search to find the optimal CFO. However, this process requires many computations and is thus difficult to implement. The projection matrix X(X H X) −1 X H is formed by the preamble signal per frame transmission. The projection matrix can be represented as follows: The main diagonals of the X H X matrix are real-valued and are denoted as R i , i = 1,…, N·L. Its adjacent diagonals parallel to the main diagonals are denoted as C i, j , i = 1,…, N·L, j = 1,…, N·L, i ≠ j, which represent conjugate symmetry on the other side parallel to the main diagonal. The matrix X(X H X) −1 X H has a similar structure. where W i , i = 1, …, K, are real-valued and E i, j , i = 1, …, K, j = 1, …, K, i ≠ j, are complex valued and are the conjugate of E j, i . By substituting (15) into (13), the estimation of ε for the i th RX antenna can be expressed as follows: To find the maximum value for the estimation of ε, (16) can be separated into multiple terms. Because the first term in (16) is a fixed real value, we only consider the second part. In the second part, because the first term is a complex value caused by the frequency offset, we can derotate its phase to the real axis; thus, the frequency offset can be compensated to achieve a maximum real value. The frequency offsets of K terms for the ith RX antenna can be obtained as follows: After finding the phase of each term of (17) associated with the maximum real value, we can average the phases of all the terms for each RX antenna and then use the calculated average values to obtain the initial CFO as follows:

Frequency-domain equalizer
The initial CFO performance degrades due to multipath interference and may not be accurate. Therefore, a frequency-domain equalizer is applied to overcome this problem and obtain a more accurate CFO. First, we utilize the initial CFO (18) to find the channel response (11). Second, we use the initial CFO to compensate the received signal, which is then transformed into the frequency domain. A frequency-domain equalizer is employed using the estimated channel response to cancel the multipath interference [30]. The aforementioned operations are expressed as follows: Ifε ¼ ε, we can rewrite the received signal (19) and transform it into the frequency domain as follows: where H i, j is a diagonal matrix that is the frequency-domain channel matrix. The term H i; j ¼ diagðh i; j Þ; h i; j ¼ ½h i; j ð0Þ; h i; j ð1Þ; …; h i; j ðK − 1Þ Η represents a K × 1 vector that includes the samples from the K-point discrete Fourier transform of the channel response, and s j is the frequency-domain signal of the jth TX antennas. We can rewrite the signal received at all antennas over the kth subcarrier as follows: whereH k ∈ℂ MÂN is the MIMO channel matrix over the kth subcarrier ands k ∈ℂ NÂ1 is the signal transmitted over the kth subcarrier. According to the signal at each subcarrier, we can equalize the compensated received signals as follows: Subsequently, we can collect K consecutive samples and transform the equalized signals back into the time domain. The derotation method is then used again to calculate the more accurate CFOε e .
In the following section, an iterative method is described. This method provides a refined estimate close to the CRB.

Small-step iterative searching
The CFO estimation is refined using an iteration method. The CFOε e derived in the previous section is used as the initial main frequency estimate. Then, adjacent frequencies are selected near the initial main frequency as candidate frequencies. The next step is to substitute the main and candidate frequencies into (16) to evaluate the CFO values. The largest CFO function value is determined to identify the search direction (refer to Fig. 4, step 1). Subsequently, the frequency selected in the previous step becomes the new main frequencyε m . The fixed-step frequency Δf is added with the new main frequencyε m to acquire its adjacent frequency, and the comparison is performed as before (refer to Fig. 4, steps 2-1 and 2-2). This iterative process continues until the function value of the main frequency in the evaluation is larger than that of its adjacent frequency (refer to Fig. 4, steps 2 and 3; e.g., the total number of iterations is 3). Finally, the maximum frequency peak ε f is achieved in the region between the last main frequency and the last adjacent frequency (refer to Fig. 4, step 3).
However, the global maximum frequency peak may not be in the search region if the setting of a small-step frequency is inappropriate. For instance, the CFO estimated with the iterative method would not be at the true frequency peak if the true global maximum frequency is in the previous range between the iterative main frequency and the iterative adjacent frequency (refer to Fig. 5, steps 2 and 3). Because the iterative adjacent frequencyε m i þ Δ f is larger than the iterative main frequencyε m i , the maximum frequency is decided to be located at the right-hand side of the iterative adjacent frequency. Therefore, the iteration continues until the result of the main frequency is larger than that of the adjacent frequency (refer to Fig. 5, steps 2-4). However, the peak of the maximum frequency is outside the final search range.
In the aforementioned circumstance, the method of reverse iterative searching is adopted. When the maximum value is at the boundary frequency, the search direction

Computational complexity and the procedure of the proposed method
We adopt the Big-Oh notation, which is a well-accepted approach for analyzing the computational complexities of algorithms. The detailed analysis of the computational complexity is as follows. Note that K is the number of subcarriers and N and M are the numbers of TX and RX antennas, respectively. The proposed scheme comprises five major steps: (1) In step 1, Eqs. (17) and (18)   Step 1 Use Eq. (17) to estimate each item of the second part of Eq. (16), and average the results to obtain the initial estimate of CFOε using Eq. (18).
Step 2 Use the result obtained in step 1 for channel estimation with Eq. (11) and design the coefficient of equalizer taps.
Step 3 Equalize the received signal by using Eq. (22) and estimate CFOε e as in step 1.
Step 4 Set the main and adjacent frequencies using Eq. (16) from step 3 and determine the frequency peak range using the iterative method.
Step 5 Search for the frequency peak ε f in the range obtained in step 4 by inputting the new main frequency into Eq. (16). Step 2 O(M (K+ KNL)) Step The procedure of the proposed algorithms and the computational complexity are summarized in Tables 1 and 2, respectively.

Adaptive mode for tracking the time variations of parameters
According to the aforementioned described process, the CFOε T ¼t and CFOε T ¼tþΔt (CFO of the next time unit) are obtained ( t can be the symbol time, slot duration, frame duration, or any other time unit in the transmission systems). The adjustment equation of the final CFO is updated with the coefficient μ as follows: By substituting (23) into (11), the channel estimation in this adaptive mode is rewritten as follows:

Simulation
This section presents the simulation results and discusses the efficacy of the proposed schemes, including the iterative mode and adaptive tracking mode. The parameters common to both the iterative mode and adaptive tracking mode were as follows: M = 2, N = 2, K = 64, and the CP size N g = 16. A preamble was inserted at the beginning of each transmitted signal per frame. The transmitted OFDM data symbols were modulated through quaternary phase-shift keying by using a three-multipath channel model, and estimations were made on the basis of the preamble. The delay profile for the multipaths is presented in Table 3. The channel gain of each multipath was randomly generated with a Gaussian distribution of variance 1. Equations (6) and (7) were used for the SNR calculation in the simulation. After the channel gain coefficients are generated, the noise power associated with the received signal power is generated for a particular SNR. A Δf value of 10 −5 was adopted in the simulation. The results of 5000 Monte Carlo simulations were averaged for each SNR. For the iterative mode, the normalized CFOs were randomly selected independently from the random variable and were uniformly distributed over [−0.5, 0.5]. The Fig. 7 The MSE performance for channel estimation compared with CRB and the ideal case  coefficients of the CIR were complex valued and generated independently from normally distributed random numbers.
Regarding the adaptive tracking mode, the coefficients of the CIR were generated using Jakes' model, with the carrier frequency equaling 3.5GHz, which is to be adopted in 5G. Various velocities (v) produced a time-varying Doppler shift. The channel variations in Jakes' model rely on the parameter of the maximum Doppler shift, which is determined by the carrier frequency and the mobile speed. Therefore, the change of the carrier frequency to scrutinize other frequencies in the 5G band is equivalent to the change of the mobile speed with a larger value for the specified carrier frequency of 3.5 GHz. Simulations were performed for various mobile speeds.
In the case of channel estimation errors, Εðkĥ − hk 2 Þ was calculated and plotted using the sample averages. To obtain benchmarks for the performance evaluation, the lower bounds for the CFO and CIR in the iterative mode were derived using the Fisher information matrix based on (4), as in [31].
where ⊗ denotes the Kronecker product. Moreover, the following equation is obtained: Similarly, the CRB was obtained for the estimatedĥ as follows: where Fig. 11 The MSE performance for CFO estimation in the adaptive mode at v = 500 km/h

Algorithm performance in the iterative mode
The mean squared error (MSE) performance of the proposed algorithm for CFO estimation in the iterative mode at different SNRs is revealed in Fig. 6. The performance of the proposed algorithm in the iterative mode was satisfactory and close to the CRB [refer to (26)]. The MSE performance of the proposed algorithm for channel estimation at different SNRs is illustrated in Fig. 7. Comparisons of the CRB [refer to (28)] and channel estimate with the estimates in the "perfect CFO estimation" condition, in which the CFOs are perfectly known, indicated that the proposed joint estimation algorithm provided satisfactory channel estimation performance.
In Fig. 8, the bit error rates (BERs) of the three schemes are illustrated as follows: (a) the "ideal channel" condition, which has perfect channel estimation and perfect CFO synchronization; (b) the "perfect CFO estimation with joint channel estimation" condition, in which CFOs are assumed to be known for estimating the CIR; and (c) "proposed joint CFO and channel estimation." The results indicate that the BER performance of the joint estimation is nearly identical to that with ideal assumptions. The proposed algorithm not only has a considerably lower computational complexity than the grid search-based method but also provides satisfactory performance, as revealed in Figs. 6, 7, and 8.

Algorithm performance in the adaptive tracking mode
The MSE performance of the proposed algorithm for CFO estimation in the adaptive mode at different SNRs is illustrated in Figs. 9, 10, and 11. To obtain intensive simulation results and the best MSE performance, the value of μ was selected as follows for various mobile speeds: (a) μ = 0.5 at v = 0 km/h in Fig. 9, (b) μ = 0.4998 at v = 60 km/h (the speed limit for automobiles in suburban areas) in Fig. 10, and (c) μ = 0.18 at v = 500 km/h (the average speed of high-speed rail, as defined for 5G) in Fig. 11. The μ values for various other mobile speeds are presented in Table 4. The adaptive mode was executed on a per frame basis. The frame duration t was 1 ms. As expected, the simulation results indicated that the coefficient μ should be decreased when the velocity v increases. The estimated CFOε T ¼tþΔt of the next frame should have more weight than the estimated CFOε T ¼t when the time variation is faster. Therefore, the weighting of the parameterized equation [refer to (23)] tends toward the second term. As illustrated in Figs. 9 and 10, the performance of the proposed algorithm in the adaptive mode was superior to its performance in the iterative mode. As depicted in Fig. 11, the MSE performance of the proposed adaptive mode began to degrade at SNR = 24 dB due to high time variations, and μ was approximately selected. This phenomenon did not occur at low mobile speeds. The MSE performance of the proposed algorithm for channel estimation in the adaptive mode is displayed in Figs. 12, 13, and 14. The value of the parameter μ is selected for various mobile speeds as: (a) μ = 0.5 at v = 0 km/h in Fig. 12, (b) μ = 0.4998 at v = 60 km/h in Fig. 13, and (c) μ = 0.18 at v = 500 km/h in Fig. 14, respectively, for the approximate best MSE performance. The results obtained with the proposed algorithm in the iteration mode and those obtained in the "perfect CFO estimation" condition were compared. The simulation results revealed that the performance obtained in the proposed adaptive mode was competitive with that obtained in the "perfect CFO estimation" condition (Figs. 9, 10, and 11). Figures 15 and 16 display the BERs of the following three schemes: (1) the "ideal channel" condition, in which perfect channel estimation and CFO estimation are assumed; (2) the "perfect CFO estimation and joint channel estimation" condition, in which the CFOs are assumed to be known to estimate the CIR; and (3) the "proposed joint CFO and channel estimation" in the adaptive mode. To obtain the best BER performance, the values of μ were selected as follows: (a) μ = 0.4998 at v = 60 km/h in Fig.   Fig. 16 The BER performance of the joint estimation in the adaptive mode at v = 500 km/h  Fig. 16. The results indicated that the proposed adaptive mode estimation offers a BER performance nearly identical to that obtained with ideal assumptions. The tracking ability of the proposed adaptive mode was examined at v = 120 km/h (Fig. 17). The estimate was compared with a real time-varying CFO. The comparison indicated that the proposed algorithm can track a time-varying CFO in the adaptive mode. For a vehicle velocity of 120 km/h, the tracking results of the h 11 are plotted in Fig. 18, where we select one of the channel gains in the illustration. The time variations of the channel were accurately tracked.

Conclusion
In this paper, we propose ML-based algorithms for joint CFO and channel estimation. The proposed methods provided fairly competitive performance to that of CRBs. Moreover, the proposed methods have a lower computational complexity than the grid search-based method does. In addition, an adaptive mode is proposed to improve the algorithm performance. The adaptive method is used to obtain the weighted CFO; thus, the adaptive mode can enhance the performance of the original proposed algorithm and provide tracking capability for time-varying parameters. The parameter μ should be adjusted according to operational environments.