1558 IEEE TRANSACTIONS ON WIRELESS COMMUNICATIONS, VOL. 5, NO.

6, JUNE 2006

Reduced Complexity Data-Aided and

Code-Aided Frequency Offset Estimation for
Flat-Fading MIMO Channels
Frederik Simoens, Student Member, IEEE, and Marc Moeneclaey, Fellow, IEEE

Abstract— This contribution deals with carrier frequency off- In [9], [10] data-aided joint ML estimation of frequency off-
set estimation for flat-fading Multiple-Input Multiple-Output set and channel gains for MIMO systems has been addressed.
(MIMO) channels. Both data-aided (using training symbols) Also an ad-hoc estimator has been proposed to synchronize
and iterative code-aided (using the unknown coded symbols)
estimation is considered. In both scenarios, maximum-likelihood the signals from the different transmit antennas separately in
(ML) frequency offset estimation involves solving a maximization [9]. In contrast to [9], the estimator proposed in the present
problem with no closed-form solution. Since numerical calcula- contribution is an ML approximation, combining information
tion of the ML estimates is computationally hard, we derive from all transmit antennas. Consequently, full diversity is
a simple closed-form approximation. Simulation results indicate exploited and the number of pilot symbols can be reduced
that the ML algorithm and the proposed reduced-complexity
algorithm operate closely to the Cramer-Rao bound (CRB). as compared to [9].
To achieve accurate synchronization without resorting to a
Index Terms— Frequency offset estimation, synchronization, large amount of pilot symbols, the frequency offset estimation
MIMO channels.
algorithm should also make use of the data portion of the
frame containing the information-bearing symbols. Recently,
I. I NTRODUCTION a great effort is being devoted to develop efficient code-aided
(CA) estimation techniques using soft information from the

U SING multiple transmit and receive antennas in wireless

fading channels has been advocated as a means to in-
crease capacity and create diversity [1]. To exploit this possible
decoding process, in order to fully exploit the code properties
during estimation. In particular the expectation-maximization
(EM) algorithm [11] is an attractive tool to carry out iterative
increase of capacity and diversity, several coding and multi- joint estimation and decoding. Most of the current technical
plexing schemes have been presented [2]–[4]. However, the literature in the field of estimation for MIMO channels has
performance of these schemes critically depends on the avail- been focused on channel estimation [12], [13], whereas, the
ability of accurate estimates of the synchronization parameters. available literature on code-aided frequency offset estimation
In particular, carrier frequency offsets can substantially affect is quite scarce. In this contribution, we investigate the EM
data detection. Moreover, the presence of a frequency offset is algorithm for joint frequency offset and channel estimation in
often detrimental in the estimation process of other parameters, a MIMO-context. We also propose a modified version, with
such as the channel gains. Hence accurate frequency offset lower computational complexity.
compensation in MIMO configurations is important to achieve This contribution is structured as follows. In section II we
a satisfactory performance. present the MIMO signal model. Section III covers data-aided
Similarly to single antenna configurations, data-aided (DA) frequency offset estimation in case the channel is known and
(using training symbols) ML frequency offset estimation in in case it is unknown. Frequency offset estimation when the
a MIMO context requires a numerical calculation, which is a channel is known is of no direct practical interest, however,
time consuming task [5], [6]. To overcome this inconvenience, its relevance will become apparent in the sequel. Code-aided
some computationally efficient estimation techniques have estimation based on the EM algorithm is addressed in section
been developed for Single-Input Single-Output (SISO) sys- IV. Finally, conclusions are drawn in section V.
tems, which result in near optimum performance [7], [8]. The
present paper deals with computationally efficient estimation II. S IGNAL MODEL
extended to a flat-fading MIMO configuration.
Let us consider a flat-fading MIMO channel with NT
Manuscript received July 28, 2004; revised February 28, 2005; accepted transmit antennas and NR receive antennas. We further assume
May 6, 2005. The editor coordinating the review of this paper and approving that all transmit/receive antenna pairs are affected by the same
it for publication is D. I. Kim. This work has been supported by the frequency offset [14]. This assumption is definitely valid for
Interuniversity Attraction Poles Program - Belgian Science Policy and the
Network of Excellence in Wireless Communications (NEWCOM), a project a frequency offset caused by an oscillator mismatch, since
funded by the European Community. the different antennas can use a common oscillator or at least
The authors are with the Department of Telecommunications and In- different collocated oscillators with a known and compensable
formation Processing, Ghent University, Gent B-9630, Belgium (e-mail:
{fsimoens,mm}@telin.ugent.be). difference. The frequency offsets induced by Doppler shifts
Digital Object Identifier 10.1109/TWC.2006.04498. will be identical if the angles-of-arrival of the signals on
1536-1276/06$20.00 
c 2006 IEEE
SIMOENS and MOENECLAEY: REDUCED COMPLEXITY DATA-AIDED AND CODE-AIDED FREQUENCY OFFSET ESTIMATION 1559

the different antennas are identical. In many communication Hence, the ML estimate maximizes the magnitude of the
systems, this is a valid assumption. Further, estimation and Fourier transform (FT) of {y(k)}. For high signal-to-noise
compensation of the mean Doppler shift, allows us to assume ratio (SNR), the mean-square error (MSE) of this estimate
the channel to be quasi static as long as the block duration closely approaches the Cramer-Rao bound (CRB), which is a
is smaller than the inverse of the remaining Doppler spread. theoretical lower bound on the MSE. Unfortunately, a simple
Sampling at symbol-rate, we can model the observation as closed-form solution to the problem of maximizing (2) does
not exist. The exact determination of the ML estimate from
yk = Hak ej2πF kT + wk k = 0...L − 1 (1) (2) would require a time consuming search over a large set
where yk = [ykn ]Tn=1...NR is a NR ×1 vector of received signal of frequency values, making ML estimation computationally
samples, H = [hnm ]n=1...NR ,m=1...NT denotes the NR × NT hard. To overcome this problem, we derive a sub-optimal
channel matrix, ak = [a m T closed-form frequency estimator.
 k ]m=1...NT is a NT × 1 vector of 2) Computationally efficient estimation: Replacing in (2)
