Blind recursive tracking of carrier frequency offset (CFO) vector in MC-CDMA systems

Transactions Papers

Blind Recursive Tracking of Carrier Frequency

Offset (CFO) Vector in MC-CDMA Systems

Feng-Tsun Chien, Member, IEEE, and C.-C. Jay Kuo, Fellow, IEEE

Abstract— A recursive algorithm for estimating and

updat-ing the effective carrier frequency offset (CFO) vector in a multicarrier code-division multiple-access (MC-CDMA) system is proposed in this work. The recursive relation is derived based on the expectation maximization (EM) algorithm with a quadratic constraint. This new approach enables the use of linear estimation theory to tackle the CFO estimation problem with or without training data, which leads to an analytic CFO estimate in closed form. Furthermore, the multiple access interference (MAI) is mitigated using the second order statistics of the interference-plus-noise vector, which is updated in a recursive manner under the EM formulation, too. When reaching a converged estimate, a fixed-norm quadratic constraint is imposed so that the final CFO estimate is robust to an imprecise covariance matrix estimate caused by insufficient data samples. It is demonstrated by computer simulation that the performance of an MC-CDMA system without the CFO information can be restored by the proposed scheme in the sense that its bit error probability (BEP) performance is close to that with perfect CFO knowledge.

Index Terms— Code division multiple access (CDMA),

multi-carrier CDMA, multi-carrier frequency offset (CFO), recursive EM algorithm, quadratic constraint.

I. INTRODUCTION

W

(CDMA) and the multicarrier modulation systems, the idea to integrate merits from both techniques has stimulated a large amount of research during the past decade. In particular, the multicarrier-CDMA (MC-CDMA) system is a potential candidate for the fourth generation (4G) wireless standard thanks to a number of promising features such as robustness to frequency selective fading channels, immunity to inter-symbol interference (ISI) by adding/removing the cyclic prefix and flexibility in allowing multiple access [1]. However, multicarrier systems, including orthogonal frequency division multiplexing (OFDM) and MC-CDMA, are sensitive to the

Manuscript received May 22, 2005; revised June 28, 2006; accepted August 26, 2006. The associate editor coordinating the review of this paper and approving it for publication was M. Uysal.

F.-T. Chien is with the Department of Electronics Engineering, National Chiao Tung University, Hsinchu, Taiwan, R.O.C. (email: ftchien@mail.nctu.edu.tw).

C.-C. Jay Kuo is with the Integrated Media Systems Center and the Department of Electrical Engineering, University of Southern California, Los Angeles, CA 90089-2564 (email: cckuo@sipi.usc.edu).

Digital Object Identifier 10.1109/TWC.2007.05380.

carrier frequency offset (CFO) effect. A prominent CFO destroys the orthogonality between sub-carriers, causing the undesirable inter-carrier interference (ICI) and dramatically deteriorating the system performance.

The range of CFO, which may come from mismatches be-tween local oscillators and/or from mobility induced Doppler shifts, in practical multicarrier applications is often too large to be acceptable. For example, the frequency accuracy 25 parts per million (ppm) specified in IEEE 802.11g can be translated

to 2.4 GHz×25 ppm = 60 KHz frequency shift. Compared

to the sub-carrier spacing 312.5 KHz, the normalized CFO is 0.192, which is well above the 0.01 threshold in maintaining a tolerable bit error probability (BEP) performance [2] and must be compensated by means of signal processing techniques in a later stage. It is expected that, when projected into 4G applications where the MC-CDMA system may come into play, there will be an even larger normalized CFO due to reduced subchannel bandwidth.

Motivated by the above observation, there have been a number of studies dedicated to the development of effective algorithms for CFO estimation in the OFDM system [2]– [8]. However, all these schemes cannot be directly applied to the MC-CDMA system due to the presence of multiple access interference (MAI). In addition to the MAI problem, the underlying nonlinear structure makes the CFO estimation a challenging issue in the receiver design for MC-CDMA. Only a limited amount of research in the literature can be found in this regard. A block-based joint CFO and channel estimation scheme for single carrier CDMA systems was considered in [9] based on the subspace projection and the polynomial root-finding algorithms. A blind channel-independent block-based approach, which is similar to [9], was developed in [10] to find a CFO estimate for MC-CDMA systems, wherein each user’s channel impulse response was implicitly assumed to be orthogonal. In [11], an iterative receiver structure for MC-CDMA was developed using the technique of generalized sidelobe cancellation, where an estimate of CFO was obtained relying on exhaustive search.

On the other hand, the application of the expectation max-imization (EM) algorithm to parameter estimation for OFDM systems has been extensively studied recently, thanks to its

power in dealing with a broad range of estimation problems with incomplete observations [12], [13]. Estimation of the symbol arrival time and the carrier phase was considered in [14] with the transmitted symbol being treated as the missing data. A blind CFO estimation scheme in MC-CDMA systems using the EM algorithm was studied in [15], where the gradient decent technique was proposed to deal with the nonlinear optimization problem in the CFO estimate. However, since there exist multiple local optima in the cost function, the solution is sensitive to the initialization of the optimiza-tion process, which could be problematic for a hill-climbing adaptive algorithm. Besides, all these schemes were operated in a block-based manner assuming that the target parameters were time invariant inside this time block. However, timely updates of system parameters would be more desirable in real time applications, which motivates the development of online algorithms to track possibly time-varying unknowns.

In this work, we address the blind CFO estimation problem in MC-CDMA systems using the sequential EM algorithm [16], [17] and the constrained optimization technique. The EM algorithm provides an iterative procedure to find the maximum likelihood (ML) estimates of target parameters with appropriate initializations. It allows us to treat the CFO estimation problem blindly, i.e. without the aid of training data. The algorithm in the training mode is a direct extension from the blind mode. Specifically, rather than finding CFO directly, we derive a recursive relation of the estimate of a time-varying vector consisting of the exponentials of the desired user’s CFO, which is referred to as the CFO vector in the sequel. By doing so, the nonlinear CFO estimation problem can be avoided. That is, the resultant cost function in the E-step is quadratic with respect to this CFO vector, which makes the M-step easier to cope with and leads to an analytically tractable estimate for the CFO vector and, hence, a closed-form solution for the CFO itself. We find that no performance loss in terms of the system bit error probability is observed by examining the CFO vector and then estimating the offset. Furthermore, the joint effect of MAI and AWGN is modeled by a colored Gaussian vector and its effect is whitened using the inverse of its covariance matrix [18]– [20]. This covariance matrix can also be updated using the recursive EM formulation. The computational cost associated with matrix inversion is reduced using the matrix inversion lemma. Finally, the performance of the proposed estimator is demonstrated by computer simulations.

