EURASIP Journal on Applied Signal Processing 2005:5, 709–717 c ○ 2005 Hindawi Publishing Corporation BER Performance for Downlink MC-CDMA Systems over Rician Fading Channels

We consider downlink multicarrier code-division multiple-access (MC-CDMA) systems using binary phase-shift keying (BPSK) modulation scheme and maximal ratio combining (MRC) in frequency-selective Rician fading channels. A time-domain method to obtain bit error rate (BER) by calculating moment generating function (MGF) of the decision variable for a tapped-delay-line channel model is proposed. This method does not require any assumption regarding the statistical or spectral distribution of multiple access interference (MAI), and it is also not necessary to assume that the fading encountered by the subcarriers is independent of each other. The analytical formula is also verified by simulations.


INTRODUCTION
MC-CDMA systems, based on the combination of code-division multiple-access (CDMA) and orthogonal frequency-division multiplexing (OFDM) techniques, were proposed in 1993 [1]. The multicarrier CDMA schemes can be categorized into two groups: MC-CDMA and MC-DS-CDMA [2]. Due to the attractive features like efficient frequency diversity and high bandwidth efficiency [3], MC-CDMA has received greater attention. Furthermore, MC-CDMA outperforms direct sequence CDMA (DS-CDMA) and MC-DS-CDMA in terms of BER performance over the downlink. Hence, MC-CDMA appears to be a suitable candidate for supporting multimedia services in mobile radio communications for the downlink.
Most of the previous papers [1,4,5], which investigated the performance of the MC-CDMA systems, assumed that the fading in different subcarriers is independent of each other, so that the variance of the interference can be approximated by using the central limit theorem. Nevertheless, the assumption is not guaranteed in practice, as the fading of the subcarriers is usually correlated due to insufficient frequency separation between the subcarriers. Also, the assumption of independent fading characteristic implies a frequency-selective fading channel at each subcarrier, since it requires sufficient independent paths uniformly distributed over the symbol duration [3], which contradicts the assumption of flat fading at each subcarrier. An exact error floor without taking into account the noise term is obtained in [6] under the assumption of exponential multipath intensity profile. In [7], a closed-form BER expression for a synchronous MC-CDMA in the uplink has been obtained assuming independent fading among subcarriers. In [8], a performance evaluation using characteristic method is proposed for MC-DS-CDMA systems. All the above-mentioned papers consider Rayleigh fading channels, and the results for other more general channel models are not available. In this paper, we propose a time-domain approach instead of usual frequency-domain approach to obtain the error performance of downlink MC-CDMA with maximal ratio combining (MRC) in correlated Rician fading channels. We can obtain an exact BER performance without having to make any assumption about the MAI distribution by calculating the moment generating function (MGF). It is also not necessary to assume that the fading of the subcarriers is independent of each other. In addition, it is not necessary to make assumption regarding the correlation property of the spreading sequence. The closed-form error performance may provide helpful insights relevant to designing the spreading sequences for MC-CDMA systems and may lead to an improved performance for MRC which has been shown to suffer severe MAI when the number of users is large.

Notation
(·) * , (·) T , and (·) H denote complex conjugate, transpose, and conjugate transpose operation, respectively. Vectors and matrices are represented in bold. E[·] denotes expectation. | · | stands for the norm of a complex variable or a vector. Re [·] and Im[·] correspond to the real and imaginary parts of a complex number. Matlab's notation FFT(·, n) and IFFT(·, n) denote n-point fast Fourier transform (FFT) and inverse FFT (IFFT). x · y stands for point multiplication between two vectors. ⊗ represents circular convolution. (m) N = mod(m, N). N(µ, σ 2 ) denotes normal distributed random variable with µ as mean and σ 2 as variance, respectively. det(·) denotes the determinant of a matrix.

Transmitter
We consider the system proposed in [1] and assume that there are Ku active users in a downlink. At the transmitter side, the spreading binary data stream is serial-to-parallel converted to N parallel substreams. All the data in N subcarriers are modulated in baseband using binary phase-shift keying (BPSK) by the means of the inverse discrete Fourier transform (IDFT). The resultant signals are then converted back into serial data. The guard interval is inserted between symbols to avoid intersymbol interference (ISI) caused by multipath fading, and finally the signal is transmitted after radio frequency (RF) upconversion. Let a k (q) and {c k,i } N i=1 denote the qth data bit and the spreading sequence of user k, respectively. The equivalent lowpass transmitted signal for user k can be expressed as [3] where E b is the power of the data bit and assumed to be the same for all users, p b (t) is a rectangular pulse defined in [0, T b ] with T b denoting the bit duration, and ∆ f = 1/T b is the minimum subcarrier separation.