2
symbols with E |am k | = Es , F is the carrier frequency the magnitude of the FT by its squared magnitude does
offset, 1/T is the symbol rate per transmit antenna, and wk not affect the value of the ML estimate. Rearranging the
is a NR × 1 noise vector. The components of the noise vector squared magnitude yields the following maximum likelihood
wk are statistically independent and complex-valued, with expression:
independent Gaussian zero-mean real and imaginary parts,  L−1
each having a variance of N0 /2; also, noise vectors at different 
−j2π F̃ mT
F̂ML;F = arg max Re RDA;F (m)e (4)
instants of time are statistically independent. The components F̃ m=1
hmn of the channel matrix are statistically independent and
complex-valued, with independent Gaussian zero-mean real where the time-correlation RDA;F (m) is defined as

 parts, each having a variance of 1/2; this yields

and imaginary L−1

E |hmn |2 = 1. RDA;F (m) = y ∗ (k − m)y(m).
It is useful to denote the NT × L space-time matrix of k=m
transmitted symbols as A = [a0 , a1 , . . . , aL−1 ] and the For further use, we decompose RDA;F (m) as the sum of a
NR × L space-time matrix of received signal samples as useful term that is independent of the additive noise wk , and
Y = [y0 , y1 , . . . , yL−1 ]. a noise term caused by wk
In this contribution, we consider both data-aided and code-
aided estimation of the parameters (F, H). In the former case, RDA;F (m) = A(m) exp(j2πF mT ) + noise (5)
the transmitted symbols are training symbols known to the
receiver, whereas in the latter case the transmitted symbols where A(m) is a positive real quantity, given by
result from channel-encoding a random information sequence. L−1
 2 2
A(m) = |Hak | |Hak−m | . (6)
III. DATA - AIDED ESTIMATION k=m

A. Frequency offset estimation with known channel Setting the derivative (with respect to F) of the function to be
maximized in (4) equal to zero, the ML-estimate satisfies the
In this section, we consider data-aided frequency offset following equation:
estimation assuming the channel matrix H to be known. This  L−1
is not a realistic scenario, but it will turn out that the resulting 
−j2πmF̂ T
algorithm is useful in the context of code-aided iterative joint Im mRDA;F (m)e = 0. (7)
estimation of the frequency offset and the channel matrix, to m=1

be investigated in section IV. Equation (7) can equivalently be written as

1) ML estimation: The problem of frequency offset esti-
mation for Additive White Gaussian Noise (AWGN) SISO- L−1
 

channels using a training sequence is well documented in m |RMF (m)| sin arg(RDA;F (m)) − 2πmF̂ T = 0.
literature [5]–[8]. When the MIMO channel matrix H is m=1
considered to be known, up to the global phase1 (H = H̃ejθ ; (8)
H̃ is known, θ remains unknown), a similar derivation as for Assuming that the training sequence is sufficiently long,
SISO yields the following ML estimate [15] the first term in (5) is strongly dominating, in which case
arg(RDA;F (m)) ∼ = 2πmF T . Hence it follows from (8) that
L−1  F̂ ∼
  = F , so that the argument of the sine function in (8) is
  small. Using in (8) the approximation sin (x) ∼
F̂ML;F = arg max  y(k)e−j2πF̃ kT  (2) = x for |x|  1
F̃   yields the following linear equation with respect to F̂
k=0

where
L−1
 

m |RDA;F (m)| arg(RDA;F (m)) − 2πmF̂ T = 0.
y (k) = aH H
k H yk (3)
m=1

1 Itis not a realistic scenario to have complete phase information but no Grouping terms appropriately and limiting the summation
frequency offset information. interval to [1, M ], with M ≤ L − 1, to reduce complexity
1560 IEEE TRANSACTIONS ON WIRELESS COMMUNICATIONS, VOL. 5, NO. 6, JUNE 2006

−3
10
Fitz (9), (11) and (12), with the CRB. The system parameters
are L = 20 and M = P = Q = L/2. We used random
Luise
DA;F
ML;F
CRB training symbols, although the performance of the different
10
−4 estimators is essentially independent of the training sequence.
We observe that the ML algorithm (2), the Fitz algorithm
(12) and the new algorithm (9) operate closely to the CRB
N =1
at moderate and large SNR, whereas the L&R algorithm (11)
MSEF

−5
10 R
suffers from a MSE floor. This MSE floor can be explained by
noticing that the derivation [7] of the L&R algorithm assumes
NR=2 A(m) in (6) to be constant; as this assumption is in general
10
−6
violated for MIMO (but not for SISO), the MIMO version
of the L&R algorithm gives rise to estimation errors in the
absence of noise. Further, we observe that the new algorithm
−7
(9) outperforms both the L&R and the Fitz algorithms. Fig. 1
10
−10 −5 0
E /N (dB)
5 10 15 also illustrates that the MSE reduces when the number of
s 0
receive antennas increases. As far as complexity is concerned,
Fig. 1. MSE-performance of different DA frequency estimation algorithms Table I indicates that the simplified algorithms are essentially
for known channel. (L = 20, NT = 4, NR = 1 − 2) equivalent for M = P = Q, as the evaluation of the time
TABLE I
correlations {RDA;F (m)} is the main computational burden.
T HE COMPLEXITY OF THE DA FREQUENCY ESTIMATION ALGORITHMS
( FOR known CHANNEL ). F OR A FAIR COMPARISON , THE PARAMETERS ARE B. Joint channel and frequency offset estimation
P = Q = M = L/2.
1) ML estimation: In this section we address the data-
aided joint maximum likelihood estimation of the channel and
complex additions complex multiplications
  frequency offset. Let us consider the log-likelihood function