The rest of this paper is organized as follows. The system model of an MC-CDMA system is provided and the problem is formulated in Sec. II. The developed recursive algorithm with a quadratic constraint is discussed in Sec. III. Numerical simulations are presented in Sec. IV. Finally, concluding remarks are given in Sec. V.

II. SYSTEMMODEL

A. Continuous-Time Signal Model

Consider an MC-CDMA system with N sub-carriers and bandwidth W . The transmitted signal of user k after perform-ing the inverse discrete Fourier transform (IDFT) and addperform-ing

the cyclic prefix can be represented as dk(t) = ∞  m=−∞ bk[m] N−1_ n=−NG pk,nψ(t− mTs− nT ), (1)

where bk[m] is the transmitted symbol of user k at the mth

symbol block, NG is the number of samples employed in the

cyclic prefix, Ts= (NG+ N )· T is the time span of the mth

transmitted symbol block with T = 1/W being the duration

of each time-domain chip pk,n, ψ(t) is the pulse-shaping

waveform assumed to be rectangular of unity amplitude and duration T in this paper and

pk,n= 1

N

N−1_ i=0

ck,iej2πNin

is the nth IDFT output of the signature sequence ck =

[c_k,0, c_k,1· · · c_k,N−1]T of user k with ck,i ∈ {−1, +1} for

0≤ i ≤ N − 1. Note that we will use smto represent b1[m]

in the sequel for notational convenience.

It is assumed that the system experiences a frequency selective fading channel and the following tapped-delay line model is adopted for the channel impulse response [21]

hk,m(τ ) = L  l=0 hk,m[l]· δ(τ − l W),

where hk,m[l]is the fading coefficient of user k on the lth path

at the mth symbol block, δ(t) is the Kronecker delta function and L is the order of the channel depending on the maximum multipath delay spread of the channel. The fading gain hk,m[l]

is modeled as a zero mean complex-valued Gaussian random process. Then, the received signal for user k in the absence of ambient noise is given by

rk(t) = dk(t)  hk(τ ) = ∞  m=−∞ bk[m] N−1_ n=−NG L  l=0 pk,nhk,m[l] × ψ(t − mTs− nT − lT ),

where  is the convolution operation. Note that the received signal rk(t)has been derived above in the absence of CFO. A

more realistic model which encompasses the CFO effect will be considered in Section II-C.

B. Discrete-Time Signal Model

At the receiver, after passing through the matched filter with a sampling rate 1/T , the discrete time signal observed at the mth symbol block and the gth chip interval for user k can be written as rk,m[g] = 1 T  mTs+(g+1)T mTs+gT rk(t)ψ(t− mTs− gT )dt, (2)

for −NG ≤ g ≤ N − 1. The above matched filtering only

involves the following integration  mTs+(g+1)T

mTs+gT

which takes nonzero values when g = n + l. This yields rk,m[g] = bk[m] 

n+l=g

pk,nhk,m[l],

which is equivalent to the result of performing

discrete-time convolution between the channel coefficient hk,m =

[hk,m[0], hk,m[1]· · · hk,m[L]]T and the time domain signature

sequence pk = [pk,0, pk,1· · · pk,N−1]T. After removing the

cyclic prefix, the received signal vector at the mth symbol block is equal to rk,m= [rk,m[0], rk,m[1]· · · rk,m[N− 1]]T = Hk,mpk· bk[m] = 1 NW H NΛk,mck· bk[m]

where Hk,m is the N × N right circular matrix with

[hk,m[0], hk,m[1]· · · hk,m[L], 0· · · 0]T being its first column and superscript H denotes the Hermitian transposition. It is

easy to see that Hk,m = _N1WHNΛk,mWN, where WN is

the standard DFT matrix with its (m, n)th element equal to e−j2πN(m−1)(n−1) _{and Λk,mis the diagonal matrix composed}

of the N-point DFT of hk,m. The discrete-time model for

the received signal with K active users in the system can be represented by rm= 1 N K  k=1 WHNΛk,mck· bk[m] + wm,

where wm is a complex additive white Gaussian noise

(AWGN) vector with zero mean and covariance matrix E[wmwHm] = σ2I. Note that, for Hk,m to be a circular matrix, the length of the cyclic prefix should be no less than the channel order. More generally, the length of the cyclic prefix is designed to absorb both the channel order as well as the residual timing delay estimation errors to avoid ISI. Here,

we take the minimum value for NG, i.e. NG = L. In other

words, we do not consider asynchronous transmission delays. Furthermore, we assume perfect channel knowledge in this work.

C. Problem Formulation

In the presence of carrier frequency offset, the received signal for user k needs to be modified by a phase shift as rk(t)ej2πΔfkt_{,where Δfk is the CFO of user k. The matched}

filter output for is therefore given by 1 T  mTs+(g+1)T mTs+gT ψ(t− mTs− nT − lT ) × ψ(t − mTs− gT ) · ej2πΔfktdt = ρk· ej2πN(m(NG+N)+12+g)εk _{for g = n + l, (3)}

where εk = Δfk/(W/N ) denotes the kth user’s normalized

CFO, which is a relative measure of the amount of CFO in Hertz compared to the sub-carrier spacing W/N and is assumed to be deterministic as well as time-invariant, and ρk = sincπ(ε_Nk), which is approximately equal to one with

a large number of sub-carriers. It is assumed that the absolute value of the normalized CFO is no larger than one half of the sub-carrier spacing, i.e. |εk| < 0.5. For a large number

of sub-carriers and the assumption|εk| < 0.5, we have scalar

ρk ≈ 1. The approximated result ej2πN(m(NG+N)+12+g)εk _in

(3) is a time-varying scalar reflecting the effect caused by the CFO, which not only rotates the symbol constellation by a different amount at different symbol timing block m, but also ruins the orthogonality between sub-carriers.

The received signal vector rm of the MC-CDMA system

