- Email: [email protected]
Blind subspace methods for code and channel estimation in Multicarrier CDMA
Signal Processing 83 (2003) 339 – 357 www.elsevier.com/locate/sigpro
Blind subspace methods for code and channel estimation in Multicarrier CDMA Carlos J. Escudero∗ , Daniel I. Iglesia, Luis Castedo Department of Electronica y Sistemas, Universidad de A Coru˜na, Campus de Elvi˜na, 15071 A Coru˜na, Spain Received 4 February 2002; received in revised form 12 August 2002
Abstract In multicarrier code-division multiple access systems the linear distortion introduced by the channel causes a perturbation in the received user codes that produces multiple access interference. This paper introduces two new techniques to combat this distortion that are based on subspace decomposition. The 6rst one consists of obtaining an approximation of the received desired user code to improve the output signal to interference noise ratio in blind linear receivers. This approximation is obtained by means of the projection of the known transmitted desired user code onto the signal subspace. The second one is a blind channel estimation method that uses the orthogonality property between the received user codes and the noise subspace. We demonstrate that the performance of channel estimation methods based on this property are severely limited by the number of simultaneous users. Nevertheless, it is shown how this limitation is overcome if user codes spanning several transmitted symbols (long codes) are used. ? 2002 Elsevier Science B.V. All rights reserved. Keywords: Multicarrier code-division multiple access systems; Blind channel estimation; Long spreading sequences; Perturbation analysis; Cramer–Rao bound
1. Introduction Multicarrier (MC) transmission methods combined with code-division multiple access (CDMA) communication systems have been proposed as an e=cient technique to combat multipath propagation and have gained an increased interest during the last years [7,3,17,4]. MC modulation and demodulation can be easily implemented using inverse fast Fourier transform (IFFT) and fast Fourier transform (FFT) algorithms, respectively. Depending on the ∗ Corresponding author. Tel.: +34-981-167150; fax: +34-981167160. E-mail address: [email protected] (C.J. Escudero).
spreading, frequency mapping and detection strategy, diDerent techniques based on this combination can be considered: MC-CDMA, MC direct sequence CDMA (MC-DS-CDMA) and multitone CDMA (MT-CDMA). This paper considers MC-CDMA schemes since they outperform the other techniques [4]. In these schemes each user is assigned to a unique identi6cation code sequence and the transmitted signal is split in diDerent subcarriers. It is assumed that the subcarrier bandwidth is smaller than the channel coherence bandwidth and, therefore, exhibits Eat fading only. As a consequence, MC-CDMA systems do not suDer from inter-symbol interference (ISI). However, the eDects of dispersive channels appear as random distortions in the amplitude and phase
0165-1684/03/$ - see front matter ? 2002 Elsevier Science B.V. All rights reserved. PII: S 0 1 6 5 - 1 6 8 4 ( 0 2 ) 0 0 4 1 9 - X
340
C. J. Escudero et al. / Signal Processing 83 (2003) 339 – 357
Nomenclature vni (k) Vni (m) ci (k) x n (k) vn;i j (k) cji (k) 2
i 2 r
kth chip of nth symbol corresponding to the ith user IDFT of vni (k) kth chip of ith user code received kth chip of the nth symbol period kth chip for the jth part of the nth block of symbols for the ith user kth chip of the jth part of the ith user long code ith user signal power noise power
xn
vector of the chips corresponding to the nth symbol period u(k : l) vector whose elements are the ones of u between the indexes k and l (both included)
of each subcarrier. This causes a perturbation in the received user codes and introduces multiple access interference (MAI) due to interfering users transmitting simultaneously with the desired one. In order to implement a multiuser detector and to reduce MAI it is necessary to characterize, implicitly or explicitly, the channel parameters. This paper introduces two new techniques to reduce MAI: the 6rst one obtains an approximation of the received user code to improve the output signal to interference noise ratio (SINR) and the second one is a blind channel estimation technique. Both techniques are based on a subspace decomposition [11] of the received autocorrelation matrix. We will derive particular algorithms for the new methods to identify the received code and channel parameters, respectively. These algorithms can be applied for both uplink and downlink scenarios. Because of the time-varying nature of mobile radio channels, adaptive 6ltering techniques have to be used in order to suppress the MAI. Conventional adaptive 6ltering techniques are based on the minimum mean square error (MMSE) criterion [10], where it is necessary to transmit a training sequence to set up the 6lter coe=cients. To avoid the transmission of
1i diag(a) F
zero vector but the ith element equal to one diagonal matrix whose elements are the elements of vector a DFT matrix
T
transpose operator transpose conjugate operator ∗ conjugate operator ∗ discrete convolution · ceiling operator · Eoor operator det(·) determinant of a square matrix Trace{·} trace of a square matrix H
training sequences, blind techniques based on the linearly constrained minimum variance (LCMV) criterion [8] have been proposed for MAI cancellation in DS-CDMA. However, these techniques are extremely sensitive to inaccuracies in the acquisition of the desired user timing and spreading code. The 6rst technique introduced in this paper is a new method to improve the performance of the blind LCMV criterion in a MC-CDMA system where we cannot obtain a perfect knowledge of the desired user code. This technique tries to approximate the received code by projecting the transmitted one onto the signal– interference subspace of the input autocorrelation matrix. The second estimation technique introduced in the paper takes into account the orthogonality property between the signal and noise subspaces to obtain an estimation of the desired channel impulse response. This technique resembles the subspace channel estimation methods already proposed for single carrier systems [11,16,9]. Nevertheless, when this technique is applied to MC-CDMA systems, it is demonstrated that, in order to obtain a good estimation of the channel, system capacity (i.e., the number of simultaneous users) is a severe constrain. This motivates the
C. J. Escudero et al. / Signal Processing 83 (2003) 339 – 357
consideration of a new structure that obtains channel estimations that are independent of the number of users and will be only limited by the number of MC-CDMA carriers. The new structure considers user codes larger than the user symbols period. In order to show the improvement achieved with the new structure, we obtain the Cramer–Rao bound expression to calculate the lower bound of the estimated parameters variance for any estimator in a MC-CDMA environment using long codes. An adaptive algorithm to obtain the channel estimation is also introduced and, using perturbation techniques, a theoretical approximation of the estimated mean-square error (MSE) is derived. Finally, the structure that can be used to detect transmitted symbols is presented. The paper is organized as follows. Section 2 introduces the basic signal model of a MC-CDMA system. Section 3 presents a new algorithm that is based on the signal subspace and it is demonstrated that it improves the output SINR. Algorithm performance is illustrated by means of some computer simulations. Section 4 introduces an alternative channel estimation technique based on the noise subspace and describes its limitations when it is used in MC-CDMA systems. We also show how to overcome these limitations by using spreading codes larger than the symbol period. The theoretical mean-squared error is obtained for this algorithm and the Cramer–Rao bound for the new signal model is also achieved. Numerical results, at this section, illustrates the performance of the algorithm and the validity of the theoretical analysis. Finally, Section 5 is devoted to the conclusions.
341
obtain the following multicarrier signal: L−1
(2=L)km 1 i vn (k)ej : L
Vni (m) = IDFT[vni (k)] =
(2)
k=0
This signal is transmitted through a dispersive channel with an impulse response hi (m); m = 0; : : : ; M − 1. At the receiver the observed signal is a superposition of the signals corresponding to N users plus an additive white Gaussian noise (AWGN). We will consider the uplink because the down link can be viewed as an special case where all user’s signals are received through the same channel. Therefore, the received signal for the nth symbol is the following: N
Xn (m) =
Vni (m) ∗ hi (m) + rn (m);
(3)
i=1
where ∗ denotes discrete convolution and rn (m) represents a white noise sequence. To recover the transmitted symbols, the receiver applies a L-DFT (discrete Fourier transform) to the received signal (3). Assuming perfect synchronization and a su=ciently large guard time between symbols, the resultant signal is x n (k) = DFT[Xn (m)] =
N
vni (k)H i (k) + n (k)
i=1
=
N
sni ci (k)H i (k) + n (k);
i=1
k = 0; : : : ; L − 1;
(4)
i
i
where H (k) and n (k) are the DFTs of h (m) and rn (m), respectively. Rewriting (4) in vector notation we obtain
2. Signal model
N
Let us consider a discrete-time baseband equivalent model of a synchronous MC-CDMA system with N users using L-chip signature codes [7]. The kth chip corresponding to the nth symbol transmitted by the ith user is given by
xn = [x n (0); : : : ; x n (L − 1)] =
vni (k) = sni ci (k);
where Ci is a diagonal matrix whose elements are the L chips of the code corresponding to the ith user, Hi = [H i (0); : : : ; H i (L − 1)]T and n = [n (0); : : : ; n (L − 1)]T . To obtain (5) we have used the relationship Hi = Fhi where F is the DFT matrix and hi = [hi (0); : : : ; hi (M − 1)]T . Note that (5)
k = 0; : : : ; L − 1;
n = 0; 1; 2; : : : ;
(1)
where ci (k) is the kth chip of the ith user code. In a MC-CDMA system the modulator computes the L-IDFT (inverse discrete Fourier transform) of (1) to
T
sni Ci Hi + n
i=1
=
N i=1
sni Ci Fhi
+ n =
N
sni c˜i + n ;
(5)
i=1
342
C. J. Escudero et al. / Signal Processing 83 (2003) 339 – 357
can be interpreted as a CDMA signal where the user codes are perturbed versions of the transmitted one, i.e., c˜i = Ci Fhi . Assuming statistical independence between users and noise, the autocorrelation matrix of the observations vector (5) can be decomposed as Rx = E[xn xnH ] =
N
∗
H
c˜i E[sni sni ]c˜i + E[n H n]
=
i2 i iH
r2
c˜ c˜ + I;
(6)
where E[ · ] is the expectation operator, superindices ∗ and H denote conjugate and conjugate transpose, 2 2 respectively, i and r are the ith user signal and noise power, respectively, and I is the identity matrix. Let us consider the eigendecomposition of Rx . There are L eigenvalues that we sort as 0 ¿ 1 ¿ · · · ¿ L−1 . It is well known [11] that the eigenvectors associated to the N most significants eigenvalues (ul ; l = 0; : : : ; N − 1) span the signal subspace where the perturbed (received) user codes, c˜i , lie. The remaining L − N eigenvectors (ul ; l = N; : : : ; L − 1) span the noise (orthogonal) subspace and their associated eigenvalues are equal 2 to the noise power, i.e., N = · · · = L−1 = r . 3. Signal subspace technique To produce the receiver output, the observations vector (5) can be processed by a linear combiner w as follows: (7)
It is well known that the SINR at the output is given by 2
d |wH c˜d |2 SINR = H ; w Ri+r w where superindex Ri+r =
N
2
d H
(8)
represents the desired user and 2
2
i c˜i c˜i + r I = Rx − d c˜d c˜d
H
w
subject to wH c˜d = 1:
(10)
It is well known [8] that the solution to this optimization problem is 2
i=1
yn = wH xn :
min E[|yn |2 ]
H
−1 d wLCMV = Rx−1 c˜d ⇒ SINR opt = d c˜d Ri+r c˜ ;
i=1 N
When considering blind adaptive MAI cancellation, the weight vector w in (7) is chosen according to the LCMV criterion [8]
(9)
i=1 i=d
is the autocorrelation matrix associated to the interferer users plus noise.
(11)
where is an arbitrary constant and SINR opt is the maximum SINR attainable. However, the weight vector (11) can only be obtained under ideal conditions where we have a perfect knowledge of the received code c˜d . In practical situations, code acquisition systems [6] provide an estimation of the transmitted code cd =[cd (0); : : : ; cd (L−1)]T . In this case, the weight vector is given by wLCMV = Rx−1 cd
(12)
that is no longer optimum since cd = c˜d . The resulting SINR is extremely low because the receiver interprets the desired signal as an interference and cancels it from the output. In order to show the degradation of the SINR when using (12), let us consider a simple environment with a single user (the desired one) and Gaussian noise. In this case the autocorrelation matrix for the interferences plus noise is reduced to a diagonal matrix 2 Ri+r = r I. Therefore, substituting (12) in (8) and using the matrix inversion lemma, the SINR can be rewritten as SINR opt cos2 ; (13) SINR = 1 + SINR opt (SINR opt + 2) sin2 where represents the angle between the vectors c˜d and cd . Note that depends on the perturbation introduced by the channel. In the absence of channel distortion, = 0, and equation in (13) corresponds with the SINR opt . However, when = 0 the SINR decreases when SINR opt increases. This degradation of the SINR can be avoided if we use in (12) any vector parallel to c˜d instead of cd . One of these vectors can be obtained, even without the knowledge of the perturbed code, by projecting the transmitted code cd onto the signal subspace since this subspace is the line spanned by the vector c˜d . Next subsection presents the projection technique for a general environment.
C. J. Escudero et al. / Signal Processing 83 (2003) 339 – 357
3.1. Projection technique When we consider a general environment with N users and Gaussian noise it is easy to determine the signal–interference subspace but not the desired signal subspace alone. Nevertheless, projecting cd onto the signal–interference subspace results in a better performance as will be shown in the sequel. This projection does not guarantee a vector parallel to the true desired user code, c˜d , but we still obtain an approximation close to it. This technique has already been used in [5] for adaptive beamforming with uncalibrated arrays. Indeed, the transmitted user code can be decomposed as the sum cd = csd + crd ; csd
(14) crd
where and are two orthogonal components that lie in the signal–interference and noise subspace, respectively. These two subspaces can be determined from the eigenvectors and eigenvalues of the autocorrelation matrix Rx [11]. Substituting (14) and (12) in (8), the SINR for the LCMV technique can be expressed as follows: H −1 H csd Rx−1 csd + crd Rx−1 crd SINR LCMV = −1 : d2 |csdH Rx−1 c˜d |2 (15) This expression can be interpreted as a generalization of (13) to a N users environment and, in general, yields to values which are considerably smaller than the optimum SINR. However, projecting the transmitted user code onto the signal plus interference subspace is equivalent to eliminating the component crd and the resulting output SINR is given by −1 H csd Rx−1 csd SINR proj = −1 : (16) d2 |csdH Rx−1 c˜d |2 It is straightforward to see that this value is higher than (15). This indicates that an improvement in performance is obtained when using the projected vector csd instead of cd . 3.2. Algorithm It is well known that the signal–interference subspace is spanned by the eigenvectors corresponding to the N most signi6cant eigenvalues of the autocorrelation matrix Rx [11]. Considering the L × N
343
unitary matrix U containing these signi6cant eigenvectors, the signal projection matrix is given by P = U(UH U)−1 UH = UUH . Therefore, the projection of cd onto the signal-interference subspace is obtained by multiplying the projection matrix with the code, i.e., Pcd = csd . If we use this projection into the linear constraint of (10), the weight vector (12) can be written as follows: wproj = Rx−1 Pcd = Rx−1 csd :
(17)
The resulting linear receiver outperforms the conventional LCMV receiver since it achieves the output SINR given by (16). In a practical implementation, the autocorrelation matrix is unknown and Ns −1 xn xnH , where Ns is estimated as Rˆ x = (1=Ns ) n=0 is the number of received symbols used to estimate the matrix. This 6nite number of data observations produce errors in the estimation that reduce the SINR similarly to situations where we do not have a perfect knowledge of the received codes. Estimation errors can be assumed as an error in the knowledge of the received codes. Therefore, for a 6xed number of data observations, the projection technique will obtain a higher output SINR. As a result, the proposed projection algorithm exhibits faster convergence even when there is no channel distortion or we have a perfect knowledge of the received code. 3.3. Numerical results To illustrate the eDectiveness of the signal subspace technique, we have carried out several computer simulations. This section presents the results of these simulations using the following common environment: channel length M = 4 and 100 averaged realizations using diDerent random channels whose coe=cients are drawn from independent Gaussian complex random variables. The remaining parameters, i.e., number of users, code length, type of code, SNR and number of subcarriers will be speci6ed in each simulation. All simulations consider an uplink scenario, i.e., the channel is assumed to be diDerent for each user. In the 6rst simulation we consider random codes with length L = 8. Here we compare the performance of the projection technique with respect to the optimum (in the sense of perfect knowledge of the desired user code) and the LCMV technique. Fig. 2
C. J. Escudero et al. / Signal Processing 83 (2003) 339 – 357
ci,0 (0)
i snp
Symbols
Serial/Parallel Converter
ci,0 (L-1) ci,j (0)
i snp+j
i vn,0 (0)
i vn,0 (L-1)
i vn,j (0)
vn,ji (L-1) ci,j (L-1) ci,p-1 (0)
i snp+p-1
ci,p-1 (L-1)
Parallel/Serial Converter
344
MC-CDMA signal
i vn,p-1 (0)
i vn,p-1 (L-1)
Fig. 1. Block diagram of the ith user MC-CDMA modulator for transmitting the nth block of p symbols using long codes of length pL, where L is the spreading gain.
12
10
8
6
SINR (dB)
4 Optimum With Projection LCMV
2
0
-2
-4
-6
-8
0
5
10 SNR (dB)
Fig. 2. SINR vs. SNR for optimum, projection and LCMV techniques.
15
C. J. Escudero et al. / Signal Processing 83 (2003) 339 – 357
345
0
10
10 -1
-2
SER
10
10 -3
10
-4
10
-5
Optimum With Projection LCMV
0
5
10
15
SNR (dB)
Fig. 3. SER vs. SNR for optimum, projection and LCMV techniques.
plots the average of the output SINR for diDerent channel realizations with respect to the input signal to noise ratio (SNR) considering three users transmitting with the same input SNR. Note that the output SINR corresponding to the LCMV reduces as the optimum increases, whereas the projection technique improves the performance by reducing the degradation of the SINR with respect to the optimum value. To illustrate the performance of the proposed projection algorithm in terms of the symbol error rate (SER), Fig. 3 plots it with respect to the input SNR for the same three users environment as before. Note that the curve corresponding to the projection method decreases with the input SNR and is close to the SER obtained when we have a perfect knowledge of the received code. Fig. 4 represents the output SINR time evolution of the projection algorithm. In this simulation we have considered the same environment as before but with the users transmitting with an input SNR = 12 dB.
Note that the projection technique presents an asymptotic convergence better than the LCMV technique. The diDerence between the convergence value of the projection technique and the optimum (in the sense of knowing the code) value depends on the channel features. Moreover, as we have said in Section 3.2 the proposed detector exhibits faster convergence than the one considering knowledge of channel impulse response (curve labelled as optimum) since the estimated channel takes into account errors in the estimated autocorrelation matrix when a small amount of symbols are available. 4. Noise subspace technique As we have seen at the signal model section, the perturbed user codes lie in the signal subspace and, therefore, are orthogonal to the noise subspace. In this case, we will try to perform an estimation of the desired user channel impulse response vector, hd .
346
C. J. Escudero et al. / Signal Processing 83 (2003) 339 – 357 15
10
SINR (dB)
5
0 Optimum With Projection LCMV -5
-10 0
100
200
300
400
500 Symbols
600
700
800
900
1000
Fig. 4. SINR time evolution for projection and LCMV techniques.
Considering the eigendecomposition of the autocorrelation matrix Rx , we have seen that the eigenvectors associated to the N most signi6cants eigenvalues (ul , l = 0; : : : ; N − 1) span the signal subspace where the perturbed (received) user codes, c˜i , lie and the remaining L − N eigenvectors (ul , l = N; : : : ; L − 1) span the noise (orthogonal) subspace. Therefore, this property can be used to obtain the channel impulse response vector for the dth user, hd , by means of the following system of equations: H
c˜d ul = (Cd Fhd )H ul = 0;
l = N; : : : ; L − 1:
(18)
Recall that this system of equations has M unknowns and L − N equations. It will be solvable if and only if the number of equations is greater or equal than the number of unknowns, that is M 6 L − N ⇒ N 6 L − M:
(19)
Note that we have obtained a limit for the system capacity that depends on the number of carriers and the channel length. This means that if we want to improve the capacity it is necessary to increase the number of
carriers, L. However, in many practical applications this is not possible, or desirable, and we have to consider another solution to increase the system capacity. 4.1. Utilization of long codes In the previous sections we have considered that the codes length is equal to the spreading gain (i.e., to the number of carriers), thus 6nding that the system capacity is limited by the number of carriers. In this section, we are going to explore the idea of increasing the system capacity by means of increasing the codes length without modifying the number of carriers, i.e., we will consider codes with a length larger than the spreading gain. To simplify the notation we consider pL-chip signature codes, where p ¿ 1 is an integer number (i.e., the code length is a multiple of the spreading gain). This means that the user codes will be repeated after p symbols. The modulator for a block of p consecutive symbols is represented in Fig. 1 where the jth symbol of this block is modulated by the jth part of
C. J. Escudero et al. / Signal Processing 83 (2003) 339 – 357
the user code to obtain the following CDMA signal: i cji (k); vn;i j (k) = snp+j
j = 0; : : : ; p − 1;
k = 0; : : : ; L − 1;
N
vn;i j (k)H i (k) + n; j (k)
i=1
=
N
i snp+j cji (k)H i (k) + n; j (k)
i=1
k = 0; : : : ; L − 1;
is a pL×p matrix. The pL×pL autocorrelation matrix of (22) can be expressed as follows: R = E[x(n)xH (n)]
n = 0; 1; 2; : : : ;
where cji (k) is the kth chip of the code associated to the jth symbol of a block for the ith user. This equation is equal to (1) for p = 1. As in Section 2, we can obtain the received signal after DFT as x n; j (k) =
(20)
where n; j (k) is the DFT of the noise sequence superposed to the jth symbol. A vector of L chips for the jth symbol of a block is
=
N
i snp+j Cij Fhi + n; j ; j = 0; : : : ; p − 1; (21)
where n; j = [n; j (0); : : : ; n; j (L − 1)]T and Cij is a L × L diagonal matrix whose elements are the L chips of the code corresponding to the jth sequence of a block for the ith user. Arranging the p consecutive vectors xn; j , corresponding to a single block, into a single vector of length pL, we obtain x(n) = [xn;T 0 ; : : : ; xn;T p−1 ]T N
Ci Hi si (np) + n ;
(22)
i=1
where Ci is a diagonal matrix containing the pL chips of the code associated to the ith user, i i si (np) = [snp ; : : : ; snp+p−1 ]T , n = [Tn; 0 ; : : : ; Tn; p−1 ]T and i 0 ··· 0 Fh .. .. i 0 . Fh . (23) Hi = .. . . .. .. . 0 0
···
0
Fhi
2
H
H
2
i Ci Hi Hi Ci + r I:
(24)
Let us consider the eigendecomposition of R and sort its eigenvalues as 0 ¿ 1 ¿ · · · ¿ pL−1 . Note that matrix Ci Hi has p independent columns and, therefore, the eigenvectors associated to the pN most signi6cants eigenvalues (ul , l = 0; : : : ; pN − 1) span the signal subspace where the diDerent user signals lie. The remaining p(L − N ) eigenvectors (ul , l = pN; : : : ; pL − 1) span the noise subspace. Taking into account the orthogonality between the signal and the noise subspace, we obtain the following system of equations for the desired user, d: H
H
(Cd Hd )H ul = 0 ⇒ hd FH Cdj ul; j = 0; l = pN; : : : ; pL − 1; j = 0; : : : ; p − 1;
i=1
=
N i=0
xn; j = [x n; j (0); : : : ; x n; j (L − 1)]T =
347
(25)
where ul; j = ul (jL + 1 : jL + L − 1) and the notation u(k : l) represents a vector containing the elements of u between the indexes k and l (both included). Recall that this system of equations still has M unknowns but p2 (L − N ) equations. It will be solvable if and only if the number of equations is greater or equal than the number of unknowns, that is, 2 p L−M 2 ; (26) p (L − N ) ¿ M ⇒ N 6 p2 where · is the Eoor operation which rounds down to the previous integer value. Therefore, the system capacity can be improved by increasing the p value, i.e., the code length. Similarly, there exists a relationship for the parameter p when there is a 6xed number of users
M p¿ ; (27) L−N where · is the ceiling operator. Finally, in order to simplify the subsequent analysis, the system of equations (25) can be rewritten as follows: H
H
hd FjH Cd ul = 0; j = 0; : : : ; p − 1;
l = pN; : : : ; pL − 1; (28)
348
where Fj = [0M ×L ; : : : ; 0M ×L , 0M ×L ]T is a pL × M matrix.
C. J. Escudero et al. / Signal Processing 83 (2003) 339 – 357
FT
, 0M ×L ; : : : ;
( j+1)th position
4.2. Cramer–Rao bound In order to show the eDect of considering long codes in a MC-CDMA environment, this section computes the channel estimation Cramer–Rao bound (CRB) for any value of p. The CRB is useful in the sense of that it allows to determine the best-possible performance (smallest variance) of any parameter estimator. Estimators whose variance is close or equal to this bound can be said to have an optimal behaviour. In our case, it is possible to compute the CRB and, therefore, determine which is the best-possible performance of the channel impulse response estimation. Moreover, by means of this CRB, we shall be able to know if the best-possible performance obtained using spreading codes equal to the symbol period would be improved by using long codes. In a few words, CRB gives a lower bound on the error covariance matrix for a parameter vector to estimate (in our case, the channel impulse response vector) [13], and it provides a benchmark against which algorithms performance can be compared. The CRB of a parameter estimator, in our case of the channel impulse response, is determined by the fact that the error covariance matrix for this parameter vector has to be greater or equal than the inverse of the Fisher information matrix (FIM). FIM is the covariance matrix of the score function which is the gradient of the log-likelihood function. The standard way of computing the CRB is explained in [9], where the FIM is obtained considering as parameters to estimate 2 2 r , i , sni and hi (m) for i = 1; : : : ; N , m = 1; : : : ; M and n = 0; : : : ; T − 1 where T is the number of observed vectors (22). Nevertheless, the complexity of this technique depends on many parameters and it is rather cumbersome. To simplify the CRB analysis. we only consider the components of the desired user channel, as parameters to estimate, that is, the remainder parameters are considered as nuisance [13]. This approximation will provide a lower variance bound, but still good to compare the performance of the new algorithm.
If we consider a 6xed number, T , of independent observed vectors, we can assume a multivariate normal model, x = (x(0); x(2); : : : ; x(T − 1)), x(n) : CN[x(n); R] where x(n) = E[x(n)]. The joint distribution of the complex random sample has density f (x) = −pLT (det(R))−T T −1 H −1 ×exp − (x(n) − x(n)) R (x(n) − x(n)) ; n=0
(29)
where det(·) represents the determinant of a square matrix. We can generalize the results in [13, pp. 235– 237] to the complex case by considering a vector T T of real parameters = [real(hd ); imag(hd )]T where real(·) and imag(·) are the real and imaginary parts of a complex vector. The FIM associated to is T @ @ ln f (x) ln f (x) : (30) J = E @ @ To compute the FIM it is easier to obtain it with respect T H to the vector h = [hd ; hd ]T and to relate it with (30) as in [1] J =
I
I
jI
−jI
Jhh
I
I
jI
−jI
T :
(31)
By doing a similar analysis to [14,1] we can obtain Jh d h d ∗ Jh d h d ; (32) Jhh = Jh d ∗ h d Jh d ∗ h d ∗ where the elements (i; j) of the four submatrices in (32) are @R @R ; Jhd hd (i; j) = T Trace R−1 d R−1 d @h (i) @h (j) @R @R −1 Jhd hd∗ (i; j) = T Trace R R @hd (i) @hd∗ (j) T −1 @xH (n) −1 @x(n) + Trace R ; @hd∗ (j) @hd (i)
−1
n=0
C. J. Escudero et al. / Signal Processing 83 (2003) 339 – 357
@R @R −1 R J (i; j) = T Trace R @hd∗ (i) @hd (j) T −1 @xH (n) @x(n) + Trace R−1 d( j) ; @hd∗ (i) @h n=0 @R @R d∗ d∗ −1 −1 Jh h (i; j)= T Trace R R ; @hd∗ (i) @hd∗ (j) (33)
hd∗ hd
−1
where Trace {·} represents the trace of a square matrix. It is interesting to note that, under our assumptions, N x(n) = i=1 Ci Hi si (np). To compute the FIM it is necessary to specify the value of the partial derivatives in (33), i.e., 2 H H @R = d C d B i H d Cd ; @hd (i)
(34)
H H @xH (n) = (sd (np))∗ Bi Cd ; ∗ d @h (i)
0
···
0
In order to solve the system of equations (28), we can consider the following equivalent system: H
H
H
j=0 l=pN
de6ned as follows: ;
H
p−1 pL−1
= arg min hd hd 2 =1
where 1i = [0; : : : ; 1 ; 0; : : : ; 0]T is a vector with the ith element equal to one and zero the others. As we saw before, by using blind estimation techniques using second-order statistics, channel impulse response can be determined up to a complex unknown factor and, therefore, FIM is rank de6cient. In order to avoid the ambiguity of the unknown factor and obtain the CRB, we can assume that one non-zero complex channel parameter, hd (i) remains constant and, as a result, its corresponding derivative vanishes. Empirically, it was observed that assuming constant the
H
p−1
= arg min hd hd 2 =1
H
FjH Cd ul ulH Cd Fj hd H
FjH Cd UU Cd Fj hd H
j=0 H
= arg min hd Qd hd ; hd 2 =1
j=0 l=pN
(35)
(36)
for l=pN; : : : ; pL−1 and j=0; : : : ; p−1. On the other hand, the solution to (36) can be found by solving the following minimization problem: p−1 pL−1 H H hˆd = arg min hd FjH Cd ul ulH Cd Fj hd hd 2 =1
F1i
ith element
4.3. Channel identi7cation algorithm
H
@x(n) = Cd Bi sd (np); @hd (i)
where B is a pL × p matrix 0 ::: 0 F1i .. .. 0 . F1i . i B = .. .. .. . . . 0
channel parameter with largest energy [1] the minimum CRB is obtained. The eDect of 6xing this complex unknown parameter [9,1] is equivalent to deleting the rows and columns corresponding with its real and imaginary parts from the FIM before computing its inverse to obtain the CRB. Section 4.6 will present numerical results of the CRB for MC-CDMA systems using long codes with diDerent values of p. We will see that systems satisfying condition (27) yield the lowest CRB.
hd FjH Cd ul 2 = hd FjH Cd ul ulH Cd Fj hd = 0
2 H H @R = d Cd Hd B i Cd ; ∗ d @h (i)
i
349
(37)
where the solution hˆd is an estimation of the channel impulse response vector, U is a pL × p(L − N ) matrix whose columns are the eigenvectors associated to the noise subspace (i.e., ul ; l = pN; : : : ; pL − 1) and p−1 H Qd = j=0 FjH Cd UUH Cd Fj . The solution can be obtained by the least-squares method and it corresponds to the eigenvector of Qd associated to its minimum eigenvalue [15]. Note that we have included the constraint hd 2 = 1 in order to avoid the null solution. In practice, we do not know a priori the autocorrelation matrix (24) but we can estimate it from the
350
C. J. Escudero et al. / Signal Processing 83 (2003) 339 – 357
sampled averaged matrix as 1 Rˆ = Ns
N s −1
x(n)xH (n);
and (38)
n=0
where Ns is the number of received pL chip blocks used to obtain the estimation. Finally, note that when using second-order statistics, the channel impulse response can be obtained up to a complex constant that should be compensated in order to remove undesirable ambiguities and analyze the algorithm performance. Towards this aim, we normalize the estimation of the impulse response vector with respect to the ith element as follows: hˆdnormalized = #hˆd ;
(39)
where # = hd (i)= hˆd (i) and hd (i) and hˆd (i) are the ith elements of the true and estimated channel impulse response vectors, respectively.
ˆ Ph −Q† PQh = −Q† (Qˆ − Q)h = −Q† Qh; (42) where Q† denotes the left pseudo-inverse of Q. The kth component of Ph is given by p−1 ˆ = −qkH Ph(k) −qkH Qh FjH CH Uˆ Uˆ H CFj h j=0 p−1 pL−1
=−
j=0
qkH FjH CH uˆl uˆH l CFj h l=pN scalar
scalar
p−1 pL−1
=−
j=0
H H H uˆH l CFj hqk Fj C uˆl l=pN
pL−1
=−
H H uˆH l CFQk C uˆl
l=pN
4.4. Mean square error analysis
H ˆ = −Trace{Uˆ H CFQH k C U}
In this section, we derive an analytical expression of the estimation MSE. 1 To clarify notation, let us denote hd = h; Qd = Q and Cd = C and consider the following identities: Qh = 0; hˆ = h + Ph;
(40)
Qˆ = Q + PQ: The perturbation terms Ph and PQ represent the estimation errors of h and Q, respectively. Thus, the channel estimation mean square error is given by E[Ph2 ]. The key point of our analysis is the utilization of a perturbation technique [12] that allows us to express Ph in terms of PQ. For a su=ciently large number of samples (Ns → ∞), Qˆ → Q, hˆ → h and Qˆ hˆ is approximately equal to the zero vector, i.e. Qˆ hˆ = (Q + PQ)(h + Ph) PQh + QPh ≈ 0; where we have consider negligible the second-order term, PQPh 0. Therefore, QPh −PQh 1
(41)
A brief description of this analysis was presented in [2].
H = −Trace{Uˆ Uˆ H CFQH k C };
(43)
where qk is the kth column of (Q† )H , Uˆ is the estimated matrix U, F = [F0 h; : : : ; Fp−1 h], Qk = [F0 qk ; : : : ; Fp−1 qk ] and Trace{·} represents the trace of a square matrix. Based on the results of [14, p. 1840, Eq. (4.11)], we can obtain the following identity: Uˆ Uˆ H CF −UUH PVVH CF;
(44)
where V is a pL × pN matrix whose columns are the eigenvectors associated to the signal subspace (i.e., ul ; l = 0; : : : ; pN − 1) and PV = Vˆ − V. Moreover, from Appendix A of [14, p. 1844, Eq. (A.2)] −1 ˆ Uˆ H PV UH RV) ;
(45) 2
2
where ) = diag(0 − r ; : : : ; pN −1 − r ) where diag(a) is a diagonal matrix whose elements are the elements of vector a. To remove the eDect of the unknown constant that we have in the estimation of the channel vector, we have to consider a normalization of the estimated vector channel. Normalization (39) is a function of the channel estimate. In order to facilitate the analysis of the mean square error approximation, we must seek a normalization independent of the error in the channel
C. J. Escudero et al. / Signal Processing 83 (2003) 339 – 357
estimate. Thus, we can approximate the normalization scalar in (39) as #= =
h(i) h(i) = ˆ h(i) − Ph(i) h(i) 1 Ph(i) 1− ; 1 + Ph(i)=h(i) h(i)
(46)
where we have used the fact 1=(1 + x) 1 − x for small x. Therefore, we can consider the following normalization: Ph(i) ˆ h−h Phnormalized = #hˆ − h = 1 − h(i) Ph(i) ˆ h h Ph − Ph(i) h(i) h(i) h1Ti Ph; (47) = I− h(i) = Ph −
4.5. Linear detector Once we have the channel estimation we can obtain the transmitted symbols corresponding to the jth sequence by using a linear combiner, i.e., yj = wjH xn; j ;
wj = %Rj−1 Cdj Fhd ;
−1 H H ˆ Ph(k) Trace{UUH RV) V CFQH kC }
SINR j =
=
1 SINR j ; p
(52)
j=0
d |wjH (Cdj Fhd )|2 wjH [Rj − d2 (Cdj Fhd hdH FH Cdj )]wj H
(53)
is the jth sequence output SINR.
−1 H H ˆ (ulH RV) V CFQH k C ul )
4.6. Numerical results ˆ l; k ; ulH Rg
(48)
l=pN H where gl; k = V)−1 VH CFQH k C ul . Finally, in Appendix A it is demonstrated that the MSE of the channel estimation algorithm is given by
E[Ph2 ] =
p−1
SINR =
2
l=pN pL−1
(51)
where % is an arbitrary constant and Rj is the L × L autocorrelation matrix for sequence xn; j . Note that Rj is a submatrix of R in (24) that considers the columns and rows between the jth and the (jL−1)th indices. In practical situations we use estimations of hd and Rj . The performance of the linear combiner will be evaluated in terms of the output SINR. We can obtain a compact expression of the SINR by using an averaged value over the p diDerent combiners, i.e.,
where
=
(50)
where wj is a coe=cients vector for the jth combiner. If we select this vector to maximize the output SINR, we obtain the following expression [8]:
where we have assumed the second-order term as negligible i.e., Ph(i)Ph 0. To include this normalization in (43) we can consider qk as the kth column of the matrix ((I − h1Ti =h(i))Q† )H . Combining (44) and (45) in (43), we obtain the following expression:
pL−1
351
2 M −1 r ˜ ˜H (Trace{UH UGH k CC Gk } Ns
k=0
2
+ r Trace{UH UGH k Gk });
(49)
where Gk = [g(pN ); k ; : : : ; g(pL−1); k ] and CQ = [ 1 C1 H1 ; : : : ; N CN HN ]. In order to obtain this expression we had to explore the fourth-order statistics of binary and Gaussian random variables.
To show the algorithm performance of this noise subspace technique, we show the results of some computer simulations using a common environment, i.e., channel length M = 4 and 100 averaged realizations using diDerent random channels whose coe=cients are drawn from independent Gaussian complex random variables. The remaining parameters, i.e., number of users, code length, type of code, SNR and number of subcarriers will be speci6ed in each simulation. All simulations consider an uplink scenario, i.e., the channel is assumed to be diDerent for each user. Next we evaluate the performance of the noise subspace technique in terms of the MSE, that it is de6ned as MSE = (hˆinormalized − hi )2 :
(54)
352
C. J. Escudero et al. / Signal Processing 83 (2003) 339 – 357 4
10
N=8 N=7 N=6 3
10
2
MSE
10
1
10
0
10
10
10
-1
-2
0
50
100
150
200 250 Symbols
300
350
400
450
500
Fig. 5. Mean square error time evolution for diDerent values of N .
Fig. 5 presents the averaged MSE time evolution for diDerent values of N . An environment using short (L=10) random codes has been considered. The users are received with a SNR = 10 dB. According to (19) a value for the number of users N = 7 is acceptable but it can be seen that for higher values, it is achieved a poor channel estimation. Fig. 6 presents the averaged MSE time evolution using long codes with diDerent values of p. We have considered the same environment as in the previous simulation but with 8 users. In this case relationship (26) reduces to p ¿ 2. It can be seen that for values below this condition (i.e., p=1), the channel estimation does not achieve a good performance. To show the eDectiveness of the proposed technique, Fig. 7 plots the MSE time evolution of the simulation and the Cramer–Rao bound (CRB) for an environment with 8 users, L=12, SNR =12 dB, random codes and p=2 (i.e., satisfying condition (27)). More-
over, it is plotted a third curve with the theoretical MSE time evolution (49). It can be seen that theoretical and simulated curves are very similar showing the validity of expression (49). On the other hand, MSE is very close to the CRB even for a small amount of received symbols. To show the bene6ts of using long codes in terms of the MSE, Fig. 8 shows the averaged CRB with the same length channel and codes as before for N = 5, L = 8, SNR = 20 dB and three diDerent values of p. It can be seen that for values satisfying (27) (i.e. p = 2 and p = 3) the resulting CRB is almost the same. However, for values not satisfying this condition the estimator yields a worse limit. Finally, Fig. 9 shows the averaged simulated and theoretical MSE vs. the SNR for an environment of N =8 of synchronous users using random codes, L=12 subcarriers and p = 1. The curves are obtained after Ns = 300 symbols. We can see that both curves are very similar even for small values of SNR.
C. J. Escudero et al. / Signal Processing 83 (2003) 339 – 357
353
4
10
p=3 p=2 p=1 3
10
2
MSE
10
1
10
0
10
10 -1
10 -2
0
50
100
150
200
250 300 Symbols
350
400
450
500
Fig. 6. Mean square error time evolution for diDerent values of p. 2
10
Simulation Theoretical CRB 1
10
0
10
-1
10
-2
10
-3
MSE
10
10-4
10
-5
0
50
100
150
200
250
Symbols
Fig. 7. CRB, theoretical MSE and simulated MSE for the channel identi6cation algorithm, vs. number of symbols (p = 2).
354
C. J. Escudero et al. / Signal Processing 83 (2003) 339 – 357 0
10
p=3 p=2 p=1 -1
10
-2
MSE
10
10-3
-4
10
-5
10
0
50
100
150 Symbols
200
250
300
Fig. 8. CRB vs. number of symbols for diDerent values of p.
5. Conclusions In this paper we have investigated subspace decomposition techniques to combat the MAI due to received user code perturbations produced in uplink, or downlink, MC-CDMA systems. This code perturbation is typically caused by the channel linear distortion. We focused our attention on blind linear adaptive receivers based on the LCMV criterion. The performance of these receivers is very sensitive to inaccuracies in the acquisition of the received user code and typically result in an extremely low input SINR. Two diDerent subspace methods have been considered. The 6rst one consists of building the linear constraint in LCMV receivers using a projection of the transmitted desired user code onto the signal subspace. This yields to a new code closer to the received one and enables us to considerably increase the output SINR. This technique also accounts for the errors in the estimation of the autocorrelation matrix and thus
exhibits faster convergence. Computer simulations illustrate the performance of the projection technique. The second method is a blind channel estimation algorithm that exploits the orthogonality property between the received user codes and the noise subspace in a way similar to single user blind channel estimation techniques [11]. Nevertheless, we demonstrate that channel estimation techniques based on this property are severely limited by the number of carriers and simultaneous users when applied to MC-CDMA systems. This limitation, however, can be overcome by using codes larger than the symbol period (long codes). This is shown by calculating the Cramer–Rao bound of any channel estimator based on this subspace property in long codes scenarios. Finally, we have derived a novel channel estimation algorithm speci6cally devised for long codes and calculated an analytical expression of the channel estimation MSE. Computer simulations have been also carried out to illustrate the performance of the proposed channel estimation
C. J. Escudero et al. / Signal Processing 83 (2003) 339 – 357 10
355
-1
Simulation Theoretical
MSE
10 -2
10
-3
10
-4
10
-5
0
2
4
6
8
10 SNR (dB)
12
14
16
18
20
Fig. 9. Simulated and theoretical MSE vs. received users SNR.
technique and to corroborate the accuracy of the obtained MSE analytical expression.
H
N
Ci Hi si (np)H n
i=1
+
Acknowledgements
H
×(sj (np))H Hj Cj + N
H H n (si (np))H Hi Ci + n H n gl; k :
i=1
This work has been supported by Ministerio de Ciencia y TechnologRSa of Spain and FEDER funds from the European Union, under grant number TIC2001-0751-C04-01.
Appendix A Taking into account the estimated autocorrelation ˆ Eq. (48) can be expressed as matrix (38), R, pL−1 N N N s −1 1 i i i Ph(k) = ulH C H s (np) Ns i=pN n=0
i=1 j=1
(A.1) ulH (Ci Hi )=0
From (25) we know that and, therefore, it is straightforward to show that pL−1 Ns −1 1 ulH Ph(k) = Ns l=pN n=0 N H i H iH iH × n (s (np)) H C + n n gl; k : i=1
Therefore, the MSE is M −1 E[Ph2 ] = E[Ph(k)Ph∗ (k)] k=0
(A.2)
356
C. J. Escudero et al. / Signal Processing 83 (2003) 339 – 357
=
M −1 k=0
H H + ulH E[n H m ]ut gt; k E[m n ]gl; k
pL−1 pL−1 Ns −1 Ns −1 1 Ns2 i=pN t=pN
4
N N
4
= r ulH ut gt;Hk gl; k '(n − m);
H E[ulH n H m ut gt; k
i=1 j=1 H
H
×Ci Hi si (mp)(sj (np))H Hj Cj gl; k ] +
pL−1 pL−1 Ns −1 Ns −1 1 2 Ns t=pN
H H × E[ulH n H n gl; k gt; k m m ut ] ;
×
H
×Hj Cj gl; k 2
2
= ulH ( r '(n − m)I)ut gt;Hk Ci Hi ( i '(i − j) H
2
2
r i ulH ut gt;Hk Ci Hi Hi
H
4
× Ci gl; k + r ulH ut gt;Hk gl; k
=
2 M −1 r ˜ ˜H (Trace{UH UGH k CC Gk } Ns
k=0
H
H i i i j H = ulH E[n H m ]ut gt; k C H E[s (mp)(s (np)) ]
N
H
H i i i j H j j E[ulH n H m ut gt; k C H s (mp)(s (np)) H C gl; k ]
H
n=0
i=1
(A.3)
H
M −1 pL−1 pL−1 Ns −1 1 Ns2 t=pN k=0 l=pN
where we have used the fact that the third order moments of a Gaussian random variable are zero. Considering statistical independence between users and noise and the user symbols i.i.d., the 6rst expectation (A.3) is
(A.6)
note that we have used the fact ulH gl; k = H ulH (V)−1 VH CFQH k C ul ) = 0. Including (A.4) and (A.6) in (A.3), it is obtained E[Ph2 ] =
n=0 m=0
l=pN
4
= r ulH gl; k gt;Hk ut + r ulH '(n − m)ut gt;Hk '(m − n)gl; k
n=0 m=0
2
+ r Trace{UH UGH k Gk });
(A.7)
where Gk = [g(pN ); k ; : : : ; g(pL−1); k ; ] and C˜ = [ 1 C1 H1 ; : : : ; N CN HN ]. That is, we have obtained the expression in (49).
H
×'(n − m)I)Hj Cj gl; k 2
2
H
References
H
= i r ulH ut gt;Hk Ci Hi Hj Cj gl; k ×'(n − m)'(i − j);
(A.4)
where '(·) is the Kronecker function. It is well known that for any four complex Gaussian variables [12], i i = 1; 2; 3; 4, E[ 1 2∗ 3 4∗ ] = E[ 1 2∗ ]E[ 3 4∗ ] + E[ 1 4∗ ]E[ 3 2∗ ]: (A.5) In the second expectation of (A.3) our four complex variables are uiH n , gl;Hk n , gt;Hk m and utH m . The expectation can be expressed as H H E[ulH n H n gl; k gt; k m m ut ] H H = ulH E[n H n ]gl; k gt; k E[m m ]ut
[1] J.L. Bapat, Partially blind identi6cation of FIR channels for QAM signals, PhD. Dissertation, Department of Electrical Engineering, Pennsylvania State University, August 1996. International Editions, Singapore, 1993. [2] C.J. Escudero, D.I. Iglesia, M.F. Bugallo, L. Castedo, Analysis of a subspace channel estimation technique for multicarrier CDMA systems, Proceedings of the IEEE Signal Processing Workshop on Statistical Signal and Array Processing, Pennsylvania-USA, August 14 –16, 2000, pp. 184 –188. [3] K. Fazel, G.P. Fettweis, Multi-Carrier Spread-Spectrum, Kluwer Academic Publishers, Dordrecht, 1997. [4] K. Fazel, S. Kaiser, Multi-Carrier Spread-Spectrum & Related Topics, Kluwer Academic Publishers, Dordrecht, 2000. [5] D.D. Feldman, L.J. Gri=ths, A projection approach for robust adaptive beamforming, IEEE Trans. Signal Process. 42 (4) (April 1994) 867–876.
C. J. Escudero et al. / Signal Processing 83 (2003) 339 – 357 [6] S. Glisic, B. Vucetic, Spread Spectrum CDMA Systems for Wireless Communications, Artech House Publishers, Norwood, MA, 1997. [7] S. Hara, R. Prasad, Overview of multicarrier CDMA, IEEE Commun. Mag. 35 (December 1997) 126 –133. [8] M.L. Honig, U. Madhow, S. VerdRu, Blind adaptive multiuser detection, IEEE Trans. Inform. Theory 41 (4) (July 1995) 944–960. [9] H. Liu, G. Xu, Closed-form blind symbol estimation in digital communications, IEEE Trans. Signal Process. 43 (11) (November 1995) 2714–2723. [10] U. Madhow, M.L. Honig, MMSE interference suppression for direct-sequence spread-spectrum CDMA, IEEE Trans. Comm. 42 (12) (December 1994) 3178–3188. [11] E. Moulines, P. Duhamel, J.F. Cardoso, S. Mayrargue, Subspace methods for the blind identi6cation of multichannel FIR 6lters, IEEE Trans. Signal Process. 43 (2) (February 1995) 516–525.
357
[12] W. Qiu, Y. Hua, Performance analysis of the subspace method for blind channel identi6cation, Signal Processing 50 (1996) 71–81. [13] L.L. Scharf, Statistical Signal Processing: Detection, Estimation, and Time Series Analysis, Addison-Wesley, Reading, MA, 1991. [14] P. Stoica, T. SYoderstrYom, Statistical analysis and subspace rotation estimates of sinusoidal frequencies, IEEE Trans. Signal Process. 39 (8) (August 1991) 1836–1847. [15] G. Strang, Linear Algebra and its Applications, 3rd Edition, Harcourt, Brace and Jovanovich, New York, 1988. [16] M. Torlak, G. Xu, Blind multi-user channel estimation in asynchronous CDMA systems, IEEE Trans. Signal Process. 45 (January 1997) 137–147. [17] N. Yee, J.P. Linnartz, G. Fettweis, Multi-carrier CDMA in indoor wireless radio networks, Proceedings of the International Symposium on Personal, Indoor and Mobile Radio Communications, Yokohama, 1993, pp. 109 –113.
Blind subspace methods for code and channel estimation in Multicarrier CDMA Carlos J. Escudero∗ , Daniel I. Iglesia, Luis Castedo Department of Electronica y Sistemas, Universidad de A Coru˜na, Campus de Elvi˜na, 15071 A Coru˜na, Spain Received 4 February 2002; received in revised form 12 August 2002
Abstract In multicarrier code-division multiple access systems the linear distortion introduced by the channel causes a perturbation in the received user codes that produces multiple access interference. This paper introduces two new techniques to combat this distortion that are based on subspace decomposition. The 6rst one consists of obtaining an approximation of the received desired user code to improve the output signal to interference noise ratio in blind linear receivers. This approximation is obtained by means of the projection of the known transmitted desired user code onto the signal subspace. The second one is a blind channel estimation method that uses the orthogonality property between the received user codes and the noise subspace. We demonstrate that the performance of channel estimation methods based on this property are severely limited by the number of simultaneous users. Nevertheless, it is shown how this limitation is overcome if user codes spanning several transmitted symbols (long codes) are used. ? 2002 Elsevier Science B.V. All rights reserved. Keywords: Multicarrier code-division multiple access systems; Blind channel estimation; Long spreading sequences; Perturbation analysis; Cramer–Rao bound
1. Introduction Multicarrier (MC) transmission methods combined with code-division multiple access (CDMA) communication systems have been proposed as an e=cient technique to combat multipath propagation and have gained an increased interest during the last years [7,3,17,4]. MC modulation and demodulation can be easily implemented using inverse fast Fourier transform (IFFT) and fast Fourier transform (FFT) algorithms, respectively. Depending on the ∗ Corresponding author. Tel.: +34-981-167150; fax: +34-981167160. E-mail address: [email protected] (C.J. Escudero).
spreading, frequency mapping and detection strategy, diDerent techniques based on this combination can be considered: MC-CDMA, MC direct sequence CDMA (MC-DS-CDMA) and multitone CDMA (MT-CDMA). This paper considers MC-CDMA schemes since they outperform the other techniques [4]. In these schemes each user is assigned to a unique identi6cation code sequence and the transmitted signal is split in diDerent subcarriers. It is assumed that the subcarrier bandwidth is smaller than the channel coherence bandwidth and, therefore, exhibits Eat fading only. As a consequence, MC-CDMA systems do not suDer from inter-symbol interference (ISI). However, the eDects of dispersive channels appear as random distortions in the amplitude and phase
0165-1684/03/$ - see front matter ? 2002 Elsevier Science B.V. All rights reserved. PII: S 0 1 6 5 - 1 6 8 4 ( 0 2 ) 0 0 4 1 9 - X
340
C. J. Escudero et al. / Signal Processing 83 (2003) 339 – 357
Nomenclature vni (k) Vni (m) ci (k) x n (k) vn;i j (k) cji (k) 2
i 2 r
kth chip of nth symbol corresponding to the ith user IDFT of vni (k) kth chip of ith user code received kth chip of the nth symbol period kth chip for the jth part of the nth block of symbols for the ith user kth chip of the jth part of the ith user long code ith user signal power noise power
xn
vector of the chips corresponding to the nth symbol period u(k : l) vector whose elements are the ones of u between the indexes k and l (both included)
of each subcarrier. This causes a perturbation in the received user codes and introduces multiple access interference (MAI) due to interfering users transmitting simultaneously with the desired one. In order to implement a multiuser detector and to reduce MAI it is necessary to characterize, implicitly or explicitly, the channel parameters. This paper introduces two new techniques to reduce MAI: the 6rst one obtains an approximation of the received user code to improve the output signal to interference noise ratio (SINR) and the second one is a blind channel estimation technique. Both techniques are based on a subspace decomposition [11] of the received autocorrelation matrix. We will derive particular algorithms for the new methods to identify the received code and channel parameters, respectively. These algorithms can be applied for both uplink and downlink scenarios. Because of the time-varying nature of mobile radio channels, adaptive 6ltering techniques have to be used in order to suppress the MAI. Conventional adaptive 6ltering techniques are based on the minimum mean square error (MMSE) criterion [10], where it is necessary to transmit a training sequence to set up the 6lter coe=cients. To avoid the transmission of
1i diag(a) F
zero vector but the ith element equal to one diagonal matrix whose elements are the elements of vector a DFT matrix
T
transpose operator transpose conjugate operator ∗ conjugate operator ∗ discrete convolution · ceiling operator · Eoor operator det(·) determinant of a square matrix Trace{·} trace of a square matrix H
training sequences, blind techniques based on the linearly constrained minimum variance (LCMV) criterion [8] have been proposed for MAI cancellation in DS-CDMA. However, these techniques are extremely sensitive to inaccuracies in the acquisition of the desired user timing and spreading code. The 6rst technique introduced in this paper is a new method to improve the performance of the blind LCMV criterion in a MC-CDMA system where we cannot obtain a perfect knowledge of the desired user code. This technique tries to approximate the received code by projecting the transmitted one onto the signal– interference subspace of the input autocorrelation matrix. The second estimation technique introduced in the paper takes into account the orthogonality property between the signal and noise subspaces to obtain an estimation of the desired channel impulse response. This technique resembles the subspace channel estimation methods already proposed for single carrier systems [11,16,9]. Nevertheless, when this technique is applied to MC-CDMA systems, it is demonstrated that, in order to obtain a good estimation of the channel, system capacity (i.e., the number of simultaneous users) is a severe constrain. This motivates the
C. J. Escudero et al. / Signal Processing 83 (2003) 339 – 357
consideration of a new structure that obtains channel estimations that are independent of the number of users and will be only limited by the number of MC-CDMA carriers. The new structure considers user codes larger than the user symbols period. In order to show the improvement achieved with the new structure, we obtain the Cramer–Rao bound expression to calculate the lower bound of the estimated parameters variance for any estimator in a MC-CDMA environment using long codes. An adaptive algorithm to obtain the channel estimation is also introduced and, using perturbation techniques, a theoretical approximation of the estimated mean-square error (MSE) is derived. Finally, the structure that can be used to detect transmitted symbols is presented. The paper is organized as follows. Section 2 introduces the basic signal model of a MC-CDMA system. Section 3 presents a new algorithm that is based on the signal subspace and it is demonstrated that it improves the output SINR. Algorithm performance is illustrated by means of some computer simulations. Section 4 introduces an alternative channel estimation technique based on the noise subspace and describes its limitations when it is used in MC-CDMA systems. We also show how to overcome these limitations by using spreading codes larger than the symbol period. The theoretical mean-squared error is obtained for this algorithm and the Cramer–Rao bound for the new signal model is also achieved. Numerical results, at this section, illustrates the performance of the algorithm and the validity of the theoretical analysis. Finally, Section 5 is devoted to the conclusions.
341
obtain the following multicarrier signal: L−1
(2=L)km 1 i vn (k)ej : L
Vni (m) = IDFT[vni (k)] =
(2)
k=0
This signal is transmitted through a dispersive channel with an impulse response hi (m); m = 0; : : : ; M − 1. At the receiver the observed signal is a superposition of the signals corresponding to N users plus an additive white Gaussian noise (AWGN). We will consider the uplink because the down link can be viewed as an special case where all user’s signals are received through the same channel. Therefore, the received signal for the nth symbol is the following: N
Xn (m) =
Vni (m) ∗ hi (m) + rn (m);
(3)
i=1
where ∗ denotes discrete convolution and rn (m) represents a white noise sequence. To recover the transmitted symbols, the receiver applies a L-DFT (discrete Fourier transform) to the received signal (3). Assuming perfect synchronization and a su=ciently large guard time between symbols, the resultant signal is x n (k) = DFT[Xn (m)] =
N
vni (k)H i (k) + n (k)
i=1
=
N
sni ci (k)H i (k) + n (k);
i=1
k = 0; : : : ; L − 1;
(4)
i
i
where H (k) and n (k) are the DFTs of h (m) and rn (m), respectively. Rewriting (4) in vector notation we obtain
2. Signal model
N
Let us consider a discrete-time baseband equivalent model of a synchronous MC-CDMA system with N users using L-chip signature codes [7]. The kth chip corresponding to the nth symbol transmitted by the ith user is given by
xn = [x n (0); : : : ; x n (L − 1)] =
vni (k) = sni ci (k);
where Ci is a diagonal matrix whose elements are the L chips of the code corresponding to the ith user, Hi = [H i (0); : : : ; H i (L − 1)]T and n = [n (0); : : : ; n (L − 1)]T . To obtain (5) we have used the relationship Hi = Fhi where F is the DFT matrix and hi = [hi (0); : : : ; hi (M − 1)]T . Note that (5)
k = 0; : : : ; L − 1;
n = 0; 1; 2; : : : ;
(1)
where ci (k) is the kth chip of the ith user code. In a MC-CDMA system the modulator computes the L-IDFT (inverse discrete Fourier transform) of (1) to
T
sni Ci Hi + n
i=1
=
N i=1
sni Ci Fhi
+ n =
N
sni c˜i + n ;
(5)
i=1
342
C. J. Escudero et al. / Signal Processing 83 (2003) 339 – 357
can be interpreted as a CDMA signal where the user codes are perturbed versions of the transmitted one, i.e., c˜i = Ci Fhi . Assuming statistical independence between users and noise, the autocorrelation matrix of the observations vector (5) can be decomposed as Rx = E[xn xnH ] =
N
∗
H
c˜i E[sni sni ]c˜i + E[n H n]
=
i2 i iH
r2
c˜ c˜ + I;
(6)
where E[ · ] is the expectation operator, superindices ∗ and H denote conjugate and conjugate transpose, 2 2 respectively, i and r are the ith user signal and noise power, respectively, and I is the identity matrix. Let us consider the eigendecomposition of Rx . There are L eigenvalues that we sort as 0 ¿ 1 ¿ · · · ¿ L−1 . It is well known [11] that the eigenvectors associated to the N most significants eigenvalues (ul ; l = 0; : : : ; N − 1) span the signal subspace where the perturbed (received) user codes, c˜i , lie. The remaining L − N eigenvectors (ul ; l = N; : : : ; L − 1) span the noise (orthogonal) subspace and their associated eigenvalues are equal 2 to the noise power, i.e., N = · · · = L−1 = r . 3. Signal subspace technique To produce the receiver output, the observations vector (5) can be processed by a linear combiner w as follows: (7)
It is well known that the SINR at the output is given by 2
d |wH c˜d |2 SINR = H ; w Ri+r w where superindex Ri+r =
N
2
d H
(8)
represents the desired user and 2
2
i c˜i c˜i + r I = Rx − d c˜d c˜d
H
w
subject to wH c˜d = 1:
(10)
It is well known [8] that the solution to this optimization problem is 2
i=1
yn = wH xn :
min E[|yn |2 ]
H
−1 d wLCMV = Rx−1 c˜d ⇒ SINR opt = d c˜d Ri+r c˜ ;
i=1 N
When considering blind adaptive MAI cancellation, the weight vector w in (7) is chosen according to the LCMV criterion [8]
(9)
i=1 i=d
is the autocorrelation matrix associated to the interferer users plus noise.
(11)
where is an arbitrary constant and SINR opt is the maximum SINR attainable. However, the weight vector (11) can only be obtained under ideal conditions where we have a perfect knowledge of the received code c˜d . In practical situations, code acquisition systems [6] provide an estimation of the transmitted code cd =[cd (0); : : : ; cd (L−1)]T . In this case, the weight vector is given by wLCMV = Rx−1 cd
(12)
that is no longer optimum since cd = c˜d . The resulting SINR is extremely low because the receiver interprets the desired signal as an interference and cancels it from the output. In order to show the degradation of the SINR when using (12), let us consider a simple environment with a single user (the desired one) and Gaussian noise. In this case the autocorrelation matrix for the interferences plus noise is reduced to a diagonal matrix 2 Ri+r = r I. Therefore, substituting (12) in (8) and using the matrix inversion lemma, the SINR can be rewritten as SINR opt cos2 ; (13) SINR = 1 + SINR opt (SINR opt + 2) sin2 where represents the angle between the vectors c˜d and cd . Note that depends on the perturbation introduced by the channel. In the absence of channel distortion, = 0, and equation in (13) corresponds with the SINR opt . However, when = 0 the SINR decreases when SINR opt increases. This degradation of the SINR can be avoided if we use in (12) any vector parallel to c˜d instead of cd . One of these vectors can be obtained, even without the knowledge of the perturbed code, by projecting the transmitted code cd onto the signal subspace since this subspace is the line spanned by the vector c˜d . Next subsection presents the projection technique for a general environment.
C. J. Escudero et al. / Signal Processing 83 (2003) 339 – 357
3.1. Projection technique When we consider a general environment with N users and Gaussian noise it is easy to determine the signal–interference subspace but not the desired signal subspace alone. Nevertheless, projecting cd onto the signal–interference subspace results in a better performance as will be shown in the sequel. This projection does not guarantee a vector parallel to the true desired user code, c˜d , but we still obtain an approximation close to it. This technique has already been used in [5] for adaptive beamforming with uncalibrated arrays. Indeed, the transmitted user code can be decomposed as the sum cd = csd + crd ; csd
(14) crd
where and are two orthogonal components that lie in the signal–interference and noise subspace, respectively. These two subspaces can be determined from the eigenvectors and eigenvalues of the autocorrelation matrix Rx [11]. Substituting (14) and (12) in (8), the SINR for the LCMV technique can be expressed as follows: H −1 H csd Rx−1 csd + crd Rx−1 crd SINR LCMV = −1 : d2 |csdH Rx−1 c˜d |2 (15) This expression can be interpreted as a generalization of (13) to a N users environment and, in general, yields to values which are considerably smaller than the optimum SINR. However, projecting the transmitted user code onto the signal plus interference subspace is equivalent to eliminating the component crd and the resulting output SINR is given by −1 H csd Rx−1 csd SINR proj = −1 : (16) d2 |csdH Rx−1 c˜d |2 It is straightforward to see that this value is higher than (15). This indicates that an improvement in performance is obtained when using the projected vector csd instead of cd . 3.2. Algorithm It is well known that the signal–interference subspace is spanned by the eigenvectors corresponding to the N most signi6cant eigenvalues of the autocorrelation matrix Rx [11]. Considering the L × N
343
unitary matrix U containing these signi6cant eigenvectors, the signal projection matrix is given by P = U(UH U)−1 UH = UUH . Therefore, the projection of cd onto the signal-interference subspace is obtained by multiplying the projection matrix with the code, i.e., Pcd = csd . If we use this projection into the linear constraint of (10), the weight vector (12) can be written as follows: wproj = Rx−1 Pcd = Rx−1 csd :
(17)
The resulting linear receiver outperforms the conventional LCMV receiver since it achieves the output SINR given by (16). In a practical implementation, the autocorrelation matrix is unknown and Ns −1 xn xnH , where Ns is estimated as Rˆ x = (1=Ns ) n=0 is the number of received symbols used to estimate the matrix. This 6nite number of data observations produce errors in the estimation that reduce the SINR similarly to situations where we do not have a perfect knowledge of the received codes. Estimation errors can be assumed as an error in the knowledge of the received codes. Therefore, for a 6xed number of data observations, the projection technique will obtain a higher output SINR. As a result, the proposed projection algorithm exhibits faster convergence even when there is no channel distortion or we have a perfect knowledge of the received code. 3.3. Numerical results To illustrate the eDectiveness of the signal subspace technique, we have carried out several computer simulations. This section presents the results of these simulations using the following common environment: channel length M = 4 and 100 averaged realizations using diDerent random channels whose coe=cients are drawn from independent Gaussian complex random variables. The remaining parameters, i.e., number of users, code length, type of code, SNR and number of subcarriers will be speci6ed in each simulation. All simulations consider an uplink scenario, i.e., the channel is assumed to be diDerent for each user. In the 6rst simulation we consider random codes with length L = 8. Here we compare the performance of the projection technique with respect to the optimum (in the sense of perfect knowledge of the desired user code) and the LCMV technique. Fig. 2
C. J. Escudero et al. / Signal Processing 83 (2003) 339 – 357
ci,0 (0)
i snp
Symbols
Serial/Parallel Converter
ci,0 (L-1) ci,j (0)
i snp+j
i vn,0 (0)
i vn,0 (L-1)
i vn,j (0)
vn,ji (L-1) ci,j (L-1) ci,p-1 (0)
i snp+p-1
ci,p-1 (L-1)
Parallel/Serial Converter
344
MC-CDMA signal
i vn,p-1 (0)
i vn,p-1 (L-1)
Fig. 1. Block diagram of the ith user MC-CDMA modulator for transmitting the nth block of p symbols using long codes of length pL, where L is the spreading gain.
12
10
8
6
SINR (dB)
4 Optimum With Projection LCMV
2
0
-2
-4
-6
-8
0
5
10 SNR (dB)
Fig. 2. SINR vs. SNR for optimum, projection and LCMV techniques.
15
C. J. Escudero et al. / Signal Processing 83 (2003) 339 – 357
345
0
10
10 -1
-2
SER
10
10 -3
10
-4
10
-5
Optimum With Projection LCMV
0
5
10
15
SNR (dB)
Fig. 3. SER vs. SNR for optimum, projection and LCMV techniques.
plots the average of the output SINR for diDerent channel realizations with respect to the input signal to noise ratio (SNR) considering three users transmitting with the same input SNR. Note that the output SINR corresponding to the LCMV reduces as the optimum increases, whereas the projection technique improves the performance by reducing the degradation of the SINR with respect to the optimum value. To illustrate the performance of the proposed projection algorithm in terms of the symbol error rate (SER), Fig. 3 plots it with respect to the input SNR for the same three users environment as before. Note that the curve corresponding to the projection method decreases with the input SNR and is close to the SER obtained when we have a perfect knowledge of the received code. Fig. 4 represents the output SINR time evolution of the projection algorithm. In this simulation we have considered the same environment as before but with the users transmitting with an input SNR = 12 dB.
Note that the projection technique presents an asymptotic convergence better than the LCMV technique. The diDerence between the convergence value of the projection technique and the optimum (in the sense of knowing the code) value depends on the channel features. Moreover, as we have said in Section 3.2 the proposed detector exhibits faster convergence than the one considering knowledge of channel impulse response (curve labelled as optimum) since the estimated channel takes into account errors in the estimated autocorrelation matrix when a small amount of symbols are available. 4. Noise subspace technique As we have seen at the signal model section, the perturbed user codes lie in the signal subspace and, therefore, are orthogonal to the noise subspace. In this case, we will try to perform an estimation of the desired user channel impulse response vector, hd .
346
C. J. Escudero et al. / Signal Processing 83 (2003) 339 – 357 15
10
SINR (dB)
5
0 Optimum With Projection LCMV -5
-10 0
100
200
300
400
500 Symbols
600
700
800
900
1000
Fig. 4. SINR time evolution for projection and LCMV techniques.
Considering the eigendecomposition of the autocorrelation matrix Rx , we have seen that the eigenvectors associated to the N most signi6cants eigenvalues (ul , l = 0; : : : ; N − 1) span the signal subspace where the perturbed (received) user codes, c˜i , lie and the remaining L − N eigenvectors (ul , l = N; : : : ; L − 1) span the noise (orthogonal) subspace. Therefore, this property can be used to obtain the channel impulse response vector for the dth user, hd , by means of the following system of equations: H
c˜d ul = (Cd Fhd )H ul = 0;
l = N; : : : ; L − 1:
(18)
Recall that this system of equations has M unknowns and L − N equations. It will be solvable if and only if the number of equations is greater or equal than the number of unknowns, that is M 6 L − N ⇒ N 6 L − M:
(19)
Note that we have obtained a limit for the system capacity that depends on the number of carriers and the channel length. This means that if we want to improve the capacity it is necessary to increase the number of
carriers, L. However, in many practical applications this is not possible, or desirable, and we have to consider another solution to increase the system capacity. 4.1. Utilization of long codes In the previous sections we have considered that the codes length is equal to the spreading gain (i.e., to the number of carriers), thus 6nding that the system capacity is limited by the number of carriers. In this section, we are going to explore the idea of increasing the system capacity by means of increasing the codes length without modifying the number of carriers, i.e., we will consider codes with a length larger than the spreading gain. To simplify the notation we consider pL-chip signature codes, where p ¿ 1 is an integer number (i.e., the code length is a multiple of the spreading gain). This means that the user codes will be repeated after p symbols. The modulator for a block of p consecutive symbols is represented in Fig. 1 where the jth symbol of this block is modulated by the jth part of
C. J. Escudero et al. / Signal Processing 83 (2003) 339 – 357
the user code to obtain the following CDMA signal: i cji (k); vn;i j (k) = snp+j
j = 0; : : : ; p − 1;
k = 0; : : : ; L − 1;
N
vn;i j (k)H i (k) + n; j (k)
i=1
=
N
i snp+j cji (k)H i (k) + n; j (k)
i=1
k = 0; : : : ; L − 1;
is a pL×p matrix. The pL×pL autocorrelation matrix of (22) can be expressed as follows: R = E[x(n)xH (n)]
n = 0; 1; 2; : : : ;
where cji (k) is the kth chip of the code associated to the jth symbol of a block for the ith user. This equation is equal to (1) for p = 1. As in Section 2, we can obtain the received signal after DFT as x n; j (k) =
(20)
where n; j (k) is the DFT of the noise sequence superposed to the jth symbol. A vector of L chips for the jth symbol of a block is
=
N
i snp+j Cij Fhi + n; j ; j = 0; : : : ; p − 1; (21)
where n; j = [n; j (0); : : : ; n; j (L − 1)]T and Cij is a L × L diagonal matrix whose elements are the L chips of the code corresponding to the jth sequence of a block for the ith user. Arranging the p consecutive vectors xn; j , corresponding to a single block, into a single vector of length pL, we obtain x(n) = [xn;T 0 ; : : : ; xn;T p−1 ]T N
Ci Hi si (np) + n ;
(22)
i=1
where Ci is a diagonal matrix containing the pL chips of the code associated to the ith user, i i si (np) = [snp ; : : : ; snp+p−1 ]T , n = [Tn; 0 ; : : : ; Tn; p−1 ]T and i 0 ··· 0 Fh .. .. i 0 . Fh . (23) Hi = .. . . .. .. . 0 0
···
0
Fhi
2
H
H
2
i Ci Hi Hi Ci + r I:
(24)
Let us consider the eigendecomposition of R and sort its eigenvalues as 0 ¿ 1 ¿ · · · ¿ pL−1 . Note that matrix Ci Hi has p independent columns and, therefore, the eigenvectors associated to the pN most signi6cants eigenvalues (ul , l = 0; : : : ; pN − 1) span the signal subspace where the diDerent user signals lie. The remaining p(L − N ) eigenvectors (ul , l = pN; : : : ; pL − 1) span the noise subspace. Taking into account the orthogonality between the signal and the noise subspace, we obtain the following system of equations for the desired user, d: H
H
(Cd Hd )H ul = 0 ⇒ hd FH Cdj ul; j = 0; l = pN; : : : ; pL − 1; j = 0; : : : ; p − 1;
i=1
=
N i=0
xn; j = [x n; j (0); : : : ; x n; j (L − 1)]T =
347
(25)
where ul; j = ul (jL + 1 : jL + L − 1) and the notation u(k : l) represents a vector containing the elements of u between the indexes k and l (both included). Recall that this system of equations still has M unknowns but p2 (L − N ) equations. It will be solvable if and only if the number of equations is greater or equal than the number of unknowns, that is, 2 p L−M 2 ; (26) p (L − N ) ¿ M ⇒ N 6 p2 where · is the Eoor operation which rounds down to the previous integer value. Therefore, the system capacity can be improved by increasing the p value, i.e., the code length. Similarly, there exists a relationship for the parameter p when there is a 6xed number of users
M p¿ ; (27) L−N where · is the ceiling operator. Finally, in order to simplify the subsequent analysis, the system of equations (25) can be rewritten as follows: H
H
hd FjH Cd ul = 0; j = 0; : : : ; p − 1;
l = pN; : : : ; pL − 1; (28)
348
where Fj = [0M ×L ; : : : ; 0M ×L , 0M ×L ]T is a pL × M matrix.
C. J. Escudero et al. / Signal Processing 83 (2003) 339 – 357
FT
, 0M ×L ; : : : ;
( j+1)th position
4.2. Cramer–Rao bound In order to show the eDect of considering long codes in a MC-CDMA environment, this section computes the channel estimation Cramer–Rao bound (CRB) for any value of p. The CRB is useful in the sense of that it allows to determine the best-possible performance (smallest variance) of any parameter estimator. Estimators whose variance is close or equal to this bound can be said to have an optimal behaviour. In our case, it is possible to compute the CRB and, therefore, determine which is the best-possible performance of the channel impulse response estimation. Moreover, by means of this CRB, we shall be able to know if the best-possible performance obtained using spreading codes equal to the symbol period would be improved by using long codes. In a few words, CRB gives a lower bound on the error covariance matrix for a parameter vector to estimate (in our case, the channel impulse response vector) [13], and it provides a benchmark against which algorithms performance can be compared. The CRB of a parameter estimator, in our case of the channel impulse response, is determined by the fact that the error covariance matrix for this parameter vector has to be greater or equal than the inverse of the Fisher information matrix (FIM). FIM is the covariance matrix of the score function which is the gradient of the log-likelihood function. The standard way of computing the CRB is explained in [9], where the FIM is obtained considering as parameters to estimate 2 2 r , i , sni and hi (m) for i = 1; : : : ; N , m = 1; : : : ; M and n = 0; : : : ; T − 1 where T is the number of observed vectors (22). Nevertheless, the complexity of this technique depends on many parameters and it is rather cumbersome. To simplify the CRB analysis. we only consider the components of the desired user channel, as parameters to estimate, that is, the remainder parameters are considered as nuisance [13]. This approximation will provide a lower variance bound, but still good to compare the performance of the new algorithm.
If we consider a 6xed number, T , of independent observed vectors, we can assume a multivariate normal model, x = (x(0); x(2); : : : ; x(T − 1)), x(n) : CN[x(n); R] where x(n) = E[x(n)]. The joint distribution of the complex random sample has density f (x) = −pLT (det(R))−T T −1 H −1 ×exp − (x(n) − x(n)) R (x(n) − x(n)) ; n=0
(29)
where det(·) represents the determinant of a square matrix. We can generalize the results in [13, pp. 235– 237] to the complex case by considering a vector T T of real parameters = [real(hd ); imag(hd )]T where real(·) and imag(·) are the real and imaginary parts of a complex vector. The FIM associated to is T @ @ ln f (x) ln f (x) : (30) J = E @ @ To compute the FIM it is easier to obtain it with respect T H to the vector h = [hd ; hd ]T and to relate it with (30) as in [1] J =
I
I
jI
−jI
Jhh
I
I
jI
−jI
T :
(31)
By doing a similar analysis to [14,1] we can obtain Jh d h d ∗ Jh d h d ; (32) Jhh = Jh d ∗ h d Jh d ∗ h d ∗ where the elements (i; j) of the four submatrices in (32) are @R @R ; Jhd hd (i; j) = T Trace R−1 d R−1 d @h (i) @h (j) @R @R −1 Jhd hd∗ (i; j) = T Trace R R @hd (i) @hd∗ (j) T −1 @xH (n) −1 @x(n) + Trace R ; @hd∗ (j) @hd (i)
−1
n=0
C. J. Escudero et al. / Signal Processing 83 (2003) 339 – 357
@R @R −1 R J (i; j) = T Trace R @hd∗ (i) @hd (j) T −1 @xH (n) @x(n) + Trace R−1 d( j) ; @hd∗ (i) @h n=0 @R @R d∗ d∗ −1 −1 Jh h (i; j)= T Trace R R ; @hd∗ (i) @hd∗ (j) (33)
hd∗ hd
−1
where Trace {·} represents the trace of a square matrix. It is interesting to note that, under our assumptions, N x(n) = i=1 Ci Hi si (np). To compute the FIM it is necessary to specify the value of the partial derivatives in (33), i.e., 2 H H @R = d C d B i H d Cd ; @hd (i)
(34)
H H @xH (n) = (sd (np))∗ Bi Cd ; ∗ d @h (i)
0
···
0
In order to solve the system of equations (28), we can consider the following equivalent system: H
H
H
j=0 l=pN
de6ned as follows: ;
H
p−1 pL−1
= arg min hd hd 2 =1
where 1i = [0; : : : ; 1 ; 0; : : : ; 0]T is a vector with the ith element equal to one and zero the others. As we saw before, by using blind estimation techniques using second-order statistics, channel impulse response can be determined up to a complex unknown factor and, therefore, FIM is rank de6cient. In order to avoid the ambiguity of the unknown factor and obtain the CRB, we can assume that one non-zero complex channel parameter, hd (i) remains constant and, as a result, its corresponding derivative vanishes. Empirically, it was observed that assuming constant the
H
p−1
= arg min hd hd 2 =1
H
FjH Cd ul ulH Cd Fj hd H
FjH Cd UU Cd Fj hd H
j=0 H
= arg min hd Qd hd ; hd 2 =1
j=0 l=pN
(35)
(36)
for l=pN; : : : ; pL−1 and j=0; : : : ; p−1. On the other hand, the solution to (36) can be found by solving the following minimization problem: p−1 pL−1 H H hˆd = arg min hd FjH Cd ul ulH Cd Fj hd hd 2 =1
F1i
ith element
4.3. Channel identi7cation algorithm
H
@x(n) = Cd Bi sd (np); @hd (i)
where B is a pL × p matrix 0 ::: 0 F1i .. .. 0 . F1i . i B = .. .. .. . . . 0
channel parameter with largest energy [1] the minimum CRB is obtained. The eDect of 6xing this complex unknown parameter [9,1] is equivalent to deleting the rows and columns corresponding with its real and imaginary parts from the FIM before computing its inverse to obtain the CRB. Section 4.6 will present numerical results of the CRB for MC-CDMA systems using long codes with diDerent values of p. We will see that systems satisfying condition (27) yield the lowest CRB.
hd FjH Cd ul 2 = hd FjH Cd ul ulH Cd Fj hd = 0
2 H H @R = d Cd Hd B i Cd ; ∗ d @h (i)
i
349
(37)
where the solution hˆd is an estimation of the channel impulse response vector, U is a pL × p(L − N ) matrix whose columns are the eigenvectors associated to the noise subspace (i.e., ul ; l = pN; : : : ; pL − 1) and p−1 H Qd = j=0 FjH Cd UUH Cd Fj . The solution can be obtained by the least-squares method and it corresponds to the eigenvector of Qd associated to its minimum eigenvalue [15]. Note that we have included the constraint hd 2 = 1 in order to avoid the null solution. In practice, we do not know a priori the autocorrelation matrix (24) but we can estimate it from the
350
C. J. Escudero et al. / Signal Processing 83 (2003) 339 – 357
sampled averaged matrix as 1 Rˆ = Ns
N s −1
x(n)xH (n);
and (38)
n=0
where Ns is the number of received pL chip blocks used to obtain the estimation. Finally, note that when using second-order statistics, the channel impulse response can be obtained up to a complex constant that should be compensated in order to remove undesirable ambiguities and analyze the algorithm performance. Towards this aim, we normalize the estimation of the impulse response vector with respect to the ith element as follows: hˆdnormalized = #hˆd ;
(39)
where # = hd (i)= hˆd (i) and hd (i) and hˆd (i) are the ith elements of the true and estimated channel impulse response vectors, respectively.
ˆ Ph −Q† PQh = −Q† (Qˆ − Q)h = −Q† Qh; (42) where Q† denotes the left pseudo-inverse of Q. The kth component of Ph is given by p−1 ˆ = −qkH Ph(k) −qkH Qh FjH CH Uˆ Uˆ H CFj h j=0 p−1 pL−1
=−
j=0
qkH FjH CH uˆl uˆH l CFj h l=pN scalar
scalar
p−1 pL−1
=−
j=0
H H H uˆH l CFj hqk Fj C uˆl l=pN
pL−1
=−
H H uˆH l CFQk C uˆl
l=pN
4.4. Mean square error analysis
H ˆ = −Trace{Uˆ H CFQH k C U}
In this section, we derive an analytical expression of the estimation MSE. 1 To clarify notation, let us denote hd = h; Qd = Q and Cd = C and consider the following identities: Qh = 0; hˆ = h + Ph;
(40)
Qˆ = Q + PQ: The perturbation terms Ph and PQ represent the estimation errors of h and Q, respectively. Thus, the channel estimation mean square error is given by E[Ph2 ]. The key point of our analysis is the utilization of a perturbation technique [12] that allows us to express Ph in terms of PQ. For a su=ciently large number of samples (Ns → ∞), Qˆ → Q, hˆ → h and Qˆ hˆ is approximately equal to the zero vector, i.e. Qˆ hˆ = (Q + PQ)(h + Ph) PQh + QPh ≈ 0; where we have consider negligible the second-order term, PQPh 0. Therefore, QPh −PQh 1
(41)
A brief description of this analysis was presented in [2].
H = −Trace{Uˆ Uˆ H CFQH k C };
(43)
where qk is the kth column of (Q† )H , Uˆ is the estimated matrix U, F = [F0 h; : : : ; Fp−1 h], Qk = [F0 qk ; : : : ; Fp−1 qk ] and Trace{·} represents the trace of a square matrix. Based on the results of [14, p. 1840, Eq. (4.11)], we can obtain the following identity: Uˆ Uˆ H CF −UUH PVVH CF;
(44)
where V is a pL × pN matrix whose columns are the eigenvectors associated to the signal subspace (i.e., ul ; l = 0; : : : ; pN − 1) and PV = Vˆ − V. Moreover, from Appendix A of [14, p. 1844, Eq. (A.2)] −1 ˆ Uˆ H PV UH RV) ;
(45) 2
2
where ) = diag(0 − r ; : : : ; pN −1 − r ) where diag(a) is a diagonal matrix whose elements are the elements of vector a. To remove the eDect of the unknown constant that we have in the estimation of the channel vector, we have to consider a normalization of the estimated vector channel. Normalization (39) is a function of the channel estimate. In order to facilitate the analysis of the mean square error approximation, we must seek a normalization independent of the error in the channel
C. J. Escudero et al. / Signal Processing 83 (2003) 339 – 357
estimate. Thus, we can approximate the normalization scalar in (39) as #= =
h(i) h(i) = ˆ h(i) − Ph(i) h(i) 1 Ph(i) 1− ; 1 + Ph(i)=h(i) h(i)
(46)
where we have used the fact 1=(1 + x) 1 − x for small x. Therefore, we can consider the following normalization: Ph(i) ˆ h−h Phnormalized = #hˆ − h = 1 − h(i) Ph(i) ˆ h h Ph − Ph(i) h(i) h(i) h1Ti Ph; (47) = I− h(i) = Ph −
4.5. Linear detector Once we have the channel estimation we can obtain the transmitted symbols corresponding to the jth sequence by using a linear combiner, i.e., yj = wjH xn; j ;
wj = %Rj−1 Cdj Fhd ;
−1 H H ˆ Ph(k) Trace{UUH RV) V CFQH kC }
SINR j =
=
1 SINR j ; p
(52)
j=0
d |wjH (Cdj Fhd )|2 wjH [Rj − d2 (Cdj Fhd hdH FH Cdj )]wj H
(53)
is the jth sequence output SINR.
−1 H H ˆ (ulH RV) V CFQH k C ul )
4.6. Numerical results ˆ l; k ; ulH Rg
(48)
l=pN H where gl; k = V)−1 VH CFQH k C ul . Finally, in Appendix A it is demonstrated that the MSE of the channel estimation algorithm is given by
E[Ph2 ] =
p−1
SINR =
2
l=pN pL−1
(51)
where % is an arbitrary constant and Rj is the L × L autocorrelation matrix for sequence xn; j . Note that Rj is a submatrix of R in (24) that considers the columns and rows between the jth and the (jL−1)th indices. In practical situations we use estimations of hd and Rj . The performance of the linear combiner will be evaluated in terms of the output SINR. We can obtain a compact expression of the SINR by using an averaged value over the p diDerent combiners, i.e.,
where
=
(50)
where wj is a coe=cients vector for the jth combiner. If we select this vector to maximize the output SINR, we obtain the following expression [8]:
where we have assumed the second-order term as negligible i.e., Ph(i)Ph 0. To include this normalization in (43) we can consider qk as the kth column of the matrix ((I − h1Ti =h(i))Q† )H . Combining (44) and (45) in (43), we obtain the following expression:
pL−1
351
2 M −1 r ˜ ˜H (Trace{UH UGH k CC Gk } Ns
k=0
2
+ r Trace{UH UGH k Gk });
(49)
where Gk = [g(pN ); k ; : : : ; g(pL−1); k ] and CQ = [ 1 C1 H1 ; : : : ; N CN HN ]. In order to obtain this expression we had to explore the fourth-order statistics of binary and Gaussian random variables.
To show the algorithm performance of this noise subspace technique, we show the results of some computer simulations using a common environment, i.e., channel length M = 4 and 100 averaged realizations using diDerent random channels whose coe=cients are drawn from independent Gaussian complex random variables. The remaining parameters, i.e., number of users, code length, type of code, SNR and number of subcarriers will be speci6ed in each simulation. All simulations consider an uplink scenario, i.e., the channel is assumed to be diDerent for each user. Next we evaluate the performance of the noise subspace technique in terms of the MSE, that it is de6ned as MSE = (hˆinormalized − hi )2 :
(54)
352
C. J. Escudero et al. / Signal Processing 83 (2003) 339 – 357 4
10
N=8 N=7 N=6 3
10
2
MSE
10
1
10
0
10
10
10
-1
-2
0
50
100
150
200 250 Symbols
300
350
400
450
500
Fig. 5. Mean square error time evolution for diDerent values of N .
Fig. 5 presents the averaged MSE time evolution for diDerent values of N . An environment using short (L=10) random codes has been considered. The users are received with a SNR = 10 dB. According to (19) a value for the number of users N = 7 is acceptable but it can be seen that for higher values, it is achieved a poor channel estimation. Fig. 6 presents the averaged MSE time evolution using long codes with diDerent values of p. We have considered the same environment as in the previous simulation but with 8 users. In this case relationship (26) reduces to p ¿ 2. It can be seen that for values below this condition (i.e., p=1), the channel estimation does not achieve a good performance. To show the eDectiveness of the proposed technique, Fig. 7 plots the MSE time evolution of the simulation and the Cramer–Rao bound (CRB) for an environment with 8 users, L=12, SNR =12 dB, random codes and p=2 (i.e., satisfying condition (27)). More-
over, it is plotted a third curve with the theoretical MSE time evolution (49). It can be seen that theoretical and simulated curves are very similar showing the validity of expression (49). On the other hand, MSE is very close to the CRB even for a small amount of received symbols. To show the bene6ts of using long codes in terms of the MSE, Fig. 8 shows the averaged CRB with the same length channel and codes as before for N = 5, L = 8, SNR = 20 dB and three diDerent values of p. It can be seen that for values satisfying (27) (i.e. p = 2 and p = 3) the resulting CRB is almost the same. However, for values not satisfying this condition the estimator yields a worse limit. Finally, Fig. 9 shows the averaged simulated and theoretical MSE vs. the SNR for an environment of N =8 of synchronous users using random codes, L=12 subcarriers and p = 1. The curves are obtained after Ns = 300 symbols. We can see that both curves are very similar even for small values of SNR.
C. J. Escudero et al. / Signal Processing 83 (2003) 339 – 357
353
4
10
p=3 p=2 p=1 3
10
2
MSE
10
1
10
0
10
10 -1
10 -2
0
50
100
150
200
250 300 Symbols
350
400
450
500
Fig. 6. Mean square error time evolution for diDerent values of p. 2
10
Simulation Theoretical CRB 1
10
0
10
-1
10
-2
10
-3
MSE
10
10-4
10
-5
0
50
100
150
200
250
Symbols
Fig. 7. CRB, theoretical MSE and simulated MSE for the channel identi6cation algorithm, vs. number of symbols (p = 2).
354
C. J. Escudero et al. / Signal Processing 83 (2003) 339 – 357 0
10
p=3 p=2 p=1 -1
10
-2
MSE
10
10-3
-4
10
-5
10
0
50
100
150 Symbols
200
250
300
Fig. 8. CRB vs. number of symbols for diDerent values of p.
5. Conclusions In this paper we have investigated subspace decomposition techniques to combat the MAI due to received user code perturbations produced in uplink, or downlink, MC-CDMA systems. This code perturbation is typically caused by the channel linear distortion. We focused our attention on blind linear adaptive receivers based on the LCMV criterion. The performance of these receivers is very sensitive to inaccuracies in the acquisition of the received user code and typically result in an extremely low input SINR. Two diDerent subspace methods have been considered. The 6rst one consists of building the linear constraint in LCMV receivers using a projection of the transmitted desired user code onto the signal subspace. This yields to a new code closer to the received one and enables us to considerably increase the output SINR. This technique also accounts for the errors in the estimation of the autocorrelation matrix and thus
exhibits faster convergence. Computer simulations illustrate the performance of the projection technique. The second method is a blind channel estimation algorithm that exploits the orthogonality property between the received user codes and the noise subspace in a way similar to single user blind channel estimation techniques [11]. Nevertheless, we demonstrate that channel estimation techniques based on this property are severely limited by the number of carriers and simultaneous users when applied to MC-CDMA systems. This limitation, however, can be overcome by using codes larger than the symbol period (long codes). This is shown by calculating the Cramer–Rao bound of any channel estimator based on this subspace property in long codes scenarios. Finally, we have derived a novel channel estimation algorithm speci6cally devised for long codes and calculated an analytical expression of the channel estimation MSE. Computer simulations have been also carried out to illustrate the performance of the proposed channel estimation
C. J. Escudero et al. / Signal Processing 83 (2003) 339 – 357 10
355
-1
Simulation Theoretical
MSE
10 -2
10
-3
10
-4
10
-5
0
2
4
6
8
10 SNR (dB)
12
14
16
18
20
Fig. 9. Simulated and theoretical MSE vs. received users SNR.
technique and to corroborate the accuracy of the obtained MSE analytical expression.
H
N
Ci Hi si (np)H n
i=1
+
Acknowledgements
H
×(sj (np))H Hj Cj + N
H H n (si (np))H Hi Ci + n H n gl; k :
i=1
This work has been supported by Ministerio de Ciencia y TechnologRSa of Spain and FEDER funds from the European Union, under grant number TIC2001-0751-C04-01.
Appendix A Taking into account the estimated autocorrelation ˆ Eq. (48) can be expressed as matrix (38), R, pL−1 N N N s −1 1 i i i Ph(k) = ulH C H s (np) Ns i=pN n=0
i=1 j=1
(A.1) ulH (Ci Hi )=0
From (25) we know that and, therefore, it is straightforward to show that pL−1 Ns −1 1 ulH Ph(k) = Ns l=pN n=0 N H i H iH iH × n (s (np)) H C + n n gl; k : i=1
Therefore, the MSE is M −1 E[Ph2 ] = E[Ph(k)Ph∗ (k)] k=0
(A.2)
356
C. J. Escudero et al. / Signal Processing 83 (2003) 339 – 357
=
M −1 k=0
H H + ulH E[n H m ]ut gt; k E[m n ]gl; k
pL−1 pL−1 Ns −1 Ns −1 1 Ns2 i=pN t=pN
4
N N
4
= r ulH ut gt;Hk gl; k '(n − m);
H E[ulH n H m ut gt; k
i=1 j=1 H
H
×Ci Hi si (mp)(sj (np))H Hj Cj gl; k ] +
pL−1 pL−1 Ns −1 Ns −1 1 2 Ns t=pN
H H × E[ulH n H n gl; k gt; k m m ut ] ;
×
H
×Hj Cj gl; k 2
2
= ulH ( r '(n − m)I)ut gt;Hk Ci Hi ( i '(i − j) H
2
2
r i ulH ut gt;Hk Ci Hi Hi
H
4
× Ci gl; k + r ulH ut gt;Hk gl; k
=
2 M −1 r ˜ ˜H (Trace{UH UGH k CC Gk } Ns
k=0
H
H i i i j H = ulH E[n H m ]ut gt; k C H E[s (mp)(s (np)) ]
N
H
H i i i j H j j E[ulH n H m ut gt; k C H s (mp)(s (np)) H C gl; k ]
H
n=0
i=1
(A.3)
H
M −1 pL−1 pL−1 Ns −1 1 Ns2 t=pN k=0 l=pN
where we have used the fact that the third order moments of a Gaussian random variable are zero. Considering statistical independence between users and noise and the user symbols i.i.d., the 6rst expectation (A.3) is
(A.6)
note that we have used the fact ulH gl; k = H ulH (V)−1 VH CFQH k C ul ) = 0. Including (A.4) and (A.6) in (A.3), it is obtained E[Ph2 ] =
n=0 m=0
l=pN
4
= r ulH gl; k gt;Hk ut + r ulH '(n − m)ut gt;Hk '(m − n)gl; k
n=0 m=0
2
+ r Trace{UH UGH k Gk });
(A.7)
where Gk = [g(pN ); k ; : : : ; g(pL−1); k ; ] and C˜ = [ 1 C1 H1 ; : : : ; N CN HN ]. That is, we have obtained the expression in (49).
H
×'(n − m)I)Hj Cj gl; k 2
2
H
References
H
= i r ulH ut gt;Hk Ci Hi Hj Cj gl; k ×'(n − m)'(i − j);
(A.4)
where '(·) is the Kronecker function. It is well known that for any four complex Gaussian variables [12], i i = 1; 2; 3; 4, E[ 1 2∗ 3 4∗ ] = E[ 1 2∗ ]E[ 3 4∗ ] + E[ 1 4∗ ]E[ 3 2∗ ]: (A.5) In the second expectation of (A.3) our four complex variables are uiH n , gl;Hk n , gt;Hk m and utH m . The expectation can be expressed as H H E[ulH n H n gl; k gt; k m m ut ] H H = ulH E[n H n ]gl; k gt; k E[m m ]ut
[1] J.L. Bapat, Partially blind identi6cation of FIR channels for QAM signals, PhD. Dissertation, Department of Electrical Engineering, Pennsylvania State University, August 1996. International Editions, Singapore, 1993. [2] C.J. Escudero, D.I. Iglesia, M.F. Bugallo, L. Castedo, Analysis of a subspace channel estimation technique for multicarrier CDMA systems, Proceedings of the IEEE Signal Processing Workshop on Statistical Signal and Array Processing, Pennsylvania-USA, August 14 –16, 2000, pp. 184 –188. [3] K. Fazel, G.P. Fettweis, Multi-Carrier Spread-Spectrum, Kluwer Academic Publishers, Dordrecht, 1997. [4] K. Fazel, S. Kaiser, Multi-Carrier Spread-Spectrum & Related Topics, Kluwer Academic Publishers, Dordrecht, 2000. [5] D.D. Feldman, L.J. Gri=ths, A projection approach for robust adaptive beamforming, IEEE Trans. Signal Process. 42 (4) (April 1994) 867–876.
C. J. Escudero et al. / Signal Processing 83 (2003) 339 – 357 [6] S. Glisic, B. Vucetic, Spread Spectrum CDMA Systems for Wireless Communications, Artech House Publishers, Norwood, MA, 1997. [7] S. Hara, R. Prasad, Overview of multicarrier CDMA, IEEE Commun. Mag. 35 (December 1997) 126 –133. [8] M.L. Honig, U. Madhow, S. VerdRu, Blind adaptive multiuser detection, IEEE Trans. Inform. Theory 41 (4) (July 1995) 944–960. [9] H. Liu, G. Xu, Closed-form blind symbol estimation in digital communications, IEEE Trans. Signal Process. 43 (11) (November 1995) 2714–2723. [10] U. Madhow, M.L. Honig, MMSE interference suppression for direct-sequence spread-spectrum CDMA, IEEE Trans. Comm. 42 (12) (December 1994) 3178–3188. [11] E. Moulines, P. Duhamel, J.F. Cardoso, S. Mayrargue, Subspace methods for the blind identi6cation of multichannel FIR 6lters, IEEE Trans. Signal Process. 43 (2) (February 1995) 516–525.
357
[12] W. Qiu, Y. Hua, Performance analysis of the subspace method for blind channel identi6cation, Signal Processing 50 (1996) 71–81. [13] L.L. Scharf, Statistical Signal Processing: Detection, Estimation, and Time Series Analysis, Addison-Wesley, Reading, MA, 1991. [14] P. Stoica, T. SYoderstrYom, Statistical analysis and subspace rotation estimates of sinusoidal frequencies, IEEE Trans. Signal Process. 39 (8) (August 1991) 1836–1847. [15] G. Strang, Linear Algebra and its Applications, 3rd Edition, Harcourt, Brace and Jovanovich, New York, 1988. [16] M. Torlak, G. Xu, Blind multi-user channel estimation in asynchronous CDMA systems, IEEE Trans. Signal Process. 45 (January 1997) 137–147. [17] N. Yee, J.P. Linnartz, G. Fettweis, Multi-carrier CDMA in indoor wireless radio networks, Proceedings of the International Symposium on Personal, Indoor and Mobile Radio Communications, Yokohama, 1993, pp. 109 –113.