3 5
L&R L 8
L + NT NR − 4
L 38 L + 2NT NR + 12 (LLF) of H and F related to the observation (1).
3 5

Fitz L 8
L + NT NR − 4
L 38 L + 2NT NR + 34
3 3

DA;F L 8
L + NT NR − 4
L 38 L + 2NT NR + 54
L−1
1   2
ln p (Y|H, F ) = − yk − Hak ej2πF kT  (13)
N0
k=0
furthermore, the new closed-form frequency estimator is ob- The joint ML-estimate corresponds to the values of H and
tained: F that maximize ln p (Y|H, F ). For a given trial value F̂ ,
the channel matrix Ĥ that maximizes (13) is readily found
M
1 m=1 m |RDA;F (m)| arg [RDA;F (m)] to be (15). Inserting (15) in (13) and rearranging yields the
F̂DA;F = M . (9) frequency offset estimate (14).
2πT 2
m=1 m |RDA;F (m)|
We should indicate that because of the presence of the arg(.)  L−1

function in (9), F̂DA;F is ambiguous when |arg (R(m))| F̂ML;H,F = arg max Re RDA;H,F (m)e −j2πmF̃ T
(14)
exceeds π. This limits the operating range of the proposed F̃ m=1
frequency offset estimation technique to the interval L−1  L−1 −1
 
|F | < [2M T ]−1 (10) ĤML;H,F = yk aH
k e
−j2π F̂M L;H,F kT
an aH
n (15)
k=0 n=0
The above algorithm (9) is similar to the straightforward
where
extensions [15] to MIMO of the Luise&Reggiannini (L&R)
L−1

[7] and Fitz [8] SISO frequency estimation algorithms. The H
corresponding MIMO algorithms are RDA;H,F (m) = yk−m yk A(k, k − m) (16)
P  k=m
1  m
L&R:F̂L&R = arg RDA;F (m) (11) and L−1 −1
π(P + 1)T m=1
L−m 
Q A(k, l) = aH
k an aH
n al .
1 m=1 m arg (RDA;F (m))
Fitz: F̂F itz = Q (12) n=0
2πT m=1 m
2
For high signal-to-noise ratio (SNR), the mean square error
where P and Q are design parameters. The L&R and Fitz (MSE) of these estimates approaches the CRB. The exact
algorithms have an operating range of |F | < [P T ]−1 and expression of the CRB for frequency offset estimation is found
|F | < [2QT ]−1 respectively [7], [8]. In a SISO context, the in [9].
L&R and Fitz algorithms both give rise to a MSE that is very As for frequency estimation with a known channel (4), a
close to the CRB for P = Q = L/2. closed-form solution to the maximization problem (14) does
Fig. 1 compares the MSE resulting from the computation- not exist. Exact determination of F̂ML;H,F requires a time
ally intensive ML algorithm (2), and the simplified algorithms consuming numerical search which makes frequency offset
SIMOENS and MOENECLAEY: REDUCED COMPLEXITY DATA-AIDED AND CODE-AIDED FREQUENCY OFFSET ESTIMATION 1561

−4
10
DA;H,F where
ML;H,F
CRB L−1

NT=1,N R=2, L=20 C(m) = aH H
k−m H Hak A(k, k − m). (19)
k=m
−5
10

It is important to note that C(m) is generally complex-

valued. Considering the phase of the time-correlation, we can
N =4, N =2, L=20
distinguish three terms:
MSEF

−6 T R
10

arg (RDA;H,F (m)) = 2πF mT + φd (m) + φn (m). (20)

10
−7 N =4, N =2, L=40
T R The first term corresponds to the phase rotation caused by
the frequency offset; the second term φd (m) = arg(C(m)) is
a self-noise term, and the last term φn (m) is introduced by
noise. It is apparent that the approximation
−8
10
−10 −5 0 5 10 15

arg (RDA;H,F (m)) ∼

E /N (dB)
= 2πF mT
s 0
(21)
Fig. 2. MSE-performance of frequency estimator (17) using randomly
that is crucial in obtaining (17) is less accurate when φd (m)+
selected training symbols compared to ML estimator and CRB (unknown
channel). φn (m) takes on larger values. Considering Fig. 2, the follow-
ing remarks are in order:
• For systems with a single transmit antenna (NT = 1), it