in the presence of CFO observed at the mth symbol block after removing the cyclic prefix of a length equal to L can be expressed by rm= 1 N K  k=1 Fk,mWHNΛk,mck· bk[m] + wm, (4)

where the effect of CFO for the kth user is modeled by the diagonal matrix

Fk,m= κk,m· diag



1, ej2πεk/N_{,· · · , e}j2π(N−1)εk/N_{, (5)}

where κk,m= ej2πN(m(NG+N)+12)εk _{is the scale factor (which}

is a function of timing index m). Apparently, the orthogonality between sub-carriers has been ruined due to the presence of Fk,m, i.e.

WNFk,mWNH= N · I,

where the resulting off-diagonal terms characterize the ICI for all sub-carrier pairs. Therefore, before facilitating the DFT operation in the receiver, this residual CFO needs to be compensated. Without loss of generality, we assume the first user to be the user of interest and its signal component can be written as dm= Gm· fm,where

Gm= diag{ 1

NW

H

NΛ1,mc1· sm}

is a diagonal matrix describing the system structure with sm= b1[m] corresponding to the scalar transmitted symbol at the

mth symbol block, and

fm= κk,m·



1, ej2πεk/N_{,· · · , e}j2π(N−1)εk/NT

is the CFO vector consisting of the diagonal entries of F1,m.

The matrix Gm collects all the parameters other than CFO,

and is generally unknown to the receiver. The signal compo-nents from all interfering users and the AWGN are modeled by a colored Gaussian vector nm=N1

K

k=2Fk,mWHNΛk,mck·

bk[m] + wm with covariance matrix R = E nmnHm

.

Therefore, we have a generic received signal model for the mth time interval

rm= Gm· fm+ nm. (6)

It is clear from (5) that the CFO vector fm satisfies the

following dynamic evolution relation

fm= E · fm−1, (7)

where E = ej2π(1+NG/N)ε1_{I. This dynamic evolution plays an}

important role in the derivation of the recursive EM algorithm for updating fm.

In this work, rather than obtaining a direct inference to

for the CFO vector fm. In a practical scenario, system

para-meters, including CFO’s, channel impulse responses and the interference-plus-noise vector correlation matrix, are generally unknown to the receiver. For the sake of clarity and brevity of the presentation, only the estimates for the CFO vector fmas

well as the interference-plus-noise vector correlation matrix R are considered. In other words, along with the assumption of synchronous transmission, we also assume perfect knowledge of the channel impulse response here. Our objective is the joint estimation of fmand correlation matrix R in a recursive manner based on all received signals up to time m without the aid of the training sequences. The maximum a posteriori (MAP) symbol detection can also be achieved at each EM iteration.

It is however worthwhile to emphasize that the estimation of possibly time-varying channel impulse responses can be incorporated in the iterative procedure using the expectation conditional maximization (ECM) algorithm [22]. Performance loss is expected due to possible phase ambiguity in the esti-mation of channel coefficients [23]. Besides, the convergence rate is likely to slow down with more unknown parameters included in the ECM iterations [22]. On the other hand, when considering timing synchronization errors, the received signal with ISI under timing inaccuracy can still be modeled by the generic representation in (6), which implies that the proposed algorithm for the CFO vector update can work as well with timing inaccuracy. Under such circumstances, system matrix

Gm needs be modified to encapsulate the effects of ISI,

channel impulse responses and timing errors [24, eq. (10)], and should be updated using the ECM algorithm.

III. CONSTRAINEDRECURSIVE ALGORITHM

In this section, we first outline the proposed recursive blind algorithm for fmwith full knowledge of the covariance matrix R. Then, an approach to determine the Lagrange multiplier for a quadratically constrained optimization problem is detailed. Finally, we relax the assumption of knowing R, and consider

joint estimation of fmand R based on the ECM algorithm.

A. Recursive EM Formulation

Let ym = [rT1rT2 · · · rTm]T be all the received signals up

to time m and θm = [f1Tf2T· · · fmT]T. Then, assuming that

covariance matrix R is available, the updated inference f_m|m(i+1) of fmbased on ymat the (i + 1)th iteration can be effectively

accomplished using the EM algorithm as stated below [16], [17]: E-step: Find Cm(θm, θ(i)_m|m) Exm  log p (ym, xm; θm)ym; θ(i)_m|m  , (8) M-step: Solve

f_m|m(i+1)= arg max

fm Cm

(θm, θ(i)_m|m), (9)

where θ_m|m(i) = [f_m|m(i)T, f_m−1|m−1(c)T · · · f_1|1(c)T]T with f_l|l(c) rep-resenting the converged estimate of fl based on yl, p(·; θ)

and E_xm[·; θ] are the probability density function (pdf) and

the expectation operator averaging over xm, respectively,

parameterized by deterministic θ, xmis the vector of missing

data which is essential to characterize the incomplete-data

likelihood function p(ym; θm) and is averaged in the

E-step. The iterative procedure of the EM algorithm guarantees convergence. With appropriate initializations, the EM iteration will converge to the ML estimate [22]. The choice of missing data is not unique in general. For the blind scenario considered in this research, xm is chosen to be the transmitted signals up to time m; i.e. xm = [s1· · · sm]T. In contrast, in the

training mode, xm belongs to the empty set. Consequently,

in the training mode, we can remove the expectation in (8). The complete-data, composed by the received data vector,

ym, and the unknown missing data vector, xm, the log

likelihood function of θmin (8) can be computed as

log p (ym, xm; θm) = log pymxm; θm + log p (xm)

where log p (xm)does not depend on θmand can be discarded

in the EM iteration. Therefore, we can equivalently state the EM iterative procedure in the following two steps:

E-step: Find Qm(θm, θ(i)_m|m) Exm  log pymxm; θm ym; θ_m|m(i)  , (10) M-step: Solve

f_m|m(i+1)= arg max

fm Qm

(θm, θ(i)_m|m). (11)

Having reached the converged estimate f_m|m(c) , the final esti-mate of the CFO vector at the mth symbol block becomes

ˆ_{fm|m= arg max}

||fm||2=NQm

(θm, θ_m|m(c) ), (12)

where the constant 2-norm constraint is imposed due to the particular structure of the CFO vector; namely, each element

of fm lies on the unit circle in the complex domain. This

quadratic constraint makes the CFO vector estimate more robust to the error of the estimated covariance matrix of the interference-plus-noise vector [25].

