- Email: [email protected]
Median equivariant adaptive separation via independence: application to communications
Neurocomputing 49 (2002) 389 – 409
www.elsevier.com/locate/neucom
Median equivariant adaptive separation via independence: application to communications Juan J. Murillo-Fuentes∗ , Francisco J. Gonz'alez-Serrano DTSC, Universidad Carlos III, Madrid, Spain Received 10 March 2001; accepted 22 October 2001
Abstract The authors propose a family of algorithms to solve the blind source separation (BSS) problem, focusing on digital communications. We 4rst include a revision of the equivariant adaptive separation via independence (EASI). These results provide a better characterization of the stability condition. Starting from this analysis we propose a new method, the Median EASI, whose stability condition presents useful properties. Next, we prove this method to be more noise-robust for some sources such as digitally modulated signals. Some experiments are included to show how the Median EASI allows phase recovering, an important issue in communications. Finally, c 2002 we propose the median EASI as a blind asynchronous CDMA multiuser detector. Elsevier Science B.V. All rights reserved. Keywords: Blind source separation; Natural gradient; Noisy signals; Maximum likelihood
1. Introduction In the linear instantaneous blind source separation (BSS) problem, m observable mixtures of n sources and the assumption of statistical independence are the only information to obtain the original signals. This may be written, in matrix form, as xt = Ast + nt ;
(1)
where A is an unknown non-singular mixture matrix m × n, t is the sample time, nt is a noise vector of variance , and s is a non-observable column vector, the sources. In the noiseless case, = 0, the aim of BSS is to identify matrix B = A−1 so that ∗
Corresponding author. E-mail address: [email protected] (J.J. Murillo-Fuentes).
c 2002 Elsevier Science B.V. All rights reserved. 0925-2312/02/$ - see front matter PII: S 0 9 2 5 - 2 3 1 2 ( 0 2 ) 0 0 5 2 3 - 4
390
J.J. Murillo-Fuentes, F.J. Gonz"alez-Serrano / Neurocomputing 49 (2002) 389 – 409
y = Bx = BAs. In [10] the author states that matrix A may be computed up to scales and permutations if no more than one Gaussian source is present in the mixture. In telecommunication applications, the entries of vector x are the transmitted signals. And in the complex case, (1) is the familiar model used in bandpass array processing. As noise is usually present, the adaptive BSS consists of updating a n × m separating matrix B to estimate the sources as yt = Bxt + Bnt = vt + Bnt = vt + rt :
(2)
Several techniques have been proposed to solve the BSS problem. Adaptive methods may optimize a contrast function of cumulants [11,16] or solve an ‘estimation equation’ [6]. This last approximation uses the natural (or relative) gradient [1,6] as learning law to minimize an objective function based on maximum likelihood (ML) [20]. In [6] the equivariant adaptive separation via independence (EASI) was presented as a BSS algorithm with a good uniform performance at low noise levels. We 4rst revise this method providing new stability conditions diJerent from those given by its authors. Then we use the Median gradient [17] along with the EASI algorithm to improve its convergence in the BSS of noisy mixtures of discrete sources such as digitally modulated signals. The new algorithm, Median-EASI (M-EASI), is applied to CDMA multiuser detection. 2. Algorithm EASI Given the statistical distributions of the sources q1 (y1 ); q2 (y2 ); : : : ; qn (yn ), the ML approximation proposes to minimize: n (B) = log |det(B)| − E log(qi (yi )) = lB (B) + E[ly (x; B)]: (3) i=1
The steepest descent technique of minimization updates the separating matrix B by a small step in a direction opposite to the gradient of the objective function, ∇ (B). The ˜ natural gradient proposes the direction ∇ (B) = ∇ (B)B T B. The stochastic version ˜ uses its instantaneous value ∇ly (x; B). The learning law yields Bt+1 = Bt + (I − ’(yt )ytT )Bt ;
(4)
where, in the ML, ’i (yi ) = −qi (yi )=qi (yi ) is the so-called score or activation function. However, the probability density functions (p.d.f.) are unknown and every author proposes each own family of score functions ’i (yi ). Some of these functions may be found in [2] for sources with positive and negative kurtoses. The EASI [6] is a ML algorithm based on the relative gradient. This method updates matrix B as follows: 1 T T T Bt+1 = Bt − yt yt − I + (’(yt )yt − yt ’(yt ) ) Bt : (5) The 4rst part of this learning law is the projection of Eq. (4) onto the space of symmetric matrices, i.e., a decorrelator. The second term corresponds to the skew-symmetric
J.J. Murillo-Fuentes, F.J. Gonz"alez-Serrano / Neurocomputing 49 (2002) 389 – 409
391
linear projection. Note that we have rewritten (5) by replacing ’(y), in [6] or [5], by (1=)’(y). We present next a novel result on the stability of the method [6,12]. 2.1. Stability of EASI The steepest descent learning law in (5) derived from the natural (relative) gradient yields @ T @ ˜ ∇(B) = B B= B; @B @X
(6)
dX = dBB−1 :
(7)
where
In the stochastic version we use the instantaneous value of the gradient: @l @ =E : @B @B
(8)
Thus, the diJerential of the instantaneous objective function dl(B) should be calculated. This function may be written as dl = yT dX y − tr(dX ) + (’(y)T dX y − yT dX ’(y)):
(9)
By developing @l=@X @l = yyT − I + (’(y)yT − y’(y)T ) @X
(10)
and multiplying it by matrix B we obtain the stochastic natural gradient @l @l = B = (yyT − I + (’(y)yT − y’(y)T ))B: @B @X
(11)
Operating as in [3,4], we consider the expected version of the gradient B˙ t = −E[@l=@X ]Bt . By linearizing it at the equilibrium point we have the variational equation [13] B˙ t = −
@(E[@l=@X ]B) B: @B
(12)
This shows that, for suOcient small step-size , only when −@(E[@l=@X ]B)=@B is de4nite-negative, the equilibrium is asymptotically stable. This operator may be evaluated in terms of dX by computing the Hessian d 2 l @l(y; B) d2 l = dbij dbkl (13) @bij @bkl as a function of X . The equilibrium is stable if and only if the expectation of the above quadratic form is positive-de4nite. We evaluated it and obtained the following results.
392
J.J. Murillo-Fuentes, F.J. Gonz"alez-Serrano / Neurocomputing 49 (2002) 389 – 409
Fig. 1. Stability conditions of the EASI in (14) (solid) and in (17) (dotted).
Theorem 1. The solution to separation of the learning law in (5) is a stable equilibrium if and only if 1 1 ((i − j ) − (j − i ))2 − 8 ((i − j ) + (j − i )) ¡ 0; (14) 2 where i = E[yi ’i (yi )];
(15)
i = E[’i (yi )]:
(16)
Proof. See Appendix A.1. The stability condition for the EASI method was given in [6] using the Lyapunov method, and also in [5], as 1 (17) (i − i + j − j ) ¿ 0: The 4rst remarkable point in the comparison of conditions in (15) and (17) is that the last one does not depends on parameter . By setting = 1 we may clearly compare both conditions to conclude that they are quite diJerent. On the other hand, it is interesting to analyze the case when → ∞, as both conditions yield the same result. Fig. 1 includes condition (14) in solid line. Points on the left or below are unstable. Condition (17) in [6] is depicted in dotted line. It can be observed in the 4gure that if sources have the same p.d.f. and the score function ’i (·) = ’j (·) i; j = 1; : : : ; n then both
J.J. Murillo-Fuentes, F.J. Gonz"alez-Serrano / Neurocomputing 49 (2002) 389 – 409
393
10 0.8 0.7 5
Stability Cond.
Perf. Index
0.6 0.5 0.4 0.3
0
_5
0.2 0.1
_ 10
0 500
(a)
1000
1500
2000
2500
500
(b)
Samples
1000
1500
2000
2500
2000
2500
Samples
10 0.02
Stability Cond.
5
Perf. Index
0.015 0
_5
_ 10 2000
0.01
0.005
2100
(c)
2200
2300
2400
500
2500
(d)
Samples
1000
1500
Samples
Fig. 2. Stability at separation using EASI for a mixture of a 2-PAM signal and noise: (a) performance index of the EASI (solid line) and a ‘stabilized version’ EASI (dotted line), (b) stability conditions in (17) (solid) and (14) (dotted), (c) instantaneous stability conditions in (17) (solid) and (14) (dotted). In (d), the performance index for a mixture of two 2-PAM.
conditions yield the same expression. However, if sources have diJerent distributions or we use diJerent score functions, then these conditions may not have the same sign. The performance index
n n n n |pij | |p | ij Q= − 1 + −1 ; (18) maxk |pik | maxk |pkj | i=1
j=1
j=1
i=1
where P = (pij ) = BA is included (solid line) in Fig. 2(a) for a mixture of a 2-PAM (pulse amplitude modulation, discrete source with possible values ±1) signal and a Gaussian random noise, score functions ’i (yi ) = |yi |4 yi , and = 2. The mixing matrix was A = I , and the sources were normalized to unit variance. In Figs. 2(b) and (c), we have depicted stability conditions (17) (solid line) and (14) (dotted line). In Fig. 2(b), these conditions were computed starting from the data available at that sample time. In Fig. 2(c), these conditions were evaluated instantaneously, i.e., for just one sample.
394
J.J. Murillo-Fuentes, F.J. Gonz"alez-Serrano / Neurocomputing 49 (2002) 389 – 409
The behavior of the performance index is not as stable as it was supposed to be. Thus, some EASI based separation systems may be unstable despite, up till now, they have been classi4ed as stable. In Fig. 2(a), we included in dotted line the performance index for a ‘stabilized’ version where the demixing matrix is not updated at instantaneous unstable points. This cannot be used directly as a stabilization algorithm [4,5] but it is useful to show how instabilities aJect separation. Fig. 2(d) depicts the performance index for a mixture of two 2-PAM signals. The score functions, mixing matrix, learning rates, etc., were those of the previous experiment. In this new example the algorithm is always stable. It is interesting to remark the diJerences in the magnitude of the index in Figs. 2(d) and (a). Once we have introduced and revised the EASI algorithm we go further and extend the concept of Median Gradient [17] to this method. We will analyze the convergence and stability of the new method, the Median-EASI (M-EASI). Then we will develop the complex normalized version of the algorithm.
3. M-EASI Algorithm 3.1. Learning law A generalized EASI learning law may be written as 1 T T T Bt+1 = Bt − h(yt )g(yt ) − I + (’(yt )(yt ) − (yt )’(yt ) ) Bt ;
(19)
where the product h(yt )g(yt )T must satisfy some conditions in order to ensure separability and convergence to a solution point. On the one hand, either E[h(yi )] = 0 or E[g(yi )] = 0 must cancel. As it will be introduced next, the learning law presented assumes symmetrically distributed sources. Thus any h(yi ) or g(yi ) being an odd function is suOcient. Besides, another restriction is that this product must be chosen in concordance with the selection of the score function ’(yi ). This point will be discussed later in this section. The median gradient introduced in [17] for symmetric sources proposes to introduce the change B ← D−1=2 B where D= diag(|y1 |; |y2 |; : : : ; |yn |) (note that matrix B may be computed up to scales), this way the stochastic median gradient yields ˆ y (B) = ∇ly (x; B)B T D−1 B ∇l
(20)
and the learning law Bt+1 = Bt − (’(yt ) sign(yt )T − I )Bt :
(21)
A geometrical perspective of the problem may be found in [18]: The sing function makes the learning law evaluate the centers of mass of the score functions at each quadrant. These centers of mass have, for symmetric sources, symmetrical properties. The median gradient exploits these properties to perform the separation. By applying
J.J. Murillo-Fuentes, F.J. Gonz"alez-Serrano / Neurocomputing 49 (2002) 389 – 409
395
these concepts to the EASI algorithm we obtain the function #i (yi ) in (19) as #i (yi ) = sign(yi ):
(22)
And in the following we will use, for mixtures of digitally modulated signals, hi (yi ) = yi ;
(23)
gi (yi ) = sign(yi ):
(24)
The M-EASI learning law yields 1 T T T Bt+1 = Bt − yt sign(yt ) − I + (’(yt ) sign(yt ) − sign(yt )’(yt ) ) Bt : (25) It is important to point out that the product h(yt )g(yt )T with hi and gi as given in (23) and (24) has been designed for mixtures of digitally modulated sources. As these functions should be chosen in concordance with the score functions ’i (yi ), the product h(yt )g(yt )T may be transposed, or just a set of its entries, for other kind of signals [4,14]. 3.2. Convergence in noisy environment The functions in (23) and (24) were selected as they provide the algorithm with a good convergence in the presence of additive noise, see (1) and (2). In the case of hi = yi and gi = yi (EASI), the expectation of the product h(yt )g(yt )T yields [9] I − E[yt ytT ] = I − BE[xt xtT ]B T + BE[nt ntT ]B T = I − E[vt vtT ] + E[rt rtT ]:
(26)
Thus, the noise contributes to the expectation of the term to minimize, at any point of the steepest descent algorithm, as a bias. If we make hi (yi ) = yi and gi (yi ) = sign(yi ), the expectation may be written, for any entry (p; q) as follows pq − E[yp sign(yq )] = pq − E[(vp + rp ) sign(vq + rq )] = pq − E&I [vp sign(vq )] − E&II [vp sign(rq )]
(27)
− E&1 [rp sign(rq )] − E&2 [rp sign(vq )];
(28)
where we have de4ned the regions &I = {vq ; rq : |vq | ¿ |rq |}; &II = {vq ; rq : |vq | ¡ |rq |}; &1 = {rq ; vq : |rq | ¿ |vq |}; &2 = {rq ; vq : |rq | ¡ |vq |}; and E& [f(y)] = & f(y)py (y) dy.
(29)
396
J.J. Murillo-Fuentes, F.J. Gonz"alez-Serrano / Neurocomputing 49 (2002) 389 – 409
The M-EASI has been designed for symmetric sources. If the noise has symmetrical p.d.f., it is easy to show that E&II [vp sign(rq )] = 0 and E&2 [rp sign(vq )] = 0. The entries of the matrix E[y sign(y)] yield I − E[y sign(y)T ] = I − E&I [v sign(v)T ] − E&1 [r sign(r)T ]:
(30)
We next analyze (26) and (30). We 4rst de4ne the following quotient: (NG =
E[r 2 ] E[v2 ]
(31)
(MG =
E&1 [|r|] : E&I [|v|]
(32)
and
If noise is present in the mixture, and whenever it is possible, the sources are designed to face its harmful eJect. This is the case of digital-modulated signals. At the receiver of a digital communication system, a hard-decision or soft-decision is a simple and eJective solution to recover the transmitted signal (the sources), avoiding the noise eJect. Note that the sign function introduced in the M-EASI is a hard-decision. Thus, as long as noise be under the distance between the transmitted constellations points, the eJect of the noise at separation is negligible. In other words, the term E&1 [r sign(r)T ] is very much lower than E&I [|v|], as the region &1 ≈ {∅} and &I ≈ {v: v ∈ [ − ∞; ∞]}. Although a detailed analysis is needed for each group of distributions, the eJect of the sign function on noisy mixtures of usual source p.d.f. such as uniforms, Laplacians, Gaussians, etc., is slightly negative or positive as they have super-Gaussian or sub-Gaussian p.d.f., respectively. It may be concluded that the whitening method based in the median gradient behaves in a similar way that the learning law based on ML-NG for usual continuous p.d.f. such as uniforms, Laplacians, Gaussians, etc. If for some of them the relation (MG slightly degradates, for the other ones it improves. Therefore, the main point is that for some discrete sources, such as digitally modulated signals, the bias reduces signi4cantly. In Fig. 3, we have depicted, in continuous line, the value (MG against the source−1 to-noise ratio SNR = E[v2 ]=E[r 2 ] = (NG , the noise being Gaussian. Notice that (NG , in dashed line, is (NG = −SNR. In Fig. 3(d), we have included a Gaussian source, either (NG or (MG are just −SNR. In Fig. 3(b), we have included these values for a uniform distribution, and in Fig. 3(c) for a Laplacian p.d.f. As already mentioned, the M-EASI is less sensitive to noise for super-Gaussian sources and more sensitive for sub-Gaussian p.d.f. The diJerence between them is no signi4cant for usual continuous distributions. Actually, this value for the uniform or the Laplacian distribution is approximately ±2 dB. Besides, the M-EASI is robust for discrete sources whose distributions have no points around the origin, usually superGaussian, as digitally modulated signals. In Fig. 3(a), we include the results for a 2-PAM.
J.J. Murillo-Fuentes, F.J. Gonz"alez-Serrano / Neurocomputing 49 (2002) 389 – 409
Uniform 10
0 _ 10
0 _ 10 MG
10
_ 20
_ 20
χ
χ
MG
2 _ PAM
_ 30
_ 30
_ 40 _ 50
_ 40 0
(a)
_ 50
40
20 SNR
40
Gaussian
10
0 _ 10
0 _ 10
_ 20 _ 30 _ 40
(c)
0
(b)
χ MG
MG
χ
20 SNR
Laplacian
10
_ 50
397
_ 20 _ 30 _ 40
0
20 SNR
40
_ 50
0
(d)
20 SNR
40
Fig. 3. Values for (MG (solid line) and (NG (dashed line) along SNR at separation for a (a) 2-PAM, (b) uniform p.d.f., (c) Laplacian p.d.f. and (d) Gaussian.
3.3. Stability As in the previous section, the 4rst step in the analysis of the stability is to evaluate dl for the learning law in (25), dl = yT dX sign(y) − tr(dX ) + (’(y)T dX sign(y) − sign(y)T dX ’(y));
(33)
where dX = dBB−1 . If we now evaluate @l=@X we have, @l = y sign(y)T − I + (’(y) sign(y)T − sign(y)’(y)T ); @X
(34)
the learning law in terms of dX . We derived this function to obtain the Hessian and the following stability condition. Theorem 2. The solution to separation provided by the learning law in (25) is a stable equilibrium if and only if 1 1 (i − j )2 − 4 (i + j ) − 4 ¡ 0; 2 where i = E[’i (yi )].
(35)
398
J.J. Murillo-Fuentes, F.J. Gonz"alez-Serrano / Neurocomputing 49 (2002) 389 – 409
Proof. See Appendix A.2. The main result of the stability condition in (35) is that it does not depend on the value of the score function but just on its derivative. In this sense and in contrast with the EASI algorithm, if all sources have identical distributions and the score functions are monotically increasing functions, the learning law is always stable. Thus, it may be stated the following. Corollary 1. If the score function ’p (·)=’q (·) ∀p; q=1; : : : ; n; this ’p (·) is a monotically increasing function; and sources are distributed with identical p.d.f.; then the M-EASI algorithm in (25) is stable. On the other hand, another interesting point is that the square of the constant is now a third term in the stability condition. This allows a better control on the stability, and the performance. We must remark that the study performed within this section has been carried out for the noiseless case, as the stability is studied at separation, i.e., assuming statistical independence at the outputs. However, these results may be extended for noisy discrete sources such as digitally modulated signals at noise levels under the distance between points of the transmitted constellations. 4. Extension to complex sources and normalization 4.1. Complex sources By assuming the sources are ‘circularly distributed’, i.e., E[si2 ] = 0, 1 6 i 6 n, the extension of the M-EASI to the complex case is immediate [6]: Bt+1 = Bt + (I − ’(yt ) sgcx(yt )H )Bt ;
(36)
sgcx(y) = sign(R(y)) + j sign(I(y)):
(37)
where Here we have replaced transposition T by conjugate transposition H. In digital communications, by forcing the correlation of ’(y) (or y) and sgcx(y) to cancel at separation, we are able to recover their phases up to multiples of .=2. This is a very important issue as communication detectors need the signal along with their phases to be recovered. The EASI algorithm itself does separate the signals, but leaves their phases indeterminate. The complex version of the M-EASI yields where
Bt+1 = Bt − QXt Bt ;
(38)
1 H H H QXt = yt sgcx(yt ) − I + (’t sgcx(yt ) − sgcx(yt )’t ) : 2
(39)
J.J. Murillo-Fuentes, F.J. Gonz"alez-Serrano / Neurocomputing 49 (2002) 389 – 409
399
4.2. Normalized M-EASI If the normalized version of the EASI algorithm [6] is yt ytH − I 1 ’(yt )ytH − yt ’(yt )H + Bt+1 = Bt − Bt ; 1 + |ytH ’(yt )| 1 + ytH yt
(40)
we propose the normalized M-EASI as follows Bt+1 = Bt − QX t Bt ; where
(41)
yt sgcx(yt )H − I 1 ’t (yt ) sgcx(yt )H − sgcx(yt )’t (yt )H QX t = + : 1 + ytH sgcx(yt ) 1 + |’H t (yt )yt |
(42)
To illustrate this method we propose the following experiment. 4.3. Separation of digitally modulated signals We generated 4ve signals. These sources were two 4-QAM signals, one 16-QAM, one GMSK and a random Gaussian noise. In [2], we found some score functions for
Perf. Index (dB)
25 20 15 10 5 0 2000
4000 Samples
6000
8000
2000
4000 Samples
6000
8000
(a)
Perf. Index (dB)
25 20 15 10 5 0 (b)
Fig. 4. Performance index in a noisy mixture of 4ve digitally modulated signals with EASI ( solid) and M-EASI (◦ solid) for (a) noise variance /2 = 25 × 10−4 , and (b) /2 = 4 × 10−2 . In (b), performance index plus its standard deviation with the EASI ( dotted) and M-EASI (◦ dotted).
400
J.J. Murillo-Fuentes, F.J. Gonz"alez-Serrano / Neurocomputing 49 (2002) 389 – 409
Fig. 5. Separated digitally modulated signals by using (a) EASI (b) M-EASI.
sources with negative kurtoses. We used here the cubic function ’i (yi ) = |yi |2 yi . The performance index in (18) was computed in dB in Fig. 4. This index may be interpreted as a signal-to-interferers ratio. Each point is the average of 1000 experiments in which the mixing matrix, A, is chosen randomly: its entries are uniformly distributed variables in the range [−1; 1]. As communication transmissions have variable attenuation, matrix A was multiplied by a Gaussian variable of unit mean and variance /2 = 104 . The learning rate was 4xed to achieve the same 4nal separation level with the normalized methods analyzed: EASI and M-EASI. These learning rates were EASI =0:96×10−3 and MEASI =0:6×10−3 . It can be observed in Fig. 4(a) that the M-EASI and EASI methods have a similar performance for low noise levels. On the contrary, it can be observed in Fig. 4(b) that if we increase the variance of the noise (/2 = 4 × 10−2 ), the EASI method deteriorates signi4cantly. However, as analyzed in Section 3.2 the convergence of the M-EASI is similar to that of Fig. 4(a). It is interesting to note that the M-EASI performance index plus its standard deviation is under the mean value of the EASI index. In addition, the mean plus the standard deviation of the EASI performance index indicates that this method does not often converge in the given range of samples. Besides, as it was previously remarked, the M-EASI provides the separated signals with a phase shift multiple of .=2, see Fig. 5(b). On the contrary, signals separated with the EASI in Fig. 5(a) have their phases undetermined at separation. The application of these solutions to beamforming [7] is straightforward (Fig. 5).
5. Application to CDMA Code division multiple access (CDMA) is a multiple access technique based on spread spectrum methods. Each user is distinguished by a unique spreading code to
J.J. Murillo-Fuentes, F.J. Gonz"alez-Serrano / Neurocomputing 49 (2002) 389 – 409
401
modulate binary antipodal data. The signal sent over the shared spectrum channel is the sum of the individual users signals. De4nitions of several key terms used in this text about CDMA systems may be found in [22]. If codes of diJerent users are synchronized the application of BSS techniques is, somehow, immediate [8,15]. On the other hand, in the asynchronous case, the delays of the users should be known to apply an ideal set of matched 4lters or a decorrelator 4lter [22]. We apply here the algorithm proposed, M-EASI, to implement a blind detector. We next describe the communication system model used and the proposed solution. 5.1. Communication system model Consider a CDMA transmitter [22] using complex-valued data symbols. The transmitted baseband signal, assuming that the propagation channel of each user is memoryless, can be represented as x(t) =
L N
(i bi ci (l − 1i )g(t − lTc − 1i );
0 6 t 6 Ts ;
(43)
l=1 i=1
where L is the length of the spreading code, N the number of active users, (i a scaling factor speci4ed by the power control loop, bi the complex symbol transmitted by user i, g(t) represents the transmitter pulse shaping 4lter, Tc the chip period, and 1i the delay of each user. The spreading code of user n is denoted by ci . Note that individual chips in ci are denoted by the elements in [ci (1); : : : ; ci (L)]T with each ci (·) ∈ ± 1. Eq. (43) usually expressed as x(t) =
N
(i bi si (t − 1i );
(44)
i=1
where si (t) is the normalized signature waveform assigned to the ith user. If the transmitted pulse-shape 4lters are Nyquist and the receiver uses a chip-rate sampled matched 4lter followed by a serial to parallel converter the matrix form of the system equation is as follows: rt = Ht−1 bt−1 + Ht bt + Ht+1 bt+1 ;
(45)
where rt is a L × 1 vector, Ht−i for i = {−1; 0; 1} is a L × N matrix representing the channel response. Note that initial assumption of the presence of the delays in reception changes the memoryless channel into a memory channel with a transfer function given by H (z) = Ht−1 z −1 + Ht + Ht+1 z:
(46)
There are several proposals in [22] for the asynchronous case. The memoryless solution is a matched 4lter given that the delays of each users are known. Consequently, both of the delays and the spreading codes should be known. We next present an alternative to this receiver.
402
J.J. Murillo-Fuentes, F.J. Gonz"alez-Serrano / Neurocomputing 49 (2002) 389 – 409
subspace separator
nt xt
_
Ht _1 z 1+ H t + H t+1 z NxL
+
rt
Θ LxN
yt
M-EASI
zt
Nx1
Fig. 6. Adaptive blind receiver for asynchronous CDMA.
5.2. Blind asynchronous receiver In the context of our work, we are considering a channel with additive Gaussian white noise and that the channel is time-invariant during the period of study. At the receiver, the knowledge of the delays and the spreading sequences would lead to a direct application of matched 4lters. We propose a blind receiver where our starting hypothesis are the number of users N transmitting in the system and the symbol timing of the user of interest, assuming no knowledge on the spreading sequences of the users. In Fig. 6, we describe the receiver scheme. The received signal is sampled at chip rate in order to provide the receiver with a L × 1 vector. This way each L samples a new vector is fed into the receiver. The subspace separator [8] projects the signal space into a N -dimensional space. Finally, this N × 1 vector is the input to a BSS block (M-EASI) (25) that provides N statistically independent outputs. Subspace separator: The aim of the block is to implement a N × L matrix, separating the signal subspace from the noise subspace. As subspace separator we used the modi4ed multi-dimensional phased-locked loop (MPLL) in [8], with associated learning constant t . This method minimizes the cost function E[yt − xˆt 2 ] that measures the energy of the 4rst N components, xˆt , of yt = r t . M-EASI: The idea of introducing BSS at this stage is based on the statistical independence of the users. We exploit the spatial diversity given by the spreading codes. Here, the M-EASI in (25) is used given its good convergence properties in noisy environments. Besides, it allows recovering signals with diJerent modulations along with their phases (Fig. 7). 5.3. Results In this section we show the performance of our proposed solution. The spreading codes used are orthogonal variable spreading factor (OVSF) sequences de4ned in [21]. The spreading factor is L = 32. In order to show the behavior of our proposal related to the near–far problem [22], we consider three users with transmitted power 0, −1:25 and −3 dB, respectively. Regarding the learning constants of the algorithms, we set = 0:001 in (25) and t = 0:8=(1 + t=200) in the subspace separator as in [8]. The comparative results in Figs. 7(a) and (b) show the success of the M-EASI separation. In Fig. 5(a) it can
J.J. Murillo-Fuentes, F.J. Gonz"alez-Serrano / Neurocomputing 49 (2002) 389 – 409
(a)
403
(b)
Fig. 7. Separatated signals for (a) synchronized matched 4lter and (b) asynchronous M-EASI based CDMA receiver.
BER
10
10
10
_2
_3
_4
10
15
20
25
Eb/No Fig. 8. BER for the matched 4lter ( ) and the M-EASI-based blind detector (◦) for the user of interest.
be observed that the matched 4lter synchronized at known delays is unable to recover the third user. On the contrary, our approach separates the three users correctly, overcoming the near–far problem, given the adverse circumstances where we do not have information at all about the spreading sequences transmitted. The results shown are for a bit energy-to-noise ratio Eb =No = 10 dB. In Fig. 8, we included the BER versus Eb =No for 4ve users with the same power, the codes and delays were chosen randomly in 200 experiments. The solid lines ( ) depicts the BER for the user of interest at the output of the matched 4lter. We have included in dotted line (◦) the BER for the proposed detector [19].
404
J.J. Murillo-Fuentes, F.J. Gonz"alez-Serrano / Neurocomputing 49 (2002) 389 – 409
6. Conclusions In this paper, the equivariant adaptive separation via independence (EASI) has been revised. We used the natural gradient formulation to provide new results on the method stability. We focused on the mixtures of digital communication signals with additive noise. The new stability conditions allow a better characterization of the performance of the method at separation. Starting from this analysis we propose to apply the Median Gradient introduced by the authors in previous works to improve the performance of the EASI algorithm. We analyze this novel algorithm, the M-EASI, to show that its stability condition has useful properties. Besides, the symmetry conditions imposed by the Median Gradient makes the M-EASI more robust against noise. We 4rst apply this algorithm to a mixture of digitally modulated signals to show how the symmetrical conditions imposed by the Median EASI allow phase recovering, an important issue in communication reception. The extension of these results to beamforming is immediate. Finally, we propose the median EASI as an asynchronous CDMA multiuser detector achieving a similar performance than that of the synchronized matched 4lter.
Acknowledgements We acknowledge 4nancial support from projects CYCIT TIC99-0210 and CAM 07T=0020=2000. The authors would like to thank Antonio Caama˜no Fern'andez for his remarks and observations, and the anonymous reviewers for their helpful comments and suggestions.
Appendix A. A.1. Stability of EASI The terms of dl in (10) may be de4ned as dlW = yT dX y;
(A.1)
dlI = tr(dX );
(A.2)
dl1 = ’(y)T dX y;
(A.3)
dl2 = yT dX ’(y):
(A.4)
Thus, dl for = 1 is the sum of these functions: d 2 l = d 2 lW − d 2 lI + (d 2 l1 − d 2 l2 ):
(A.5)
J.J. Murillo-Fuentes, F.J. Gonz"alez-Serrano / Neurocomputing 49 (2002) 389 – 409
405
We evaluate these terms. The function d 2 lW in terms of dX yields d 2 lW = yT dX T dX y + yT dX dX y
(A.6)
and its expectation is evaluated as E[yT dX T dX y] = E[yi d xji d xjk yk ] = (d xji )2 + (d xii )2 ; ij
E[yT dX dX y] =
E[yi d xij d xjk yk ] =
(A.7)
i
(d xij )(d xji ):
(A.8)
ij
Starting from the terms in (A.7) and (A.8), and for each pair of outputs (i; j), we have ij the quadratic form qW = (d xij + d xji )2 so that d 2 lW may be written as ij d 2 lW = qW + 2(d xii )2 : (A.9) i
i=j
On the other hand, the term d 2 lI cancels. We next evaluate the skew-symmetric part. Its 4rst term yields d 2 l1 = yT dX T ’(y) ˙ dX y + ’(y)T dX dX y;
(A.10)
where ’(y) ˙ is a diagonal matrix whose entries are ’i (yi ). Its expectation is given as follows: E[yT dX T ’(y) ˙ dX y] = E[yi d xji ’j (yj ) d xjk yk ] =
j (d xji )2 +
=
(A.11)
i
j=i
E[’(y)T dX dX y] =
mi (d xii )2 ;
E[’i (yi ) d xij d xjk yk ] i d xij d xji +
i=j
i (d xii )2 ;
(A.12)
i
where mi = E[(yi )2 ’i (yi )], and i and i were de4ned in (15) and (16). The quadratic form of this term may be written as q1ij = j (d xji )2 + i (d xij )2 + ( i + j ) d xij d xji and the term d 2 l1 yields ij d 2 l1 = q1 + (mi + i )(d xii )2 : i=j
(A.13)
(A.14)
i
The last term in (A.5) is evaluated as follows: d 2 l2 = yT dX ’(y) ˙ dX y + ’(y)T dX T dX y:
(A.15)
406
J.J. Murillo-Fuentes, F.J. Gonz"alez-Serrano / Neurocomputing 49 (2002) 389 – 409
By evaluating its expectation we have E[yT dX ’(y) ˙ dX y] = E[yi d xij ’˙ j (yj ) d xjk yk ] =
j d xij d xji +
=
(A.16)
i
j=i
E[’(y)T dX T dX y] =
mi (d xii )2 ;
E[’i (yi ) d xji d xjk yk ] i (d xji )2 +
i=j
i (d xii )2 ;
(A.17)
i
the quadratic form for this term is written as q2ij = (j + i ) d xji d xij + j (d xij )2 + i (d xji )2 and its Hessian ij d 2 l2 = q2 + (mi + i )(d xii )2 : i=j
(A.19)
i
We may rewrite d 2 l in terms of the quadratic forms above as ij d2 l = qEASI + 2(d xii )2 ; ij
(A.18)
(A.20)
i
where ij ij = qW + (q1ij − q2ij ) qEASI
(A.21)
=(i − j + 1)(d xij )2 + (j − i + 1)(d xji )2 +( i + j − j − i + 2) d xji d xij :
(A.22)
ij Finally, qEASI is positive if and only if the following conditions hold:
i − j + 1 ¿ 0;
(A.23)
j − i + 1 ¿ 0;
(A.24)
((i − j ) − (j − i ))2 − 8((i − j ) + (j − i )) ¡ 0:
(A.25)
It can be shown that conditions (A.23) and (A.24), derived from the Sylvester criterion for positive-de4nite quadratic forms, are satis4ed whenever (A.25) is true. Thus, it follows that condition (14) should be met ∀. A.2. Stability of M-EASI The terms of dl in (34) may be de4ned as dlW = yT dX sign(y);
(A.26)
J.J. Murillo-Fuentes, F.J. Gonz"alez-Serrano / Neurocomputing 49 (2002) 389 – 409
407
dlI = tr(dX );
(A.27)
dl1 = ’(y)T dX sign(y);
(A.28)
dl2 = sign(y)T dX ’(y)
(A.29)
so that the Hessian d 2 l, for = 1, yields d 2 l = d 2 lW − d 2 lI + (d 2 l1 − d 2 l2 ):
(A.30)
The pair of terms d 2 lW − d 2 lI were evaluated in [17] for dlMG = ’(y)T dX sign(y). If we set ’(y) = y we have d 2 lW , ij qW + d xii2 ; (A.31) E[d 2 lW − d 2 lI ] = i
i=j
ij where qW = d xij2 + d xji2 is the quadratic form in terms of (d xij ; d xji ) j = i. Regarding the skew-symmetric part, its 4rst term yields
˙ dX y; d 2 l1 = sign(y)T dX T ’(y)
(A.32)
where ’(y) ˙ is a diagonal matrix whose entries are ’i (yi ). Its expectation is evaluated as follows: ˙ dX y] = E[sign(yi ) d xji ’j (yj ) d xjk yk ] E[sign(y)T dX T ’(y) =
j (d xji )2 +
mˆ i (d xii )2 ;
(A.33)
i
j=i
where mˆ i = E[|yi |’i (yi )], and i and i were de4ned in (15) and (16). Hence, ij q1 + mˆ i (d xii )2 ; (A.34) d 2 l1 = i=j
i
where it has been de4ned the quadratic form q1ij = j (d xji )2 + i (d xij )2 . The last term may be evaluated in a similar way, ˙ dX y: d 2 l2 = sign(y)T dX ’(y) Its expectation yields ˙ dX y] = E[sign(y)T dX ’(y) =
j=i
(A.35) E[sign(yi ) d xij ’j (yj ) d xjk yk ] j d xij d xji +
mˆ i (d xii )2
(A.36)
i
and the quadratic form for this term is q2ij = (j + i ) d xji d xij . Thus, d 2 l2 yields ij q2 + mˆ i (d xii )2 : (A.37) d 2 l2 = i=j
i
408
J.J. Murillo-Fuentes, F.J. Gonz"alez-Serrano / Neurocomputing 49 (2002) 389 – 409
Finally, we have the condition ij qM-EASI + 2(d xii )2 ¿ 0; d2 l = ij
(A.38)
i
where ij ij ij ij 2qM -EASI = qW + (q1 − q2 )
= (i + 1)(d xij )2 + (j + 1)(d xji )2 − (j + i ) d xji d xij :
(A.39) (A.40)
Thus, the Hessian is de4nite positive, ∀, if the condition in (35) is true. References [1] S.I. Amari, Neural learning in structured parameter spaces—natural Riemannian gradient, in: M.C. Mozer, M.I. Jordan, T. Petsche (Eds.), Advances in Neural Information Processing Systems, Vol. 9, The MIT Press, Cambridge, 1997, p. 127. [2] S.I. Amari, Natural gradient works eOciently in learning, Neural Comput. 10 (2) (1998) 251–276. [3] S. Amari, SupereOciency in blind source separation, IEEE Trans. Signal Process. 47 (4) (1999) 936–944. [4] S.I. Amari, T. Chen, A. Cichochi, Stability analysis of learning algorithms for blind source separation, Neural Networks 10 (8) (1997) 1345–1351. [5] J.-F. Cardoso, On the stability of source separation algorithms, J. VLSI Signal Process. Systems 26 (1=2) (2000) 7–14. [6] J.F. Cardoso, B.H. Laheld, Equivariant adaptive source separation, IEEE Trans. Signal Process. 44 (12) (1996) 3017–3030. [7] J.F. Cardoso, A. Souloumiac, Blind beamforming for non-Gaussian signals, Proc. IEE F 140 (6) (1993) 362–370. [8] R.T. Causey, J.R. Barry, Blind multiuser detection using linear prediction, IEEE J. Selected Areas Commun. 16 (9) (1998) 1702–1710. [9] A. Cichocki, S. Douglas, S. Amari, Robust techniques for independent component analysis (ICA) with noisy data, Neurocomputing (22) (1998) 113–129. [10] P. Comon, Independent component analysis, a new concept? Signal Process. 36 (3) (1994) 287–314. [11] N. Defolsse, P. Loubaton, Adaptive separation of independent sources: a deWation approach, in: ICASSP, 1994. [12] S.C. Douglas, Self-stabilized gradient algorithms for blind source separation with orthogonality constraints, IEEE Trans. Neural Networks 11 (06) (2000) 1490. [13] I. Gelfand, S. Fomin, Calculus of Variations, Prentice-Hall Inc., Englewood CliJs, NJ, 1963. [14] T.P.V. HoJ, A.G. Lindgren, A. Kaelin, Transpose properties in the stability and performance of the classic adaptive algorithms for blind source separation and deconvolution, Signal Process. (80) (2000) 1807–1822. [15] J. Joutsensalo, T. Ristaniemi, Learning algorithms for blind multiuser detection in CDMA downlink, in: PIMRC, 1998. [16] C. Jutten, J. Herault, Blind separation of sources I. An adaptive algorithm based on neuromimetic architecture, Signal Process. 24 (1) (1991) 1–10. [17] J. Murillo-Fuentes, F. Gonz'alez-Serrano, Higher order moments algorithms for blind signal separation, in: International Workshop on Independent Component Analysis and Blind Signal Separation, Helsinki, Finland, 2000, pp. 345 –350. [18] J. Murillo-Fuentes, F. Gonz'alez-Serrano, Improving stability in blind source separation with the stochastic median gradient, Electron. Lett. 36 (19) (2000) 1662–1663. [19] J.J. Murillo-Fuentes, M. S'anchez-Fern'andez, A. Caamano-Fern'andez, F.J. Gonz'alez-Serrano, Adaptive blind joint source-phase separation in digital communications, in: IEEE International Conference on Communications, Helsinki, Finland, 2001.
J.J. Murillo-Fuentes, F.J. Gonz"alez-Serrano / Neurocomputing 49 (2002) 389 – 409
409
[20] P.G.D.T. Pahm, C. Jutten, Separation of a mixture of independent sources through a ml approach, in: EUSIPCO, 1992, pp. 771, 774. [21] G.T. Ran, Spreading and modulation (FDD), Technical Report, 3GPP, 2000. [22] S. Verd'u, Multiuser Detection, Cambridge University Press, Cambridge, 1998. Juan Jos+e Murillo Fuentes received the telecommunication engineering degree from the University of Seville, Spain, in 1996. He is now pursuing the Ph.D. degree from the University Carlos III de Madrid, Spain. Since 1999, he has been an Assistant Professor at the University Carlos III. His main research interests include neural networks and their application to communications, array processing with emphasis on algorithms for blind separation of sources and multiuser detection, and image processing. Francisco J. Gonz+alez Serrano was born in Betanzos, Spain, on October 10, 1968. He received the telecommunication engineering and Ph.D. degrees from the University of Vigo, Spain, in 1992 and 1997, respectively. From 1992 to 1993 he was with DSP research group, University of Vigo, as a Research Associate. From 1993 to 1998, he was an Assistant Professor with the University Carlos III, Madrid, Spain. His research interest are in the areas of digital signal processing, wireless communications and multimedia communications.
www.elsevier.com/locate/neucom
Median equivariant adaptive separation via independence: application to communications Juan J. Murillo-Fuentes∗ , Francisco J. Gonz'alez-Serrano DTSC, Universidad Carlos III, Madrid, Spain Received 10 March 2001; accepted 22 October 2001
Abstract The authors propose a family of algorithms to solve the blind source separation (BSS) problem, focusing on digital communications. We 4rst include a revision of the equivariant adaptive separation via independence (EASI). These results provide a better characterization of the stability condition. Starting from this analysis we propose a new method, the Median EASI, whose stability condition presents useful properties. Next, we prove this method to be more noise-robust for some sources such as digitally modulated signals. Some experiments are included to show how the Median EASI allows phase recovering, an important issue in communications. Finally, c 2002 we propose the median EASI as a blind asynchronous CDMA multiuser detector. Elsevier Science B.V. All rights reserved. Keywords: Blind source separation; Natural gradient; Noisy signals; Maximum likelihood
1. Introduction In the linear instantaneous blind source separation (BSS) problem, m observable mixtures of n sources and the assumption of statistical independence are the only information to obtain the original signals. This may be written, in matrix form, as xt = Ast + nt ;
(1)
where A is an unknown non-singular mixture matrix m × n, t is the sample time, nt is a noise vector of variance , and s is a non-observable column vector, the sources. In the noiseless case, = 0, the aim of BSS is to identify matrix B = A−1 so that ∗
Corresponding author. E-mail address: [email protected] (J.J. Murillo-Fuentes).
c 2002 Elsevier Science B.V. All rights reserved. 0925-2312/02/$ - see front matter PII: S 0 9 2 5 - 2 3 1 2 ( 0 2 ) 0 0 5 2 3 - 4
390
J.J. Murillo-Fuentes, F.J. Gonz"alez-Serrano / Neurocomputing 49 (2002) 389 – 409
y = Bx = BAs. In [10] the author states that matrix A may be computed up to scales and permutations if no more than one Gaussian source is present in the mixture. In telecommunication applications, the entries of vector x are the transmitted signals. And in the complex case, (1) is the familiar model used in bandpass array processing. As noise is usually present, the adaptive BSS consists of updating a n × m separating matrix B to estimate the sources as yt = Bxt + Bnt = vt + Bnt = vt + rt :
(2)
Several techniques have been proposed to solve the BSS problem. Adaptive methods may optimize a contrast function of cumulants [11,16] or solve an ‘estimation equation’ [6]. This last approximation uses the natural (or relative) gradient [1,6] as learning law to minimize an objective function based on maximum likelihood (ML) [20]. In [6] the equivariant adaptive separation via independence (EASI) was presented as a BSS algorithm with a good uniform performance at low noise levels. We 4rst revise this method providing new stability conditions diJerent from those given by its authors. Then we use the Median gradient [17] along with the EASI algorithm to improve its convergence in the BSS of noisy mixtures of discrete sources such as digitally modulated signals. The new algorithm, Median-EASI (M-EASI), is applied to CDMA multiuser detection. 2. Algorithm EASI Given the statistical distributions of the sources q1 (y1 ); q2 (y2 ); : : : ; qn (yn ), the ML approximation proposes to minimize: n (B) = log |det(B)| − E log(qi (yi )) = lB (B) + E[ly (x; B)]: (3) i=1
The steepest descent technique of minimization updates the separating matrix B by a small step in a direction opposite to the gradient of the objective function, ∇ (B). The ˜ natural gradient proposes the direction ∇ (B) = ∇ (B)B T B. The stochastic version ˜ uses its instantaneous value ∇ly (x; B). The learning law yields Bt+1 = Bt + (I − ’(yt )ytT )Bt ;
(4)
where, in the ML, ’i (yi ) = −qi (yi )=qi (yi ) is the so-called score or activation function. However, the probability density functions (p.d.f.) are unknown and every author proposes each own family of score functions ’i (yi ). Some of these functions may be found in [2] for sources with positive and negative kurtoses. The EASI [6] is a ML algorithm based on the relative gradient. This method updates matrix B as follows: 1 T T T Bt+1 = Bt − yt yt − I + (’(yt )yt − yt ’(yt ) ) Bt : (5) The 4rst part of this learning law is the projection of Eq. (4) onto the space of symmetric matrices, i.e., a decorrelator. The second term corresponds to the skew-symmetric
J.J. Murillo-Fuentes, F.J. Gonz"alez-Serrano / Neurocomputing 49 (2002) 389 – 409
391
linear projection. Note that we have rewritten (5) by replacing ’(y), in [6] or [5], by (1=)’(y). We present next a novel result on the stability of the method [6,12]. 2.1. Stability of EASI The steepest descent learning law in (5) derived from the natural (relative) gradient yields @ T @ ˜ ∇(B) = B B= B; @B @X
(6)
dX = dBB−1 :
(7)
where
In the stochastic version we use the instantaneous value of the gradient: @l @ =E : @B @B
(8)
Thus, the diJerential of the instantaneous objective function dl(B) should be calculated. This function may be written as dl = yT dX y − tr(dX ) + (’(y)T dX y − yT dX ’(y)):
(9)
By developing @l=@X @l = yyT − I + (’(y)yT − y’(y)T ) @X
(10)
and multiplying it by matrix B we obtain the stochastic natural gradient @l @l = B = (yyT − I + (’(y)yT − y’(y)T ))B: @B @X
(11)
Operating as in [3,4], we consider the expected version of the gradient B˙ t = −E[@l=@X ]Bt . By linearizing it at the equilibrium point we have the variational equation [13] B˙ t = −
@(E[@l=@X ]B) B: @B
(12)
This shows that, for suOcient small step-size , only when −@(E[@l=@X ]B)=@B is de4nite-negative, the equilibrium is asymptotically stable. This operator may be evaluated in terms of dX by computing the Hessian d 2 l @l(y; B) d2 l = dbij dbkl (13) @bij @bkl as a function of X . The equilibrium is stable if and only if the expectation of the above quadratic form is positive-de4nite. We evaluated it and obtained the following results.
392
J.J. Murillo-Fuentes, F.J. Gonz"alez-Serrano / Neurocomputing 49 (2002) 389 – 409
Fig. 1. Stability conditions of the EASI in (14) (solid) and in (17) (dotted).
Theorem 1. The solution to separation of the learning law in (5) is a stable equilibrium if and only if 1 1 ((i − j ) − (j − i ))2 − 8 ((i − j ) + (j − i )) ¡ 0; (14) 2 where i = E[yi ’i (yi )];
(15)
i = E[’i (yi )]:
(16)
Proof. See Appendix A.1. The stability condition for the EASI method was given in [6] using the Lyapunov method, and also in [5], as 1 (17) (i − i + j − j ) ¿ 0: The 4rst remarkable point in the comparison of conditions in (15) and (17) is that the last one does not depends on parameter . By setting = 1 we may clearly compare both conditions to conclude that they are quite diJerent. On the other hand, it is interesting to analyze the case when → ∞, as both conditions yield the same result. Fig. 1 includes condition (14) in solid line. Points on the left or below are unstable. Condition (17) in [6] is depicted in dotted line. It can be observed in the 4gure that if sources have the same p.d.f. and the score function ’i (·) = ’j (·) i; j = 1; : : : ; n then both
J.J. Murillo-Fuentes, F.J. Gonz"alez-Serrano / Neurocomputing 49 (2002) 389 – 409
393
10 0.8 0.7 5
Stability Cond.
Perf. Index
0.6 0.5 0.4 0.3
0
_5
0.2 0.1
_ 10
0 500
(a)
1000
1500
2000
2500
500
(b)
Samples
1000
1500
2000
2500
2000
2500
Samples
10 0.02
Stability Cond.
5
Perf. Index
0.015 0
_5
_ 10 2000
0.01
0.005
2100
(c)
2200
2300
2400
500
2500
(d)
Samples
1000
1500
Samples
Fig. 2. Stability at separation using EASI for a mixture of a 2-PAM signal and noise: (a) performance index of the EASI (solid line) and a ‘stabilized version’ EASI (dotted line), (b) stability conditions in (17) (solid) and (14) (dotted), (c) instantaneous stability conditions in (17) (solid) and (14) (dotted). In (d), the performance index for a mixture of two 2-PAM.
conditions yield the same expression. However, if sources have diJerent distributions or we use diJerent score functions, then these conditions may not have the same sign. The performance index
n n n n |pij | |p | ij Q= − 1 + −1 ; (18) maxk |pik | maxk |pkj | i=1
j=1
j=1
i=1
where P = (pij ) = BA is included (solid line) in Fig. 2(a) for a mixture of a 2-PAM (pulse amplitude modulation, discrete source with possible values ±1) signal and a Gaussian random noise, score functions ’i (yi ) = |yi |4 yi , and = 2. The mixing matrix was A = I , and the sources were normalized to unit variance. In Figs. 2(b) and (c), we have depicted stability conditions (17) (solid line) and (14) (dotted line). In Fig. 2(b), these conditions were computed starting from the data available at that sample time. In Fig. 2(c), these conditions were evaluated instantaneously, i.e., for just one sample.
394
J.J. Murillo-Fuentes, F.J. Gonz"alez-Serrano / Neurocomputing 49 (2002) 389 – 409
The behavior of the performance index is not as stable as it was supposed to be. Thus, some EASI based separation systems may be unstable despite, up till now, they have been classi4ed as stable. In Fig. 2(a), we included in dotted line the performance index for a ‘stabilized’ version where the demixing matrix is not updated at instantaneous unstable points. This cannot be used directly as a stabilization algorithm [4,5] but it is useful to show how instabilities aJect separation. Fig. 2(d) depicts the performance index for a mixture of two 2-PAM signals. The score functions, mixing matrix, learning rates, etc., were those of the previous experiment. In this new example the algorithm is always stable. It is interesting to remark the diJerences in the magnitude of the index in Figs. 2(d) and (a). Once we have introduced and revised the EASI algorithm we go further and extend the concept of Median Gradient [17] to this method. We will analyze the convergence and stability of the new method, the Median-EASI (M-EASI). Then we will develop the complex normalized version of the algorithm.
3. M-EASI Algorithm 3.1. Learning law A generalized EASI learning law may be written as 1 T T T Bt+1 = Bt − h(yt )g(yt ) − I + (’(yt )(yt ) − (yt )’(yt ) ) Bt ;
(19)
where the product h(yt )g(yt )T must satisfy some conditions in order to ensure separability and convergence to a solution point. On the one hand, either E[h(yi )] = 0 or E[g(yi )] = 0 must cancel. As it will be introduced next, the learning law presented assumes symmetrically distributed sources. Thus any h(yi ) or g(yi ) being an odd function is suOcient. Besides, another restriction is that this product must be chosen in concordance with the selection of the score function ’(yi ). This point will be discussed later in this section. The median gradient introduced in [17] for symmetric sources proposes to introduce the change B ← D−1=2 B where D= diag(|y1 |; |y2 |; : : : ; |yn |) (note that matrix B may be computed up to scales), this way the stochastic median gradient yields ˆ y (B) = ∇ly (x; B)B T D−1 B ∇l
(20)
and the learning law Bt+1 = Bt − (’(yt ) sign(yt )T − I )Bt :
(21)
A geometrical perspective of the problem may be found in [18]: The sing function makes the learning law evaluate the centers of mass of the score functions at each quadrant. These centers of mass have, for symmetric sources, symmetrical properties. The median gradient exploits these properties to perform the separation. By applying
J.J. Murillo-Fuentes, F.J. Gonz"alez-Serrano / Neurocomputing 49 (2002) 389 – 409
395
these concepts to the EASI algorithm we obtain the function #i (yi ) in (19) as #i (yi ) = sign(yi ):
(22)
And in the following we will use, for mixtures of digitally modulated signals, hi (yi ) = yi ;
(23)
gi (yi ) = sign(yi ):
(24)
The M-EASI learning law yields 1 T T T Bt+1 = Bt − yt sign(yt ) − I + (’(yt ) sign(yt ) − sign(yt )’(yt ) ) Bt : (25) It is important to point out that the product h(yt )g(yt )T with hi and gi as given in (23) and (24) has been designed for mixtures of digitally modulated sources. As these functions should be chosen in concordance with the score functions ’i (yi ), the product h(yt )g(yt )T may be transposed, or just a set of its entries, for other kind of signals [4,14]. 3.2. Convergence in noisy environment The functions in (23) and (24) were selected as they provide the algorithm with a good convergence in the presence of additive noise, see (1) and (2). In the case of hi = yi and gi = yi (EASI), the expectation of the product h(yt )g(yt )T yields [9] I − E[yt ytT ] = I − BE[xt xtT ]B T + BE[nt ntT ]B T = I − E[vt vtT ] + E[rt rtT ]:
(26)
Thus, the noise contributes to the expectation of the term to minimize, at any point of the steepest descent algorithm, as a bias. If we make hi (yi ) = yi and gi (yi ) = sign(yi ), the expectation may be written, for any entry (p; q) as follows pq − E[yp sign(yq )] = pq − E[(vp + rp ) sign(vq + rq )] = pq − E&I [vp sign(vq )] − E&II [vp sign(rq )]
(27)
− E&1 [rp sign(rq )] − E&2 [rp sign(vq )];
(28)
where we have de4ned the regions &I = {vq ; rq : |vq | ¿ |rq |}; &II = {vq ; rq : |vq | ¡ |rq |}; &1 = {rq ; vq : |rq | ¿ |vq |}; &2 = {rq ; vq : |rq | ¡ |vq |}; and E& [f(y)] = & f(y)py (y) dy.
(29)
396
J.J. Murillo-Fuentes, F.J. Gonz"alez-Serrano / Neurocomputing 49 (2002) 389 – 409
The M-EASI has been designed for symmetric sources. If the noise has symmetrical p.d.f., it is easy to show that E&II [vp sign(rq )] = 0 and E&2 [rp sign(vq )] = 0. The entries of the matrix E[y sign(y)] yield I − E[y sign(y)T ] = I − E&I [v sign(v)T ] − E&1 [r sign(r)T ]:
(30)
We next analyze (26) and (30). We 4rst de4ne the following quotient: (NG =
E[r 2 ] E[v2 ]
(31)
(MG =
E&1 [|r|] : E&I [|v|]
(32)
and
If noise is present in the mixture, and whenever it is possible, the sources are designed to face its harmful eJect. This is the case of digital-modulated signals. At the receiver of a digital communication system, a hard-decision or soft-decision is a simple and eJective solution to recover the transmitted signal (the sources), avoiding the noise eJect. Note that the sign function introduced in the M-EASI is a hard-decision. Thus, as long as noise be under the distance between the transmitted constellations points, the eJect of the noise at separation is negligible. In other words, the term E&1 [r sign(r)T ] is very much lower than E&I [|v|], as the region &1 ≈ {∅} and &I ≈ {v: v ∈ [ − ∞; ∞]}. Although a detailed analysis is needed for each group of distributions, the eJect of the sign function on noisy mixtures of usual source p.d.f. such as uniforms, Laplacians, Gaussians, etc., is slightly negative or positive as they have super-Gaussian or sub-Gaussian p.d.f., respectively. It may be concluded that the whitening method based in the median gradient behaves in a similar way that the learning law based on ML-NG for usual continuous p.d.f. such as uniforms, Laplacians, Gaussians, etc. If for some of them the relation (MG slightly degradates, for the other ones it improves. Therefore, the main point is that for some discrete sources, such as digitally modulated signals, the bias reduces signi4cantly. In Fig. 3, we have depicted, in continuous line, the value (MG against the source−1 to-noise ratio SNR = E[v2 ]=E[r 2 ] = (NG , the noise being Gaussian. Notice that (NG , in dashed line, is (NG = −SNR. In Fig. 3(d), we have included a Gaussian source, either (NG or (MG are just −SNR. In Fig. 3(b), we have included these values for a uniform distribution, and in Fig. 3(c) for a Laplacian p.d.f. As already mentioned, the M-EASI is less sensitive to noise for super-Gaussian sources and more sensitive for sub-Gaussian p.d.f. The diJerence between them is no signi4cant for usual continuous distributions. Actually, this value for the uniform or the Laplacian distribution is approximately ±2 dB. Besides, the M-EASI is robust for discrete sources whose distributions have no points around the origin, usually superGaussian, as digitally modulated signals. In Fig. 3(a), we include the results for a 2-PAM.
J.J. Murillo-Fuentes, F.J. Gonz"alez-Serrano / Neurocomputing 49 (2002) 389 – 409
Uniform 10
0 _ 10
0 _ 10 MG
10
_ 20
_ 20
χ
χ
MG
2 _ PAM
_ 30
_ 30
_ 40 _ 50
_ 40 0
(a)
_ 50
40
20 SNR
40
Gaussian
10
0 _ 10
0 _ 10
_ 20 _ 30 _ 40
(c)
0
(b)
χ MG
MG
χ
20 SNR
Laplacian
10
_ 50
397
_ 20 _ 30 _ 40
0
20 SNR
40
_ 50
0
(d)
20 SNR
40
Fig. 3. Values for (MG (solid line) and (NG (dashed line) along SNR at separation for a (a) 2-PAM, (b) uniform p.d.f., (c) Laplacian p.d.f. and (d) Gaussian.
3.3. Stability As in the previous section, the 4rst step in the analysis of the stability is to evaluate dl for the learning law in (25), dl = yT dX sign(y) − tr(dX ) + (’(y)T dX sign(y) − sign(y)T dX ’(y));
(33)
where dX = dBB−1 . If we now evaluate @l=@X we have, @l = y sign(y)T − I + (’(y) sign(y)T − sign(y)’(y)T ); @X
(34)
the learning law in terms of dX . We derived this function to obtain the Hessian and the following stability condition. Theorem 2. The solution to separation provided by the learning law in (25) is a stable equilibrium if and only if 1 1 (i − j )2 − 4 (i + j ) − 4 ¡ 0; 2 where i = E[’i (yi )].
(35)
398
J.J. Murillo-Fuentes, F.J. Gonz"alez-Serrano / Neurocomputing 49 (2002) 389 – 409
Proof. See Appendix A.2. The main result of the stability condition in (35) is that it does not depend on the value of the score function but just on its derivative. In this sense and in contrast with the EASI algorithm, if all sources have identical distributions and the score functions are monotically increasing functions, the learning law is always stable. Thus, it may be stated the following. Corollary 1. If the score function ’p (·)=’q (·) ∀p; q=1; : : : ; n; this ’p (·) is a monotically increasing function; and sources are distributed with identical p.d.f.; then the M-EASI algorithm in (25) is stable. On the other hand, another interesting point is that the square of the constant is now a third term in the stability condition. This allows a better control on the stability, and the performance. We must remark that the study performed within this section has been carried out for the noiseless case, as the stability is studied at separation, i.e., assuming statistical independence at the outputs. However, these results may be extended for noisy discrete sources such as digitally modulated signals at noise levels under the distance between points of the transmitted constellations. 4. Extension to complex sources and normalization 4.1. Complex sources By assuming the sources are ‘circularly distributed’, i.e., E[si2 ] = 0, 1 6 i 6 n, the extension of the M-EASI to the complex case is immediate [6]: Bt+1 = Bt + (I − ’(yt ) sgcx(yt )H )Bt ;
(36)
sgcx(y) = sign(R(y)) + j sign(I(y)):
(37)
where Here we have replaced transposition T by conjugate transposition H. In digital communications, by forcing the correlation of ’(y) (or y) and sgcx(y) to cancel at separation, we are able to recover their phases up to multiples of .=2. This is a very important issue as communication detectors need the signal along with their phases to be recovered. The EASI algorithm itself does separate the signals, but leaves their phases indeterminate. The complex version of the M-EASI yields where
Bt+1 = Bt − QXt Bt ;
(38)
1 H H H QXt = yt sgcx(yt ) − I + (’t sgcx(yt ) − sgcx(yt )’t ) : 2
(39)
J.J. Murillo-Fuentes, F.J. Gonz"alez-Serrano / Neurocomputing 49 (2002) 389 – 409
399
4.2. Normalized M-EASI If the normalized version of the EASI algorithm [6] is yt ytH − I 1 ’(yt )ytH − yt ’(yt )H + Bt+1 = Bt − Bt ; 1 + |ytH ’(yt )| 1 + ytH yt
(40)
we propose the normalized M-EASI as follows Bt+1 = Bt − QX t Bt ; where
(41)
yt sgcx(yt )H − I 1 ’t (yt ) sgcx(yt )H − sgcx(yt )’t (yt )H QX t = + : 1 + ytH sgcx(yt ) 1 + |’H t (yt )yt |
(42)
To illustrate this method we propose the following experiment. 4.3. Separation of digitally modulated signals We generated 4ve signals. These sources were two 4-QAM signals, one 16-QAM, one GMSK and a random Gaussian noise. In [2], we found some score functions for
Perf. Index (dB)
25 20 15 10 5 0 2000
4000 Samples
6000
8000
2000
4000 Samples
6000
8000
(a)
Perf. Index (dB)
25 20 15 10 5 0 (b)
Fig. 4. Performance index in a noisy mixture of 4ve digitally modulated signals with EASI ( solid) and M-EASI (◦ solid) for (a) noise variance /2 = 25 × 10−4 , and (b) /2 = 4 × 10−2 . In (b), performance index plus its standard deviation with the EASI ( dotted) and M-EASI (◦ dotted).
400
J.J. Murillo-Fuentes, F.J. Gonz"alez-Serrano / Neurocomputing 49 (2002) 389 – 409
Fig. 5. Separated digitally modulated signals by using (a) EASI (b) M-EASI.
sources with negative kurtoses. We used here the cubic function ’i (yi ) = |yi |2 yi . The performance index in (18) was computed in dB in Fig. 4. This index may be interpreted as a signal-to-interferers ratio. Each point is the average of 1000 experiments in which the mixing matrix, A, is chosen randomly: its entries are uniformly distributed variables in the range [−1; 1]. As communication transmissions have variable attenuation, matrix A was multiplied by a Gaussian variable of unit mean and variance /2 = 104 . The learning rate was 4xed to achieve the same 4nal separation level with the normalized methods analyzed: EASI and M-EASI. These learning rates were EASI =0:96×10−3 and MEASI =0:6×10−3 . It can be observed in Fig. 4(a) that the M-EASI and EASI methods have a similar performance for low noise levels. On the contrary, it can be observed in Fig. 4(b) that if we increase the variance of the noise (/2 = 4 × 10−2 ), the EASI method deteriorates signi4cantly. However, as analyzed in Section 3.2 the convergence of the M-EASI is similar to that of Fig. 4(a). It is interesting to note that the M-EASI performance index plus its standard deviation is under the mean value of the EASI index. In addition, the mean plus the standard deviation of the EASI performance index indicates that this method does not often converge in the given range of samples. Besides, as it was previously remarked, the M-EASI provides the separated signals with a phase shift multiple of .=2, see Fig. 5(b). On the contrary, signals separated with the EASI in Fig. 5(a) have their phases undetermined at separation. The application of these solutions to beamforming [7] is straightforward (Fig. 5).
5. Application to CDMA Code division multiple access (CDMA) is a multiple access technique based on spread spectrum methods. Each user is distinguished by a unique spreading code to
J.J. Murillo-Fuentes, F.J. Gonz"alez-Serrano / Neurocomputing 49 (2002) 389 – 409
401
modulate binary antipodal data. The signal sent over the shared spectrum channel is the sum of the individual users signals. De4nitions of several key terms used in this text about CDMA systems may be found in [22]. If codes of diJerent users are synchronized the application of BSS techniques is, somehow, immediate [8,15]. On the other hand, in the asynchronous case, the delays of the users should be known to apply an ideal set of matched 4lters or a decorrelator 4lter [22]. We apply here the algorithm proposed, M-EASI, to implement a blind detector. We next describe the communication system model used and the proposed solution. 5.1. Communication system model Consider a CDMA transmitter [22] using complex-valued data symbols. The transmitted baseband signal, assuming that the propagation channel of each user is memoryless, can be represented as x(t) =
L N
(i bi ci (l − 1i )g(t − lTc − 1i );
0 6 t 6 Ts ;
(43)
l=1 i=1
where L is the length of the spreading code, N the number of active users, (i a scaling factor speci4ed by the power control loop, bi the complex symbol transmitted by user i, g(t) represents the transmitter pulse shaping 4lter, Tc the chip period, and 1i the delay of each user. The spreading code of user n is denoted by ci . Note that individual chips in ci are denoted by the elements in [ci (1); : : : ; ci (L)]T with each ci (·) ∈ ± 1. Eq. (43) usually expressed as x(t) =
N
(i bi si (t − 1i );
(44)
i=1
where si (t) is the normalized signature waveform assigned to the ith user. If the transmitted pulse-shape 4lters are Nyquist and the receiver uses a chip-rate sampled matched 4lter followed by a serial to parallel converter the matrix form of the system equation is as follows: rt = Ht−1 bt−1 + Ht bt + Ht+1 bt+1 ;
(45)
where rt is a L × 1 vector, Ht−i for i = {−1; 0; 1} is a L × N matrix representing the channel response. Note that initial assumption of the presence of the delays in reception changes the memoryless channel into a memory channel with a transfer function given by H (z) = Ht−1 z −1 + Ht + Ht+1 z:
(46)
There are several proposals in [22] for the asynchronous case. The memoryless solution is a matched 4lter given that the delays of each users are known. Consequently, both of the delays and the spreading codes should be known. We next present an alternative to this receiver.
402
J.J. Murillo-Fuentes, F.J. Gonz"alez-Serrano / Neurocomputing 49 (2002) 389 – 409
subspace separator
nt xt
_
Ht _1 z 1+ H t + H t+1 z NxL
+
rt
Θ LxN
yt
M-EASI
zt
Nx1
Fig. 6. Adaptive blind receiver for asynchronous CDMA.
5.2. Blind asynchronous receiver In the context of our work, we are considering a channel with additive Gaussian white noise and that the channel is time-invariant during the period of study. At the receiver, the knowledge of the delays and the spreading sequences would lead to a direct application of matched 4lters. We propose a blind receiver where our starting hypothesis are the number of users N transmitting in the system and the symbol timing of the user of interest, assuming no knowledge on the spreading sequences of the users. In Fig. 6, we describe the receiver scheme. The received signal is sampled at chip rate in order to provide the receiver with a L × 1 vector. This way each L samples a new vector is fed into the receiver. The subspace separator [8] projects the signal space into a N -dimensional space. Finally, this N × 1 vector is the input to a BSS block (M-EASI) (25) that provides N statistically independent outputs. Subspace separator: The aim of the block is to implement a N × L matrix, separating the signal subspace from the noise subspace. As subspace separator we used the modi4ed multi-dimensional phased-locked loop (MPLL) in [8], with associated learning constant t . This method minimizes the cost function E[yt − xˆt 2 ] that measures the energy of the 4rst N components, xˆt , of yt = r t . M-EASI: The idea of introducing BSS at this stage is based on the statistical independence of the users. We exploit the spatial diversity given by the spreading codes. Here, the M-EASI in (25) is used given its good convergence properties in noisy environments. Besides, it allows recovering signals with diJerent modulations along with their phases (Fig. 7). 5.3. Results In this section we show the performance of our proposed solution. The spreading codes used are orthogonal variable spreading factor (OVSF) sequences de4ned in [21]. The spreading factor is L = 32. In order to show the behavior of our proposal related to the near–far problem [22], we consider three users with transmitted power 0, −1:25 and −3 dB, respectively. Regarding the learning constants of the algorithms, we set = 0:001 in (25) and t = 0:8=(1 + t=200) in the subspace separator as in [8]. The comparative results in Figs. 7(a) and (b) show the success of the M-EASI separation. In Fig. 5(a) it can
J.J. Murillo-Fuentes, F.J. Gonz"alez-Serrano / Neurocomputing 49 (2002) 389 – 409
(a)
403
(b)
Fig. 7. Separatated signals for (a) synchronized matched 4lter and (b) asynchronous M-EASI based CDMA receiver.
BER
10
10
10
_2
_3
_4
10
15
20
25
Eb/No Fig. 8. BER for the matched 4lter ( ) and the M-EASI-based blind detector (◦) for the user of interest.
be observed that the matched 4lter synchronized at known delays is unable to recover the third user. On the contrary, our approach separates the three users correctly, overcoming the near–far problem, given the adverse circumstances where we do not have information at all about the spreading sequences transmitted. The results shown are for a bit energy-to-noise ratio Eb =No = 10 dB. In Fig. 8, we included the BER versus Eb =No for 4ve users with the same power, the codes and delays were chosen randomly in 200 experiments. The solid lines ( ) depicts the BER for the user of interest at the output of the matched 4lter. We have included in dotted line (◦) the BER for the proposed detector [19].
404
J.J. Murillo-Fuentes, F.J. Gonz"alez-Serrano / Neurocomputing 49 (2002) 389 – 409
6. Conclusions In this paper, the equivariant adaptive separation via independence (EASI) has been revised. We used the natural gradient formulation to provide new results on the method stability. We focused on the mixtures of digital communication signals with additive noise. The new stability conditions allow a better characterization of the performance of the method at separation. Starting from this analysis we propose to apply the Median Gradient introduced by the authors in previous works to improve the performance of the EASI algorithm. We analyze this novel algorithm, the M-EASI, to show that its stability condition has useful properties. Besides, the symmetry conditions imposed by the Median Gradient makes the M-EASI more robust against noise. We 4rst apply this algorithm to a mixture of digitally modulated signals to show how the symmetrical conditions imposed by the Median EASI allow phase recovering, an important issue in communication reception. The extension of these results to beamforming is immediate. Finally, we propose the median EASI as an asynchronous CDMA multiuser detector achieving a similar performance than that of the synchronized matched 4lter.
Acknowledgements We acknowledge 4nancial support from projects CYCIT TIC99-0210 and CAM 07T=0020=2000. The authors would like to thank Antonio Caama˜no Fern'andez for his remarks and observations, and the anonymous reviewers for their helpful comments and suggestions.
Appendix A. A.1. Stability of EASI The terms of dl in (10) may be de4ned as dlW = yT dX y;
(A.1)
dlI = tr(dX );
(A.2)
dl1 = ’(y)T dX y;
(A.3)
dl2 = yT dX ’(y):
(A.4)
Thus, dl for = 1 is the sum of these functions: d 2 l = d 2 lW − d 2 lI + (d 2 l1 − d 2 l2 ):
(A.5)
J.J. Murillo-Fuentes, F.J. Gonz"alez-Serrano / Neurocomputing 49 (2002) 389 – 409
405
We evaluate these terms. The function d 2 lW in terms of dX yields d 2 lW = yT dX T dX y + yT dX dX y
(A.6)
and its expectation is evaluated as E[yT dX T dX y] = E[yi d xji d xjk yk ] = (d xji )2 + (d xii )2 ; ij
E[yT dX dX y] =
E[yi d xij d xjk yk ] =
(A.7)
i
(d xij )(d xji ):
(A.8)
ij
Starting from the terms in (A.7) and (A.8), and for each pair of outputs (i; j), we have ij the quadratic form qW = (d xij + d xji )2 so that d 2 lW may be written as ij d 2 lW = qW + 2(d xii )2 : (A.9) i
i=j
On the other hand, the term d 2 lI cancels. We next evaluate the skew-symmetric part. Its 4rst term yields d 2 l1 = yT dX T ’(y) ˙ dX y + ’(y)T dX dX y;
(A.10)
where ’(y) ˙ is a diagonal matrix whose entries are ’i (yi ). Its expectation is given as follows: E[yT dX T ’(y) ˙ dX y] = E[yi d xji ’j (yj ) d xjk yk ] =
j (d xji )2 +
=
(A.11)
i
j=i
E[’(y)T dX dX y] =
mi (d xii )2 ;
E[’i (yi ) d xij d xjk yk ] i d xij d xji +
i=j
i (d xii )2 ;
(A.12)
i
where mi = E[(yi )2 ’i (yi )], and i and i were de4ned in (15) and (16). The quadratic form of this term may be written as q1ij = j (d xji )2 + i (d xij )2 + ( i + j ) d xij d xji and the term d 2 l1 yields ij d 2 l1 = q1 + (mi + i )(d xii )2 : i=j
(A.13)
(A.14)
i
The last term in (A.5) is evaluated as follows: d 2 l2 = yT dX ’(y) ˙ dX y + ’(y)T dX T dX y:
(A.15)
406
J.J. Murillo-Fuentes, F.J. Gonz"alez-Serrano / Neurocomputing 49 (2002) 389 – 409
By evaluating its expectation we have E[yT dX ’(y) ˙ dX y] = E[yi d xij ’˙ j (yj ) d xjk yk ] =
j d xij d xji +
=
(A.16)
i
j=i
E[’(y)T dX T dX y] =
mi (d xii )2 ;
E[’i (yi ) d xji d xjk yk ] i (d xji )2 +
i=j
i (d xii )2 ;
(A.17)
i
the quadratic form for this term is written as q2ij = (j + i ) d xji d xij + j (d xij )2 + i (d xji )2 and its Hessian ij d 2 l2 = q2 + (mi + i )(d xii )2 : i=j
(A.19)
i
We may rewrite d 2 l in terms of the quadratic forms above as ij d2 l = qEASI + 2(d xii )2 ; ij
(A.18)
(A.20)
i
where ij ij = qW + (q1ij − q2ij ) qEASI
(A.21)
=(i − j + 1)(d xij )2 + (j − i + 1)(d xji )2 +( i + j − j − i + 2) d xji d xij :
(A.22)
ij Finally, qEASI is positive if and only if the following conditions hold:
i − j + 1 ¿ 0;
(A.23)
j − i + 1 ¿ 0;
(A.24)
((i − j ) − (j − i ))2 − 8((i − j ) + (j − i )) ¡ 0:
(A.25)
It can be shown that conditions (A.23) and (A.24), derived from the Sylvester criterion for positive-de4nite quadratic forms, are satis4ed whenever (A.25) is true. Thus, it follows that condition (14) should be met ∀. A.2. Stability of M-EASI The terms of dl in (34) may be de4ned as dlW = yT dX sign(y);
(A.26)
J.J. Murillo-Fuentes, F.J. Gonz"alez-Serrano / Neurocomputing 49 (2002) 389 – 409
407
dlI = tr(dX );
(A.27)
dl1 = ’(y)T dX sign(y);
(A.28)
dl2 = sign(y)T dX ’(y)
(A.29)
so that the Hessian d 2 l, for = 1, yields d 2 l = d 2 lW − d 2 lI + (d 2 l1 − d 2 l2 ):
(A.30)
The pair of terms d 2 lW − d 2 lI were evaluated in [17] for dlMG = ’(y)T dX sign(y). If we set ’(y) = y we have d 2 lW , ij qW + d xii2 ; (A.31) E[d 2 lW − d 2 lI ] = i
i=j
ij where qW = d xij2 + d xji2 is the quadratic form in terms of (d xij ; d xji ) j = i. Regarding the skew-symmetric part, its 4rst term yields
˙ dX y; d 2 l1 = sign(y)T dX T ’(y)
(A.32)
where ’(y) ˙ is a diagonal matrix whose entries are ’i (yi ). Its expectation is evaluated as follows: ˙ dX y] = E[sign(yi ) d xji ’j (yj ) d xjk yk ] E[sign(y)T dX T ’(y) =
j (d xji )2 +
mˆ i (d xii )2 ;
(A.33)
i
j=i
where mˆ i = E[|yi |’i (yi )], and i and i were de4ned in (15) and (16). Hence, ij q1 + mˆ i (d xii )2 ; (A.34) d 2 l1 = i=j
i
where it has been de4ned the quadratic form q1ij = j (d xji )2 + i (d xij )2 . The last term may be evaluated in a similar way, ˙ dX y: d 2 l2 = sign(y)T dX ’(y) Its expectation yields ˙ dX y] = E[sign(y)T dX ’(y) =
j=i
(A.35) E[sign(yi ) d xij ’j (yj ) d xjk yk ] j d xij d xji +
mˆ i (d xii )2
(A.36)
i
and the quadratic form for this term is q2ij = (j + i ) d xji d xij . Thus, d 2 l2 yields ij q2 + mˆ i (d xii )2 : (A.37) d 2 l2 = i=j
i
408
J.J. Murillo-Fuentes, F.J. Gonz"alez-Serrano / Neurocomputing 49 (2002) 389 – 409
Finally, we have the condition ij qM-EASI + 2(d xii )2 ¿ 0; d2 l = ij
(A.38)
i
where ij ij ij ij 2qM -EASI = qW + (q1 − q2 )
= (i + 1)(d xij )2 + (j + 1)(d xji )2 − (j + i ) d xji d xij :
(A.39) (A.40)
Thus, the Hessian is de4nite positive, ∀, if the condition in (35) is true. References [1] S.I. Amari, Neural learning in structured parameter spaces—natural Riemannian gradient, in: M.C. Mozer, M.I. Jordan, T. Petsche (Eds.), Advances in Neural Information Processing Systems, Vol. 9, The MIT Press, Cambridge, 1997, p. 127. [2] S.I. Amari, Natural gradient works eOciently in learning, Neural Comput. 10 (2) (1998) 251–276. [3] S. Amari, SupereOciency in blind source separation, IEEE Trans. Signal Process. 47 (4) (1999) 936–944. [4] S.I. Amari, T. Chen, A. Cichochi, Stability analysis of learning algorithms for blind source separation, Neural Networks 10 (8) (1997) 1345–1351. [5] J.-F. Cardoso, On the stability of source separation algorithms, J. VLSI Signal Process. Systems 26 (1=2) (2000) 7–14. [6] J.F. Cardoso, B.H. Laheld, Equivariant adaptive source separation, IEEE Trans. Signal Process. 44 (12) (1996) 3017–3030. [7] J.F. Cardoso, A. Souloumiac, Blind beamforming for non-Gaussian signals, Proc. IEE F 140 (6) (1993) 362–370. [8] R.T. Causey, J.R. Barry, Blind multiuser detection using linear prediction, IEEE J. Selected Areas Commun. 16 (9) (1998) 1702–1710. [9] A. Cichocki, S. Douglas, S. Amari, Robust techniques for independent component analysis (ICA) with noisy data, Neurocomputing (22) (1998) 113–129. [10] P. Comon, Independent component analysis, a new concept? Signal Process. 36 (3) (1994) 287–314. [11] N. Defolsse, P. Loubaton, Adaptive separation of independent sources: a deWation approach, in: ICASSP, 1994. [12] S.C. Douglas, Self-stabilized gradient algorithms for blind source separation with orthogonality constraints, IEEE Trans. Neural Networks 11 (06) (2000) 1490. [13] I. Gelfand, S. Fomin, Calculus of Variations, Prentice-Hall Inc., Englewood CliJs, NJ, 1963. [14] T.P.V. HoJ, A.G. Lindgren, A. Kaelin, Transpose properties in the stability and performance of the classic adaptive algorithms for blind source separation and deconvolution, Signal Process. (80) (2000) 1807–1822. [15] J. Joutsensalo, T. Ristaniemi, Learning algorithms for blind multiuser detection in CDMA downlink, in: PIMRC, 1998. [16] C. Jutten, J. Herault, Blind separation of sources I. An adaptive algorithm based on neuromimetic architecture, Signal Process. 24 (1) (1991) 1–10. [17] J. Murillo-Fuentes, F. Gonz'alez-Serrano, Higher order moments algorithms for blind signal separation, in: International Workshop on Independent Component Analysis and Blind Signal Separation, Helsinki, Finland, 2000, pp. 345 –350. [18] J. Murillo-Fuentes, F. Gonz'alez-Serrano, Improving stability in blind source separation with the stochastic median gradient, Electron. Lett. 36 (19) (2000) 1662–1663. [19] J.J. Murillo-Fuentes, M. S'anchez-Fern'andez, A. Caamano-Fern'andez, F.J. Gonz'alez-Serrano, Adaptive blind joint source-phase separation in digital communications, in: IEEE International Conference on Communications, Helsinki, Finland, 2001.
J.J. Murillo-Fuentes, F.J. Gonz"alez-Serrano / Neurocomputing 49 (2002) 389 – 409
409
[20] P.G.D.T. Pahm, C. Jutten, Separation of a mixture of independent sources through a ml approach, in: EUSIPCO, 1992, pp. 771, 774. [21] G.T. Ran, Spreading and modulation (FDD), Technical Report, 3GPP, 2000. [22] S. Verd'u, Multiuser Detection, Cambridge University Press, Cambridge, 1998. Juan Jos+e Murillo Fuentes received the telecommunication engineering degree from the University of Seville, Spain, in 1996. He is now pursuing the Ph.D. degree from the University Carlos III de Madrid, Spain. Since 1999, he has been an Assistant Professor at the University Carlos III. His main research interests include neural networks and their application to communications, array processing with emphasis on algorithms for blind separation of sources and multiuser detection, and image processing. Francisco J. Gonz+alez Serrano was born in Betanzos, Spain, on October 10, 1968. He received the telecommunication engineering and Ph.D. degrees from the University of Vigo, Spain, in 1992 and 1997, respectively. From 1992 to 1993 he was with DSP research group, University of Vigo, as a Research Associate. From 1993 to 1998, he was an Assistant Professor with the University Carlos III, Madrid, Spain. His research interest are in the areas of digital signal processing, wireless communications and multimedia communications.