Semiblind linear multiuser interference cancellation: a maximum likelihood approach

Signal Processing 81 (2001) 2041–2057

www.elsevier.com/locate/sigpro

Semiblind linear multiuser interference cancellation: a maximum likelihood approach M'onica F. Bugallo ∗ , Joaqu'-n M'-guez, Luis Castedo Departamento de Electronica e Sistemas, Universidade da Coru˜na, Facultade de Informatica (lab. 1.2), Campus de Elvi˜na s=n, 15071 A Coru˜na, Spain Received 3 April 2000; received in revised form 9 April 2001

Abstract This paper introduces a linear interference suppression algorithm based on the maximum likelihood principle. The method is termed semiblind because it simultaneously exploits the transmission of training sequences and the statistics of the unknown transmitted symbols. We propose both iterative and adaptive implementations of the receiver using the expectation maximization algorithm and the matrix inversion lemma, respectively. Computer simulations show that the resulting multiuser receivers attain practically the same performance as the linear minimum mean square error multiuser detector. ? 2001 Elsevier Science B.V. All rights reserved. Keywords: Maximum likelihood; Multiuser detection; Semiblind receivers; Interference suppression; CDMA

1. Introduction Code division multiple access (CDMA) has been chosen as the basic multiple access scheme of third generation wideband mobile communication systems. The most outstanding advantages of CDMA over time division multiple access (TDMA) schemes employed in second generation systems are its higher spectral e=ciency and superior >exibility in the radio interface [1,3–5]. Nevertheless, important limitations still exist for the capacity of practical CDMA systems. More speciBcally, multiple access interference (MAI), due to the diCerent users simultaneously using the channel, and inter-symbol interference (ISI), due to multipath propagation, are major problems when high data rates must be attained. In order to mitigate both MAI and ISI, sophisticated signal processing capabilities have to be incorporated into wireless receivers. Due to the prohibitive computational and informative requirements of the optimum CDMA multiuser receiver [24] alternative approaches based on linear Bltering have been investigated. The decorrelating receiver [24] requires a perfect knowledge of the desired user signature waveform. This signature is not likely to be available because it depends on the channel distortion, which is usually unknown. To avoid this drawback, minimum mean square error (MMSE) detectors [16,24] resort to the transmission ∗