Although it is feasible in finding the newly update of (11) and the final estimate (12) at each symbol block interval m with batch processing, a recursive structure is more desirable for lower computational cost and real time implementations. By applying the Taylor series expansion to Qm(θm, θ(i)_m|m)at

f_m|m(i) , we can derive a recursive relation for the newly update f_m|m(i+1) in (11) as [16, eq. (3.20)], [26, theorem 2]

f_m|m(i+1)= f_m|m(i) − ⎛ ⎝∂2Qm(θm, θ(i)m|m) ∂fm2    fm=f_m|m(i) ⎞ ⎠ −1 × ⎛ ⎝∂Qm(θm, θm|m(i) ) ∂fm∗    fm=f_m|m(i) ⎞ ⎠ , (13)

where the partial derivative relative to complex variable z is defined as _∂z∂  1₂  ∂ ∂(z) − j∂(z)∂  with (z) and (z)

representing the real and the imaginary parts of z, respec-tively, and the 2nd order derivative is defined as ∂2/∂fm2 

∂2/∂fm∗∂fmT with the superscript * denoting the complex

conjugate [17]. Note that (13) is an equality instead of an approximation sinceQm(θm, θ(i)_m|m)is quadratic with respect to fm.

B. Blind Recursive EM Algorithm

We describe the procedure to obtain the final estimate in (12) and the recursive update in (13) in this subsection. We focus on the development of the recursive algorithm without the aid of training symbols. The training data based approach can be straightforwardly deduced by removing the expectation operator in (10) since the set of missing data is empty. Let

Pm ∂2Qm(θm, θ(i)_m|m) ∂fm2    fm=f_m|m(i) ,

where the evaluation at fm= f_m|m(i) actually has no effect on

the result of the matrix Pmdue to the quadratic structure of

Qm(θm, θ(i)_m|m)with respect to fm. It is shown in the appendix

that Pm has the following recursive structure

Pm= Pm−1− GHmR−1Gm. (14)

Then, carrying out the derivatives in (13) yields

f_m|m(i+1)= f_m|m(i) − P−1 m · Esm ¢ GH mR−1   rm− Gmf_m|m(i) ¡¬ ¬r m; f_m|m(i) £ , (15)

which is of the standard prediction and correction form. The above recursion is initialized by f_m|m(0) , which is the maximum likelihood prediction of fmbased on ym−1, the received signal

vector up to time m− 1. It can be shown from the dynamic

evolution given in (7) that

f_m|m(0) = E · f_m−1|m−1(c) , (16)

where transition matrix E is generally unknown to the receiver and should be estimated separately. The estimate of E, or equivalently the estimate of ej2π(1+L/N)ε1_{, at the mth time}

instance can be obtained according to the previously updated estimate f_m−1|m−1(c) through all its elements, as will be detailed later in subsection III-E.

Having reached the converged estimate f_m|m(c) during the mth symbol block time interval, we employ the constrained opti-mization in (12) to find the final estimate. With a replacement ofQm(θm, θ(i)_m|m)in (13) by

Jm(λm, fm) =Qm(θm, θ_m|m(c) ) + λm(||fm||2− N), (17)

which is the cost function of the constrained optimization problem using the Lagrangian multiplier technique and λmis the Lagrange multiplier at the mth symbol block time interval, we have the following prediction-correction form

ˆ fm|m= fm|m(c) −  Pm+ λmI −1 (18) ×Es_m[GHmR−1  rm− Gmf_m|m(c) rm; f_m|m(c) ] + λmf_m|m(c)  , as in (15). After some algebraic manipulations, (18) is reduced to ˆ fm|m=  Pm+ λmI −1 · Pm−1f_m|m(c) − Es_m[GmHR−1rmrm; f_m|m(c) ]  .(19) Note that, through the manipulations from (18) to (19), the

Lagrange multiplier λm is contained only in the inverse in

(19). This rearrangement is particularly useful in determining the Lagrange multiplier, as will be detailed in the next sub-section. We relegate the derivation of (19) to the appendix. C. Lagrange Multiplier

We discuss the determination of the Lagrange multiplier for the constrained optimization in (12) for each symbol block time interval in this subsection. More specifically, we need to determine the Lagrange multiplier λmin (19) in order to reach the final estimate ˆfm|mthat satisfies the constraint||ˆfm|m||2=

N.

It is clear from (14) that Pmis Hermitian. Thus, it can be

decomposed as

Pm= UmΓmUHm, (20)

where Um is an unitary matrix consisting of eigenvectors of

Pmand Γm= diag(γm,1· · · γm,N)with γm,1≤ · · · ≤ γm,N

being the corresponding eigenvalues. The Lagrange multiplier λmcan be determined as a solution to the equation||ˆfm|m||2=

N. Substituting (20) into (19) and taking the square of its norm gives g(λ) (Γm+ λI)−1qm 2 = N  l=1 |qm,l|2 (λ + γm,l)2 , (21) where qm= UHm  Pm−1fm|m(c) −Esm GHmR−1rmrm; f_m|m(c)

with qm,l being its lth element. It was shown in [27] that

the solution to the Lagrange multiplier of the constrained optimization problem in least squares formulation is equal to the largest real root of g(λ) = N , and can be resorted to root-finding algorithms for its solution. Here, we employ the Newton method and derive the upper and lower bounds for this optimal root. The obtained upper and lower bounds provide insights into the determination of a good initial guess for the Newton method.

Observing the relation in (21), it is clear that g(λ) is monotonically decreasing in (−γ_m,1,∞) since g(λ) < 0 in

this interval with g(−γ_m,1)→ ∞ and g(∞) → 0. Therefore,

there exists an unique solution ¯λ in this interval such that

g(¯λ) = N. Apparently, this is the largest real root of g(λ) =

N and is the solution to the Lagrange multiplier, which may

take the form

¯

λ =−γm,1+,

where is a positive real number. Replacing γm,l in (21) by

γm,1 and γm,N, we can obtain the following upper and lower bounds max  ||qm|| √ N − (γm,N− γm,1), 0  ≤ ≤ ||q√m|| N .

In the above interval lies a good initial guess for and, hence, ¯

λ. The upper bound is not a good candidate of the initial

guess λ1, since we have λ1 = −γ_m,1 +||qm||/√N at this

point, where g(λ1)has a mild slope, which is likely to result