estimation a computationally hard problem. On the other is readily seen from (19) that C(m) is real-valued, and
hand, the channel estimation algorithm (15) requires only a no self-noise term φd (m) appears. This illustrates that
small computational effort once the frequency offset estimate the self-noise is induced by the multiple transmit antenna
has been found. Hence it is of practical interest to search set-up.
• For systems with more than one transmit antenna, C(m)
for a frequency offset estimation algorithm with reduced
complexity, avoiding the exhaustive search implied in (14). is generally not real, resulting in a self-noise term φd (m).
Consequently, even in the absence of noise, the approx-
2) Computationally efficient estimation: Here, we derive imation (21) is not valid, and an error-floor is observed.
a new closed-form joint frequency offset and channel matrix • When the number L of randomly selected pilot symbols
estimator. As will be apparent in the sequel, the multi-antenna increases, the error floor reduces. This effect is readily
set-up gives rise to self-noise when the training sequence is understood since, for large L, C(m) converges to its
arbitrary. In section III-B.3 we show how proper design of the expectation Ea [C(m)], which is real-valued.
training sequence can avoid this problem. • The ML-estimator (14) does not suffer from an er-
Considering the similarity between (4) and (14), our closed- ror floor, which might be unexpected since the time-
form approximation to (14) is correlation in (14) contains the self-noise term φd (m).
However, as shown in appendix I, F̂ML;H,F = F in the
absence of noise. This is clearly not true for F̂DA;H,F ,
L−1 indicating that the self-noise term φd (m) is too large for
1 m=1 m|R DA;H,F (m)| arg(RDA;H,F (m))
F̂DA;H,F = L−1 2 . the approximation (21) to be accurate.
2πT m=1 m |RDA;H,F (m)|
(17) 3) Training sequence design: We will show how a
The corresponding channel matrix estimate is given by (15), proper selection of the training symbols can eliminate
with F̂ML;H,F replaced by F̂DA;H,F . In contrast with the the self-noise and the resulting MSE floor. As explained
algorithm (9) for the known channel case, the performance of L−1 above, the MSE floor arises from the fact that C(m) =
H H
the approximated closed-form solution (17) does not approach a
k=m k−m H Ha k A(k, k − m) is complex-valued. The
the CRB at high SNR for an arbitrary training sequence. In following proposition constitutes an approach to make C(m)
Fig. 2 we observe an error floor in the MSE performance of the real-valued.
the closed-form estimator (17), whereas the ML estimator (14) proposition 1: C(m) is real-valued ∀m if every sym-
does not suffer from such floor. It is interesting to note that bolvector ak is chosen from a complete orthogonal set, i.e.
a similar flooring effect occurs with closed-form estimators a k ∈ {α 0 , α1 , . . . , αNT −1 } ∀k, with αH i αj = NT Es δi−j ,
designed for SISO frequency selective channels [16]. and every vector α j is chosen at least once.
To explain this discrepancy between the performance of the proof: see appendix II.
sub-optimal solution (17) and the CRB or ML performance, Accordingly, the self-noise term φd (m) in (20) equals zero
we examine (17) and the time-correlation (16). Separating and the error-floor disappears. We point out that the number
signal from noise terms, the latter can be written as of symbolvectors in the orthogonal set is equal to the number
of transmit antennas NT .
In order to optimize the training sequence furthermore, we
take the training sequence constraint, which minimizes the
RDA;H,F (m) = C(m) exp (j2πF mT ) + noise (18) MSE of the channel estimate, into account. In [17], the training
1562 IEEE TRANSACTIONS ON WIRELESS COMMUNICATIONS, VOL. 5, NO. 6, JUNE 2006

−4
10
sequence that minimizes the MSE related to the ML estimation DA;H,F

of the channel matrix H was found to be orthogonal across

ML;H,F
CRB

all transmit antennas, i.e. the different rows of the matrix

[a0 , a1 , . . . , aL−1 ] should be orthogonal, which means that 10
−5


ak aH
k = LEs I, (22) N =4, N =2, L=20
T R
k
with I denoting the identity matrix. Choosing every vector ak

MSEF
−6
10

from an orthogonal set, it is easy to see that (22) holds, as

long as each orthogonal vector αj is chosen an equal number
of times (see Remark in appendix II). 10
−7

Let us now consider a periodic pilot sequence consisting of NT=4, N R=2, L=40

K systematic repetitions of the NT orthogonal symbolvectors

(L = KNT )
−8
10
aiNT +j = αj , i = 0 . . . K − 1, j = 0 . . . NT − 1 (23) −10 −5 0
Es/N0 (dB)
5 10 15

With this definition, no error floor will occur and (22) is

Fig. 3. MSE-performance of DA frequency estimator (24) using orthogonal
satisfied. Due to the periodicity, RDA;H,F (m) is zero except training compared to ML estimator and CRB (unknown channel).
for indices m = kNT , k ∈ N, yielding a significant
complexity reduction since only these non-zero terms have to
be taken into account in the computation of (17). To reduce Fig. 3 shows the MSE performance of estimator (24) using
complexity furthermore, we limit the summation interval in orthogonal training symbolvectors. We notice a substantial
(17) to [1, M NT ], with M < K. We have verified that this performance improvement compared to the case of randomly
does not result in a significant performance degradation, as selected pilot symbols, shown in Fig. 2. In Table II, we
long as M ≥ K/2. Assuming that the training sequence illustrate the significant reduction in complexity of (24) com-
satisfies (23), (17) can be simplified yielding the following pared to the fast-fourier transformation (FFT) implementation
closed-form frequency offset estimator: of the ML algorithm. The complexity of the ML estimator
corresponds to the ν-point FFT calculation; the parameter
M
 ν depends on the required precision, but is typically much
F̂DA;H,F = Bn arg(r(n)), (24)
greater than L. We notice a reduction in complexity of at least
n=1
a factor 43 NT , but the reduction in complexity is even greater
where
L−1
 for ν  L. Furthermore, the actual maximization search is
H
r(n) = yk−nN T
yk not considered in the complexity of the ML estimator.
k=nNT

and IV. C ODE - AIDED ESTIMATION