Channel model
We consider a slowly varying frequency-selective Rician fading channel. The channel is modeled as a tapped delay line (TDL) having the following baseband equivalent impulse response [9]: where δ(·) is the Dirac delta function, L 1 is the number of resolvable paths of the channel, T c is the chip duration, b(l) is the lth path gain which is a complex Gaussian random process with zero mean (Rayleigh) or nonzero mean (Rician) and are mutually independent for different l. The Rician probability density function (PDF) is obtained as the PDF of |b(l)| = B 2 lx + B 2 ly , where B lx ∼ N(µ lx , σ 2 l ), B ly ∼ N(µ ly , σ 2 l ), and B lx , B ly are independent. The Rician K-factor is defined as the ratio of signal power in dominant component over the (local-mean) scattered power. We have Since only the nonzero paths need to be considered, we define {v l } L−1 l=0 as the propagation delay for the nonzero paths normalized by the sample interval, and σ 2 l as the variance, L is the number of the nonzero paths.
In the frequency domain, all subcarriers are assumed to experience flat but correlated fading. The channel gain for the ith subcarrier is It has been shown in [10] that The received signal can be written as where n(t) is additive white Gaussian noise (AWGN) having a double-sided power spectrum density of N 0 /2 for both real and imaginary components.

Receiver
At the receiver, the RF signal is first converted to baseband signal. After the portion of the signal corresponding to the prefix is removed, DFT is performed on the signal samples. A coherent correlation receiver with MRC is then used. As the focus of this paper is on the evaluation of BER, we assume perfect subcarrier synchronization with no frequency offset and no nonlinear distortion and also assume perfect subcarrier amplitude/phase estimation. Assuming user u is the desired user, the decision variable of the zeroth data bit is given by

PERFORMANCE ANALYSIS
We aim at obtaining a concise form of the decision variable for further processing. With the knowledge that there are usually only a few active taps compared to the number of subcarriers, the idea is to transform the calculation of the decision variable from the frequency domain to the time domain by applying discrete Parseval's theorem, so that the MGF of the decision variable can be derived. The exact BER is then obtained from the MGF by the Laplace inversion integral or numerical methods. Assuming the length of the cyclic prefix is longer than the length of the channel impulse response, so that there is no ISI. Here we discuss real spreading codes, which are commonly used in practice, however, extension to complexvalued codes is straightforward. The power loss due to the cyclic prefix is not considered in the analysis. The decision variable of the zeroth data bit in (6) with proper sampling time can be written as where , n i corresponds to the complex additive Gaussian noise at the ith subcarrier.
In order to normalize the AWGN noise, a factor equal to 2/N 0 is multiplied to (7). We obtain After normalization, theñ i in (8) is zero-mean, complex Gaussian random variable with a variance of σ 2 n = 1 for the real and imaginary components, respectively, and the ratio γ b = E b /N 0 represents the signal-to-noise ratio (SNR) per bit.
Since the noise is AWGN, the multiplication of the noise by spreading codes of the desired user will not change the distribution of the noise. The AWGN noise in the frequency domain is the DFT of the normalized AWGN noise in the time domain η(i) multiplied by 1/ √ N, where η(i) is actually n(t) sampled at the rate of 1/T c : Assume y is the corresponding signal of β u in the time domain given the value y = IFFT(β u , N), x is the signal corresponding to β u · h in the time domain, that is x = IFFT( β u · h , N). The multiplication of the DFT of two sequences is equivalent to the circular convolution of the two sequences in the time domain [11], the relationship between x and y can be thus expressed as Applying discrete Parseval's theorem [11], the summation in the frequency domain in (8) can be converted to the summation in the time domain. Omitting a factor equal to √ N, the decision variable in (8) can be rewritten in the time domain as From (10) and (11), the decision variable in the time domain can be obtained as For In order to obtain an expression of exact error performance, we calculate the Laplace transform of , that is, the MGF of the decision variable Z u (0). Usually it is not an easy task to calculate the MGF for a random vector since it is by definition the mean of an exponential function of the random variables involved. This generally requires the calculation of multidimensional integral and the knowledge of the joint probability density of the random variables. In this case, since the calculation is converted to the time domain, the channel coefficients {b(v i )} L−1 i=0 for different paths are statistically independent Gaussian random variables [9] and the exponential of the function turns out to be a noncentral Gaussian quadratic form. The MGF can be calculated as (see the appendix for details) where λ m , for m = 1, . . . , M, are distinct eigenvalues of A defined in (A.10) where M is the total number of distinct eigenvalues, w m is the multiplicity of the eigenvalue λ m , b m is given by (A.14). It is clear that