in a newly update smaller than −γ_m,1 and converge to an

undesired root. Therefore, it is advised to initialize the Newton algorithm with the lower bound if it is not zero. Otherwise, the Newton algorithm is initialized by

λ1=−γm,1+ ||qm||

with a proper choice of positive integer p to set the initial

guess between the desired root and−γ_m,1, where the slope

of g(λ) is sharper.

D. Interference-plus-Noise Correlation Matrix

In Sec. III-A, it was assumed that correlation matrix R of the interference-plus-noise vector nm is known in the

devel-opment of the recursive algorithm for updating fm. However,

this is in general not true in practice. Here, we will show that R can also be jointly estimated in a recursive manner based on the ECM algorithm, which is a generalized EM algo-rithm aiming at finding the ML estimates of multidimensional parameters that do not have a simultaneous closed form to maximize the obtained cost function in the E-step of the EM algorithm [22], [28].

The ECM algorithm states that, conditioned on the previous estimate of one parameter, it is much easier to obtain a closed form representation of the update for another parameter. Likewise, conditioned on this new update, we can go back to update the original parameter. This procedure yields an analytically tractable way to iteratively update all parameters one by one and guarantees an increase of the likelihood as the iteration goes on.

For the problem under consideration, the set of parameters contains R and fm, and the ECM algorithm is applicable for

their joint estimation. More specifically, conditioned on the ith update of R(i)_|m and f_m|m(i) , the new update f_m|m(i+1)for the CFO vector can be derived using (13) with associated R’s being replaced by R(i)_|m. Then, having obtained f_m|m(i+1), the update of the correlation matrix in a sequential form is given by

R(i+1)_|m =  1− 1 m  R(c)_|m−1 + 1 mEsm 

n(i)_|mn(i)_|mHrm; f_m|m(i) , R(i)_|m 

, (22)

where R(c)_|m−1is the converged estimate of R based on all re-ceived signals up to time m−1 and n(i+1)_|m = rm−Gmf_m|m(i+1).

This recursive estimate can be initialized either by a warm-up period of length Np, or by an identity matrix multiplied with

a small scalar. With the recursive and iterative update for the correlation matrix R in (22), the calculation of matrix Pmin (14) should be modified as

P(i+1)m = P(c)m−1− GHmR(i+1)_|m −1

Gm. (23)

A final remark can be made that, by using the matrix inversion lemma repeatedly, the inversion of R(i+1)_|m only involves a pre-viously calculated matrix inverse and scalar inversions, which greatly reduces the complexity of the proposed algorithm. The computational cost mostly comes from calculating the inverse of P(i+1)m and the eigenvalue decomposition of P(c)m required at every time instance.

E. Symbol Detection

After arriving at the convergence stage with a final estimated parameter set f_m|m(c) , the maximum a posteriori probability

...

...

...

...

....

....

r

f

r

f

f

f

f

s

R

P

P

R

f

f

f

f

f

s

Fig. 1. The update order of estimates in the proposed recursive algorithm at each symbol blockm.

(MAP) detection for sm has already been realized in the

E-step via [29] ˆ

sm= arg max

sm∈A

P [sm|rm; ˆf_m|m(c) , R(c)_|m].

However, it should be noted that the EM algorithm does not guarantee convergence to the global optimum. This effect may lead to phase ambiguity in the CFO vector estimate in the complex domain and results in a total erroneous decision for symbol detection.

One way to resolve this problem is to employ the differential encoding/decoding scheme at the cost of 3 dB performance loss. Alternatively, thanks to the particular structure of the CFO vector, the phase ambiguity can be resolved without invoking differential detection, observing that the CFO vector

fm is actually a geometric sequence with a common factor

ej2πNε at every symbol block m. For N ≥ 2, this common

factor can be extracted by comparing adjacent entries of the updated CFO vector, regardless of the presence of the phase ambiguity. Therefore, the scalar CFO can be estimated via

ˆ ε_|m= 1 N− 1 N−1_ l=1 N 2πarg  (ˆfm|m(l))∗· ˆfm|m(l + 1)  , (24)

from which the CFO vector can be reconstructed without the phase ambiguity, and the transition matrix can be inferred, too. Then, the MAP symbol detection can be realized using

ˆ

sm= arg max

sm∈AP [sm|rm; ˆfm, R

(c)

|m], (25)

0 10 20 30 40 50 60 70 80 90 100 −2 −1 0 1 2

Tracking of the real part of the CFO vector

Estimate True value 0 10 20 30 40 50 60 70 80 90 100 −1.5 −1 −0.5 0 0.5 1 1.5 Time (m)

Tracking of the imag part of the CFO vector

Estimate True value

Fig. 2. The tracking behavior of the real and the imaginary parts of the first element off_m|m(1) in the training mode.

0 20 40 60 80 100 0 0.02 0.04 0.06 0.08 0.1 0.12 time (m) CFO ( ε1 ) K=5 (with MAI) K=1 (single user)

Fig. 3. The normalized CFO vector estimates as a function of tracking time in the training mode.

F. Summary of the Algorithm

The proposed recursive algorithm is summarized below and the order of estimates update is depicted in Fig. 1.

1) Initialize f_m|m(0) using (16) at the beginning of each mth symbol block for m > 1.

2) Iteratively update R(i+1)_|m , P(i+1)m , and f_m|m(i+1)using (22),

(23), and (15), respectively.

3) Impose the quadratic constraint and find ˆfm|m in (19).

4) Calculate ˆε_|musing (24) and reconstruct the CFO vector ˆ

fm using ˆε_|m.

5) Perform MAP symbol detection in (25).

IV. SIMULATIONRESULTS

In this section, we verify the effectiveness of the proposed recursive algorithm for the CFO vector estimation in the MC-CDMA systems with computer simulations. All numerical ex-periments were conducted under the environment of the BPSK transmission, the equal power assumption for all users, and a

0 50 100 150 200 250 300 350 400 450 500 10−1 100 101 102 Time (m) MSE No constraint, N_p=100 With constraint, N_p=100 No constraint, N_p=10 With constraint, Np=10

Fig. 4. The MSE of the CFO vector estimate as a function of time using the recursive EM algorithm in the training mode with and without the quadratic constraint.

