Optimal combination of fourth-order cumulant based contrasts for blind separation of noncircular signals

Signal Processing 93 (2013) 842–855

Contents lists available at SciVerse ScienceDirect

Signal Processing journal homepage: www.elsevier.com/locate/sigpro

Optimal combination of fourth-order cumulant based contrasts for blind separation of noncircular signals$ Christophe De Luigi a,b,n, Eric Moreau a,b a b

Aix Marseille Universite´, CNRS, ENSAM, LSIS, UMR 7296, 13397 Marseille, France Universite´ de Toulon, CNRS, LSIS, UMR 7296, 83957 La Garde, France

a r t i c l e in f o

abstract

Article history: Received 20 December 2011 Received in revised form 27 September 2012 Accepted 8 October 2012 Available online 17 October 2012

In this paper, we have considered the problem of blind source separation of noncircular signals. We have proposed a new Jacobi-like algorithm that achieves optimization of the optimal combination of complex fourth-order cumulant based contrasts. We have investigated the application to separation of non-Gaussian sources using fourth order cumulants. And computer simulations applied to noncircular signals illustrate the proposed approach. & 2012 Elsevier B.V. All rights reserved.

Keywords: High order statistics Contrast Jacobi-like algorithm Blind source separation Noncircular signal

1. Introduction The problem of the blind source separation (BSS) appears in many signal processing application such as telecommunications, RADAR and the bio-medical field. The goal is to restore a set of source signals through observation values of a mixture of these source signals without using explicit knowledge. The BSS problem has given rise to numerous theoretical and applied works, see [1] for a review of recent results and applications. The instantaneous linear mixture model, where no time delays occur in the propagation from sources to sensors, may be very useful and can be accurate in many of BSS applications.

$ Nota bene: Part of this work has been presented at the conference ICA’2007, see [24]. n Corresponding author at: Universite´ de Toulon, CNRS, LSIS, UMR7296, 83957 La Garde, France. Tel.: þ 33 494142128. E-mail addresses: [email protected] (C. De Luigi), [email protected] (E. Moreau).

0165-1684/$ - see front matter & 2012 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.sigpro.2012.10.006

In this paper, we have considered an instantaneous mixing system for which some solutions were proposed in [2–4]. And, we have focused our attention on contrast functions based on high-order statistics. In this field, it was shown that maximizing contrast solves the source separation problem e.g. [5–8]. A previous whitening process helps to reduce the number of unknown parameters and it results in a set of normalized uncorrelated components in relation with the sources through a unitary transformation. The resulting problem is then to identify the unitary matrix. One way to optimize contrast to estimate the unknown unitary mixing matrix is to use Jacobi optimization techniques because of their favorable stability, satisfactory practical convergence and computational parallelism, see [9,10]. Moreover, the analytic solution of the optimal Jacobi angles can be done in the case of two source signals. The Jacobi technique was spread out to the joint approximate diagonalization of matrix sets. It has yielded the two popular algorithms JADE [5] and SOBI [11] while it has led later to the algorithm STOTD [12] and the ones proposed in [6] for tensor sets. A generalization

C. De Luigi, E. Moreau / Signal Processing 93 (2013) 842–855

of these algorithms has been done in [7]. It was shown that optimal combination of contrasts, proposed in several works see [8,7,13,28], allows to improve blind source separation. In [13], in the real valued case, a generalized ICA algorithm, based on Jacobi technique, shows very efficient results optimizing an optimal combination of contrasts based on fourth-order cumulants. In the twosource case, we have derived the large-sample approximations of the mean square error (MSE) in order to reach an optimal weight. Noncircular signals have attracted a high interest these last 20 years, e.g. in wireless telecommunication applications. In this context, more statistics may be used than we can use for circular random variables. See the initial works of Picinbono on circularity [14] and on the second-order statistics of complex signals described by the covariance and the relation functions (the last one is also named pseudo-covariance matrix) [15]. In the case of complex signals, and for a given order, cumulants depend on the considered relative number of complex conjugates. Hence, in [16,17] two fourth-order cumulants based on different statistics are exploited in order to identify the mixing matrix. This is done in the context of blind separation of noncircular sources whatever the number of sources and sensors (within a limited range). In [18,19], it was pointed out the interest of second order optimal widely linear filters for cyclo-stationarity signals under non-circularity conditions. In [20], the covariance and the relation matrices are used in order to perform the MUSIClike estimation of direction of arrival for noncircular sources. Algorithms were proposed in the case of parametric entropy rate estimator using a widely linear autoregressive model [21] and in the case of second order statistics for blind separation of noncircular sources [22]. In [23], the NC-FastICA algorithm was proposed in order to improve the separation of noncircular sources. This algorithm, built under unitary constraint, is based on fixed-point update for the constrained optimization problem. The idea is to take into account the pseudocovariance matrix, that is the second-order noncircular information, which is nonzero if the sources are secondorder noncircular. Several functionals were proposed in the cost function. One of them is even motivated by kurtosis. Nevertheless it does not offer a cost function depending on all the fourth-order cumulants. A comparison to the NC-FastICA algorithm is done in the simulation experiment section. In this paper, we have focused on fourth-order cumulants with complex valued sources. In this context, there are three different fourth-order cumulants without taking into account the conjugate cumulants. The idea is, that with noncircular signal, at least one of the two nonclassical cumulants is nonzero. So its use should lead to the improvement of the performance of the BSS. We have generalized the contrast proposed in [12] by using noncircular statistics. The purpose of this paper is to provide

 a generalized contrast which exploits all the statistics available through the R-order cumulants in the context of blind noncircular signal separation,

843

 an original and optimal combination of the different



statistics obtained through the minimization of the asymptotical error variance on the Jacobi parameters in the context of two sources, an adapted Jacobi-like algorithm in the context of noncircular sources (more than two sources) under unitary constraint.

The paper is organized as follows. First, in Section 2, recalling the BSS model, we have proposed to associate a contrast for each statistic available in the context of complex source signals. Thus, a new contrast is defined based on a combination of the different contrasts previously proposed. In Section 3, we have proposed, for each contrast, a Jacobi-like algorithm in order to optimize the corresponding criterion. And we have presented a generalized ICA algorithm based on the Givens rotation in the case of more than two complex valued sources. In Section 4, a statistical local study of the Jacobi angles is presented in the case of two sources in order to estimate an optimal combination of the available contrasts. It allows us to compare, in Section 5, different algorithms thanks to computer simulations. The latter illustrate the usefulness of considering an optimal combination of the different statistics in fourth-order cumulant based contrasts in the context of noncircular signals even in the case of more than two sources. 2. Blind source separation Let z½n an observed signal vector follows the linear model: z½n ¼ As½n þ n½n,

ð1Þ

where n 2 Z is the discrete time index; s½n the ðN,1Þ vector of N unobservable complex sources si ½n, i 2 f1, . . . ,Ng; z½n the ðM,1Þ vector of observations zi ½n, i 2 f1, . . . ,Mg; n½n the ðM,1Þ vector of additive noises which follow a zero-mean Gaussian random law and A the (M,N) full rank mixing matrix where it is assumed that M Z N. The signals s½n and n½n are assumed statistically mutually independent. The above model is called the instantaneous mixture [2]. We consider that the sources si ½n, with i 2 f1, . . . ,Ng, are zero-mean, unit power, and statistically mutually independent. Here the source signals are assumed complex values. For each source si ½n at least one of the R-th order cumulants, where the optional complex conjugates ðnÞk (k 2 f1, . . . ,Rg) are considered, ðnÞ

ðnÞ

Cumfsi 1 , . . . ,si R g ¼ CðRnÞ fsi g |fflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflffl}

ð2Þ

R terms

is nonzero. The blind source separation problem consists in estimating a matrix H in such a way that the vector Hz½n restores one of the different sources on each of its different components. At this step, a pre-whitening stage is considered, hence the system (1) becomes x½n ¼ Qs½n þ nw ½n,

ð3Þ

844

C. De Luigi, E. Moreau / Signal Processing 93 (2013) 842–855

where the unknown matrix Q is unitary and the vector x½n and nw ½n are obtained from the pre-whitening stage, for a more detailed description see [6,5]. Then, we define the vector y½n by

Using (4), we can see that C4,I ðyÞ, I ¼ 0,1,2, depends on the matrix U. So, the aim is to estimate this unitary matrix U thanks to the CðyðUÞÞ ¼ CðUÞ contrast optimization.

y½n ¼ Ux½n,

3. Jacobi like algorithms

ð4Þ

where U is the unknown unitary matrix that is estimated so that y½n restores the N sources. It is useful to define the global matrix S as S ¼ UQ :

ð5Þ

We classically have proposed to do the (11) contrast optimization through a Jacobi like procedure. In Jacobi like algorithm, optimization is done through a sequence of plane (or Givens) rotations of the unitary matrix U as Y U¼ Ha,b :

Because sources are assumed unobservable, there are some inherent indeterminations in their restitution. That is, why in general, we cannot identify the power and the order of each source. Hence, the separation is done when the matrix S reads

The 12 NðN1Þ matrices Ha,b are identity matrices whose four components in position (a,a), (b,b), (a,b) and (b,a) are respectively replaced by

S ¼ DP,

ðHa,b Þa,a ¼ ðHa,b Þb,b ¼ cosðaa,b Þ,

ð6Þ

where D is an invertible diagonal matrix and P a permutation matrix. Perhaps one of the most useful ways to solve the separation problem is to use contrast functions [6]. They correspond to objective functions which depend on the outputs of the separating system and they have to be maximized to get a separating solution. We then have proposed the following result whose proof is reported in Appendix A. Proposition 2.1. Let R be an integer such that R 43, using the notation ðnÞ1

CðRnÞ fy,i,jg ¼ Cumfyi

ðnÞ2

,yi

ðnÞ3

,yi

ðnÞ4

,yj

1

ðnÞR

, . . . ,yj

R3

g,

ð7Þ

where the optional complex conjugations ðnÞk , k 2 f1, . . . ,Rg, are considered and j ¼ ðj1 , . . . ,jR3 Þ. The function N X

CR ðyÞ ¼

2

9CRðnÞ fy,i,jg9

ð8Þ

i,j1 ,...,jR3 ¼ 1