Closed-form expression
If an exact calculation of P(Z u (0) < 0) is sought, a common starting point is the Laplace inversion formula [9] where c is a small positive constant between zero and the smallest positive pole of Φ(s)/s. This complex contour integral might be evaluated exactly by calculating the residues of Φ(s)/s [12] over the poles in the right-hand side of the complex s-plane. Since Φ(s) contains essential singularities which make the residue evaluation very difficult, power series expansion is applied to solve the problem instead. Assuming that among the M distinct poles of Φ(s) (13), M 1 of them are positive. Following the method described in [13], we obtain a closed-form expression for (14) as where In (16), g (n) k (λ) is given by

Numerical method
Equation (15) gives a general closed-form BER expression for MC-CDMA systems over arbitrary Rician multipath fading channels. Since (15) is derived using the fast-convergent power series expansion of the MGF about its positive poles, it is expected to converge rapidly with respect to the index n. However, as shown in [13], the series (15) converges more slowly as the Rician K-factor increases, therefore it is more suitable for analysis of Rician fading channels with a relatively small K-factor, for example, K ≤ 7 dB. To achieve a more convenient solution to this problem, a different approach towards the exact calculation of (14) is used instead [14].
Since the left-hand side of (14) is a real quantity, we can write The change of variable ω = c Finally, using a Gauss-Chebyshev quadrature rule with an even number v of nodes, we have where τ k = tan((2k − 1)π/(2v)) and the error term E v vanishes as v → ∞. To achieve the desired degree of accuracy, it is practical increasing values of v and accepting the result when it does not change significantly. In general, 32 to 64 nodes are sufficient to obtain an accuracy better than 10 −5 . The value of c affects the number of nodes necessary to achieve a preassigned accuracy. The detailed discussion concerning the selection of c can be found in [14]. In our numerical results, we assume v = 64, and c is set equal to one half the smallest real part of the poles of Φ(s). However, even with a suboptimum choice of c, the value of v does not grow large enough as to become unmanageable. In the above section, a closed-form expression and a numerical method are given to obtain the error probability conditioned on the user data based on the MGF. For coherent detection, we have By averaging the conditioned BER (15) or (20) over {a k (0)} Ku k=1, k =u which are independent and identically distributed binary random variables, we obtain the BER of the user u: The average BER of all users is given by As the computational complexity to evaluate (22) increases exponentially with the number of users, the proposed approach is appropriate for systems with relatively small number of users, say Ku < 20. When the number of users is large, we may resort to traditional Gaussian approximation. In the following section, we propose a simple evaluation method to remedy this problem.

Gaussian approximation method for large number of users
From (7), the decision variable of the desired user can be rewritten as where The term D u is the desired output. I u corresponds to the MAI component. N u is the noise term contributed by AWGN. Due to the equal power, equally likely, antipodal data modulation a k (0) in (26), we apply the central limit theorem and approximate g u,i = Ku k=1, k =u a k (0)c u,i c k,i (i = 0, 1, . . . , N − 1) as a zero mean Gaussian random variable. Note that g u,i is correlated between subcarriers. The covariance matrix M u is an N × N matrix with each element given by M u (i1, i2) = Ku k=1, k =u c u,i1 c k,i1 c u,i2 c k,i2 .
The decision variable conditioned on {h i } N−1 i=0 is a Gaussian random variable with Equation (29) is derived from [15] for the variance of the linear combination of correlated random variables. The probability of error conditioned on {ρ i } N−1 i=0 is simply given by  The BER for the zeroth data stream is obtained via averaging P[e|ρ] over ρ: (32) Equation (32) can be evaluated by Monte Carlo integration [4,16].