frequency-selective fading channel model with order L = 5. Each user’s signature sequence was randomly generated, with a total K = 5 users in the system. The number of sub-carriers and the desired user’s CFO were set to N = 8 and ε₁= 0.1, respectively. All simulation tests adopted a common set of parameters, including the channel state information, CFO and thus the covariance matrix of MAI, all assuming time-invariant. The Monte Carlo simulation technique was employed to plot all simulation curves averaged over 250 realizations.

First, we demonstrate the effectiveness of the proposed algorithm in the training mode. Fig. 2 shows the tracking behavior of the first element of the estimated CFO vector. It is observed that the CFO vector is successfully tracked as the number of received signals grows. And, the behavior of the CFO estimate itself is shown in Fig. 3, where the estimate of the CFO itself is obtained via (24). In the presence of MAI, we see that the algorithm requires a larger number of samples for establishing the statistics of the interference-plus-noise correlation matrix until a reliable estimation of the CFO can be reached.

The effect of imposing the quadratic constraint in the

training mode with K = 5 and Eb/No = 0 dB is presented

in Fig. 4. For the two curves with a warm-up period of

Np = 10, the number of data samples is not sufficient

in the beginning stage to support a reliable estimate for R and the recursive algorithm with the quadratic constraint provides a more accurate estimate. This indicates that the proposed algorithm is more robust to imprecise estimation of the correlation matrix. A similar trend is observed in the case when the warm-up period is 100. However, as the number of samples grows, this discrepancy narrows down.

Next, we show the performance of the algorithm in the blind scenario. Fig. 5 demonstrates the mean-squared error (MSE) between the updated and true CFO vectors with respect to the symbol block time index m in the blind case. It is observed that there appears multiple high MSE time instances in the figure. This is due to possible erroneous symbol detection at

0 100 200 300 400 500 600 700 800 900 1000 10−3 10−2 10−1 100 101 Time (m)

MSE of the CFO vector (f

m| m (1 )) E b/No=20 dB E b/No=5 dB E_b/N_o=0 dB

Fig. 5. The MSE of the estimate ˆf_m|m as a function of time for different Eb/No in the blind scenario.

0 2 4 6 8 10 10−2 10−1 100 E b/No (dB) BEP Known CFO No CFO estimation With 3 iterations Steady state, 3 iters. With 1 iteration

Fig. 6. The evolution of the BEP performance of the MC-CDMA system in the blind scenario withε1= 0.1, N = 8, K = 5 and L = 5.

those time instances occurring in any of the 250 realizations. In other words, a wrong decision of a symbol results in a common sign change to all the components of the CFO vector estimate, which in turn affects the accuracy of the estimate of it. Nevertheless, this effect on the estimate of CFO vector doesn’t influence the estimate of the CFO itself, because the CFO estimate is obtained by comparing adjacent components of the updated CFO vector and the phase ambiguity can be resolved in the correlation operation.

The evolution of the BEP performance of the system using the proposed joint blind recursive estimation and detection algorithm is presented in Fig. 6. When no CFO estimation and compensation are employed, the system performance is too bad to be useful. Using the proposed algorithm, we can see from the figure that the BEP performance improves significantly even with only 1 iteration. After 3 iterations, the BEP performance can be recovered within a comparable level to that of an MC-CDMA system with perfect knowledge of CFO. 0 1 2 3 4 5 6 7 8 9 10 10−2 10−1 100 E b/No BEP

Tureli −Liu, block size 200 Li −Liu, block size 50 Li −Liu, block size 200 Li −Liu, block size 500 Proposed, 3 iter.

Fig. 7. The BEP performance of the MC-CDMA system in the blind scenario employing the proposed recursive algorithm, Li-Liu’s scheme [9] and Toreli-Liu’s scheme [10] withε₁= 0.1, N = 8, K = 5, and L = 5.

Next, we compare the BEP performance for MC-CDMA systems employing the proposed recursive algorithm, Tureli-Liu’s block-based scheme [10] and Li-Tureli-Liu’s block subspace scheme [9] in Fig. 7, where the number of samples M simulated at each point is at least 104 for each realization. We assume that the channel is known to the receiver in all schemes for fair comparison. Tureli-Liu’s block-based scheme fails to provide an accurate CFO estimate due to its unrealistic assumption on the orthogonality between channel impulse responses of different users. Consequently, we see an unacceptable level of BEP performance in Fig. 7. Li-Liu’s block subspace scheme with the minimum mean squared-error (MMSE) detection1has performance comparable to that of the proposed recursive algorithm when the number of samples B in each observation block is 500, as seen in Fig. 7. When processing delay and complexity are not of primary concern, increasing the block size in Li-Liu’s scheme will improve the system performance. However, there is an implicit assumption in this scheme. That is, the channel impulse responses of all users have to be static for the entire observation block. This requirement can be relaxed in our proposed algorithm. In addition, the computational cost required by the singular value decomposition (SVD) in Li-Liu’s scheme with N = 8

and block size B = 500 is about 4N2B +13B3≈ 1.625×109

floating point operations (flops) using the R-SVD algorithm [30, p.254]. In contrast, the eigenvalue decomposition in the

proposed algorithm using B samples demands about B(4N3+

13N3) ≈ 4.352 × 106 flops. It is clear that the size of

observation samples significantly impedes the efficiency of Li-Liu’s scheme in terms of detection delay and computational complexity.

On the other hand, the proposed recursive algorithm does need a sufficient number of samples to get statistics for ac-curate estimation for both the CFO vector and the correlation 1_{Due to the multiplication factors by the block size in the phase terms in}

[9, eq. (26)], differential MMSE detection is employed in the simulation to avoid magnifying the compensated residual phase.

0 2 4 6 8 10 10−2 10−1 100 E b/No BEP Proposed, M=50, 3 iter. Proposed, M=200, 3 iter. Proposed, M=103, 3 iter. Proposed, M=104, 3 iter.

Li −Liu, block size 200

Fig. 8. The BEP performance of the MC-CDMA system in the blind scenario employing the proposed recursive algorithm and Li-Liu’s scheme [9] with ε1= 0.1, N = 8, K = 5, and L = 5.

matrix. Fig. 8 shows that, when the total number of samples

M simulated in each realization is 200 (as opposed to 104

samples simulated in Fig. 7), the proposed recursive algorithm does not perform as well as Li-Liu’s scheme with a block size of 200 symbols. Nevertheless, we see a trend of performance improvement in the proposed algorithm as the number of the received samples increases. For practical applications such as the IEEE 802.16 standard, the time duration for one symbol