Considering the particular case of the fourth-order cumulant (R¼ 4), we can write three cumulants for which none of these cumulants is the conjugate of the two others. By convenience, we have used the following notation instead of the one used in (7): C4,0 fy,i,jg ¼ Cumfyi ,yi ,yi ,yj g, C4,1 fy,i,jg ¼ Cumfyi ,yni ,yi ,yj g, C4,2 fy,i,jg ¼ Cumfyi ,yni ,yni ,yj g:

ð9Þ

For each of these three cumulants, a contrast from (8) is associated as N X

2

9C4,I fy,i,jg9 ,

I ¼ 0,1,2:

ð10Þ

i,j ¼ 1

From these three contrasts, we define a new contrast CðyÞ ¼ ðl1 l2 ÞC4,2 ðyÞ þ l1 ð1l2 ÞC4,1 ðyÞ þ ð1l1 ÞC4,0 ðyÞ, T

2

ðHa,b Þa,b ¼ ðHa,b Þnb,a ¼ eıfa,b sin ðaa,b Þ, where ı2 ¼ 1 and aa,b and fa,b are two real-valued angles associated to the matrix Ha,b , 1 r a ob r N. 3.1. The two-source case Using such decomposition reduces the N-dimensional problem into a product of two-dimensional problems. One of the main advantages is that the two-dimensional problem is often simpler. More particularly when it admits an analytical solution. In the two-dimensional problem, only one rotation has to be determined (a ¼ 1, b ¼ 2). It is written in the following as ! cos a eıf sin a U¼ : ð12Þ eıf sin a cos a We have proposed the following result whose proof is reported in Appendix B. Proposition 3.1. When N ¼2, the criteria C4,I ðUÞ, I ¼ 0,1,2, can be written as

is a contrast for white vectors y.

C4,I ðyÞ ¼

1raobrN

ð11Þ

where k ¼ ðl1 l2 Þ 2 ½0 1 . This parametrization has been chosen in order to impose each coefficient to be in ½0 1. The sum of the three coefficients is equal to one.

C4,I ðuÞ ¼ uT B4,I u,

ð13Þ

where uT ¼ ðcos 2a sin 2a sinf sin 2a cos fÞ

ð14Þ

and B4,I , I ¼ 0,1,2, are three real symmetric matrices. The components of the matrices B4,0 and B4,1 can be found in Appendix B and the ones of the matrix B4,2 are defined in [12]. In this proposition, the result concerning the contrast C4,2 ðUÞ was already done in [12]. The maximization of the criteria C4,I ðUÞ, with I ¼ 0,1,2 defined in (13) can be easily found by computing the eigenvector of the matrices B4,I associated with the largest eigenvalue. Then the angles can be obtained e.g. in using qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi cos a ¼ 1 þ 12u1 , ejf sin a ¼

1 ðu3 þ ju2 Þ, 2 cosðaÞ

ð15Þ

with a 2 ½p=4, p=4 and u1, u2 and u3 are the components of vector u.

C. De Luigi, E. Moreau / Signal Processing 93 (2013) 842–855

Let us now consider the contrast CðUÞ defined in (11). Hence in the case of N ¼2, using (13) we directly have CðuÞ ¼ uT Bu,

ð16Þ

where B ¼ ðl1 l2 ÞB4,2 þ l1 ð1l2 ÞB4,1 þð1l1 ÞB4,0 :

ð17Þ

Thus the combination of the three criteria can be considered very easily. 3.2. The more than two-source case In a Jacobi like algorithm, due to its flexibility, the twosource approach may be easily extended to the N-source approach, see [5,6]. For each of the NðN1Þ=2 pairs, the algorithm computes the Givens angles, using (15), from the eigenvector (associated to the largest eigenvalue) of the matrix B defined in (17). The parameter k may be chosen arbitrarily or calculated in an optimal way. For this latter choice, in the next section, an optimal weight is proposed in the two-source case. For more than two sources, and for each pair of sources, it is easy to compute the parameter k through the statistics of the sources. If the statistics are unknown, we do a numerical optimization w.r.t. parameter k (classical and basic numerical method). Several strategies are possible for the first stage. In this paper, we propose to compute at first only the B4,2 matrix that is kð1Þ ¼ ð1 1ÞT . At the end of this first stage, we obtain an estimation of the mixing matrix U and therefore that of the sources and their statistics. At the second stage, the optimal parameter kð2Þ can be computed from these estimated statistics obtained at the end of the first stage.

Ne X b p ½ng ¼ 1 Efy yp ½k, i Ne k ¼ 1 i

ð19Þ

where Ne is the number of available data. The error in the estimation of h~ o resulting from the use b hÞ is denoted e and defined as of the empirical Cð b h~ o , e ¼ ðea ,ef ÞT ¼ h

ð20Þ

where b hÞ: hb ¼ arg max Cð h

So, we have to determine the asymptotical error variance Efe2 g which yields a quality measure of the estimation. It is common to have this variance depending on the statistics of the sources. For the following, we assume that the sources are independent and identically b obtained from (17) depends on distributed. The matrix B the observation vector x½n. Then, using x½n ¼ U1 y½n b will depend on the source from (4), and y½n ¼ s½n, B signals and the statistics of the source signals. b hÞ may be studied through the The maximization of Cð solution of the equation b Þ ¼ 0: bh rCð

ð21Þ

b Þ, bh Using a first order Taylor series expansion of rCð we have ð22Þ

b h~ o Þ and rCð b h~ o Þ are respectively the Hessian and where HCð the gradient of the empirical criterion calculated in h~ o . b to be in the neighborhood of h~ o , we obtain Assuming h

4.1. Asymptotical error variance The general statistical study of the parameter estimation based on the optimization of the contrasts presented in the hereabove section seems to be difficult to achieve. For this reason, we consider a simpler problem: a local asymptotical analysis in the case of two complex and noncircular sources, where the constraints are taken into account through the parametrization of the separating matrix as defined in (12). In this context, we are able to propose a value of the parameter k which allows an optimal combination of the three contrasts C4,I ðuÞ, I ¼ 0,1,2. Let us define the following vector parameter:

h ¼ ða, fÞT ,

now focus on the asymptotical variance property. For this purpose, expectations in the different matrices B40 , B41 , B42 are replaced by sample averages leading to their b 40 , B b 41 , B b 42 . The empirical empirical versions denoted B b hÞ when the different expecversion of CðhÞ is denoted Cð tations are estimated by the classical formula

b h~ o Þ, b Þ ¼ rCð b h~ o Þðh bh b h~ o Þ þ HCð rCð

4. Small error analysis

845

ð18Þ

where a and f are the two angles from (12). We have to estimate the two angles in h according to the maximization of CðuðhÞÞ as a function of h. We consider here that h~ o corresponds to the following case of outputs: y½n ¼ s½n, where s½n is the source vector from (1). That is to say, we consider no noise (or no error) in the separation model (1) in order to simplify the following calculations. The statistical properties of the estimator of h~ o do not depend on a particular optimization method. Hence, we

b h~ o Þe  rCð b h~ o Þ: HCð

ð23Þ

Using Eq. (16), we can give the following expressions of the Hessian and the gradient of the empirical criterion w.r.t. h~ o noted as @uT b Bu, @h T 2 b @u þ 2uT B b @ u , b h~ o Þ ¼ 2 @u B ½HCð i,j @yi @yj @yi @yj

b h~ o Þ ¼ 2 rCð

ð24Þ

with @uT =@h and @2 u=@yi @yj ði ¼ 1,2 and j ¼ 1,2Þ given in Appendix C. 4.2. A local analysis In order to simplify these calculations, we have chosen to do the local analysis around the value y~ o ¼ ð0,0Þ. We can notice that it should be possible to choose any value of y~ o . This choice can be justified by the fact that the estimator of h is obtained by optimizing under unitary constraint a contrast with Ux½n as entry. So, it satisfies the orthogonal invariance property and then the fully

846

C. De Luigi, E. Moreau / Signal Processing 93 (2013) 842–855

invariance property, for that see [25,26]. In the absence of noise, it was shown in [25] that the performance of a fully invariant estimator does not depend on mixing matrix and therefore not on the values of h. It follows the values of the gradient and the Hessian in this local analysis: !T b b h~ o Þ ¼ 4B 1,3 rCð ð25Þ 0

4.3. Optimal coefficients This variance depends on the statistics of the sources and the parameters l1 and l2 . Looking for the minimum of Efe2a g w.r.t. l1 and l2 , it follows a system of two equations to solve: 8 2 > > < l1 l2 JA1 þ2l1 l2 JA2 þ l1 JA3 þ l2 JA4 þJA5 ¼ Efz2b gððJ 4,2 J 4,1 Þl2 þ J4,1 J 4,0 Þ, ð32Þ > > : l1 l2 JA þ l1 JA þ JA ¼ Efz2 gðJ J Þ: 1

and 0 b h~ o Þ ¼ @ HCð

b 1,1 B b 3,3 Þ 8ðB b 1,2 4B

b 1,2 4B 0

1 A:

ð26Þ

For the empirical Hessian, using the large number law, it b 1,2 converges to zero (the sources is easy to show that B b 1,1 B b 3,3 converges are independent and zero-mean) and B to J where J ¼ 34ðl1 l2 J 4,2 þ l1 ð1l2 ÞJ 4,1 þð1l1 ÞJ4,0 Þ,

ð27Þ

with 2

2

J 4,I ¼ 9C4,I fs,1,1g9 þ 9C4,I fs,2,2g9 ,

I ¼ 0,1,2:

ð28Þ

So, we obtain that the empirical Hessian converges to ð8J 0 b 1,3 is given by For the gradient, the expression of B Ne X b 1,3 ¼ 1 B z ½k, Ne k ¼ 1 b

0 0Þ.

ð29Þ

where zb ½k is given in Appendix D. It is easy to see that Efzb g ¼ 0 and thus the central limit pffiffiffiffiffiffi P e theorem ensures that ð1= Ne Þ N k ¼ 1 zb ½k converges (in distribution) to a normal random variable with zero mean and a variance equal to 2 2

2

2

Efz2b g ¼ l1 l2 A1 þ 2l1 l2 A2 þ l1 A3 þ 2l1 l2 A4 þ2l1 A5 þ A6 ,

ð30Þ