1 n(K − n) In this section we address code-aided (CA) frequency offset
Bn = 
2πNT T M
n=1 (K − n)n
2 and channel estimation. We assume no knowledge of the
and where we have carried out the following approximation, transmitted symbols at receiver side, yet information about
valid for small noise: the code structure is available. As we will point out, however,
a limited number of pilot symbols is necessary to initialize
NT −1 the code-aided estimation scheme. We resort to the iterative
K −n  EM algorithm [11] to solve the estimation problem. In section
|RDA;H,F (nNT )| = aH H
k H Hak .
K IV-A we briefly outline the theory behind the EM algorithm.
k=0
Application of the EM algorithm to the specific frequency
As the summation over k in the above expression does not
offset and channel estimation problem is considered in section
depend on n, it cancels from the numerator and denominator
IV-B. Similarly to data-aided estimation considered in section
in (17). Again, the operating range of the algorithm (24) is
III, we propose in section IV-C some modifications to reduce
limited, because of the presence of the arg(.) function in
complexity. Computer simulations compare the performance
(24), and is given by |F | < [2M NT T ]−1 . Note that only
of the estimator resulting from the EM algorithm and the
the correlations between two time-instants corresponding to
approximated estimators.
identical symbolvectors are considered in the calculation of
r(n).
The estimator ’s implementation is reasonably simple, be- A. ML estimation through the EM algorithm
cause the windowing function Bn can be precalculated for We consider the estimation of a parameter vector b from an
a given set-up. Once the frequency offset estimate has been observation r in the presence of a so-called nuisance vector
obtained, the channel estimate is found with little additional x, with distribution p(x), which does not depend on b. In the
computation, using (22), context of the present synchronization problem, x denotes the
L−1 coded symbols A, whereas b contains the synchronization
1  parameters (frequency offset and channel gains), i.e. b 
Ĥ = yk aH
k e
−j2π F̂ kT
. (25)
LEs {F, H}.
k=0
SIMOENS and MOENECLAEY: REDUCED COMPLEXITY DATA-AIDED AND CODE-AIDED FREQUENCY OFFSET ESTIMATION 1563

TABLE II
T HE COMPLEXITY OF M L; H, F FREQUENCY ESTIMATOR (14) COMPARED TO THE DA; H, F FREQUENCY ESTIMATOR (24). (M IS CHOSEN M = K/2)

complex additions complex multiplications

L(L−1) ν L(L−1)
ML;H,F ν log ν + 2
(NR + NT − 2) 2
log ν + 2
(NR + NT + 1)
L L
DA;H,F 8NT
((NR − 1)(3L + 2NT − 4) + 4) 8NT
(NR (3L + 2NT ) + 8)

The ML estimate b̂ML of b maximizes the log-likelihood 2) Maximization step: The maximization of (31) is similar
function (LLF): to the maximization of (14) for the data-aided estimation

  problem. Hence, the following estimate update equations are

b̂ML = arg max ln p r b̃ (26) obtained:
b̃
where  L−1


  
  (i+1)
F̂EM = arg max Re
(i)
RCA;H,F (m)e−j2πF̃ mT (34)
 
p r b̃ = p r x, b̃ p (x) dx (27) F̃ m=1
x

 
 L−1  L−1 −1
denotes the likelihood function. Often p r b̃ is very diffi-  (i+1)  (i)
(i+1) (i)H
cult to calculate, since the number of possible coded symbol ĤEM = yk āk e−j2πF̂EM kT · H
an an
combinations grows exponential with the block length. The k=0 n=0
(35)
EM algorithm is a method that iteratively solves (26). The
where
EM algorithm breaks up in two parts: the Expectation part
(28) and the Maximization part (29): L−1

(i) H
 RCA;H,F (m) = yk−m yk A(i) (k, k − m)


 
 k=m
Q b̃, b̂(i) = ln p(r|x, b̃)p x r, b̂(i) dx (28)
x

 
(i) −1

 and A(i) (k, l) = āk
(i)H H al (i) . For large L,
n an an
b̂(i+1) = arg max Q b̃, b̂(i) . (29) (i)
b̃ RCA;H,F (m) can be simplified, using the approximation
It has been shown that b̂(n) converges to a stationary point of L−1 H
(i)
∼
n=0 an an = LEs I.
the LLF under fairly general conditions [11].

B. Code-aided frequency offset and channel estimation C. Reduced complexity estimation

In this section we apply the EM algorithm to the synchro- The EM estimation update equations (35) and (34) are very
nization of a MIMO system with b  {F, H} the parameters similar to the data aided ML estimation described in section
to be estimated. The nuisance parameter x and the observation III-B.1. Again computational complexity is dominated by the
r correspond to the coded symbols A and the received signal frequency offset estimation part, since the frequency update
Y, respectively. requires a numerical search (see (34)), while the channel
1) Expectation step: The conditional LLF (conditioning is update is closed-form (see (35)). In this section we propose
on A) of H and F is given by the use of closed-form frequency offset estimation techniques
to avoid the numerical complexity associated with (34).
L−1 Note the resemblance of the update equation (34) to the
1   2
ln p (Y|H, F, A) = − yk − Hak ej2πF kT  (30) data-aided ML equation (14). A very similar derivation yields
N0 the following code-aided (CA) frequency offset estimator as
k=0
After some straightforward calculation, the E-step in the EM an alternative for (34).
algorithm can be written as (31) (see top of next page) where    
we defined M  (i)  (i)
1 m=1 m RCA;H,F (m)  arg R CA;H,F (m)
(i+1)
F̂CA1 = M  
   2πT 2  (i) 
(i)
 m R (m) 
ak aH
k = Ea ak aHk Y, b̂
(i)
(32) m=1 CA;H,F
  (36)
ak (i) = Ea ak | Y, b̂(i) . (33) As demonstrated in section III-B.2, this estimator is far less
complex then the estimator (34). On the other hand, the
and tr(.) denotes the trace of the matrix. We dropped terms reduced complexity DA estimator (17) experiences a flooring
not depending on the parameters to be estimated. effect for randomly selected pilot symbols. To avoid this
The computation of the expectations
in (32) and (33) behavior, an orthogonal training sequence was proposed. One
requires the posterior probabilities p ak |Y, b̂(i) , ∀k of the can expect that (36) will suffer from a similar flooring effect,
symbolvectors. For uncoded systems these are easily obtained since coded symbols are very unlikely to be orthogonal.
from the observation, whereas for coded system these are Hence, we introduce another ad hoc method to reduce com-
provided by a MAP-decoder. plexity of the frequency update (34). In fact it is possible to use
1564 IEEE TRANSACTIONS ON WIRELESS COMMUNICATIONS, VOL. 5, NO. 6, JUNE 2006