block is 100.8 μs. It takes 1.008s to collect 104 samples,

by which a reasonably acceptable BEP performance can be achieved as shown in Fig. 8.

V. CONCLUSION

A blind recursive CFO estimation and tracking technique in an MC-CDMA system was proposed in this research. We treated the problem from a new viewpoint so that techniques from linear estimation theory can be used. A recursive relation for the CFO vector was developed based on the EM algorithm with a quadratic constraint. The recursive update for the interference-plus-noise correlation matrix was also derived. It was shown by computer simulation that the BEP of an MC-CDMA system without CFO information can be restored to a level comparable to that with perfect CFO knowledge using the proposed CFO estimation and tracking algorithm.

APPENDIX

In the appendix, we provide the derivation of (14) and (19). First, we derive the result in (14). From (10) and the fact

that each received sample rl, 1 ≤ l ≤ m, is independent

with each other at different time instants, the cost function Qm(θm, θ(i)_m|m)can be decomposed as

Qm(θm, θ_m|m(i) ) = m−1_ l=1 Esl  log p (rl| sl; θl)ym; θ(i)_m|m  +Esm  log p(rm| sm; fm)ym; θ(i)_m|m  . (26)

The first expectation on the right hand side of (26) can be furthermore carried out as

Es_l  log p(rl| sl; fl)ym; θ(i)_m|m  = |A|  n=1 p(sl= ξn | ym; θ_m|m(i) ) log p(rl| sl= ξn; fl) = |A|  n=1 p(ym| sl= ξn; θ(i)_m|m)p(sl= ξn)

p(ym; θ_m|m(i) ) log p(rl| sl= ξn; fl) = |A|  n=1 p(ym−1| sl= ξn; θ(c)_m−1|m−1)p(rm; f_m|m(i) )p(sl= ξn) p(ym−1; θ_m−1|m−1(c) )p(rm; f_m|m(i) ) × log p(rl| sl= ξn; fl)  since l= m  =Esl  log p(rl| sl; fl)ym−1; θ(c)m−1|m−1 

For the second expectation in (26), the condition on ymcan be

reduced to rm. Therefore, we can represent the decomposition

in (26) as Qm(θm, θ_m|m(i) ) =Qm−1(θm−1, θ_m−1|m−1(c) ) +Esm  log p(rm| sm; fm)rm; θ(i)_m|m  . With the definition of Pm, we take the 2nd order derivative

of the above and obtain Pm=  ∂fm−1H ∂fm∗  Pm−1  ∂fm−1 ∂fmT  − Esm  GHmR−1Gm|rm; f_m|m(i)  ,

where the expectation can be further removed since sm has

been cancelled out in GmR−1Gm. Also, from the dynamic

evolution in (7), it follows that ∂fm−1/∂fmT = EH and

∂fm−1H /∂fm∗ = E, both of which are equivalent to mutually

conjugate complex scalars and can also be cleared. Thus, we

prove the relationship in (14). Note that Pm and Pm−1 do

not depend on any specific realizations of fm, because the

cost function Qm(θm, θ(i)_m|m) and Qm−1(θm−1, θ_m−1|m−1(c) )

are quadratic with respect to fm and fm−1, respectively.

Next, we show the derivation of (19). It can be easily shown that ˆ_f m|m= fm|m(c) −  Pm+ λmI −1 (27) ×Es_m[GHmR−1  rm− Gmf_m|m(c) rm; f_m|m(c) ] + λmf_m|m(c)  , which is the result from replacing Qm(θm, θ_m|m(i) ) in (13) by Jm(θm, λm). Gathering all factors associated with f_m|m(c)

yields  I −Pm+ λmI −1− GmR−1Gm+ λmI f_m|m(c) =  I −Pm+ λmI −1Pm− Pm−1+ λmI f_m|m(c) =  Pm+ λmI −1 Pm−1f_m|m(c) , (28)

where the relation in (14) is used in the derivation. Then, plugging (28) into (27), we obtain the result in (19).

REFERENCES

[1] N. Yee, J. P. Linnartz, and G. Fettweis, “Multi-carrier CDMA in indoor wireless radio networks,” in Proc. IEEE PIMRC 1993, pp. 109–113. [2] L. Wei and C. Schlegel, “Synchronization requirements for multi-user

OFDM on satellite mobile and two-path Rayleigh fading channels,” IEEE Trans. Commun., vol. 43, pp. 887–895, Feb./Mar./Apr. 1995. [3] P. H. Moose, “A technique for orthogonal frequency division

multiplex-ing frequency offset correction,” IEEE Trans. Commun., vol. 42, pp. 2908–2914, Oct. 1994.

[4] T. Pollet and M. Moeneclaey, “The effect of carrier frequency offset on the performance of band limited single carrier and OFDM signals,” in Proc. IEEE GLOBECOM’96, vol. 1, pp. 719–723.

[5] M. Luise and R. Reggiannini, “Carrier frequency acquisition and track-ing for OFDM systems,” IEEE Trans. Commun., vol. 44, pp. 1590–1598, Nov. 1996.

[6] J.-J. van de Beek, M. Sandell, and P.-O. Borjesson, “ML estimation of timing and frequency offset in OFDM Systems,” IEEE Trans. Signal Processing, vol. 45, pp. 1800–1805, July 1997.

[7] T. M. Schmidl and D. C. Cox, “Robust frequency and timing synchro-nization for OFDM,” IEEE Trans. Commun., vol. 45, pp. 1613–1621, Dec. 1997.

[8] H. Bolcskei, “Blind estimation of symbol timing and carrier frequency offset in wireless OFDM systems,” IEEE Trans. Commun., vol. 49, no. 6, pp. 988–999, 2001.

[9] K. Li and H. Liu, “Joint channel and carrier offset estimation in CDMA communications,” IEEE Trans. Signal Processing, vol. 47, no. 7, pp. 1811–1822, July 1999.

[10] U. Tureli, D. Kivanc, and H. Liu, “MC-CDMA uplink-blind carrier frequency offset estimation,” in Proc. 34th Asilomar Conference 2000, vol. 1, pp. 241–245.