where the coefficient AI, I ¼ 0, . . . ,6, are given in Appendix E. In this local analysis, the problem comes to onedimension with a as unknown parameter. It follows that the asymptotical error variance on the parameter a reads 1 1 Efz2b g : Efe2a g  Ne 4 J2

ð31Þ

2

b

4,2

4,1

This system (32) can be rewritten in a system of two polynomial equations in l1 and l2 : ( a1 ðl2 Þl1 þ a0 ðl2 Þ ¼ 0, ð33Þ 2 b2 ðl2 Þl1 þb1 ðl2 Þl1 þ b0 ðl2 Þ ¼ 0, with 8 2 > > > a1 ðl2 Þ ¼ a12 l2 þ a11 l2 þ a10 , > > > > < a0 ðl2 Þ ¼ a01 l2 þ a00 , b2 ðl2 Þ ¼ b21 l2 þ b20 , > > > b1 ðl2 Þ ¼ b11 l2 þ b10 , > > > > : b0 ðl2 Þ ¼ b00 ,

ð34Þ

where the coefficients of (34) are given in Appendix F. When a1 ðl2 Þa0, we substitute l1 from the first equation in the second one from the system (33). We obtain, after some reduction, a polynomial of degree three in the single variable l2 : 3

2

p3 l2 þp2 l2 þp1 l2 þp0 ¼ 0,

ð35Þ

with 8 2 2 > p3 ¼ b00 b21 b11 a00 þb11 b00 a11 b11 b00 b10 , > > > > 2 > > > p2 ¼ 2a00 b00 b21 þb20 b00 a00 b10 b11 a00 a11 b11 > > < a b b þb a2 þ a b b , 11 00 10

00 11

10 00 11

> p1 ¼ a200 b21 þ2a00 b00 b20 a10 a00 b11 a11 a00 b10 > > > > > a10 b00 b10 þ2a11 a10 b00 , > > > > : p0 ¼ a2 b20 þ a2 b00 a10 a00 b10 : 00

ð36Þ

10

The roots of this cubic give l^ 2 and we obtain l^ 1 by

l^ 1 ¼ a0 ðl2 Þ=a1 ðl2 Þ: Remark. In [27], in the circular case, the CRLB for parameters a and f was calculated. It was shown that these two angles are decoupled. To calculate the CRLB in the noncircular case in order to compare the two CRLB is clearly out of the scope of this paper even if it should be certainly interesting. Unfortunately, the calculus of the CRLB in the noncircular case is very difficult to obtain. In the simulation section, we will show that the asymptotic variance obtained for the optimal combination of the three cumulants outperforms the one obtained in the circular case (for which only the classical C4,2 cumulant in nonzero). This simulation allows the reader to better appreciate the potential benefits of exploiting non-circularity.

4

ð37Þ

When a1 ðl2 Þ ¼ 0, we have l2 ¼ a00 =a01 and we obtain l1 by the second equation of the system (33). We can notice that all the coefficients in (36) depend only on the statistics of the sources. In order to illustrate these results, some computer simulations will be presented in the next section. Remark. In some cases, it might occur that the roots of the system do not fall in the permissible range ½0 12 . When it happens, one efficient strategy is to reduce the combination of three contrasts to two contrasts. In the next subsection, the optimal coefficient of the combination of two contrasts is directly calculated (without passing by any root of a system).

C. De Luigi, E. Moreau / Signal Processing 93 (2013) 842–855

4.4. Two particular cases We close this section with two particular cases where two contrasts are combined. In these cases, a theoretical result of the optimal coefficient can be easily obtained. When l1 ¼ 1, we obtain J ¼ l2 J 4,2 þ ð1l2 ÞJ 4,1 :

ð38Þ

Looking for the minimum of E½e2a , we obtain an equation in l2 to be solved. It follows the optimal coefficient l^ 2 which only depends on the statistics of the sources, that is: J 4,2 G1 J 4,1 G2,1 , J 4,1 ðG2 G2,1 Þ þ J 4,2 ðG1 G2,1 Þ

ð39Þ

where the coefficients G1 , G2 and G2,1 are given in Appendix G. Similarly when l2 ¼ 1, we have 2

Efz2b g ¼ l1 ðA1 þ2A2 þA3 Þ þ 2l1 ðA4 þA5 Þ þ A6 , J ¼ l1 J 4,2 þ ð1l1 ÞJ 4,0 :

ð40Þ

Hence, we obtain l^ 1

l^ 1 ¼

J 4,2 G0 J4,0 G2,0 , J 4,0 ðG2 G2,0 Þ þ J 4,2 ðG0 G2,0 Þ

where the coefficients G0 Appendix G.

and G2,0

ð41Þ are given in

5. Computer simulations In our experiments, we consider the following complex noncircular source signals.

 S1: 5-states modulation: { 1 þ1ı, 1 1ı, 0, bð11ıÞ,  

This non-negative index is zero if S satisfies (6). A small value indicates the proximity to the desired solutions. The performances of the presented algorithms are evaluated by Monte Carlo simulations in which we average over 500 runs. 5.1. Numerical analysis of the optimal parameter

2

Efz2b g ¼ l2 A1 þ 2l2 ðA2 þ A4 Þ þ ðA3 þ2A5 þA6 Þ,

l^ 2 ¼

847

bð1 þ 1ıÞ}, with the probabilities f1=2ð1 þ bÞ,1=2 ð1þ bÞ, ðb1Þ=b,1=2bð1 þ bÞ,1=2bð1 þ bÞg: S2: 8-states modulation: f1þ 1i,11i,1 þ1i,1 1i, 3 þ1i,31i,3 þ 1i,31ig, with uniform probabilities. S3: QPSK modulation: f1 þ1i,11i,1 þ 1i,11ig, with uniform probabilities.

For S1, the parameter b may be chosen so that 9C4,0 ðsÞ9 4 9C2,2 ðsÞ9 (we choose b ¼ 1:4). Four algorithms will be compared:

 G-STOTD algorithm which corresponds to the optimization of CðuÞ as defined by (16) and (17).

 STOTD algorithm from [12] (case with l1 ¼ l2 ¼ 1Þ.  NC-STOTD algorithm (case with l1 ¼ 1 and l2 ¼ 0Þ.  NC-FastICA from [23]. The performances of the algorithms are associated with the following error index [7] defined on the global matrix S as 0 0 1 2 N N X X 9ðSÞi,j 9 1 @ @ IðSÞ ¼ 1A NðN1Þ i ¼ 1 j ¼ 1 max‘ 9ðSÞ 92 i,‘ !1 2 N N X X 9ðSÞi,j 9 ð42Þ þ 1 A: 2 j ¼ 1 i ¼ 1 max‘ 9ðSÞ‘,j 9

The objective is to emphasize the existence of a parameter for the contrast as defined in (11) which allows an optimal performance of the general algorithm. First, we have studied the two particular situations mentioned in Section 4.4. In these two cases, in order to numerically find the optimal coefficient, we sample l1 in the first case and l2 in the second case from 0 to 1 with a 1 100 step. Finally, for each case, we find a numerical value of the parameter that gives the minimum error index. And, we compare this value with the theoretical optimal coefficient found by (39) and (41). At each process, Ne ¼ 5  104 , we keep the same unitary mixing (no pre-whitening step is needed) with N ¼ 2, a ¼ p=3, f ¼ 2p=3. The two sources follow the 5 state modulation S1. In Figs. 1 and 2 we respectively plot the error index for the G-STOTD in cases l1 ¼ 1 and l2 ¼ 1. We have added an arrow on these figures whose head is positioned to indicate the theoretical optimal value obtained respectively from (39) and (41). We can see in the two cases that there exists k^ in CðuÞ so that the G-STOTD performs the STOTD [12]. Moreover, when we compare the coefficients obtained numerically and theoretically, they are very close. Second, we have studied the two dimensional optimization problem (16) and (17). To numerically find the optimal coefficient, we have sampled l1 and l2 from 0 to 1 1 with a 100 step. In Fig. 3, we have plotted the curves of the error index versus l1 and l2 for one Monte Carlo run. We have added an arrow on the figure whose head is positioned to indicate the theoretical optimal value obtained from the system (33). Regarding the average on 500 runs, we have obtained the numerical solution whose value is ð0:83,0:75Þ. The theoretical solution from the system (33) has the value ð0:83,0:76Þ that is to say very close to the numerical solution value. 5.2. Asymptotical optimal variance The aim is to show that the asymptotical variances obtained in Section 4 are good enough estimations. Hence, we emphasize the benefits taking into account to noncircular fourth-order statistics. At each process, we have kept the same unitary mixing (no pre-whitening step is needed) with N ¼ 2, a ¼ p=3, f ¼ 2p=3. The two sources follow the 5 state modulation S1. For each simulation, we have given the empirical variance, obtained by Monte-Carlo simulations, and the asymptotical variance, obtained from Efe2a g in (31), for each criterion. The number of samples goes from 5000 to 500 000 in order to show the effect of Ne . In Fig. 4, we have shown the results: the performances obtained by the G-STOTD algorithm are better than the one obtained by the STOTD criterion. In both scenarios, a

848

C. De Luigi, E. Moreau / Signal Processing 93 (2013) 842–855

Optimal Combination of two criteria 10−4

10−5

Error Index

λ2 = 0,78 10−6

10−7

10−8

10−9 0

0.1

0.2

0.3

0.4

0.5 λ2

0.6

0.7

0.8

0.9

1

Fig. 1. Performance of the G-STOTD algorithm versus the sampling of l for S1 when l1 ¼ 1.

Optimal Combination of two criteria

Error Index

10−4.8132

λ1 = 0.6 10−4.8133

10−4.8134 0.4

0.5

0.6

0.7

0.8

0.9

1

λ1 Fig. 2. Performance of the G-STOTD algorithm versus the sampling of l for S1 when l2 ¼ 1.

better estimation of a is obtained with the optimal algorithm. So, combining the different statistics in an optimal way allows to improve the results of the classical algorithms in case of noncircular sources. 5.3. The two-source case Here, we have compared the performances of the three algorithms STOTD, NC-STOTD and G-STOTD. For G-STOTD, parameters l2 and l1 were chosen respectively by the roots of (35) and by (37). We have considered successively two cases of mixture. In the first one, the mixing matrix is a unitary with N ¼2 sources and a ¼ p=3, f ¼ 2p=3. In the second case, the dimensions of the mixing matrix are N ¼2 and M¼4. Its components follow a Gaussian law with zero-mean and

unitary power. Hence, a pre-whitening step is done. The two sources follow the 5 state modulation S1. First, we have compared the three algorithms w.r.t. Ne in a noiseless context. In the two cases of mixture, the results are given in Tables 1 and 2. The best performances are obtained with G-STOTD whatever Ne. Without any surprise, the unitary mixing scenario shows better results than in the nonunitary mixture scenario. Second, we compare the three algorithms, with Ne ¼ 5  104 samples, in a context of added Gaussian noise on the mixture signal. The power of noise, noted pw, takes the values 0.125, 0.25, 0.5, 0.75 and 1. The results are given in Tables 3 and 4. Once again, it is shown that the best performances are obtained with G-STOTD whatever the power of noise, especially in the unitary mixing scenario.

C. De Luigi, E. Moreau / Signal Processing 93 (2013) 842–855

849

x 10−8 10

Optimal Combination of three criteria

9

0.95 (0.76,0.83)

8

0.9

7 6

0.85 λ1

5 0.8

4 3

0.75

2

0.7

1 0.65 0.1

0.2

0.3

0.4

0.5 λ2

0.6

0.7

0.8

0.9

1

Fig. 3. Performance of the G-STOTD algorithm versus the sampling of l1 and l2 for S1 .

10−4

Comparison of the different algorithms for the variance of the parameter alpha

10−5

Variance of alpha

10−6 10−7 10−8 10−9 10−10 10−11 10−12

Algo NC−STOTD Asymp NC−STOTD Algo STOTD Asymp STOTD Algo G−STOTD Asymp G−STOTD

104

105 Number of samples

Fig. 4. Asymptotical and empirical variance of the different ICA criteria versus the number of samples. Table 1 Performances of the algorithms w.r.t. the number of samples for a fixed unitary mixture in the two-source case.

Table 2 Performances of the algorithms w.r.t. the number of samples for a nonunitary mixture in the two-source case.

Algorithms

STOTD

NC-STOTD

G-STOTD

Algorithms

STOTD

NC-STOTD

G-STOTD

N e ¼ 2  103

4:2  104

3:6  105

3:2  107

N e ¼ 2  103

6:5  104

2:8  104

2:4  104

N e ¼ 1  104

7:7  105

6:3  106

1:5  108

N e ¼ 1  104

1:2  104

5:7  105

5:1  105

N e ¼ 4  104

1:8  105

1:5  106

9:2  1010

N e ¼ 4  104

3:2  105

1:5  105

1:3  105

N e ¼ 6  104

1:3  105

1:0  106

4:2  1010

N e ¼ 6  :104

2:0  105

9:1  106

8:0  106

N e ¼ 1  105

7:3  106

5:8  107

1:5  1010

N e ¼ 1  105

1:2  105

5:3  106

4:7  106

5.4. The more than two-source case In this experiment we have compared the performances of the three algorithms STOTD, G-STOTD and

NC-FastICA versus the number of sources in a noiseless context. The mixing matrix is unitary and random, the number of samples is Ne ¼ 105 . We have tested the performance of the algorithms in three contexts of

850

C. De Luigi, E. Moreau / Signal Processing 93 (2013) 842–855

Table 3 Performances of the algorithms w.r.t. the power of noise for a fixed unitary mixture in the two-source case. Algorithms

STOTD

NC-STOTD

G-STOTD

pw ¼ 0:125

1:9  105

2:3  106

2:4  107

pw ¼ 0:250

2:9  10

5

6

5:6  10

1:1  10

6

pw ¼ 0:500

9:2  105

2:3  105

5:1  106

pw ¼ 0:750

2:2  104

6:5  105

1:6  105

pw ¼ 1:000

4

4

4:3  105

5:7  10

1:5  10

Table 4 Performances of the algorithms w.r.t. the power of noise for a nonunitary mixture in the two-source case. Algorithms

STOTD

NC-STOTD

G-STOTD

pw ¼ 0:125

2:6  105

8:2  106

3:0  106

pw ¼ 0:250

4:0  105

1:2  105

4:4  106

pw ¼ 0:500

1:2  104

3:4  105

1:1  105

pw ¼ 0:750

4:3  104

8:5  105

2:6  105

pw ¼ 1:000

7:7  104

2:6  104

7:3  105

different source modulation in order to show the influence of the characteristics of the sources on the performance of the algorithms.

 For S1 with b ¼ 1:4: 9C4,0 ðsÞ9 ¼ 1:56, 9C2,2 ðsÞ9 ¼ 0:44,  

9C3,1 ðsÞ9 ¼ 0, and EfS21 g ¼ 0. For S2: 9C4,0 ðsÞ9 ¼ 9C2,2 ðsÞ9 ¼ 36, EfS22 ga0. For S3: 9C4,0 ðsÞ9 ¼ 9C2,2 ðsÞ9 ¼ 4, EfS23 g ¼ 0.

9C3,1 ðsÞ9 ¼ 32,

and

9C3,1 ðsÞ9 ¼ 0,

and

For each source, the stability, described in [23], is guaranteed. The results are given in Tables 5, 6 and 7 for each source. First, we can see that the proposed algorithm G-STOTD outperforms STOTD and NC-FastICA whatever the source. For S1 and S3, the difference of performance becomes important enough. It may come from the asymmetric values of the three contrasts. Indeed, for S2, since all the cumulants have approximatively the same absolute value, G-STOTD and STOTD have the same order of performance (106 for both). 6. Conclusion

Table 5 Performances of the algorithms w.r.t. the number of sources for a random unitary mixture with S1 and N e ¼ 105 samples. Algorithms

STOTD

NC-FastICA

G-STOTD

N¼2

1:4  105

1:9  105

6:7  1010

N¼4

1:1  10

5

5

1:6  10

7:6  1010

N¼6

1:1  105

1:6  105

9:2  1010

N¼8

9:2  10

6

5

1:2  10

6:0  1010

N ¼ 10

1:0  105

1:2  105

8:2  1010

The optimization, via a Jacobi-like algorithm, of an optimal combination of complex fourth-order cumulant based contrasts has been presented. In the specific case of two noncircular sources, a local asymptotical analysis of the estimated parameter has been derived. The error variance was shown to depend on the statistics of the sources. Minimization of this variance has allowed us to derive the optimal coefficient. Finally, computer simulations have shown the interest of considering optimal combination of all fourth-order statistics in the context of noncircular sources even when there were more than two sources. Appendix A. Proof of Proposition 2.1

Table 6 Performances of the algorithms w.r.t. the number of sources for a random unitary mixture with S2 and N e ¼ 105 samples. Algorithms

STOTD

NC-FastICA

G-STOTD

N¼2 N¼4

2:7  106

7:1  106

2:1  106

6

6

2:3  106

6

6

3:1  10

7:4  10

N¼6

3:0  10

8:2  10

2:2  10

N¼8

4:0  106

7:4  106

2:2  106

N ¼ 10

3:8  106

8:5  106

2:3  106

6

In this proof, we consider R¼4 for simplicity. One easily has for all permutation matrix P, C4 ðPyÞ ¼ C4 ðyÞ and for all orthonormal diagonal matrix D, C4 ðDyÞ ¼ C4 ðyÞ. Let the two following propositions to help the proof. Proposition A.1. For any statistically independent random vector a, any orthonormal matrix S, we have C4 ðSaÞ rC4 ðaÞ:

ð43Þ

Proof. With S ¼ ðSi,j Þ, the multi-linearity of cumulants and the independence of sources, we have Table 7 Performances of the algorithms w.r.t. the number of sources for a

Cð4nÞ fS,i,jg ¼

N X

ðnÞ

ðnÞ

ðnÞ

ðnÞ4

Sil 1 ,Sil 2 ,Sil 3 ,Sjl

C4,c fa‘ g,

‘¼1

5

random unitary mixture with S3 and N e ¼ 10 samples. Algorithms

STOTD

N¼2 N¼4

where

NC-FastICA

G-STOTD

1:8  1010

5:2  106

3:0  1014

10

6

4:3  1014

6

14

2:2  10

5:4  10

N¼6

2:4  10

4:9  10

4:5  10

N¼8

2:5  1010

5:1  106

5:2  1014

N ¼ 10

2:6  1010

5:4  106

5:8  1014

10

C4,c fa‘ g ¼ Cumfan‘ , . . . ,an‘ ,a‘ , . . . ,a‘ g: |fflfflfflfflfflffl{zfflfflfflfflfflffl} |fflfflfflfflffl{zfflfflfflfflffl} c terms

4-c terms

Now because S is a unitary matrix then 8‘1 ,‘2 ,

N X m¼1

Sm,‘1 Snm,‘2 ¼ d‘1 ,‘2 ,

C. De Luigi, E. Moreau / Signal Processing 93 (2013) 842–855

P n 0 Using j ujk4 ujk04 ¼ dk4 ,k04 , we have k4 ¼ k4 and we replace k4 and k04 by j, we obtain X aj , C4,0 ðyÞ ¼

where di,j ¼ 1 if i¼j and 0 otherwise. Thus ! N N X X 4 2 9Si,‘ 9 9C4,c fa‘ g9 C4,c ðyÞ ¼ ‘¼1

i¼1

j

and because 8‘1 ,

N X

4

9Si1 ,‘1 9 r

i1 ¼ 1

851

N X

where X aj ¼

2

9Si1 ,‘1 9 ¼ 1,

X

uik1 uik2 uik3 unik0 unik0 unik0 1

i k1 ,k2 ,k3 ,k01 ,k02 ,k03

i1 ¼ 1

n

C4,c ðyÞ r

n

3

n

Cumfxk1 ,xk2 ,xk3 ,xj gCumfxk0 ,xk0 ,xk0 ,xj g:

then we have N X

2

n

1

2

3

Using the values of the matrix U from the Jacobi model (12), we regroup the 64 cumulants terms as

2

9C4,c fa‘ g9 ¼ C4,c ðaÞ:

‘¼1

2

where C4,c ðyÞ ¼ C4 ðyÞ with c complex conjugates.

2

þsin6 ðaÞÞ þ 2:25ð9Cumfx1 ,x1 ,x2 ,xj g9

Proposition A.2. For any statistically independent random vector a having at most one null cumulant of 4-order, C4 ðSaÞ ¼ C4 ðaÞ if and only if S ¼ DP where P is a permutation 2 and D ¼ diagðd1 , . . . ,dN Þ such that 8i, 9di 9 ¼ 1. Proof. The equality in (43) requires the equalities PN 4 ðnÞ i ¼ 1 9sij 9 ¼ 1 for all j such that 9C4 fa,i,jg9a0. Since 8j, PN 2 i ¼ 1 9sij 9 ¼ 1, this is possible if and only if columns j of S have only one nonzero component of modulus 1. That is for at least N1 columns because we assume that a has at most one null cumulant of 4-order. Because S is orthonormal then S ¼ DP where P is a permutation and 2

D ¼ diagðd1 , . . . ,dN Þ such that 8i, 9di 9 ¼ 1.

2

aj ¼ ð9Cumfx1 ,x1 ,x1 ,xj g9 þ 9Cumfx2 ,x2 ,x2 ,xj g9 Þðcos6 ðaÞ

&

&

2

2

þ9Cumfx2 ,x2 ,x1 ,xj g9 Þsin ð2aÞ þ3 RefCumfx2 ,x2 ,x2 ,xj gCumfx1 ,x2 ,x2 ,xj gn eıf Cumfx1 ,x1 ,x1 ,xj gCumfx1 ,x1 ,x2 ,xj gn eıf gsinð2aÞcosð2aÞ þ1:5RefCumfx1 ,x1 ,x1 ,xj gCumfx1 ,x2 ,x2 ,xj gn eı2f þCumfx2 ,x2 ,x2 ,xj gCumfx1 ,x1 ,x2 ,xj gn eı2f g sin2 ð2aÞ, which can be written as aj ¼ uT Bð0Þ ðjÞu, where u is a real ð3  1Þ vector defined in (14) and such that uT u ¼ 1. Bð0Þ ðjÞ is a real symmetric matrix ð3  3Þ for each j and defined by 2

When N ¼2, using (10) and (9), the criteria C4,0 ðyÞ can be written as C4,0 ðyÞ ¼

2 X 2 X

2 X 2 X

k¼1

k¼1

X

n

n

n

1

C4,0 ðyÞ ¼

2

1

2

X

2

3

4

uik1 uik2 uik3 unik0 unik0 unik0 1

i k1 ,k2 ,k3 ,k4 ,k01 ,k02 ,k03 ,k04

2

3

X Cumfxk1 ,xk2 ,xk3 ,xk4 gCumfxnk0 ,xnk0 ,xnk0 ,xnk0 g ujk4 unjk0 : 1

2

3

4

j

4

Cumfx2 ,x2 ,x2 ,xj gCumfx1 ,x1 ,x2 ,xj gn g, 2

2

ð0Þ B3,3 ðjÞ ¼ 0:25ð9Cumfx1 ,x1 ,x1 ,xj g9 þ 9Cumfx2 ,x2 ,x2 ,xj g9 Þ 2

2

þ 2:25ð9Cumfx1 ,x1 ,x2 ,xj g9 þ 9Cumfx1 ,x2 ,x2 ,xj g9 Þ þ 1:5 RefCumfx1 ,x1 ,x1 ,xj gCumfx1 ,x2 ,x2 ,xj gn þ Cumfx2 ,x2 ,x2 ,xj gCumfx1 ,x1 ,x2 ,xj gn g:

n

uik1 uik2 uik3 ujk4 uik0 uik0 uik0 ujk0

j k1 ,k2 ,k3 ,k4 ,k01 ,k02 ,k03 ,k04

Cumfxk1 ,xk2 ,xk3 ,xk4 gCumfxnk0 ,xnk0 ,xnk0 ,xnk0 g, X

2

þ 2:25ð9Cumfx1 ,x1 ,x2 ,xj g9 þ 9Cumfx1 ,x2 ,x2 ,xj g9 Þ

ð0Þ ðjÞ ¼ 1:5 ImfCumfx1 ,x1 ,x1 ,xj gCumfx1 ,x2 ,x2 ,xj gn B2,3

k¼1

With the multi-linearity of the cumulants and k1 ,k2 ,k3 ,k4 , k01 ,k02 ,k03 ,k04 2 f1,2g we can write  2  XX X    u u u u Cumfx ,x ,x ,x g C4,0 ðyÞ ¼ ik1 ik2 ik3 jk4 k1 k2 k3 k4  ,   i j k1 ,k2 ,k3 ,k4

i

2

þ Cumfx2 ,x2 ,x2 ,xj gCumfx1 ,x1 ,x2 ,xj gn g,

 ( )2   2 2 2 2 X X X X   uik xk , uik xk , uik xk , ujk xk  : Cum  

C4,0 ðyÞ ¼

2

ð0Þ ðjÞ ¼ 0:25ð9Cumfx1 ,x1 ,x1 ,xj g9 þ 9Cumfx2 ,x2 ,x2 ,xj g9 Þ B2,2

1:5 RefCumfx1 ,x1 ,x1 ,xj gCumfx1 ,x2 ,x2 ,xj gn

i¼1j¼1

XX

þ Cumfx1 ,x1 ,x1 ,xj gCumfx1 ,x1 ,x2 ,xj gn g,

Cumfx1 ,x1 ,x1 ,xj gCumfx1 ,x1 ,x2 ,xj gn g,

Using y½n ¼ Ux½n with the Jacobi model (12), that is P yi ½n ¼ 2k ¼ 1 uik xk ½n, we obtain

k¼1

ð0Þ B1,2 ðjÞ ¼ 1:5 ImfCumfx2 ,x2 ,x2 ,xj gCumfx1 ,x2 ,x2 ,xj gn

ð0Þ B1,3 ðjÞ ¼ 1:5 RefCumfx2 ,x2 ,x2 ,xj gCumfx1 ,x2 ,x2 ,xj gn

2

9Cumfyi ,yi ,yi ,yj g9 :

i¼1j¼1

C4,0 ðyÞ ¼

2

ð0Þ ðjÞ ¼ 9Cumfx1 ,x1 ,x1 ,xj g9 þ 9Cumfx2 ,x2 ,x2 ,xj g9 , B1,1

Appendix B. Proof of Proposition 3.1

3

4

So, the maximization of the contrast C4,0 is obtained by C4,0 ¼ uT B4,0 u where the real matrix B4,0 is defined by P B4,0 ¼ 2j ¼ 1 Bð0Þ ðjÞ: Using the same reasoning, it is easy to show that C4,1 ðyÞ ¼ uT B4,1 u, where the real matrix B4,1 is defined P by B4,1 ¼ 2j ¼ 1 Bð1Þ ðjÞ with Bð1Þ ðjÞ a real symmetric matrix ð3  3Þ for each j defined by 2

2

n n Bð1Þ 1,1 ðjÞ ¼ 9Cumfx1 ,x1 ,x1 ,xj g9 þ 9Cumfx2 ,x2 ,x2 ,xj g9 ,

852

C. De Luigi, E. Moreau / Signal Processing 93 (2013) 842–855

Appendix D. Expression of B^ 13

n n n Bð1Þ 1,2 ðjÞ ¼ ImfCumfx1 ,x2 ,x2 ,xj gCumfx2 ,x2 ,x2 ,xj g , n

n

n

0:5 Cumfx2 ,x1 ,x2 ,xj gCumfx2 ,x2 ,x2 ,xj g g

b 1,3 , we find In order to obtain the empirical version of B the empirical statistics thanks to some approximations:

þ Imf0:5 Cumfx1 ,xn1 ,x1 ,xj gCumfx1 ,xn2 ,x1 ,xj gn Cumfx1 ,xn1 ,x1 ,xj gCumfx2 ,xn1 ,x1 ,xj gn g,

Ne Ne 1 X 1 X z~ p ½k y~ pi ½k Ne k ¼ 1 Ne k ¼ 1 i

n n n Bð1Þ 1,3 ðjÞ ¼ RefCumfx1 ,x2 ,x2 ,xj gCumfx2 ,x2 ,x2 ,xj g

is approached by

þ 0:5 Cumfx2 ,xn1 ,x2 ,xj gCumfx2 ,xn2 ,x2 ,xj gn g n

n

n

 Ref0:5 Cumfx1 ,x1 ,x1 ,xj gCumfx1 ,x2 ,x1 ,xj g n

n

Efy~ pi g

n

þ Cumfx1 ,x1 ,x1 ,xj gCumfx2 ,x1 ,x1 ,xj g g, 2

For example, the statistic Cumfsi1 ,si2 ,si3 ,sj g is approached by

2

n n Bð1Þ 2,2 ðjÞ ¼ 0:25ð9Cumfx1 ,x1 ,x1 ,xj g9 þ 9Cumfx2 ,x2 ,x2 ,xj g9

2

Ne Ne 1 X 1 X si1 ½ksi2 ½ksi3 ½ksj ½kEfsi1 si2 g s ½ksj ½k Ne k ¼ 1 N e k ¼ 1 i3

2

þ 9Cumfx1 ,xn2 ,x1 ,xj g9 þ9Cumfx2 ,xn1 ,x2 ,xj g9 Þ 2

2

þ 9Cumfx2 ,xn1 ,x1 ,xj g9 þ9Cumfx1 ,xn2 ,x2 ,xj g9 n

n

þ RefCumfx1 ,x1 ,x1 ,xj gCumfx1 ,x2 ,x2 ,xj g

Efsi1 si3 g

n

0:5 Cumfx1 ,xn1 ,x1 ,xj gCumfx2 ,xn1 ,x2 ,xj gn g 0:5 Cumfx1 ,xn2 ,x1 ,xj gCumfx2 ,xn2 ,x2 ,xj gn g

zb ½k ¼ ðl1 l2 Þ z2 ½k þ l1 ð1l2 Þ z1 ½k þ ð1l1 Þ z0 ½k

RefCumfx1 ,xn1 ,x2 ,xj gCumfx1 ,xn2 ,x1 ,xj gn

and the coefficients zI ½k, I ¼ 0,1,2, are given by

þ Cumfx2 ,xn1 ,x2 ,xj gCumfx2 ,xn2 ,x1 ,xj gn g,

zI ½k ¼ RefwI ½kg,

n n n Bð1Þ 2,3 ðjÞ ¼ ImfCumfx2 ,x1 ,x2 ,xj gCumfx2 ,x2 ,x1 ,xj g

þ 0:5 Cumfx1 ,x1 ,x1 ,xj gCumfx2 ,x1 ,x2 ,xj g n

n

w0 ½k ¼ 32C4,0 fs,2,2gðs1 ½ks32 ½k3Efs22 gs1 ½ks2 ½kÞn

n

32C4,0 fs,1,1gðs2 ½ks31 ½k3Efs21 gs1 ½ks2 ½kÞn ,

2

w1 ½k ¼ C4,1 fsn ,2,2gð12 sn1 ½ks32 ½k32Efs22 gsn1 ½ks2 ½k 2

Efs22 gs1 ½ksn2 ½k þ s1 ½ksn2 ½ks22 ½k2Ef9s22 9gs1 ½ks2 ½kÞ

0:25ð9Cumfx1 ,x1 ,x1 ,xj g9 þ 9Cumfx2 ,x2 ,x2 ,xj g9 n

n

2

2

C4,1 fsn ,1,1gð12 sn2 ½ks31 ½k32Efs21 gs1 ½ksn2 ½k

þ 9Cumfx2 ,xn1 ,x2 ,xj g9 þ 9Cumfx1 ,xn2 ,x1 ,xj g9 Þ 2

Efs21 gsn1 ½ks2 ½k þ s21 ½ksn1 ½ks2 ½k2Ef9s21 9gs1 ½ks2 ½kÞ

2

þ 9Cumfx2 ,xn1 ,x1 ,xj g9 þ9Cumfx1 ,xn2 ,x2 ,xj g9

ð46Þ

þ RefCumfx1 ,xn1 ,x1 ,xj gCumfx1 ,xn2 ,x2 ,xj gn and

þ 0:5 Cumfx1 ,xn1 ,x1 ,xj gCumfx2 ,xn1 ,x2 ,xj gn g n

n

þ RefCumfx2 ,x1 ,x1 ,xj gCumfx2 ,x2 ,x2 ,xj g

w2 ½k ¼ C4,2 fs,2,2gðs1 ½ks2 ½ksn2 2½k2Ef9s22 9gs1 ½ksn2 ½k

n

þ 0:5 Cumfx1 ,x2 ,x1 ,xj gCumfx2 ,x2 ,x2 ,xj g g

En fs22 gs1 ½ks2 ½k þ 12sn1 ½ksn2 ½ks22 ½k

þ RefCumfx1 ,xn1 ,x2 ,xj gCumfx1 ,xn2 ,x1 ,xj gn

Ef9s22 9gsn1 ½ks2 ½k12Efs22 gsn1 ½ksn2 ½kÞ

n

n

n

ð45Þ

n

 Cumfx1 ,x1 ,x2 ,xj gCumfx1 ,x2 ,x1 ,xj g g, Bð1Þ 3,3 ðjÞ ¼

ð44Þ

with

þ 0:5 Cumfx1 ,xn2 ,x1 ,xj gCumfx2 ,xn2 ,x2 ,xj gn n

Ne Ne 1 X 1 X si2 ½ksj ½kEfsi1 sj g s ½ksi3 ½k: Ne k ¼ 1 N e k ¼ 1 i2

b 1,3 ¼ ð1=Ne Þ PNe zb ½k where For recall (29) gives: B k¼1 the signal zb ½k is defined as

þ RefCumfx2 ,xn1 ,x1 ,xj gCumfx2 ,xn2 ,x2 ,xj gn

n

Ne 1 X p z~ ½k: Ne k ¼ 1 i

n

n

C4,2 fs,1,1gðs1 ½ks2 ½ksn1 2½k2Ef9s21 9gs2 ½ksn1 ½k

n

þ Cumfx2 ,x1 ,x2 ,xj gCumfx2 ,x2 ,x1 ,xj g g:

En fs21 gs1 ½ks2 ½k þ 12 sn1 ½ksn2 ½ks21 ½k Ef9s21 9gs1 ½ksn2 ½k12Efs21 gsn1 ½ksn2 ½kÞ:

ð47Þ

Appendix C. Gradient and Hessian of u We have uT ¼ ðcos 2a sin 2a sin f sin 2a cos fÞ and

h ¼ ða fÞT , so the gradient is @uT ¼ @h

2 sin 2a

2 cos 2a sin f

2 cos 2a cos f

0

sin 2a cos f

sin 2a sin f

!

T

@ u ¼ @h@y1

2 2 Efz2b g ¼ l1 l2

2

2

A1 þ 2l1 l2 A2 þ l1 A3 þ2l1 l2 A4 þ 2l1 A5 þ A6 ,

where the coefficients AI, I ¼ 0, . . . ,6, are given by

4 cos 2a

4 sin 2a sin f

4 sin 2a cos f

0

2 cos 2a cos f

2 cos 2a sin f

and for j¼ 2, it is @2 u T ¼ @h@y2

From Appendix D, it is easy to obtain

:

For j ¼1, the Hessian is 2

Appendix E. Expressions of the coefficients in (30)

0

2 cos 2a cos f

2cos 2a sin f

0

sin 2a sin f

sin 2a cos f

! :

! ,

A1 ¼ Efz21 g þ Efz22 g2Efz2 z1 g, A2 ¼ Efz21 g þEfz2 z1 gEfz0 z2 g þ Efz1 z0 g, A3 ¼ Efz21 g þ Efz20 g2Efz1 z0 g, A4 ¼ Efz0 z2 gEfz1 z0 g, A5 ¼ Efz20 g þ Efz1 z0 g, A6 ¼ Efz20 g:

C. De Luigi, E. Moreau / Signal Processing 93 (2013) 842–855

The expected values, Efz2I g, I ¼ 0,1,2, can be obtained from the definition of zI ½k that is using (44)–(47) and the following result: 2

Efz2I g ¼ EfRfwI g2 g ¼ 12ðRefEfw2I gg þEf9wI 9 gÞ:

ð48Þ

with

k1,1 ¼ 12C3,1 fsn ,2,2g, k1,2 ¼ 32C3,1 fsn ,1,1gEfs21 gC3,1 fsn ,2,2gEfs22 g 2

k1,3 ¼ 2C3,1 fsn ,1,1gEf9s1 92 g2C3,1 fsn ,2,2gEf9s2 92 g, k1,4 ¼ 12C3,1 fsn ,1,1g,

Efw20 g ¼ k20,1 Efs21 gEfs62 g þ k20,2 Efs21 gEfs22 g þ k20,3 Efs61 gEfs22 g þ 2k0,1 k

k k

2 4 0,2 Efs1 gEfs2 g

þ 2k0,2 k 2

2

2

For Efz22 g, we have

6

2

2

4

6

4

3

þ k22,4 Efs21 sn1 gEfs22 g þ 2k2,1 k2,2 Efs21 gEfs2 sn2 g

2

þ9k0,3 9 Ef9s1 9 gEf9s2 9 g

4

4

4

2

4

2

2

2

þ 2k2,1 k2,3 Efs21 gEf9s2 9 g þ2k2,1 k2,4 Ef9s1 9 gEf9s2 9 g

3

þ2Refk0,1 kn0,3 Efs1 sn1 gEfsn2 s32 g

þ 2k2,2 k2,3 Efs21 gEf9s2 9 g þ2k2,2 k2,4 Ef9s1 9 gEf9s2 9 g

2

þ k0,1 kn0,2 Ef9s1 9 gEfsn2 s32 g 2

3

2

Efw22 g ¼ k22,1 Efs21 gEfs22 sn2 g þ k22,2 Efs21 gEfsn2 g þ k22,3 Efs21 gEfs22 g

2

Ef9w0 9 g ¼ 9k0,1 9 Ef9s1 9 gEf9s2 9 g þ 9k0,2 9 Ef9s1 9 gEf9s2 9 g 2

2

þ C3,1 fs,1,1gEfsn1 g32C3,1 fs,2,2gEfsn2 g,

For Efz20 g, we have 4 4 0,3 Efs1 gEfs2 g þ 2 0,1 4 2 0,3 Efs1 gEfs2 g,

853

4

þ k0,2 kn0,3 Efs1 sn1 gEf9s2 9 gg,

þ 2k2,3 k2,4 Ef9s1 9 gEfs22 g,

with

2

2

2

6

2

Ef9w2 9 g ¼ 9k2,1 9 Ef9s1 9 gEf9s2 9 g þ ð9k2,2 9 þ 9k2,3 9 ÞEf9s1 9 g,

n 3 2C4,0 fs ,2,2g,

k0,1 ¼ k0,2 ¼ 92 C4,0 fsn ,1,1gEfs21 g92C4,0 fsn ,2,2gEfs22 g, k0,3 ¼ 32C4,0 fsn ,1,1g:

2

2

6

2

Ef9s2 9 g þ 9k2,4 9 f9s1 9 gEf9s2 9 g 2

4

þ2Ref k2,1 kn2,2 Ef9s1 9 gEf9s2 9 g 2

3

3

þ k2,1 kn2,3 Ef9s1 9 gEfs2 sn2 g þ k2,1 kn2,4 Efs31 sn1 gEfs2 sn2 g

For Efz21 g, we have 2

2

2

þ k2,3 kn2,4 Efs31 sn1 gEf9s2 9 gg, with

4

2

þ 4k21,4 Efs41 sn1 gEfs22 g þ 4k1,1 k1,2 Efs21 gEf9s2 9 g

k2,1 ¼ C2,2 fs,2,2g þ 12C2,2 fsn ,2,2g,

2

þ 4k1,1 k1,3 Efs21 gEfs32 sn2 g þ 4k21,1 Ef9s1 9 gEfs52 sn2 g 4

þ 4k1,1 k1,4 Efs41 gEf9s2 9 g þ10k1,1 k1,4 Efs31 sn1 gEfs32 sn2 g 2

2

2

2

þ k21,1 Efsn1 gEfs62 g þ k21,4 Efs61 gEfsn2 g

2

þ k2,2 kn2,3 Ef9s1 9 gEfsn2 g þ k2,2 kn2,4 Efs31 sn1 gEfsn2 g

2

Efw21 g ¼ 4k21,1 Efs21 gEfs42 sn2 g þ k21,2 Efs21 gEfsn2 g þ k21,3 Efs21 gEfs22 g

2

þ 2k1,2 k1,3 Efs21 gEf9s2 9 g þ2k1,2 k1,1 Ef9s1 9 gEfs32 sn2 g

k2,2 ¼ Ef9s1 92 gðC2,2 fs,1,1g þ 2C2,2 fsn ,1,1gÞ 2

Ef9s2 9 gð2C2,2 fs,2,2g þ C2,2 fsn ,2,2gÞ, 2

k2,3 ¼ Efsn1 gðC2,2 fs,1,1g þ 12C2,2 fsn ,1,1gÞ

2

2

þ 2k1,2 k1,4 Efs41 gEfsn2 g þ 4k1,2 k1,4 Efs31 sn1 gEf9s2 9 g 2

2

Efsn2 gðC2,2 fs,2,2g þ 12C2,2 fsn ,2,2gÞ,

2

þ 2k1,3 k1,1 Ef9s1 9 gEfs42 g þ2k1,3 k1,4 Efs41 gEf9s2 9 g

k2,4 ¼ C2,2 fs,1,1g12C2,2 fsn ,1,1g:

4

þ 4k1,3 k1,4 Efs31 sn1 gEfs22 g þ 4k1,1 k1,4 Ef9s1 9 gEfs42 g 2

þ 4k21,4 Efs51 sn1 gEf9s2 9 g, 2

2

2

6

2

6

2

Ef9w1 9 g ¼ 59k1,1 9 Ef9s1 9 gEf9s2 9 g þ 59k1,4 9 f9s1 9 gEf9s2 9 g 2

2

2

EfzI zJ g ¼ EfRfwI gRfwJ gg ¼ 12RefEfwI wJ g þ EfwI wnJ gg:

2

þð9k1,2 9 þ9k1,3 9 ÞEf9s1 9 gEf9s2 9 g

For Efz0 z1 g, we have

2

þ2Ref2k1,1 kn1,2 Ef9s1 9 gEfsn2 s32 g 2

4 2

þ2k0,2 k1,1 Efs21 gEfs32 sn2 g þ k0,2 k1,2 Efs21 gEf9s2 9 g

4

2

2

þ k0,2 k1,3 Efs21 gEfs22 g þ k0,2 k1,1 Ef9s1 9 gEfs42 g

4

þ k1,2 kn1,3 Ef9s1 9 gEfsn2 g þ k1,2 kn1,1 Efs21 gEfsn2 g

2

þ k0,2 k1,4 Efs41 gEf9s2 9 g þ 2k0,2 k1,4 Efs31 gEfs22 g

2

þ k1,2 kn1,4 Efs1 sn1 gEf9s2 9 g 4

2

n2

þ2k1,2 k1,4 Ef9s1 9 gEfs2 g þ k1,3 k

þ2k0,3 k1,1 Efs41 gEfs32 sn2 g þ k0,3 k1,2 Efs41 gEf9s2 9 g

2 n3 1,1 Efs1 gEfs2 s2 g n

4

3

þ k0,3 k1,3 Efs41 gEfs22 g þ k0,3 k1,1 Efs31 sn1 gEfs42 g 2

þ k1,3 kn1,4 Efs1 sn1 gEfs22 g þ 2k1,3 kn1,4 Ef9s1 9 gEf9s2 9 g 4

2

2

þ4k1,1 kn1,4 Ef9s1 9 gEf9s2 9 g

n

3

þ k0,1 k1,4 Efs41 gEf9s2 9 g þ 2k0,1 k1,4 Efs31 gEf9s2 9 g

3

þ4k1,1 kn1,4 Efs1 sn1 gEfs32 sn2 g

3

2

þ k0,1 k1,3 Efs21 gEf9s2 9 g þ k0,1 k1,1 Ef9s1 9 gEfs2 sn2 g

2

þ29k1,1 9 Efs21 gEfs2 n 4s22 g

2

4

3

Efw0 w1 g ¼ 2k0,1 k1,1 Efs21 gEfs2 sn2 g þ k0,1 k1,2 Efs21 gEf9s2 9 g

4

þ2k1,1 kn1,3 Ef9s1 9 gEf9s2 9 g

4

The expected values, EfzI zJ g, ðI,JÞ 2 f0,1,2g2 with J 4 I, can be obtained from the definition of zI ½k that is using (44)–(47) and the following result:

2

2

2

þ k1,1 kn1,4 Efsn1 gEfs42 g þ 29k1,4 9 Efs41 sn1 gEfsn2 gg,

2

þ k0,3 k1,4 Efs61 gEf9s2 9 g þ 2k0,3 k1,4 Efs51 sn1 gEfs22 g, 2

2

2

Efw0 wn1 g ¼ 2k0,1 kn1,1 Ef9s1 9 gEfs42 sn2 g þ k0,1 kn1,2 Ef9s1 9 gEfs42 g

854

C. De Luigi, E. Moreau / Signal Processing 93 (2013) 842–855 2

6

2

4

3

4

2

2

2

þ k1,3 kn2,1 Ef9s1 9 gEfs32 sn2 g þ k1,3 kn2,2 Ef9s1 9 gEfs22 g

þ 2k0,2 kn1,1 Ef9s1 9 gEf9s2 9 g þ k0,2 kn1,2 Ef9s1 9 gEfs22 g 2

2

2

þ k0,1 kn1,4 Efs1 sn1 gEfs42 g þ 2k0,1 kn1,4 Ef9s1 9 gEfs32 sn2 g 2

2

þ k1,2 kn2,3 Ef9s1 9 gEfsn2 g þ k1,2 kn2,4 Efs31 sn1 gEfsn2 g

þ k0,1 kn1,3 Ef9s1 9 gEfs32 sn2 g þ k0,1 kn1,1 Efs21 gEf9s2 9 g

2

2

2

þ k1,3 kn2,3 Ef9s1 9 gEf9s2 9 g þ k1,3 kn2,4 Efs31 sn1 gEf9s2 9 g 2

2

þ k1,1 kn2,1 Efsn1 gEfs52 g þ k1,1 kn2,2 Efsn1 gEfs42 g

3

þ k0,2 kn1,3 Ef9s1 9 gEf9s2 9 g þ k0,2 kn1,1 Efs21 gEfs2 sn2 g

4

2

4

3

2

þ k0,2 kn1,4 Efs1 sn1 gEfs22 g þ 2k0,2 kn1,4 Ef9s1 9 gEf9s2 9 g þ 2k0,3 k

4 3 n 1,1 Efs1 s1 gEf9s2 9 g þ n

3 n 2 1,2 Efs1 s1 gEfs2 g

k0,3 k

2

6

þ k0,3 k1,4 Ef9s1 9

4

2

k k

2

þ k1,4 kn2,3 Efs31 sn1 gEfsn2 g þ k1,4 kn2,4 Efs51 sn1 gEfsn2 g

3

gEfs22 g þ 2 0,3

2

þ k1,4 kn2,1 Efs31 sn1 gEf9s2 9 g þ k1,4 kn2,2 Efs31 sn1 gEf9s2 9 g

n

4

þ k0,3 kn1,3 Efs31 sn1 gEf9s2 9 g þ k0,3 kn1,1 Efs41 gEfs2 sn2 g n

þ k1,1 kn2,3 Efsn1 gEfs32 sn2 g þ k1,1 kn2,4 Ef9s1 9 gEfs32 sn2 g

4

þ 2k1,4 kn2,1 Ef9s1 9 gEfs32 sn2 gþ 2k1,4 kn2,2 Ef9s1 9 gEfs22 g

2 4 n2 1,4 Efs1 s1 gEf9s2 9 g: n

4

2

2

þ 2k1,4 kn2,3 Ef9s1 9 gEf9s2 9 g þ 2k1,4 kn2,4 Efs41 gEf9s2 9 g:

For Efz0 z2 g, we have 2

Efw0 w2 g ¼ k0,1 k2,1 Efs21 gEfs42 sn2 g þ k0,1 k2,2 Efs21 gEfs32 sn2 g 4

Appendix F. Expressions of the coefficients in system (34)

þ k0,1 k2,3 Efs21 gEfs42 g þ k0,1 k2,4 Ef9s1 9 gEfs42 g 4

2

þ k0,2 k2,1 Efs21 gEf9s2 9 g þ k0,2 k2,2 Efs21 gEf9s2 9 g 4

þ k0,2 k2,3 Efs21 gEfs22 g þ k0,2 k2,4 Ef9s1 9 gEfs22 g 4 2,1 Efs1 gEf9s2 9

þ k0,3 k

4

2 4 2,2 Efs1 gEf9s2 9 g

g þ k0,3 k

2

þ k0,3 k2,3 Efs41 gEfs22 g þ k0,3 k2,4 Efs41 sn1 gEfs22 g, 2

2

Efw0 wn2 g ¼ k0,1 kn2,1 Ef9s1 9 gEfs52 sn2 g þ k0,1 kn2,2 Ef9s1 9 gEfs42 g 2

þ k0,1 kn2,3 Ef9s1 9 gEfs32 sn2 g þ k0,1 kn2,4 Efs31 sn1 gEfs32 sn2 g 2

2

þ k0,2 kn2,1 Ef9s1 9 gEfs32 sn2 g þ k0,2 kn2,2 Ef9s1 9 gEfs22 g 2

2

2

þ k0,2 kn2,3 Ef9s1 9 gEf9s2 9 g þ k0,2 kn2,4 Efs31 sn1 gEf9s2 9 g 3 n 3 n 2,1 Efs1 s1 gEfs2 s2 g þ 2 n 3 n 2,3 Efs1 s1 gEf9s2 9 g

þ k0,3 k

n

3 n 2 2,2 Efs1 s1 gEfs2 g

k0,3 k

n

þ k0,3 k

2

þ k0,3 kn2,4 Efs51 sn1 gEf9s2 9 g: For Efz1 z2 g, we have 6

4

Efw1 w2 g ¼ 2k1,1 k2,1 Efs21 gEf9s2 9 g þ 2k1,1 k2,2 Efs21 gEf9s2 9 g 4

þ 2k1,1 k2,3 Efs21 gEfs32 sn2 g þ 2k1,1 k2,4 Ef9s1 9 gEfs32 sn2 g 3

2

þ k1,2 k2,1 Efs21 gEfs2 sn2 g þ k1,2 k2,2 Efs21 gEfsn2 g 2

4

2

þ k1,2 k2,3 Efs21 gEf9s2 9 g þ k1,2 k2,4 Ef9s1 9 gEf9s2 9 g 4

2

þ k1,3 k2,1 Efs21 gEf9s2 9 g þ k1,3 k2,2 Efs21 gEf9s2 9 g 4

þ k1,3 k2,3 Efs21 gEfs22 g þ k1,3 k2,4 Ef9s1 9 gEfs22 g 2

8 > a12 ¼ J 4,0 ðEfz21 g þ Efz22 g2Efz2 z1 gÞ > > > > > þ ðJ 4,1 J 4,2 ÞðEfz0 z2 gEfz1 z0 gÞ, > > > > 2 > ¼ J a > 4,0 ð2Efz2 z1 g2Efz1 gEfz0 z2 g þ Efz1 z0 gÞ > < 11 þ J 4,1 ð2Efz1 z0 gEfz20 gEfz0 z2 gÞ þ J4,2 ðEfz20 gEfz1 z0 gÞ, > > > > a10 ¼ J 4,0 ðEfz21 g Efz1 z0 gÞ þ J 4,1 ðEfz20 g Efz1 z0 gÞ, > > > > > a01 ¼ J 4,0 ðEfz0 z2 gEfz1 z0 gÞ þ ðJ 4,1 J 4,2 ÞEfz20 g, > > > > : a00 ¼ J 4,0 Efz1 z0 gJ4,1 Efz2 g 0 and 8 > b ¼ J4,0 ð2Efz2 z1 gEfz21 gEfz22 gÞ > > 21 > > > þ J 4,1 ðEfz22 gEfz2 z1 gEfz0 z2 g þ Efz1 z0 gÞ > > > > > > þ J 4,2 ðEfz21 gEfz2 z1 g þ Efz0 z2 gEfz1 z0 gÞ, > > > > > b20 ¼ J 4,0 ðEfz21 gEfz2 z1 g þ Efz0 z2 gEfz1 z0 gÞ > > > > > > þ J 4,1 ðEfz20 g þEfz2 z1 gEfz0 z2 gEfz1 z0 gÞ > > > > < þ J 4,2 ð2Efz1 z0 gEfz21 gEfz20 gÞ, > b11 ¼ J4,0 ðEfz21 g þ Efz22 g2Efz2 z1 gÞ > > > > > þ ðJ4,1 J4,2 ÞðEfz0 z2 gEfz1 z0 gÞ, > > > > > b10 ¼ J ðEfz2 g þ Efz2 z1 g2Efz0 z2 g þ 2Efz1 z0 gÞ > 4,0 > 1 > > > > þ J 4,1 ðEfz0 z2 g þ Efz1 z0 g2Efz20 gÞ > > > > > > þ 2J 4,2 ðEfz20 gEfz1 z0 gÞ, > > > > : b00 ¼ J 4,0 ðEfz0 z2 gEfz1 z0 gÞ þ ðJ4,1 J4,2 ÞEfz20 g:

2

2

þ k1,1 k2,1 Ef9s1 9 gEfs42 sn2 g þ k1,1 k2,2 Ef9s1 9 gEfs32 sn2 g 2

All these expected values can be found in Appendix E.

3

þ k1,1 k2,3 Ef9s1 9 gEfs42 g þ k1,1 k2,4 Efs1 sn1 gEfs42 g 3

Appendix G. Expressions of the coefficients in (39) and (41)

2

þ k1,4 k2,1 Efs41 gEfs2 sn2 g þ k1,4 k2,2 Efs41 gEfsn2 g 2

3

þ k1,4 k2,3 Efs41 gEf9s2 9 g þ k1,4 k2,4 Efs41 gEfs2 sn2 g 4

þ 2k1,4 k2,1 Efs31 sn1 gEf9s2 9 g

G1 ¼ Efz21 g, G2 ¼ Efz22 g, G2,1 ¼ Efz2 z1 g,

2 3 n 2,2 Efs1 s1 gEf9s2 9 g

þ 2k1,4 k

G0 ¼ Efz20 g, G2,0 ¼ Efz2 z0 g:

6

þ 2k1,4 k2,3 Efs31 sn1 gEfs22 g þ 2k1,4 k2,4 Ef9s1 9 gEfs22 g, 2

ð49Þ

All these expected values can be found in Appendix E.

2

n2

Efw1 wn2 g ¼ 2k1,1 kn2,1 Ef9s1 9 gEfs42 s2 g þ 2k1,1 kn2,2 Ef9s1 9 gEfs32 sn2 g 2

4

References 4

þ 2k1,1 kn2,3 Ef9s1 9 gEf9s2 9 g þ 2k1,1 kn2,4 Efs31 sn1 gEf9s2 9 g 2

4

2

2

þ k1,2 kn2,1 Ef9s1 9 gEf9s2 9 g þ k1,2 kn2,2 Ef9s1 9 gEf9s2 9 g

[1] P. Comon, C. Jutten editors. Handbook of Blind Source Separation, Independent Component Analysis and Applications, Academic Press, 2010.

C. De Luigi, E. Moreau / Signal Processing 93 (2013) 842–855

[2] J. He´rault, C. Jutten, B. Ans, De´tection de grandeurs primitives dans un message composite par une architecture de calcul neuromime´tique en apprentissage non supervise´, in: Proceedings of the Xeme Colloque GRETSI, Nice, France, 20–24 Mai 1985, pp. 1017–1022. [3] C. Jutten, J. He´rault, Blind separation of source, Part i: an adaptative algorithm based on neuromimetic architecture, Signal Processing 24 (1991) 1–10. [4] E. Moreau, O. Macchi, High order contrasts for self-adaptative source separation, International Journal of Adaptive Control and Signal Processing 10 (January (1)) (1996) 19–46. [5] J.-F. Cardoso, A. Souloumiac, Blind beamforming for non Gaussian signals, IEE Proceedings—F 40 (1993) 362–370. [6] P. Comon, Independent component analysis, a new concept? Signal Processing 36 (1994) 287–314. [7] E. Moreau, A generalization of joint-diagonalization criteria for source separation, IEEE Transactions on Signal Processing 49 (March (3)) (2001) 530–541. [8] E. Moreau, N. Thirion-Moreau, Nonsymmetrical contrasts for sources separation, IEEE Transactions on Signal Processing 47 (August (8)) (1999) 2241–2252. [9] A. Bunse-Gerstner, R. Byers, V. Mehrmann, Numerical methods for simultaneous diagonalization, SIAM Journal of Matrix Analysis and Applications 14 (4) (1993) 927–949. [10] J.-F. Cardoso, A. Souloumiac, Jacobi angles for simultaneous diagonalization, SIAM Journal of Matrix Analysis and Applications 17 (January (1)) (1996) 161–164. [11] A. Belouchrani, K. Abed-Meraim, J.-F. Cardoso, E. Moulines, A blind source separation technique using second order statistics, IEEE Transactions on Signal Processing 45 (February (2)) (1997) 434–444. [12] L. De Lathauwer, B. De Moor, J. Vanderwalle, Independent component analysis and (simultaneous) third-order tensor diagonalization, IEEE Transactions on Signal Processing 49 (2001) 2262–2271. [13] V. Zarzoso, J.J. Murillo-Fuentes, R. Boloix-Tortosa, A.K. Nandi, Optimal pairwise fourth-order independent component analysis, IEEE Transactions on Signal Processing 54 (August (8)) (2006) 3049–3063. [14] B. Picinbono, On circularity, IEEE Transactions on Signal Processing 42 (December) (1994) 3473–3482. [15] B. Picinbono, P. Bondon, Second-order statistics of complex signals, IEEE Transactions on Signal Processing 45 (2) (1997) 411–419. [16] L. De Lathauwer, B. De Moor, J. Vanderwalle, ICA techniques for more sources than sensors, in: Proceeding of the IEEE Signal

[17]

[18]

[19]

[20]

[21]

[22]

[23]

[24]

[25]

[26]

[27]

[28]

855

Processing Workshop on Higher-Order Statistics (HOS ’99), Caesarea, Israel, June 1999, pp. 121–124. L. De Lathauwer, B. De Moor, On the blind separation of noncircular sources, in: Proceeding of EUSIPCO-02, vol. II, September 2002, Toulouse, France, pp. 99–102. P. Chevalier, Optimal time invariant and widely linear spatial filtering for radiocommunications, in: Proceedings of the EUSIPCO’96, Trieste, Italy, September 1996, pp. 559–562. P. Chevalier, Optimal array processing for non stationary signals, in: Proceedings of the ICASSP’96, Atlanta, USA, May 1996, pp. 2868– 2871. H. Abeida, J.P. Delmas, MUSIC-like estimation of direction of arrival for non-circular sources, IEEE Transactions on Signal Processing 54 (July (7)) (2006) 2678–2690. X.-L. Li, T. Adali, Blind separation of noncircular correlated sources using Gaussian entropy rate, IEEE Transactions on Signal Processing 59 (June (6)) (2011) 2969–2975. X.-L. Li, T. Adali, M. Anderson, Noncircular principal component analysis and its application to model selection, IEEE Transactions on Signal Processing 59 (October (10)) (2011) 4516–4528. M. Novey, T. Adali, On extending the complex FastICA algorithm to noncircular sources, IEEE Transactions on Signal Processing 56 (5) (2008) 2148–2154. C. De Luigi, E. Moreau, Optimal joint diagonalization of complex symmetric third-order tensors. Application to separation of non circular signals, in: Proceedings of the ICA’07, London, UK, September 2007, pp. 25–32. J.-F. Cardoso, On the performance of orthogonal source separation algorithms, in: Proceedings of the EUSIPCO-94 VII European Signal Processing Conference, Edinburgh, UK, September 13–17, 1994, pp. 776–779. J.-F. Cardoso, B.H. Laheld, Equivariant adaptive source separation, IEEE Transactions on Signal Processing 44 (December (12)) (1996) 3017–3030. V. Zarzoso, F. Herrmann, A.K. Nandi, Weighted closed-form estimators for blind source separation, in: Proceedings of the SSP-2001, 11th IEEE Workshop on Statistical Signal Processing, Singapore, August 6–8, 2001, pp. 456–459. E. Moreau, J.-C. Pesquet, N. Thirion-Moreau, Convolutive blind signal separation based on asymmetrical contrast functions, IEEE Transactions on Signal Processing 55 (January (1)) (2007) 356–371.