Q b| b̂(i)
 
= EA ln p ( Y| A, b)| Y, b̂(i)
L−1 L−1 
 (i)  (i)
= −tr H ak aH
k HH + 2 yk e−j2πF kT aH
k HH . (31)
k=0 k=0

TABLE III
T HE COMPLEXITY OF THE CODE - AIDED ESTIMATORS (34), (36) AND (37) WITHIN EACH ITERATION .

complex additions complex multiplications

L(L−1) ν L(L−1)
EM ν log ν + (NR + NT − 2) log ν + (NR + NT + 1)
 2  2  2 
M −1
CA1 M L− 2
(N T + NR − 2) M L − M2−1 (NT + NR + 1)
   
CA2 M L − M2−1 + L(NR − 1)(NT − 1) M L − M2−1 + 2LNR NT

the frequency offset estimator (9) for a known channel, since a

channel estimate is available from the previous EM iteration.
Estimator (9) does not suffer from a MSE flooring effect for B it DE CODE R
DE T E CTOR Π− 1
randomly selected pilot symbols, hence we do not expect a B CJR

MSE-floor for coded symbols either. In each EM iteration, the

frequency estimate is updated using the channel estimate from E M-iterations
previous iteration. The channel is consecutively updated using
this new frequency offset. The resulting frequency estimator Π
has the following form: B it Prob
E ST IMAT ION
M     to
 (i)  (i) Ĥ , Fˆ
1 m=1 m R CA;F (m)  arg RCA;F (m) PIL OT
Symbol Prob
(i+1)
F̂CA2 = M  
2πT 2  (i) 
m=1 m RCA;F (m) 
(37) Fig. 4. EM based receiver.
where
L−1
 and the current estimates of H and F . These likelihoods are
(i)
RCA;F (m) = y (i) (k − m)∗ y (i) (k) (38) fed to a MAP decoder. Using the BCJR algorithm [18], the
k=m decoder computes extrinsic bit probabilities of the transmitted

H bits, which are converted to symbolvector probabilities. Fi-
(i)
and y (i) (k) = Ĥ(i) āk yk . The channel estimates are nally, new parameter estimates, computed from the posterior
obtained according to (35). The number of arithmetic oper- probabilities, are provided to the detector. Note that the
ations required by estimators CA1 and CA2 is indicated in synchronization algorithms are also fit for other (space-time
Table III. We note that CA2 requires less computations2 for coded) MIMO set-ups as long as the computation of the a
large L (reduction of a factor NT + NR − 2). The performance posteriori probabilities is feasible.
of the estimators is discussed in next section. We have carried out computer simulations for a rate 1/2
recursive systematic binary convolutional coded system with
D. Simulation results octal generators (7, 5)8 . We applied BP SK signaling on a
In this section we illustrate the performance of the algo- NT = 4 , NR = 2 antenna set-up. Each block consists of
rithms presented in section IV-B and IV-C for a particular set- 32 orthogonal pilot symbolvectors and 128 coded symbols
up. We briefly describe the transmitter and receiver structure. per transmit antenna, and our simulations involved 500.000
A detailed outline of the configuration can be found in [4]. At blocks. All code-aided estimators are initialized with the DA
the transmitter, a block of information bits is channel-encoded, estimator (24).
bit-interleaved by a random interleaver and mapped onto a Fig. 5 and Fig. 6 show performance plots of the DA
normalized Mc -point signaling constellation Ω. The resulting algorithm (24), the EM algorithm using (34) and the modified
frame, is applied to a serial-to-parallel (S/P) converter, which versions CA1 (36) and CA2 (37), where we chose M =
L
feeds the NT transmit antennas. 4 = 32. The MSE-performance obtained after 5 EM-iterations
Fig. 4 illustrates the receiver structure. Each EM iteration is compared to the Modified Cramer Rao Bound (MCRB),
entails the following steps: The detector determines the mar- which is a lower bound on the estimation error of unbiased
ginal likelihood of the coded bits, based on the received signal estimators in presence of nuisance parameters [19]. We note
that estimator CA1 experiences a flooring effect in the MSE-
2 Parameter M is typically proportional to L. performance because of the randomness of the coded symbols
SIMOENS and MOENECLAEY: REDUCED COMPLEXITY DATA-AIDED AND CODE-AIDED FREQUENCY OFFSET ESTIMATION 1565

−5
10
DA
the channel is unknown, however, it is possible to achieve ML
EM CA1
EM CA2 performance provided that the training sequence is carefully
EM
−6
MCRB selected.
10
Furthermore, we have derived an iterative EM-based code-
aided estimation scheme for the joint frequency offset and
−7
10
channel estimation. To reduce the complexity in the compu-
tation of the frequency offset update estimate, two alternative
MSE F

algorithms CA1 and CA2 were presented. The resulting

−8
10 performance of both estimators at low to moderate SNR is
close to the EM performance. At high SNR, however, we
observe a minor performance degradation for CA1 that is
caused by self-noise, while the performance of CA2 does not
−9
10

degrade at all.
−10
10
0 5 10 15
E /N (dB)
b 0
ACKNOWLEDGMENT
The first author gratefully acknowledges the support
Fig. 5. MSE performance of the data-aided estimation algorithm (24) from the Fund for Scientific Research in Flanders (FWO-
using orthogonal pilot symbols compared to the iterative EM-algorithm and
approximations (36) and (37). (L = 128, NT = 4, NR = 2, ν = 1024)
Vlaanderen).