[11] J.-H. Deng and T.-S. Lee, “An iterative maximum SINR receiver for multicarrier CDMA systems over a multipath fading channel with frequency offset,” IEEE Trans. Wireless Commun., vol. 2, no. 3, pp. 560–569, May 2003.

[12] Y. Xie and C. N. Georghiades, “Two EM-type channel estimation algorithms for OFDM with transmitter diversity,” IEEE Trans. Commun., vol. 51, pp. 106–115, Jan. 2003.

[13] W. Mo, Z. Wang, and A. Dogandˆzi´c, “Iterative channel estimation and decoding for coded MIMO system in unknown spatially correlated noise,” in Proc. Allerton Conf. 2003, pp. 228–232.

[14] E. Panayirci and C. N. Georghiades, “Joint ML timing and phase estimation in OFDM systems using the EM algorithm,” in Proc. IEEE ICASSP 2000, vol. 1, pp. 2949–2952.

[15] F.-T. Chien and C.-C. J. Kuo, “Joint symbol detection and channel estimation for MC-CDMA systems in the presence of carrier frequency offset,” in Proc. IEEE VTC’04 Fall.

[16] V. Krishnamurthy and J. B. Moore, “On-line estimation of hidden Markov model parameters based on the Kullback-Leibler information measure,” IEEE Trans. Signal Processing, vol. 4, no. 8, pp. 2557–2573, Aug. 1993.

[17] H. Zamiri-Jafarian and S. Pasupathy, “EM-based recursive estimation of channel parameters,” IEEE Trans. Commun., vol. 47, pp. 1297–1302, Sept. 1999.

[18] S. E. Bensley and B. Aazhang, “Maximum likelihood synchronization of a single-user for code-division multiple-access communication systems,” IEEE Trans. Commun., vol. 46, pp. 392–399, Mar. 1998.

[19] E. Ertin, U. Mitra, and S. Siwamogsatham, “Maximum-likelihood-based multipath channel estimation for code-division multiple-access systems,” IEEE Trans. Commun., vol. 49, pp. 290–302, Feb. 2001.

[20] A. A. D’Amico and M. Morelli, “Frequency estimation and timing ac-quisition in the uplink of a DS-CDMA system,” IEEE Trans. Commun., vol. 52, no. 10, pp. 1809–1819, Oct. 2004.

[21] J. G. Proakis, Digital Communications, Fourth Edition. New York: McGraw-Hill, 2000.

[22] X. Meng and D. B. Rubin, “Maximum likelihood estimation via the ECM algorithm: a general framework,” Biometrika, vol. 80, no. 2, pp. 263–278, June 1993.

[23] S. H. Wu, U. Mitra, and C.-C. J. Kuo, “Graph representation for joint channel estimation and symbol detection,” in Proc. IEEE GLOBE-COM’04, vol. 1, pp. 2329–2333.

[24] M. Morelli, “Timing and frequency synchronization for the uplink of an OFDMA systems,” IEEE Trans. Commun., vol. 52, no. 2, pp. 296–306, Feb. 2004.

[25] H. Cox, R. Zeskind, and M. Owen, “Robust adaptive beamforming,” IEEE Trans. Acoustic, Speech, Signal Processing, vol. 35, no. 10, pp. 1365–1376, Oct. 1987.

[26] D. M. Titterington, “Recursive parameter estimation using incomplete data,” J. Roy. Statist. Soc. B, vol. 46, no. 2, pp. 257–267, 1984. [27] W. Gandar, “Least squares with a quadratic constraint,” Numer. Math.,

vol. 36, pp. 291–307, 1981.

[28] L. B. Nelson and H. V. Poor, “Iterative multiuser receivers for CDMA channels: an EM based approach,” IEEE Trans. Commun., vol. 44, no. 12, pp. 1700–1710, Dec. 1996.

[29] G. K. Kaleh and R. Vallet, “Joint parameter estimation and symbol detection for linear or nonlinear unknown channels,” IEEE Trans. Commun., vol. 42, no. 7, pp. 2406–2413, July 1994.

[30] G. H. Golub and C. F. V. Loan, Matrix Computations, Third Edition The Johns Hopkins University Press, 1996.

Feng-Tsun Chien (S’02-M’05) received the B.S.

degree from the National Tsing Hua University, Hsinchu, in 1995, the M.S. degree from the National Taiwan University, Taipei, in 1997, and the Ph.D. degree from the University of Southern California, Los Angeles, in 2004, respectively, all in Electrical Engineering.

He joined the Department of Electronics En-gineering of the National Chiao Tung University, Hsinchu, in July 2005, as an Assistant Professor. His current research interests include signal process-ing aspects on communications, cross-layer considerations for OFDM and OFDMA systems, multicarrier CDMA, and MIMO-OFDM systems.

C.-C. Jay Kuo (S’83-M’86-SM’92-F’99) received

the B.S. degree from the National Taiwan University, Taipei, in 1980 and the M.S. and Ph.D. degrees from the Massachusetts Institute of Technology, Cam-bridge, in 1985 and 1987, respectively, all in Elec-trical Engineering. Dr. Kuo was Computational and Applied Mathematics (CAM) Research Assistant Professor in the Department of Mathematics at the University of California, Los Angeles, from October 1987 to December 1988. Since January 1989, he has been with the University of Southern California, where he is currently Professor of Electrical Engineering, Computer Science and Mathematics and Director of the Signal and Image Processing Institute. His research interests are in the areas of digital signal and image processing, multimedia compression, communication and networking technologies. Dr. Kuo has guided about 80 students to their Ph.D. degrees and supervised 15 postdoctoral research fellows. He is co-author of seven books and about 800 technical publications in international conferences and journals.

Dr. Kuo is a Fellow of IEEE and SPIE and a member of ACM. He is Editor-in-Chief for the Journal of Visual Communication and Image Representation, and Editor for the Journal of Information Science and Engineering and the EURASIP Journal of Applied Signal Processing. He was on the Editorial Board of the IEEE Signal Processing Magazine in 2003-2004. He served as Associate Editor for IEEE Transactions on Image Processing in 1995-98, IEEE Transactions on Circuits and Systems for Video Technology in 1995-1997 and IEEE Transactions on Speech and Audio Processing in 2001-2003. He received the National Science Foundation Young Investigator Award (NYI) and Presidential Faculty Fellow (PFF) Award in 1992 and 1993, respectively.