NUMERICAL RESULTS
In this section, we present numerical results. To calculate the BER, it is assumed that the number of subcarriers is 32, and two channel delay profiles are considered, a simple tworay multipath delay profile often encountered in urban and hilly areas [14], and a three-path channel with linear delay power profile [17]. For the two-ray channel, we assume v 0 = 0, v 1 = 14. Furthermore, Walsh Hadamard codes are used as the spreading codes.
For the channels and the signals described above, we can compute the BER based on the closed-form expression given by (15) and (22) and the numerical method (20). Shown in Figure 1, the results obtained by the two methods match very well as long as the value of n is sufficiently large. To check how quickly the infinite series in (15) converges with respect to n, we give the BER values based on the truncated series for the index n summing up from 0 to n max in Tables 1 and  2 for two-ray and three-path channels. The BER for n max is obtained by setting n max large enough in (15), so that the first  nine significant digits for each BER value converge. It is also shown in the tables that the results of the closed-form expression conform exactly with those computed by the numerical method. The results in Tables 1 and 2 also show that the series in (15) converges very rapidly for a relatively small K-factor, and n max = 5-10 may be accurate enough for practical interest. For large K-factors, it requires n max = 15-30 to get an accurate approximation. Since the numerical method appears more computationally efficient compared with the closed-form expression especially when the K-factor is large, in Figures 2-7, the analytical results are obtained using the numerical method. The results of our proposed time-domain method described herein are verified by comparing with the Monte Carlo simulations. In Figures 2 and 3, the K-factors for different taps are equal, while in Figure 4, the K-factors for each tap are different. In all these cases, the simulation results agree very well with the theoretical results obtained by the technique proposed in this paper. The simulation results are based on the calculation of the decision variable for each bit, and averaging the results over large number of bits (say 100 000 bits). Figure 5 shows the accuracy of our lowcomplexity method proposed in Section 3.3. The results indicate that our low-complexity method gives quite a good approximation as compared with the accurate result, and the accuracy improves when the number of users is high.
Next, the effect of the Rician K-factor on the BER performance is investigated using the analytical formula derived in Section 3. The Rician K-factor is defined as the ratio of signal power in dominant component over the (local-mean) scattered power which corresponds to a deterministic strong wave received. It is natural to expect a better performance for larger values of K. From the analytical results shown in (Figure 6), it can be seen that the performance improvement due to the increase in K is evident when the number of user is small (i.e., MAI is not dominant), while for 20 users     system performance. Simulation results are also given to validate the analytical results. Figure 7 shows the BER comparison for various K when SNR equals 10 dB for 1, 10, and 20 users. A distinct variation can be seen when the number of users is small, while for 20 users, there is hardly any variation in the BER as MAI dominates the performance in this case.

CONCLUSION
In this paper, we have proposed an accurate and simple method to calculate the performance of MC-CDMA systems with MRC combining scheme over downlink correlated Rician fading channels. A general algorithm is proposed to calculate the moment generating function by expressing the decision variables in Gaussian quadratic forms. Based on the moment generating function, an exact BER is obtained in the time domain. However, the complexity of this method increases exponentially with the number of users. To alleviate this problem, we have also proposed a low-complexity approximate BER evaluation method by using the Monte Carlo integration. The results obtained by the analytical formula have been thoroughly verified by computer simulations.

DERIVATION OF MGF
In this appendix, we show the detailed derivation of (13) starting from (12). The MGF of the decision variable for frequency-selective channels under Rician distribution is derived. The decision variable Z u (0) in (12) can be rewritten as If the decision variable in (A.1) is written in vector form, the derivation of the MGF can be expressed more concisely. The vector form of the decision variable for the desired user is given by where v(0) With the assumption that the channel coefficients b(l) for different paths are statistically independent of each other and the AWGN noise term is independent of the channel impulse response, a simple expression of the covariance matrix of v(0) is achieved: where Λ(s) = (I − sΛ v ) −1 − I is a diagonal matrix.
In the case of repeated eigenvalues, the MGF of Z u (0) is given by (13). In (13) b 2 m = k∈κm μ k 2 , (A.14) where κ m denotes the set of k indices associated with the mth distinct eigenvalue.
(A. 16) The region of convergence is the vertical strip enclosing the jω axis bounded by the closest poles on either side. The BER performance can be obtained through the method described above.