A PPENDIX I
(as explained in sections III-B.2 and IV-C). Such behavior In this appendix, we prove that the ML-estimate F̂ML;H,F
is not observed for estimator CA2. The BER performance equals the actual frequency offset F in the absence of noise.
resulting from the various estimators is shown as a function of Proof: In the absence of noise, the received signal has the
the iteration number in Fig. 6. Above a certain SNR threshold, following form
the CA2 estimator yields a BER performance that is very
close to the BER corresponding to known H and F . For high yk = Hak ej2πF kT (39)
SNR, estimator CA1 yields a BER floor, that is caused by the
self-noise in the estimates. As explained in section III-B.2, Setting the derivative of the function to be maximized in (14)
this error floor drops as the number of symbols increases. equal to zero, and substituting (39), we obtain for F̂ = F
Furthermore, since this flooring effect only occurs at very low L−1

values of the BER, estimator CA1 will yield a satisfactory (k − l)aH H
k H Hal A(l, k) = U − U
H

performance in many scenarios. k,l=0

As explained above, we can avoid the flooring effect in any where
case by resorting to estimator CA2. Estimator CA2, which
uses channel information from the previous EM iteration, does
L−1
L−1  L−1 −1
not suffer from convergence problems when using a sufficient   
amount of training symbols to initialize the EM algorithm. The U = kaH H
k H H al aH
l an aH
n ak
k=0 l=0 n=0
latter is not considered as a shortcoming since a reasonable
L−1

amount of pilot symbols is needed anyway (also when using
= kaH H
k H Hak
the optimal EM estimator (34)) to prevent the EM algorithm
k=0
from converging to a local maximum of the LLF (27). In
fact, the LLF with respect to the frequency offset has many It is readily seen that U − U H = 0, confirming that F solves
local maxima, which can only be avoided through proper (14) and F̂ML;H,F = F in the absence of noise. 
initialization.
We conclude that, while interesting to see that the CA2 A PPENDIX II
algorithm can achieve optimal MSE-performance, the actual Proof of proposition 1:
BER-performance obtained by algorithm CA1 will be suffi- if all symbolvectors are chosen from an orthogonal set:
cient for most practical implementations. ak ∈ {α0 , α1 , . . . , αNT −1 }, where αHi αj = NT Es δi−j , and
every vector αl is chosen at least once, then arg(C(m)) is
V. C ONCLUSIONS equal to zero, where C(m) is defined  in section
 III-B.2.
H −1
We have investigated data-aided and code-aided frequency Proof: First we prove that αH k a
i i ia αl = 0 for
offset estimation in a MIMO context. We have indicated that k
= l. It is readily seen that, since all vectors ak belong to
ML estimation is computationally hard. To overcome this the orthogonal set {α0 , α1 , . . . , αNT −1 },

inconvenience, we have developed new reduced complexity 
estimation algorithms. ai aH
i αl = NT Es γl αl
Considering data-aided estimation, we have shown that a i

closed-form estimator can achieve ML performance irrespec- where γl ∈ N is the multiplicity of αl , i.e. the num-
tive of the training sequence if the channel is known. When ber of times αl appears in the transmitted sequence
1566 IEEE TRANSACTIONS ON WIRELESS COMMUNICATIONS, VOL. 5, NO. 6, JUNE 2006

0 0
10 10
DA DA
EM CA1 iter1 EM CA2 iter1
EM CA1 iter2 EM CA2 iter2
−1 EM CA1 iter3 −1
10 EM CA2 iter3
10
EM CA1 iter4 EM CA2 iter4
EM CA1 iter5 EM CA2 iter5
perfect synchro perfect synchro
−2 −2
10 10

−3 −3
10 10

BER
BER

−4 −4
10 10

−5 −5
10 10

−6 −6
10 10

−7 −7
10 10
0 5 10 15 0 5 10 15
E /N (dB) Eb/N0 (dB)
b 0

Fig. 6. BER performance of the data-aided estimation algorithm (24) using orthogonal pilot symbols compared to the iterative EM-algorithm with
approximations (36) and (37). (L = 128, NT = 4, NR = 2, ν = 1024)

{ak , k = 0, . . . , L − 1}. Note that γl ≥ 1.

This illustrates
 R EFERENCES
H
that αl is an eigenvector of the matrix i ai a
i with
H
corresponding eigenvalue NT Es γl . We represent i ai ai
[1] G. J. Foshini and M. J. Gans, “On limits of wireless communication
in a fading environment when using multiple amtennas,” Wireless Pers.
according to its eigenvalue decomposition, Commun., vol. 6, no. 3, pp. 311–335, Mar. 1998.
 [2] A. Stefanov and T. M. Duman, “Turbo-coded modulation for systems
ai aH
i = DΛD
−1
with transmit and receive antenna diversity over block fading channels:
i system model, decoding approaches, and practical considerations,” IEEE
J. Select. Areas Commun., vol. 19, no. 5, pp. 958–968, May 2001.
where D = [α0 , α1 , . . . , αNT −1 ], and Λ is a diagonal matrix [3] V. Tarokh, N. Seshadri, and A. R. Calderbank, “Space-time codes for
high data rate wireless communications: Performance criterion and code
containing the eigenvalues. It is now easy to calculate the construction,” IEEE Trans. Commun., vol. 44, no. 2, pp. 744–765, Mar.
inverse 1998.
[4] J. J. Boutros, F. Boixadera, and C. Lamy, “Bit-interleaved coded
−1 modulations for multiple-input multiple-output channels,” in Proc. IEEE
 1 6th Int. Symp. on Spread-Spectrum Tech. & Appli., New Jersey, USA,