Corresponding author. Tel.: +34-981-167000; fax: +34-981-167160. E-mail addresses: [email protected] (M.F. Bugallo), [email protected] (J. M'-guez), [email protected] (L. Castedo). 0165-1684/01/$ - see front matter ? 2001 Elsevier Science B.V. All rights reserved. PII: S 0 1 6 5 - 1 6 8 4 ( 0 1 ) 0 0 0 8 5 - 8

2042

M.F. Bugallo et al. / Signal Processing 81 (2001) 2041–2057

Nomenclature bi (n) b(n) bt D g(n) x(n) w wo y fu (·) fu; (·) fu|v; (·) fu|v; (·) Eu [ · ] Eu;J [ · ] Eu|v;J [ · ] Eu|v; z;J [ · ] L(wo ) L(wo ) | bt T (·)

JT

ith user transmitted symbol transmitted symbols vector training sequence received code matrix noise vector observation vector Blter coe=cient vector optimum Blter coe=cient vector soft estimate probability density function (p.d.f.) of the random vector u p.d.f. of the random vector u that depends on the parameter vector  p.d.f. of the random vector u conditioned to the vector of values v. The conditioned random vector u | v depends on the parameter vector  p.d.f. of the random vector u conditioned to the vectors of values v and u. The conditioned random vector u | v; z depends on the parameter vector  statistical expectation with respect to (w.r.t.) the random vector u statistical expectation w.r.t. the random vector u depending on the parameter vector  statistical expectation w.r.t. the random vector u conditioned to the vector values v. The conditioned random vector u | v depends on the parameter vector  statistical expectation w.r.t. the random vector u conditioned to the vectors values v and z. The conditioned random vector u | v; z depends on the parameter vector  log-likelihood of wo w.r.t. the observed data y log-likelihood of wo w.r.t. the observed data y conditioned to the training sequence bt linear invertible transformation of (·) Jacobian of the transformation T (·)

of training sequences. Such sequences, however, are not always available in many practical applications. Therefore, alternative blind implementations have been studied [6,9,15,18,19,23,25]. Most blind schemes based on the linearly constrained minimum variance (LCMV) criterion [9,15] require a precise knowledge of the desired user signature which is not likely to be available in practice. This drawback is overcome in [23], which only employs the transmitted spreading code. Nevertheless, all LCMV multiuser receivers suCer from a slow convergence rate and steady-state performance degradation [7]. Constant modulus (CM) receivers [18,19] are more robust to code mismatches but do not improve the convergence speed. Subspace techniques with a somewhat faster convergence rate have also been suggested [2,6,14,25] but they are limited by their high computational complexity and their poor performance in the low signal to noise ratio (SNR) region. In third generation communication systems information symbols are transmitted in a burst way. Data frames are constructed according to communication protocols that typically add a priori known synchronization bits and training sequences to the unknown information symbols for many diCerent purposes. In this paper, we introduce a novel semiblind technique for linear multiuser interference cancellation that exploits both the a priori known transmitted symbols and the statistical features of the unknown symbols transmitted by the desired user. No knowledge of the user codes or the channel response is assumed. Our approach consists in applying the maximum likelihood (ML) principle in order to estimate the coe=cients of the linear receiver that suppresses both MAI and ISI. We show how this estimation

M.F. Bugallo et al. / Signal Processing 81 (2001) 2041–2057

2043

Fig. 1. Baseband discrete-time equivalent model of a DS CDMA system with time dispersive channels.

problem can be solved using the expectation maximization (EM) method. We also derive a computationally e=cient adaptive implementation of the EM algorithm that attains fast convergence. It is interesting to note that since the proposed receiver jointly exploits the knowledge of the training sequences and the statistical properties of the unknown symbols, it achieves a signiBcant performance gain with respect to (w.r.t.) conventional linear detectors that rely exclusively on the training sequence and blind linear receivers that ignore the known symbols in the data frames. The remainder of this paper is organized as follows. Section 2 introduces the baseband discrete-time equivalent signal model of an asynchronous CDMA communication system with time-dispersive channels. In Section 2 we introduce the ML linear multiuser receiver. Section 4 describes the optimization algorithms used to compute the Blter coe=cients. Section 5 presents simulation results and Section 6 is devoted to the conclusions. 2. Signal model Fig. 1 depicts the signal >ow diagram of a discrete-time baseband equivalent model of an asynchronous time-dispersive direct sequence (DS) CDMA communication system. The ith user nth transmitted symbol, bi (n), is multiplied by a unique binary-valued spreading sequence with L chips per symbol, ci (k); k = 0; : : : ; L − 1. The resulting transmitted signal passes through a linear time-dispersive channel, hi (k); k = 0; : : : ; P − 1, that accounts for the continuous-time channel response, the relative time delays of the diCerent users and the transmitter and receiver front-end Blters. The received signal is the superposition of the signals from the N users plus the additive white Gaussian noise (AWGN) sequence g(k), i.e., x(n; k) =

N


di (k)bi (n) + g(k);

k = 0; : : : ; L + P − 2:

(1)

i=1

 The sequence di (k) = ci (k) ∗ hi (k) = P−1 p=0 hi (p)ci (k − p), where ∗ denotes convolution, has length L + P − 1 and will be termed the ith user signature waveform. Due to the channel time scattering eCect, when a symbol stream is transmitted the symbol bi (n) interferes with bi (n − 1); : : : ; bi (n − m + 1), where m = (L + P − 1)=L is the channel memory size (· denotes upper integer part). Therefore, the overall received sequence during the nth symbol period is

x(n; k) = x(nL + k) =

N m−1

i=1 r=0

di (r; k)bi (n − r) + g(k);

k = 0; : : : ; L − 1;

(2)

2044

M.F. Bugallo et al. / Signal Processing 81 (2001) 2041–2057

Fig. 2. Linear multiuser receiver.

where di (r; k) = di (rL + k) is the rth segment of length L in di (k). Bringing together all the observations that involve the symbol of interest, bi (n), we obtain an expression for the overall received sequence using the vector notation x(n) = Db(n) + g(n);

(3) T

where x(n) = [x(n; 0); : : : ; x(n; L − 1); : : : ; x(n + m − 1; 0); : : : ; x(n + m − 1; L − 1)] is the Lm × 1 observation vector,   D m−1 D m−2 · · · D0 0 ··· 0    0 D m−1 D m−2 · · · D0 ··· 0     D=  ..  . . .. . . . . . .  . . . ··· . .  .   · · · D m−1 D m−2 · · · D 0 0 0 is the Lm × N (2m − 1) signature waveform matrix, composed of the submatrices D r = [d1 (r); : : : ; dN (r)]T , with di (r) = [di (r; 0); : : : ; di (r; L − 1)]T . The transmitted symbol vector is b(n) = [b1 (n − m + 1); : : : ; bN (n − m+1); : : : ; b1 (n); : : : ; bN (n); : : : ; b1 (n+m − 1); : : : ; bN (n+m − 1)]T and g(n) = [g(n; 0); : : : ; g(n+m − 1; L − 1)]T is a vector of independent and identically distributed (i.i.d.) Gaussian random variables with zero mean and variance g2 . The linear multiuser receiver consists of a Bnite impulse response (FIR) Blter, w = [w1 ; : : : ; wLm ]T , followed by a threshold detector as shown in Fig. 2. The soft estimate corresponding to the nth symbol period can be written as y(n) = wH x(n);

(4)

where the superindex H denotes Hermitian transposition. 3. Selection of the receiver coe"cients In this section we derive a novel statistical criterion to select the receiver coe=cients in order to obtain MAI and ISI free estimates of the desired user symbols. The criterion is based on the fact that, when the MAI and the ISI are totally suppressed, the soft symbol estimate, y(n), consists of just two components:

M.F. Bugallo et al. / Signal Processing 81 (2001) 2041–2057

2045

the desired user symbol, b1 (n), and an additive Gaussian noise term, gf (n). Indeed, let wo denote the optimum value of the Blter coe=cients that completely eliminates the MAI and the ISI. Then, we can write y(n) = woH x(n) = b1 (n) + gf (n);

(5)

where b1 (n) is the desired user symbol and gf (n) is a complex Gaussian random variable with zero mean and variance f2 = g2 woH wo . In this section, the Bltered noise variance, f2 , will be considered a constant for the sake of simplicity. In Section 4, however, we will also propose an easy-to-implement updating rule for f2 . It is straightforward to show that the probability density function (p.d.f.) of y(n) = b1 (n) + gf (n) is given by fy;wo (y(n)) = fb+gf (y(n)) =

2 2 1 Eb [e−|y(n)−b| =f ]; f2

(6)

where Eb [ · ] denotes statistical expectation w.r.t. the desired user transmitted symbols. In digital communications, the transmitted symbols are usually modelled as discrete i.i.d. random variables with known p.d.f. and Bnite alphabet. As a consequence, the soft estimates obtained with the optimum Blter wo can also be considered as i.i.d. random variables and, when a block of K observation vectors is available, the joint p.d.f. of the resulting data frame y = [y(0); : : : ; y(K − 1)]T is K−1



1 fy;wo (y(n)) = fy;wo (y) = f2 n=0

K K−1

2

2

Eb [e−|y(n)−b| =f ]

(7)

n=0

that can be analytically calculated because a statistical expectation w.r.t. a Bnite alphabet discrete random variable reduces to a simple summation. Note that the p.d.f. of y given by Eq. (7) depends on the unknown parameter wo . Thus, the ML estimate of wo corresponds to

wˆ = arg max L(wo ) = log wo

K−1

−|y(n)−b|2 =f2

Eb [e

n=0

]=

K−1


−|y(n)−b|2 =f2

log Eb [e

] ;

(8)

n=0

where L(wo ) is the log-likelihood of wo w.r.t. the observed data y = [y(0); : : : ; y(K − 1)]. Note that the expression of the log-likelihood, L(wo ), involves averaging over the p.d.f. of the unknown transmitted symbols. Therefore, the proposed approach can be classiBed within the family of stochastic ML methods [22]. Nevertheless, such classiBcation only applies to the computation of the detector coe=cients, since the receiver itself is a linear processor. As a consequence, it is suboptimum in terms of bit error rate (BER) but it presents a very reduced computational complexity when compared to true ML detectors [24]. The log-likelihood L(wo ) is a nonquadratic function with several maxima. In particular, the solutions to problem (8) guarantee that the soft estimates, y(n), have a p.d.f. close to fb+gf (·) but this is not enough to ensure that the desired user is extracted. Since in CDMA all users transmit symbols with the same modulation format, the p.d.f. of the ith interference at the receiver is fb+gf (·) which does not diCer from the desired user p.d.f.. Therefore, solving the optimization problem (8) may lead to the capture of the ith interference. In order to alleviate this limitation we propose to exploit the transmission of a short training sequence of M ¡ K symbols. Such sequences are usually employed for diCerent purposes in currently

2046

M.F. Bugallo et al. / Signal Processing 81 (2001) 2041–2057

standardized mobile communication systems. Indeed, let us assume that the Brst M symbols (i.e., bt = [b1 (0); : : : ; b1 (M − 1)]T ) are known a priori by the receiver. Conditioning the expectations in (8) w.r.t. the known symbols, bt , we arrive at a semiblind receiver where the Blter coe=cients are computed as the solution to

wˆ = arg max L(wo )|bt = wo

K−1


log Eb1 (n)|bt [e

−|y(n)−b1 (n)|2 =f2

]

n=0

M −1

K−1

2 −|y(n)−b|2 =f2 = arg min |y(n) − b1 (n)| − log Eb [e ] : wo

n=0

(9)

n=M

The computer simulations in Section 5 show that relatively short training sequences are enough to avoid the capture problem. This is because the Brst term in (9) is a purely quadratic form with a single maximum that corresponds to the rough extraction of the desired user. Note that, since wˆ is still computed according to the ML principle, all the available statistical information is employed to obtain the Blter coe=cients, and hence the proposed semiblind ML linear receiver outperforms conventional MMSE multiuser detectors [16,24] that only exploit the training sequence bt . It is interesting to remark that the semiblind formulation of the linear receiver in (9) is just a particular case of the more general ML approach that led to problem (8) (which is a totally blind formulation). On the contrary, most semiblind receivers proposed so far in the literature [11–13] result from an ad hoc approach that consists of combining two diCerent criteria, one for the training sequence and another one for the unknown symbols.

3.1. Existence of wo According to Eq. (5), wo is the decorrelating solution for user 1, i.e., it is the Brst row in the left pseudoinverse of the distorted code matrix, D, that we denote as D† . It is apparent that matrix D† only exists when N; L and m verify the inequality N (2m − 1) 6 Lm:

(10)

Nevertheless, even when (10) does not hold, the existence of the Wiener solution [8] wW = (E[x(n)xH (n)])−1 E[x(n)b∗1 (n)]

(11)

is guaranteed. The soft estimates resulting from Wiener Bltering can be written as H y(n) = wW x(n) = b1 (n) + gf (n) + r(n);

(12)

where r(n) is the residual MAI and ISI. Since the interference power, E[r H (n)r(n)], is very low (specially for medium and high values of the SNR) the associated p.d.f. fy;wW (·) is very close to Eq. (6). Indeed, it can be shown (see Appendix A) that the solution to problem (9) has the same structure as the Wiener Blter.

M.F. Bugallo et al. / Signal Processing 81 (2001) 2041–2057

2047

4. Implementation 4.1. Iterative-block implementation Unfortunately, it is not possible to Bnd a closed form solution to problem (9) and, therefore, some optimization algorithm must be used to obtain the Blter coe=cients. In this work, we propose to compute wo using the EM algorithm [20] that provides an iterative procedure to perform ML estimation when direct maximization of the likelihood is not feasible. The EM approach postulates the existence of some missing (unobserved) data that, if known, would aid in the estimation problem. The algorithm consists of a two-step iteration: use the incomplete (observed) data and the current parameter estimates to compute su=cient statistics of the complete data (E-step) and re-estimate the parameters using the computed complete data su=cient statistics (M-step). The sequence of estimates thus obtained exhibits the desirable property of being monotonically nondecreasing in likelihood. In our problem, the incomplete-data set is given by the soft estimates, y(n); n = 0; : : : ; K − 1, whereas the complete-data set is given by the extended vectors ye (n) = [y(n) b1 (n)]T ; n = 0; : : : ; K − 1. Let us build the complete-data block Ye = [ye (0); : : : ; ye (K − 1)], with joint p.d.f. fYe ;wo (·). It is easy to decompose fYe ;wo (·) as fYe ;wo (Ye ) = fYe |y;wo (Ye ) · fy;wo (y):

(13)

Taking logarithms and conditional expectations on both sides of (13) yields EYe |y;bt ;w(i) ˆ [log fy;wo (y)] = EYe |y;bt ;w(i) ˆ [log (fYe ;wo (Ye )) − log (fYe |y;wo (Ye ))]

(14)

and it is straightforward to realize that EYe |y;bt ;w(i) ˆ [log fy;wo (y)] = log fy;wo (y):

(15)

Let us now introduce the following notation: ˆ = EYe |y;bt ;w(i) U (wo ; w(i)) ˆ [log (fYe ;wo (Ye ))];

(16)

ˆ V (wo ; w(i)) = EYe |y;bt ;w(i) ˆ [log (fYe |y;wo (Ye ))]:

(17)

Substituting (15) – (17) into (14) we arrive at the relationship ˆ ˆ − V (wo ; w(i)); log fy;wo (y) = U (wo ; w(i))

(18)

which provides a tractable form of the log-likelihood function log fy;wo (·). Indeed, an application of Jensen’s inequality allows us to show that [17] ˆ ˆ 6 V (ˆ w(i); w(i)) V (wo ; w(i))

(19)

for any value of wo and, as a consequence, the sequence of estimates ˆ ˆ + 1) = arg max U (wo ; w(i)) w(i wo

(20)

is nondecreasing in likelihood since ˆ ˆ log fy;w(i+1) (y) = U (ˆ w(i + 1); w(i)) − V (ˆ w(i + 1); w(i)) ˆ ˆ ˆ ¿ U (ˆ w(i); w(i)) − V (ˆ w(i); w(i)) = log fy;w(i) ˆ (y): Substituting the joint p.d.f. K−1 1 2 2 e−|y(n)−b1 (n)| =f · fb (b1 (n)) fYe ;wo (Ye ) = 2 f n=0

(21)

(22)

2048

M.F. Bugallo et al. / Signal Processing 81 (2001) 2041–2057

(see Appendix B) into (20) leads to the single iterative rule ˆ ˆ + 1) = arg max{U (wo ; w(i)) } w(i wo

M −1

K−1

2 2 = arg min |y(n) − b1 (n)| + Eye (n)|y(n);w(i) ˆ [|y(n) − b1 (n)| |] ; wo

(23)

n=M

n=0

where we have neglected all terms that are constant w.r.t. wo . At this point, we have cast problem (9), which does not have a closed form solution, into a sequence of quadratic problems that can be analytically solved. It can be demonstrated that the solution to (23) is given by ˆ + 1) = w(i

K−1


−1 M −1  K−1

∗ x(n)xH (n) · x(n)b∗1 (n) + Eb1 (n)|y(n);w(i) ˆ [b1 (n)]x(n) ;

n=0

n=0

(24)

n=M

where we have also used the fact that the only random part in ye (n) | y(n) is b1 (n). Notice that, for M = K, the previous algorithm reduces to the least squares (LS) Blter. Appendix C shows the derivation of the conditioned expectation Eb1 (n)|y(n);w(i) ˆ [ · ] (n ¿ M ) using the Bayes theorem. Finally, note that the Bltered noise variance parameter, f2 , is required in order to compute Eb1 (n)|y(n);w(i) ˆ [ · ], n ¿ M , in Eq. (24). Since f2 = woH wo g2 is not known a priori, it must be estimated. A very simple estimation method consists of iteratively updating f2 using the estimates of wo obtained from (24), i.e., ˆ ˆ2f (i): ˆ2f (i + 1) =wˆ H (i)w(i)

(25)

ˆ + 1). According to our computer Thus, the updated value of ˆ2f (i + 1) can be used to compute w(i 2 2 simulations, an adequate initialization of (25) is ˆf = ˆg , where ˆ2g is a rough estimate of the AWGN variance. 4.2. Adaptive implementation ˆ ˆ − 1) according to (24), a system of Lm linear equations with Lm In order to compute w(i) from w(i unknowns has to be solved. In practice, however, it may be desirable to avoid this operation, specially ˆ when the spreading code length, L, is very large. We would also like to compute the estimate w(i) recursively in order to minimize the number of operations per iteration. Both objectives can be realized by using the matrix inversion lemma [21]. Applying this result in (24) and denoting P(0) = ILm ;  i −1
H P(i) = x(n)x (n) = P(i − 1) − k(i)xH (i)P(i − 1);

i¿0

(26)

n=0

and k(i) =

P(i − 1)x(i) ; 1 + xH (i)P(i − 1)x(i)

(27)

M.F. Bugallo et al. / Signal Processing 81 (2001) 2041–2057

2049

Table 1 Adaptive EM receiver Initialize the algorithm by setting P(0) = ILm f2 (0) ≈ g2 For each available observation vector i = 1; 2; : : : ; compute P(i−1)x(i) 1+xH (i)P(i−1)x(i)

k(i) =

ˆ − 1); (i ¿ M ) f2 (i) = f2 (i − 1)wˆ H (i − 1)w(i

d(i) =

b∗ (i);

06i6M − 1 ∗

Eb(i)|y(i);w(i) ˆ b (i); M 6 i 6 K − 1

ˆ =w(i ˆ − 1) + k(i)[d(i) − xH (i)w(i − 1)] w(i) P(i) = P(i − 1) − k(i)xH (i)P(i − 1)

where ILm is the Lm × Lm identity matrix, we get the desired equation for updating the Blter coe=cients, ˆ =w(i ˆ − 1) + k(i)[d(i) − xH (i)w(i − 1)]; w(i)

(28)

where

d(i) =

b∗ (i);

0 6 i 6 M − 1; ∗

Eb(i)|y(i);w(i) ˆ b (i);

M 6 i 6 K − 1:

(29)

We would like to remark that the adaptive EM algorithm (28) is a generalization of the well-known recursive least squares (RLS) method [8]. Actually, (28) reduces to the RLS rule when M = K. Table 1 summarizes the proposed adaptive algorithm, including the update of the Bltered noise variance estimate, ˆ2f .

5. Computer simulations Finally, we present computer simulations that illustrate the validity of our approach. We have considered an asynchronous DS CDMA communication system with N = 4 users transmitting QPSK symbols and length L = 8 random binary spreading codes,   −1 1 1 1 −1 −1 −1 −1    −1 −1 −1 1 1 −1 −1 −1    T C = :  1 1 1 1 1 1 −1 −1    1 1 1 1 1 −1 1 1

2050

M.F. Bugallo et al. / Signal Processing 81 (2001) 2041–2057

Fig. 3. MSE for several values of the SNR in a dispersive asynchronous DS CDMA system with N = 4 users, length L = 8 random binary spreading codes and length P = 16 discrete channels. Left: K = 300, M = 15. Right: K = 600; M = 30.

The user signals propagate through a randomly chosen unknown channel with P = 16 complex coe=cients per user. The resulting multiuser channel matrix is   0:23 + 0:26i 0:01 − 0:24i 0:25 − 0:36i −0:00 − 0:17i    0:63 + 0:09i 0:53 − 0:26i −0:01 − 0:13i 0:44 − 0:06i       −0:33 + 0:15i 0:57 − 0:13i −0:30 − 0:94i 1:14 − 0:39i       −0:64 − 0:08i −0:00 + 0:54i −0:06 − 0:47i −0:12 + 0:03i       0:17 − 0:02i −0:50 + 0:08i 0:41 + 0:94i 0:23 − 0:60i       −0:34 + 0:38i 0:28 + 0:75i −0:09 + 0:10i 0:17 − 0:57i         −0:18 + 0:61i −0:21 + 0:22i −0:18 + 0:74i −0:64 − 0:21i    h1 (0) · · · hN (0)      0:71 + 0:36i −0:43 + 0:54i 1:08 + 0:25i −0:31 + 0:59i      .. . . .. .. ;  = .    0:46 − 1:10i 0:39 + 0:07i −0:93 − 0:18i 0:24 − 0:01i      h1 (P − 1) · · · hN (P − 1)  0:12 + 0:45i −0:27 − 0:09i 0:05 − 0:16i 0:12 − 0:55i       −0:37 + 0:30i 0:26 − 0:47i −0:48 + 0:23i 0:70 − 0:38i       −0:19 − 0:80i 0:22 + 0:62i 0:58 − 0:86i 0:29 + 0:68i       −0:90 − 1:03i −0:01 − 0:63i −0:12 + 0:15i −0:75 − 1:36i       −1:06 − 0:47i 0:04 + 0:20i 0:42 + 0:04i 0:38 − 0:19i       0:63 + 0:00i 1:35 + 0:17i 0:68 + 0:51i −0:43 + 0:64i    0:72 − 0:13i

0:34 + 0:48i

0:39 + 1:10i

0:62 − 0:62i

√ where i = −1. Fig. 3 (left) plots the mean square error (MSE) attained by the proposed iterative EM receiver (labelled EM) for several values of the input SNR, deBned as

SNR = 10 log10

E[|b1 |2 ]d1H d1 (dB) g2

(30)

M.F. Bugallo et al. / Signal Processing 81 (2001) 2041–2057

2051

Fig. 4. MSE for several values of the SNR in a dispersive asynchronous DS CDMA system with N = 4 users, length L = 8 random binary spreading codes, length P = 16 discrete channels and SIR j = − 5 dB ∀j. Left: K = 300; M = 15. Right: K = 600; M = 30.

when the number of observation vectors available to estimate the receiver coe=cients is K = 300, including a training sequence of length M = 15. It is apparent that the proposed iterative-block algorithm performs close to the theoretical MMSE receiver (labelled MMSE) constructed with perfect knowledge of all user codes and channel coe=cients. In this Bgure we have also plotted the MSE achieved by the conventional LS receiver (labelled LS) implemented with the RLS algorithm (as deBned in [8]) using the training sequence, bt , alone, as it is done in conventional supervised multiuser receivers. It can be seen that this detector performs considerably worse than the theoretical MMSE because the training sequence is too short to obtain an adequate estimation of the Blter coe=cients. Fig. 3 (right) illustrates the performance of the adaptive EM receiver (labelled AEM) in terms of MSE using K = 600 observation vectors and a training sequence of M = 30 symbols. The reduction in computational complexity of the AEM algorithm w.r.t. its iterative counterpart is achieved at the expense of an increase in the number of observation vectors required to attain an adequate performance. The AEM receiver is compared, in terms of MSE, with the conventional RLS receiver that exploits only the training sequence, bt . It is clear that, for practical values of the SNR, the proposed adaptive detector performs very close to the theoretical MMSE and achieves a clear improvement w.r.t. the conventional RLS detector. Fig. 4 (left) depicts the MSE of the theoretical MMSE, LS and the iterative EM linear receivers in a near-far environment (with K = 300 and M = 15). Let us deBne the signal to interference ratio (SIR) of the desired user w.r.t. the jth interference as dH d1 E[|b1 |2 ] (31) SIR j = 10 log10 1H dj dj E[|bj |2 ] and let us choose the value of E[|bj |2 ] so that SIR j = − 5 dB ∀j. Simulation results show that all detectors suCer a moderate degradation in performance but the proposed iterative receiver still approaches the theoretical limit for moderate and high SNR values. Similar results are obtained with the AEM receiver, as shown in Fig. 4 (right) with parameters K = 600 and M = 30. Another important measure of the receiver performance is the MSE achieved for diCerent system loads. Fig. 5 (left) plots the MSE for several values of the number of users, N , with K = 300; M = 15; SNR = 12 dB and SIR j = 0 dB ∀j. The resulting curve shows that the performance degradation of the iterative EM receiver with increasing system load is the same one suCered by the theoretical MMSE receiver, whereas the conventional LS receiver performance is considerably worse. The MSE of the iterative receiver becomes considerably higher than the theoretical MMSE only when the system is heavily loaded. Fig. 5 (right) shows the corresponding results for the AEM receiver (with K = 600; M = 30).

2052

M.F. Bugallo et al. / Signal Processing 81 (2001) 2041–2057

Fig. 5. MSE vs. number of users in a dispersive asynchronous DS CDMA system with length L = 8 random binary spreading codes, length P = 16 discrete channels, SNR = 12 dB and SIR j = 0 dB ∀j. Left: K = 300, M = 15. Right: K = 600; M = 30.

Fig. 6. SER for several values of the SNR in a dispersive asynchronous DS CDMA system with N = 4 users, length L = 8 random binary spreading codes, length P = 16 discrete channels and SIR j = 0 dB ∀j. Left: K = 300; M = 15. Right: K = 600; M = 30.

Fig. 6 (left) depicts the symbol error rate (SER) attained by the proposed iterative EM algorithm for several values of the input SNR with system parameters N = 4; K = 300; M = 15 and SIR j = 0 dB ∀j. The resulting curves show that the proposed receiver performs close to the theoretical MMSE detector, whereas the performance of the LS receiver is considerably worse. The same conclusion is drawn from the simulation results shown in Fig. 6 (right) with K = 600 and M = 30. Fig. 7 illustrates the fast convergence speed of the proposed algorithms. With the same simulation parameters as in the previous experiments and SNR = 12 dB, we have plotted the MSE at the receiver output as a function of the number of iterations (left) and the number of observation vectors used in the estimation of the Blter coe=cients (right) for the iterative and adaptive algorithms, respectively. It is apparent that the iterative receiver converges in less than 20 iterations and the recursive detector attains the theoretical limit with a moderate number of observation vectors. This may be an important advantage when time or computational load constraints have to be fullBlled. Finally, we have carried out computer simulations in order to evaluate the relative merit of the proposed EM-based receivers when compared with the other blind and semiblind multiuser detectors. In particular, we have considered the following receivers: • A block-iterative decision directed LS (DD LS) receiver. • The semiblind CM (SB CMA) receiver proposed in [11] implemented with an oC-line gradient algorithm.

M.F. Bugallo et al. / Signal Processing 81 (2001) 2041–2057

2053

Fig. 7. Convergence speed of the proposed algorithms in a dispersive asynchronous DS CDMA system with N = 4 users, length L = 8 random binary spreading codes, length P = 16 discrete channels, SNR = 12 dB and SIR j = 0 dB: ∀j Left: MSE vs. no. of iterations (K = 300 and M = 15). Right: MSE vs. no. of observation vectors (K = 600 and M = 30).

Fig. 8. Comparison with other methods: MSE for several values of the SNR in a dispersive asynchronous DS CDMA system with N = 4 users, length L = 8 random binary spreading codes, length P = 16 discrete channels and SIR j = 0 dB: ∀j. Left: K = 300; M = 15. Right: K = 600; M = 30.

• A closed-form solution implementation of the blind LCMV (B LCMV) multiuser detector [9] which requires perfect knowledge of the desired user channel coe=cients, h1 (0); : : : ; h1 (P − 1). • The blind Subspace Algorithm (BSA) proposed by Wang and Poor [26] with a smoothing factor

of 8.

• A decision directed implementation (hard decisions are employed) of the RLS (DD RLS)

algorithm [8]. The DDLS, SB CMA, B LCMV and BSA approaches are compared with the iterative EM receiver in Fig. 8 (left). It is apparent that the proposed method practically matches the theoretical MMSE whereas the remaining receivers present a signiBcantly poorer performance. The simulation parameters are K = 300 and M = 15. The AEM algorithm is compared with the DD RLS approach in Fig. 8 (right). The simulation parameters are K = 600 and M = 30. In this case, there is only a slight performance advantage of the AEM algorithm w.r.t. the conventional DD RLS approach. This result is actually a rather intuitive one, since the AEM algorithm can also be seen as a decision directed approach where symbol estimates are obtained according to the nonlinear mean-squared method [21] (see Eq. (29)).

2054

M.F. Bugallo et al. / Signal Processing 81 (2001) 2041–2057

6. Conclusions We have introduced a new semiblind approach to linear interference cancellation in DS CDMA. The proposed receiver relies on the ML principle to estimate the coe=cients of a linear FIR Blter that suppresses both MAI and ISI. The method exploits both the transmission of short training sequences and the knowledge of the transmitted symbols and noise p.d.f. Two possible implementations are suggested: an iterative-block one based on the EM algorithm and an adaptive one that uses the matrix inversion lemma in order to reduce the computational complexity. Finally, computer simulations show that, using relatively short training sequences and small blocks of observations, the resulting semiblind receivers attain the same performance as the theoretical linear MMSE multiuser detector. Acknowledgements This work has been supported by FEDER funds (grant number 1FD97-0082) and Xunta de Galicia (PGIDT00PXI10504PR). We would also like to thank the anonymous reviewers for their comments (specially reviewer A) that have helped to signiBcantly improve the Bnal version of the paper. Appendix A. Relationship between the proposed semiblind receiver and the linear Wiener solution Let us consider the linear MMSE multiuser receiver, i.e., wW = arg min{Ex(n) [|wH x(n) − b1 (n)|2 ]}:

(A.1)

w

It is straightforward to show that the solution to the above problem is wW = Rx−1 p;

(A.2)

where Rx = Ex(n) [x(n)xH (n)] and p = Ex(n) [x(n)b∗1 (n)]. In this appendix we will show the close relationship between the proposed receiver and the Wiener solution given by (A.2). With this aim, let us characterize the local maxima of the log-likelihood function L(w) | bt =

K−1


2

2

log Eb|bt [e−|y(n)−b| =f ]

(A.3)

n=0

w.r.t. the Blter coe=cients w. The stationary points of L(w) are found on calculating the gradient ∇w L and equalling it to zero: ∇w L =

K−1


x(n)Eb1 (n)|y(n);w [y∗ (n) − b∗1 (n)] = 0:

(A.4)

n=0

Taking into account that y(n) = wH x(n), the previous equation can be elaborated to yield K−1
n=0

x(n)xH (n)w =

K−1


x(n)d(n);

(A.5)

n=0

where



d(n) = Eb1 (n)|y(n);bt ;w [b∗1 (n)] =

b∗ (n);

0 6 n 6 M − 1;

Eb(n)|y(n);w b∗ (n);

M 6n6K − 1

(A.6)

M.F. Bugallo et al. / Signal Processing 81 (2001) 2041–2057

is the nonlinear mean-squared estimate of b∗1 (n) [21]. Solving for w we arrive at ˆ −1 ˆ w =R x p;

2055

(A.7)

where ˆ −1 R x =

K−1


x(n)xH (n)

(A.8)

is the empirical autocorrelation matrix and K−1
pˆ = x(n)d(n)

(A.9)

n=0

n=0

is an empirical cross-correlation vector where the transmitted symbols are substituted by their meansquared estimates (when n ¿ M ). It is apparent that (A.7) converges to the Wiener solution (A.2) when the block size, K, is large enough. Note that Eq. (A.7) is not a useful result from a practical point of view since it is not an explicit expression of w because p depends on w. This coe=cient vector must be known in order to compute the mean-squared estimates of the symbols. The EM algorithm proposed in this paper is actually an iterative method to numerically approximate (A.7). Appendix B. Derivation of fye ;wo (·) When w is adequately chosen, i.e, w = wo , y(n) = woH x(n) = b1 (n) + gf (n) the extended vector ye (n) = [y(n) b1 (n)]T is easily obtained through a linear invertible transformation of the extended symbol vector be (n) = [gf (n) b1 (n)]T as  b1 (n) + gf (n) : (B.1) ye (n) = T (be (n)) = b1 (n) It is well known that the p.d.f.’s of ye (n) and be (n) are related by the following expression [10]: fb (T −1 (ye (n))) fye ;wo (ye (n)) = e ; (B.2) |JT | where JT is the Jacobian of the transformation and | · | denotes absolute value. It is straightforward to show that  @y(n) @b (n)  1    @gf   1 0  @g f  (B.3) JT = det  =1  @y(n) @b (n)  =  1 1 1 @b1 (n) @b1 (n) and, assuming gf (n) is statistically independent of b1 (n), fye ;wo (ye (n)) = fbe (T −1 (ye (n))) 1 −|y(n)−b1 (n)|2 =f2 e : (B.4) f2 Since the transmitted symbols are i.i.d., the joint p.d.f. of the K extended vectors, Ye = [ye (0); : : : ; ye (K − 1)], can be written as K−1 1 2 2 fYe ;wo (Ye ) = e−|y(n)−b1 (n)| =f · fb (b1 (n)): (B.5) 2  f n=0 = fb (b1 (n)) · fgf (y(n) − b1 (n)) = fb (b1 (n)) ·

2056

M.F. Bugallo et al. / Signal Processing 81 (2001) 2041–2057

Appendix C. Derivation of the expectation Eb|y;w(i) ˆ [·] Let us consider the general expression for Eb|y;w(i) ˆ [ · ],
Eb|y;w(i) fb|y;w(i) ˆ [%(b)] = ˆ (b) · %(b);

(C.1)

b

is the p.d.f. of symbol b where the summation is over all possible values of a single symbol, fb|y;w(i) ˆ ˆ conditioned to the output y obtained with Blter coe=cients w(i) and %(·) is an arbitrary function. Using the Bayes theorem we can express the a posteriori density fb|y;w(i) ˆ (·), fy|b;w(i) ˆ (·) in terms of fy;w(i) ˆ (·) and fb (·) as 2

fb|y;w(i) ˆ (b(n)) =

2

2

2

fy|b;w(i) fb (b(n))e−|y−b(n)| =f e−|y−b(n)| =f ˆ (y) · fb (b(n)) = ; =  2 2 −|y−b|2 =f2 fy;w(i) Eb [e−|y−b| =f ] ˆ (y) be

(C.2)

where fy|b;w(i) ˆ (·) by simply dropping the expectation (see Eq. (6)). Substituting ˆ (·) is obtained from fy;w(i) (C.2) into (C.1), the following relationship is obtained:  (−1=2 )|y−b|2 2 2 f %(b) Eb [e−|y−b| =f %(b)] be : (C.3) Eb|y;w(i) =  ˆ [%(b)] = 2 2 −|y−b| =f (−1=f2 )|y−b|2 Eb [e ] be References [1] F. Adachi, M. Sawahashi, H. Suda, Wideband CDMA for next-generation mobile communications systems, IEEE Commun. Mag. (1998) 56–59. [2] S.E. Bensley, B. Aazhang, Subspace-based channel estimation for code division multiple access communication systems, IEEE Trans. Commun. 44 (8) (1996) 1009–1020. [3] P. Chaudhury, W. Mohr, S. Onoe, The 3GPP proposal for IMT-2000, IEEE Commun. Mag. 37 (12) (1999) 72–81. [4] B. Chen, P. Willet, The theoretical advantage of CDMA over FDMA in a Gaussian MAC, IEEE Trans. Inform. Theory 45 (6) (1999) 2046–2053. [5] E. Dahlman, B. Gudmundson, M. Nilsson, UMTS=IMT-2000 Based on Wideband CDMA, IEEE Commun. Mag. (1998) 70–80. [6] C.J. Escudero, U. Mitra, D. Slock, A toeplitz displacement method for blind multipath estimation for long code DS=CDMA signals, IEEE Trans. Signal Process., to appear. [7] D.D. Feldman, L.J. Gri=ths, A projection approach for robust adaptive beamforming, IEEE Trans. Signal Process. 42 (4) (1994) 867–876. [8] S. Haykin, in: Adaptive Filter Theory, 3rd Edition, Information and System Sciences Series, Prentice-Hall, Englewood CliCs, NJ, 1996. [9] M.L. Honig, U. Madhow, S. Verd'u, Blind adaptive multiuser detection, IEEE Trans. Inform. Theory 41 (4) (1995) 944–960. [10] F. Kai-Tai, Z. Yao-Ting, Generalized Multivariate Analysis, Science Press, Beijing, Springer, Berlin, 1990. [11] A.M. Kuzminskiy, Finite amount of data eCects in spatio-temporal Bltering for equalization and interference rejection in short burst wireless communications, Signal Processing (80) (2000) 1987–1997. [12] A.M. Kuzminskiy, D. Hatzinakos, Semi-blind spatio-temporal processing with temporal scanning for short burst sdma systems, Signal Processing (80) (2000) 2063–2073. [13] G. Li, Z. Ding, Semi-blind channel identiBcation for individual data bursts in gsm wireless systems, Signal Processing (80) (2000) 2017–2031. [14] H. Liu, G. Xu, Closed-form blind symbol estimation in digital communications, IEEE Trans. Signal Process. 43 (12) (1995) 2714–2723. [15] U. Madhow, Blind adaptive interference suppression for the near-far resistant acquisition and demodulation of direct-sequence CDMA signals, IEEE Trans. Signal Process. 45 (1) (1997) 124–136. [16] U. Madhow, M.L. Honig, MMSE interference suppression for direct-sequence spread-spectrum CDMA, IEEE Trans. Commun. 42 (12) (1994) 3178–3188. [17] G.J. McLachlan, T. Krishnan, The EM Algorithm and Extensions, Wiley Series in Probability and Statistics, Wiley, New York, 1997.

M.F. Bugallo et al. / Signal Processing 81 (2001) 2041–2057

2057

[18] J. M'-guez, L. Castedo, A constant modulus blind adaptive receiver for multiuser interference suppression, Signal Processing 71 (issue 1) (1998) 15–27. [19] J. M'-guez, L. Castedo, A linearly constrained constant modulus approach to blind adaptive multiuser interference suppression, IEEE Commun. Lett. 2 (8) (1998) 217–219. [20] L.B. Nelson, V. Poor, Iterative multiuser receivers for CDMA channels: An EM-based approach, IEEE Trans. Commun. 44 (12) (1996) 1700–1710. [21] C.W. Therrien, Discrete Random Signals and Statistical Signal Processing, Prentice-Hall, Englewood CliCs, NJ, USA, 1992. [22] L. Tong, S. Perreau, Multichannel blind identiBcation: From subspace to maximum likelihood methods, IEEE Trans. Signal Process. 45 (5) (1998) 1343–1348. [23] M.K. Tsatsanis, Inverse Bltering criteria for CDMA systems, IEEE Trans. Signal Process. 45 (1) (1997) 102–112. [24] S. Verd'u, Multiuser Detection, Cambridge University Press, Cambridge, UK, 1998. [25] X. Wang, H.V. Poor, Blind equalization and multiuser detection in dispersive CDMA channels, IEEE Trans. Commun. 46 (1) (1998) 91–103. [26] X. Wang, H.V. Poor, Blind multiuser detection: A subspace approach, IEEE Trans. Inform. Theory 44 (2) (1998) 677–690.