ai aH
i = DΛ−1 D−1 = DΛ−1 DH Sept. 2000, pp. 123–126.
i
NT Es [5] D. C. Rife and R. R. Boorstyn, “Single tone parameter estimation from
discrete-time observations,” IEEE Trans. Inform. Theory, vol. 20, pp.
Again making use of αH i αj = NT Es δi−j , we obtain after
591–598, Sept. 1974.
[6] H. Meyr, M. Moeneclaey, and S. A. Fechtel, Synchronization, Channel
some straightforward calculation. Estimation, and Signal Processing, vol. 2 of Digital Communication
Receivers. John Wiley & Sons, 1997.
−1
 [7] M. Luise and R. Reggiannini, “Carrier frequency recovery in all-digital
αH
k ai aH
i αl = γl−1 δk−l modems for burst-mode transmissions,” IEEE Trans. Commun., vol. 43,
pp. 1169–1178, Feb./Mar./Apr. 1995.
i [8] M. P. Fitz, “Further results in the Fast Estimation of a Single Frequency,”
 
H −1
IEEE Trans. Commun., vol. 42, pp. 862–864, Feb./Mar./Apr. 1994.
indicating that αHk i ai ai αl = 0 when k
= l, [9] O. Besson and P. Stoica, “On parameter estimation of MIMO flat-fading
allowing us to simplify (19) channels with frequency offsets,” IEEE Trans. Signal Process., vol. 51,
pp. 602–613, Mar. 2003.
[10] Y.-D. Kim, J. K. Lim, C.H. Suh, E.-R. Jeong, and Y. H. Lee, “Carrier
L−1
 frequency estimation for transmissions with antenna diversity,” in IEEE
C(m) = aH H
k−m H Hak A(k, k − m)
Veh. Technol. Conf., Birmingham, 2002, pp. 1569–1573.
[11] A. P. Dempster, N. M. Laird, and D. B. Rubin, “Maximum likelihood
k=m
 from incomplete data via the EM algorithm,” Journal of the Royal
−1
= aH H
k H Hak γk
Statistical Society, vol. 39, no. 1, pp. 1–38, 1977, Series B.
[12] Z. Baranski, A. M. Haimovich, and J. Garcia-Frias, “EM-based iterative
k: ak =ak−m receiver for space-time coded modulation with noise variance estima-
tion,” in Global Telecommunications Conference, Nov. 2002, pp. 355–
since A(k, k − m) = 0, for ak
= ak−m . It is readily seen that 359.
C(m)H = C(m), and thus arg(C(m)) = 0.  [13] C. Cozzo and B. L. Hughes, “Joint channel estimation and data detection
Remark: If the vectors
 α l have the same multiplicity, i.e. in space-time communications,” IEEE Trans. Commun., vol. 51, no. 8,
pp. 1266–1270, Aug. 2003.
γl = K, ∀l, the matrix i ai aH i becomes a diagonal matrix, [14] D. K. Hong, Y.-J. Lee, D. Hong, and C. E. Kang, “Robust frequency
since offset estimation for pilot assisted packet CDMA with MIMO antenna
systems,” IEEE Commun. Lett., vol. 6, no. 6, pp. 262–264, June 2002.
[15] F. Simoens and M. Moeneclaey, “A reduced complexity frequency offset

= NT Es K DD−1 = LEs I.
estimation technique for flat-fading MIMO channels,” in IEEE CAS
ai aH
i = DΛD
−1
Symposium on Emerging Technologies, Shanghai, May 2004, pp. 705–
i 708.
SIMOENS and MOENECLAEY: REDUCED COMPLEXITY DATA-AIDED AND CODE-AIDED FREQUENCY OFFSET ESTIMATION 1567

[16] M. Morelli and U. Mengali, “Carier-frequency estimation for transmis- Marc Moeneclaey received the Diploma and the
sions over selective channels,” IEEE Trans. Commun., vol. 48, no. 9, Ph.D. Degree, both in electrical engineering, from
pp. 1580–1589, Sept. 2000. Ghent University, Gent, Belgium, in 1978 and 1983,
[17] T. L. Marzetta, “BLAST training: estimating channel characteristics for respectively. He is currently a Professor in the
high capacity space-time wireless,” in Ann. Allerton Conf., Monticello, Department of Telecommunications and Informa-
pp. 958–966, Sept. 1999. tion Processing at Ghent University. His main re-
[18] L. R. Bahl, J. Cocke, F. Jelinek, and J. Raviv, “Optimal decoding of search interest are in statistical communication the-
linear codes for minimising symbol error rate,” IEEE Trans. Inform. ory, carrier and symbol synchronization, bandwidth-
Theory, vol. 20, pp. 284–287, Mar. 1974. efficient modulation and coding, spread spectrum,
[19] A. N. D’Andrea, U. Mengali and R. Reggiannini, “The modified cramer- satellite and mobile communication. He is the author
Rao bound and its application to synchronization parameters,” IEEE of more than 250 scientific papers in international
Trans. Commun., vol. 42, pp. 1391–1399, Feb./Mar./Apr. 1994. journals and conference proceedings. Together with H. Meyr (RWTH Aachen)
and S. Fechtel (Siemens AG), he is coauthor of the book ”Digital Communica-
Frederik Simoens received the diploma of electrical tion Receivers - Synchronization, Channel estimation, and Signal Processing”
engineering in 2003 from Ghent University, Gent, (New York: Wiley, 1998)
Belgium, where he is currently working toward the
Ph.D. degree in the Department of Telecommunica-
tions and Information Processing. His main research
interests include parameter estimation, modulation
and coding for wireless digital communications.