Bernard Xerri -Bruno Borloz

An Iterative Method Using Conditional Second-Order Statistics Applied to the Blind Source Separation Problem

Keywords: blind source separation, conditional statistics, iterative method, performance index, contrast functions

This paper is concerned with the problem of blind separation of instantaneous mixture of sources (BSS) which has been addressed in many ways. When power spectral densities of the sources are different, methods using second-order statistics are sufficient to solve this problem. Otherwise these methods fail and others (higher-order statistics,...) must be used.

In this paper, we propose an iterative method to process the case of sources with the same power spectral density. This method is based on an evaluation of conditional first and second-order statistics only. Restrictions on characteristics of sources are given to reach a solution, and proofs of convergence of the algorithm are provided for particular cases of probability density functions. Robustness of this algorithm with respect to the number of sources is shown through computer simulations.

A particular case of sources which have a probability density function with unbounded domain of definition is described; here the algorithm does not lead directly to a separation state, but to an a priori known mixture state. Finally, prospects of links with contrast functions are mentioned, with a possible generalization of them based on results obtained with particular sources.

I. Introduction

Blind source separation problem is currently receiving an increased interest [START_REF] Comon | Independent component analysis, a new concept ?[END_REF][6] [START_REF] Jutten | Blind separation of sources, Part I : An adaptative algorithm based on neuromimatic architecture[END_REF][9] [START_REF] Pham | Blind separation of instantaneous mixture of sources via an independent component analysis[END_REF] in numerous engineering applications. This problem consists in restoring N unknown, statistically independent random sources from N available observations which are linear combinations of these sources.

Methods based on second-order statistics (AMUSE [START_REF] Tong | AMUSE: A new blind identification algorithm[END_REF], SOBI [START_REF] Belouchrani | A blind source separation technique using second order statistics[END_REF] or IMISO [START_REF] Cavassilas | Séparation autodidacte de sources temporellement corrélées[END_REF]) have been extensively used to solve BSS problem. These methods perform well when power spectral densities (PSD) of sources are different. Otherwise High-Order Statistics (HOS) methods are required [START_REF] Lacoume | Statistiques d'ordres supérieurs pour le traitement du signal[END_REF] [START_REF] Tong | Indeterminacy and identificability of blind identification[END_REF][8] [START_REF] Pham | Blind separation of instantaneous mixture of sources via order statistics[END_REF] [START_REF] Ynouye | Cumulant-based blind identification of linear multi-input-multi-output systems driven by colored inputs[END_REF], subjected to delicate estimation: for example iterative methods using fourth-order cumulants have been early studied [START_REF] Cardoso | Iterative techniques for blind source separation using only fourth-order cumulants[END_REF].

In this paper, we present a new iterative method based on the evaluation of only first and secondorder statistics estimated from extracted series of observations (we talk of conditional statistics), avoiding calculation of high-order quantities. The main advantage of this method is that, like usual HOS techniques, it allows to process certain cases where PSD's of sources are the same; hence, it can doubtless be interpreted as an HOS method. However, in this paper, we deliberately do not compare HOS methods to ours from the theoretical point of view: for example, we do not approach the theoretical problems of robustness with respect to the number of sources. But numerical comparisons with a classical HOS method (JADE) are provided in this paper to draw conditions for which our method brings an improvement. Indeed this method is sensitive to probability density functions (PDF) of sources, and conditions on these PDF's are necessary to reach separation: these conditions are not very restrictive and usual domains of application as digital communications, etc... respect them. We will see that this method restores sources of same or different PSD's as well. We prove that our method converges. Classically HOS methods use contrast functions and impose a convergence to a separation state [START_REF] Moreau | High-order contrasts for self-adaptative source separation[END_REF][27] [START_REF] Pham | Contrastes pour la séparation aveugle de sources[END_REF]. We will prove that our method converges to a state which can be known according to PDF's of sources (then of engineering [START_REF] Pham | Blind separation of instantaneous mixture of sources via an independent component analysis[END_REF] DRAFT applications), and which is not necessarily a separation state: this method cannot be linked to a contrast function but nevertheless yields the separation. Some particular types of PDF's are more detailed.

To solve BSS problem for an instantaneous mixture, there exist different methods: HOS methods (algorithm JADE for example), geometrical methods [START_REF] Barrere | Chabriel The Archimedes principle applied to separation of uniformly distributed sources, SEE PSIP[END_REF][25] [START_REF] Puntonet | A competitive neural network for blind separation of sources based on geometry properties, biological and artificial computation: from neuroscience to technology[END_REF] or neural network methods [3][15]. Severe restrictive conditions must be verified for geometrical methods to work: PDF's of sources must have a bounded domain of definition and results are very sensible to the density on the edges.

Some authors use the property of non-stationarity of signals (for example: speech signals) to solve BSS problem: in the case of instantaneous mixtures [START_REF] Matsuoka | A Neural Net for Blind Separation of Non-stationary Signals[END_REF] or of convolutive mixtures [START_REF] Parra | Convolutive Blind Separation of Non-Stationary Sources[END_REF]: this assumption of non-stationarity is crucial for their algorithms. For instance, this assumption allows [START_REF] Matsuoka | A Neural Net for Blind Separation of Non-stationary Signals[END_REF] to use only second-order moments and an iterative procedure.

Concerning the present paper, in the mind of authors, the main idea of the method proposed is that it uses conditional statistics (restricted up to second-order ones) through an iterative algorithm. We are aware that HOS are efficient to solve the BSS problem for the instantaneous mixture model provided that cumulants are not null (non-gaussian signals). We are also aware that the convolutive mixture model is presently the main field of researches for the BSS problem, and that a method which works for instantaneous mixture could seem to bring little. However, as far as we know, it appears to us that the approach proposed is original enough; in a first time, this method needed to be tried and justified within the framework of the more simple model of instantaneous mixture; proofs are already not obvious to bring in this more simple model. Adaptation of this method for the more general convolutive mixtures model (reverberations for example) will be treated in a forthcoming paper. Hence, in a first step, this paper restrains deliberately its topic to the instantaneous mixture model for which proofs can be provided.

Formulation of the problem

For blind source separation problem, we consider N unknown and statistically independent sources X,

and N available observations Y = M X where M is an unknown deterministic nonsingular constant matrix called mixture matrix. The BSS problem consists in restoring sources X by applying to Y an unknown linear transformation S, so that SM = ΛΠ: this is a product of a nonsingular diagonal matrix and a permutation matrix. Thus S is a separating matrix. The algorithm proposed calculates at each step a matrix K i so that the resulting matrix K n K n-1 ...K 1 associated with the a priori knowledge of the state at the end of convergence tends to a separating matrix when n increases.

Organization of the paper

This paper is organized as follows. Section II formulates the mathematical model and introduces tools for study. Section III summarizes a few classical methods using only second order moments and underlines their limitations; the aim of the method proposed is to reduce them. After a detailed presentation of the reasons which lead us to the method proposed in section IV, our algorithm is explained in section V; its limitations and methods used to quantify its performances are described. Section VI is dedicated to analysis of convergence of this algorithm. Calculation of conditional moments in different interesting April 11, 2003 DRAFT cases is done. Proof of convergence of the algorithm is given according to PDF's of sources. Monte-Carlo simulation results are presented in section VII; robustness of this algorithm with respect to the number of sources is illustrated for some specific cases.

II. Basic assumptions and problem formulation

A. Basic assumptions

Let us consider N zero-mean random sources assumed to be unknown, stationary and statistically independent:

X(t) = X 1 (t) ... X N (t) T .
Through a set of receivers, N observations are available:

Y (t) = Y 1 (t) ... Y N (t) T ,
which are linearly linked to sources by:

Y (t) = M X(t), ( 1
)
where M is an unknown deterministic nonsingular constant matrix called mixture matrix. This relation can also be written:

Y i (t) = N j=1 α ij X j (t) ∀i = 1, ..., N.
We will note R Y (τ ) the spatial covariance matrix of Y (t) for a lag τ : R Y (τ ) E Y (t)Y T (t + τ ) . By assumption, R X (τ ) E X(t)X T (t + τ ) is a diagonal matrix.

B. Problem formulation

The principle of separation, i.e. restoration of sources, consists in applying a linear transformation S to Y to generate N signals z i statistically independent and proportional to sources x i :

Z SY = SM X.
There are two inherent indeterminacies concerning the order i and the power σ 2 i of each source. Hence, the matrix S is valid if it verifies: SM = ΛΠ with Λ nonsingular diagonal matrix and Π permutation matrix. Furthermore, the indeterminacy concerning σ 2 i allows, without loss of generality, to consider σ 2 i = 1; thus the spatial covariance matrix of sources at lag 0 is

R X (0) = E X(t)X T (t) = I N
where I N denotes the N -by-N identity matrix. [START_REF] Pham | Blind separation of instantaneous mixture of sources via an independent component analysis[END_REF] DRAFT

C. Measure index of closeness to ΛΠ

For the needs of our algorithm, we have to define an index which measures the closeness of a matrix A to ΛΠ. To quantify such a proximity, we will use an adapted version of the performance index proposed by E. Moreau [START_REF] Moreau | High-order contrasts for self-adaptative source separation[END_REF]; this one, which indicates how a matrix A is close to ΛΠ, is noted ind 0 (M ) and is defined as follows:

ind 0 (A) 1 2N (N -1)   N i=1   N j=1 |a ij | 2 max k |a ik | 2 -1   + • • • • • • N j=1   N i=1 |a ij | 2 max k |a kj | 2 -1     ,
and has the following properties:

0 ind 0 (A) 1 ind 0 (A) = ind 0 (ΠA) ind 0 (A) = 0 ⇔ A = ΛΠ if α is a nonzero scalar, ind 0 (αA) = ind 0 (A)
Notwithstanding its properties, this performance index does not suit for our needs. We want an index ind(.) which verifies also the additional property: ind (ΛA) = ind (A). A way to complete this requirement is to apply to A the following transformation:

A = [a ij ] u -→ u(A) =      a11 δ1 • • • a 1N δ1 . . . • • • . . . a N 1 δ N • • • a N N δ N     
with δ i = a 2 i1 + ... + a 2 iN , and then to apply the classical performance index to u(A) in lieu of A:

ind(A) ind 0 (u(A)).
Note first that if U is a unitary matrix, u(U ) = U and then ind(U ) = ind 0 (U ), and secondly that a normalization of rows of A by δ i = max j (a ij ) would be effective as well. From here, we will call 'performance index of A' the value ind(A).

III. Classical second order methods

Methods using second-order statistics have been successfully proposed by several authors for the case of sources with different PSD's. Some of them are presented hereafter.

A. AMUSE

The AMUSE (Algorithm for Multiple Unknown Signals Extraction) algorithm proposed by L. Tong and al. [START_REF] Tong | AMUSE: A new blind identification algorithm[END_REF] processes the more general following model:

Y (t) = M X(t) + n(t)

B. SOBI

The SOBI (Second Order Blind Identification) algorithm has been proposed by A. Belouchrani and al. [START_REF] Belouchrani | A blind source separation technique using second order statistics[END_REF].

A matrix W can be found to ensure spatial decorrelation of the process Y w (t) W Y (t). For a lag τ = 0, the covariance matrix of whitened observations

Y w (t) is R Y w (τ ) W R Y (τ )W T which can be written R Y w (τ ) = W M R X (τ )M T W T = U R X (τ )U T ,
where U is a unitary matrix to be evaluated.

Matrices R X (τ ) are diagonal and then U is the eigenmatrix for all matrices R Y (τ ), ∀τ . If a value τ p of τ for which R Y (τ p ) has distinct eigenvalues can be found, U can be completely determined. An interesting option is to choose a set of values of τ say {τ 1 , ..., τ k } and to perform a simultaneous diagonalization

of the set of matrices [R Yw (τ 1 ), R Yw (τ 2 ), ..., R Yw (τ k )]
. This technique which has been introduced in [START_REF] Cardoso | Jacobi angles for simultaneous diagonalization[END_REF],

consists in a generalization of Jacobi's method to diagonalize an hermitian matrix. Knowing U , sources are restored by X(t) = U T W Y (t).

C. IMISO

The IMISO (Identification de Mélanges Instantanés au Second Ordre -Second Order Identification of Instantaneous Mixtures) algorithm proposed by J.F. Cavassilas and al. [START_REF] Cavassilas | Séparation autodidacte de sources temporellement corrélées[END_REF][13] is a natural extension of previous one and can be summarized as follows: the covariance matrix of observations is

R Y (τ ) = M R X (τ )M T .
For commonly encountered signals (band-limited), it is reasonable to postulate that R Y (τ ) is indefinitely differentiable with respect to τ . Then

R (2k) Y (0)M -T = R Y (0)M -T R -1 X (0)R (2k) X (0).
This is a generalized eigen-problem:

R Y (0)V = R (2k) Y (0)V D, where D = R -1 X (0)R (2k)
X (0) is obviously a diagonal matrix, and V = M -T (except for a diagonal and a permutation matrix). Then, applying V T = M -1 to Y leads to X (except for a diagonal and a permutation matrix).

At least two conditions are required to assert such a result: (0) must be proportional to a row of R Y (0); such a thing would correspond to the case where second-order series expansion of normalized autocorrelation functions of sources are identical or very close in vicinity of 0; this assumption is generally not realistic. This case is involved if

R -1 Y (0)R (2k) Y (0)
has a multiple eigenvalue, i.e. there is an indeterminacy within the eigen-subspace concerning the choice of eigenvectors to take into account. If PSD's (therefore autocorrelation functions) of sources are identical, there is one eigenvalue of multiplicity N : indeterminacy is complete: it is not possible to extract privileged directions.

Experimentally a value k=1 (second order derivative) is then often sufficient, but a value k=2 (or more) may be necessary to solve the problem.

A generalized method, including AMUSE and IMISO as particular cases, has been recently proposed by J. Barrère and G. Chabriel [START_REF] Barrere | Chabriel A compact array for blind separation of sources[END_REF].

D. Limitations of these methods

There is unfortunately a severe limitation to these methods: they work in one step only with secondorder statistics. This is why it is not possible to use them to restore sources which have the same PSD:

sources must be differently colored. In such a case, people usually try to solve the problem by using higher-order statistics: the main disadvantage of such an approach lies in the estimation of higher-order moments which becomes quickly difficult when N increases.

For sources with identical PSD, the way of using first and second-order conditional statistics through an iterative formulation seems attractive. Nevertheless, the idea remains to find a linear transformation that makes independent observations.

IV. Use of conditional statistics

A. Introduction

The idea of previous algorithms is always to find two functions of Y = M X, f 1 and f 2 verifying the property of linearity

f i (Y ) = M f i (X) ∀i
and such as the two matrices

R i = E f i (Y )f T i (Y ) = M E f i (X)f T i (X) M T = M D i M T verify that R -1 1 R 2 = M -T D -1 1 D 2 M T , ( 2
)
is an eigenvalue problem (this is the case provided that D -1 1 D 2 is a diagonal matrix) with non-multiple eigenvalue. Then, applying the transpose of an eigenvectors matrix to Y leads to X (except for a diagonal and a permutation matrix). In the case of identical PSD, there is a unique eigenvalue of multiplicity N , and this approach (classical second-order moments) is hopeless.

B. Conditional statistics

Let us introduce conditional moments. Our problem is very specific since practically only linear combinations of sources are available (1), say the observations Y i . For example, the condition used can be:

Y i > 0.
Conditional first and second-order moments are E (Y /Y i > 0) and E Y Y T /Y i > 0 . The second one is naturally one out of the two matrices needed, say

R 1 . Denoting ν E(Y /Y i > 0), it is tempting to define R 2 = E (Y -ν)(Y -ν) T /Y i > 0 . Then R -1 1 R 2 = R -1 1 (R 1 -νν T ) = I N -R -1 1 νν T .
Noting conditional mean of sources (N -dimension vector)

µ = E (X/Y i > 0) , (3)
and conditional covariance matrix of sources

D = E XX T /Y i > 0 , (4)
it comes

R -1 1 R 2 = M -T I N -D -1 µµ T M T .
Generally I N -D -1 µµ T is not a diagonal matrix. However, this must be the case if separation is already reached: Y = X (M = I N ). In such a case, µ is a vector whose all components are null except the i th which is strictly positive. The term µµ T is already a diagonal matrix (with only one non-zero term).

Matrix D remains to be studied.

If we note D i the domain of definition of source X i , the hyper-volume described by all sources is then

Σ = D 1 × ... × D N ⊂ R N . Let us consider one observation Y i = N k=1 α ik X k . The equation "Y i = 0" is those of an hyperplane Π passing through O; Π divides Σ into two parts Σ i (Y i = N k=1 α ik X k > 0) and Σ i (Y i < 0). Using condition "Y i > 0" describes the hyper-volume Σ i .
Making appear the joint conditional PDF of X under the condition Y i > 0: p X/Y i >0 (x 1 , ..., x N ), it comes that the j th component of µ is given by:

µ j = R N x j p X/Yi>0 (x 1 , ..., x N )dx 1 ...dx N ,
and element of row j and column k of D is [START_REF] Pham | Blind separation of instantaneous mixture of sources via an independent component analysis[END_REF] DRAFT Like in the case of expected values under a discrete random variable condition, we can calculate these expressions as following: noting P (Y i > 0) = Σ i p X (x 1 , ..., x N )dx 1 ...dx N the probability that the i th observation is positive, they are both calculated by integrating first or second-order moments in Σ i :

D jk = R N x j x k p X/Y i >0 (x 1 , ..., x N )dx 1 ...dx N .
µ j = 1 P (Y i > 0) Σi x j p X (x 1 , ..., x N )dx 1 ...dx N , (5)
and

D jk = 1 P (Y i > 0) Σ i x j x k p X (x 1 , ..., x N )dx 1 ...dx N . ( 6
)
Of course as sources are independent, p X (x 1 , ..., x N ) = p 1 (x 1 )...p N (x N ).

D is not diagonal in the general case. However, as the condition used is linear passing through O, it can be proved that if PDF's of sources are symmetrical (

x i ∈ D i ⇒ -x i ∈ D i and p i (x i ) = p i (-x i )), then P (Y i > 0) = 1 2 and matrix D is diagonal of terms D jj = Dj x 2 p j (x)dx = σ 2 j .
We can note that this term does not depend on the hyperplane Π chosen but only on PDF's of sources.

proof Let us state that the PDF's p i (x) of sources are symmetrical. Σ i and Σ i are symmetrical: effectively,

x = (x 1 , x 2 , ...x N ) ∈ Σ i ⇒ -x = (-x 1 , -x 2 , ... -x N ) ∈ Σ i . (7)
Then defining f nm ij (x) = x n i x m j p 1 (x 1 )...p N (x N ), which verifies the property

f nm ij (-x) = (-1) n+m f nm ij (x) ∀x ∈ Σ, ( 8
) Σ f nm ij (x)dx = Σ i f nm ij (x)dx + Σ i f nm ij (x)dx + Π f nm ij (x)dx.
f nm ij (x) being a continuous function with 'good properties', the last term is null, and

Σ f nm ij (x)dx = Σ i f nm ij (x)dx + Σ i f nm ij (x)dx.
Furthermore, from ( 7) and ( 8),

Σ i f nm ij (x)dx = (-1) n+m Σi f nm ij (x)dx. If n + m is even Σi f nm ij (x)dx = 1 2 Σ f nm ij (x)dx = 1 2 E(X n i X m j ).
This term depends no more on the condition chosen but only on PDF's of X i and X j .

As X i and X j are independent, if i = j, this expression becomes 1 2 E(X n i )E(X m j ); in this case, if n or m is odd, it becomes null. If i = j, then we obtain the (n + m) th -moment of X j .

If n = m = 0 we obtain P (Y i > 0) = 1 2 . If n = m = 1 we obtain the elements of the matrix D which is diagonal of term D jj = σ 2 j . Obviously if PDF's are identical, this matrix becomes scalar. [START_REF] Pham | Blind separation of instantaneous mixture of sources via an independent component analysis[END_REF] DRAFT

C. Conclusion

In conclusion, if we want to respect the condition "D is a diagonal matrix", two sufficient (but not necessary) conditions are:

• PDF's of sources are symmetrical,

• the condition to take into account is the equation of an hyperplane passing through O.

First condition is a very acceptable one: in the field of digital communications, signals are usually uniformly distributed (with a limited or unlimited number of states), ...

The second one is in fact imposed by the problem: available quantities are observations, i.e. linear combinations of sources.

For reasons explained in this section and others that will be detailed later, those both conditions will be stated to be verified in next sections.

V. Presentation of the algorithm

This section is devoted to a general presentation of the algorithm proposed. In the light of previous section, this algorithm relies on an iterative application of an adaptation of the generalized algorithm described in section IV with conditional statistics: a linear condition (based on observations) is used to make appear those conditional statistics.

A. Principle of the algorithm

Let us expose the mechanism to restore sources. As said in previous section, we first extract from observations a conditional subsequence of them, called conditional observations and noted

Y c i = (Y /Y i > 0),
or:

Y c i = Y ci,1 ... Y ci,N T .
This is a highly nonlinear operation and

Y ci,i > 0. Obviously E (Y ci ) = E(Y /Y i > 0) = 0, and we can define Ỹc i Y c i -q c i where q c i = E (Y c i ).
From these extracted observations, we calculate the N -by-N covariance matrices of Y c i and Ỹc i noted:

R c i E Y c i Y T ci = E Y Y T /Y i > 0 and Rci E Ỹci Ỹ T ci = R ci -E (Y /Y i > 0) E (Y /Y i > 0) T .

Now let us define

Γ c i R -1 ci Rc i ,
and let us process the eigen-decomposition of Γ c i :

Γ c i P = P Q. April 11, 2003 DRAFT
If we write (1) as Y (1) = M (1) X, denoting the first iteration of our algorithm, and then apply P T to Y (1) (this is natural according to the comment linked to equation ( 2) in section IV), we obtain a new combination of observations called "new observations"

Y (2) = P T Y (1) = P T M (1) X = M (2) X,
which appears as a new mixture of sources; M (2) is the new mixture matrix:

M (2) = P T M (1) . ( 9
)
Here, the first step of our algorithm is finished. If properties of M (2) are satisfying, this step can be iterated again on new observations Y (2) .

P is determined except for a permutation and a power. In other words, P ΛΠ (or P ΠΛ) are also eigenvectors of Γ c1 : those both indeterminacies involved in Λ and Π correspond to the inherent indeterminacies aforementioned regarding power of sources and their order. P can be normalized readily. When applying P T to observations Y , we apply in fact the transformation ΛΠP T (or ΠΛP T ). According to section II-C,

ind(ΛΠP T M (1) ) = ind(P T M (1) ) = ind(M (2)
).

If P is a diagonal matrix, no improvement has been brought by the transformation: this observation leads to determine the criterion to stop convergence.

test for end of convergence As M is unknown, mixture matrix M (n) at step n is never reachable.

Even if we use it for proofs of convergence, we cannot use it practically during convergence to test the closeness of our solution to the good one.

Nevertheless, matrix P being available at each step, we stop iterative process for an ind(P ) smaller than a value ε defined beforehand, i.e. when the transformation matrix P is close to the product ΛΠ, in other words when convergence is achieved.

We will study in section VI the state at the end of convergence.

general methodology of algorithm

Y (0) = Y .
Let's consider modified observations available at step n:

Y (n) = Y (n) 1 ... Y (n) N T . 1) Create extracted observations Y ci = (Y /Y i > 0). 2) Compute zero-mean extracted observations Ỹci = Y ci -E (Y ci ).
3) Estimate R c i and Rc i .

4) Calculate Γ c i = R -1 c i Rc i . 5) Calculate the eigenvectors P of Γ c i .

6) Apply the transformation

P T to Y (n) : Y (n+1) = P T Y (n) = M (n+1) X.

7)

If ind(P ) < ε stop iterations. [START_REF] Pham | Blind separation of instantaneous mixture of sources via an independent component analysis[END_REF] DRAFT remark: Γ ci can also be written

Γ ci = I N -R -1 c i q ci q T c i . Eigenvectors of Γ ci are then R -1 c i q ci (associated with eigenvalue 1 -q T c i R -1 c i q c i )
and N -1 vectors orthogonal to q c i associated with the same eigenvalue 1, say (q ⊥ c i ,1 , ...q ⊥ c i ,N -1 ). Then P = R -1 ci q ci q ⊥ ci,1

... q ⊥ ci,N -1 . We note that there is an eigenvalue with multiplicity N -1 which corresponds to a subspace of dimension N -1 in which there is an indeterminacy for the choice of N -1 eigenvectors (this indeterminacy is cancelled if N = 2). In fact practically we observe that the algorithm applied stricto sensu extracts one out of the sources. If one source is found, it is theoretically possible to extract it and to reduce the problem to a (N -1)-dimension problem. Practically however, a significant deterioration appears quickly, and such a way is not conceivable. The way of a simultaneous extraction of all the sources is greatly preferable. The algorithm needs to be modified: as mixture matrix is nonsingular, it is reasonable to suppose that the same test (of series extraction) applied to each observation will lead to N satisfactory vectors, i.e. linearly independent. Consequently the algorithm is modified as follows:

Summarized basic algorithm at step n

Y (0) = Y . Let's consider modified observations available at step n: Y (n) = Y (n) 1 ... Y (n) N T . 1) For i=1 to N , 1-a) Create extracted observations Y c i = (Y /Y i > 0). 1-b) Compute q c i E (Y c i ). 1-c) Estimate R c i E Y c i Y T c i . 1-d) Calculate p i R -1 ci q c i . 2) Create P = R -1 c 1 q c 1 ... R -1 c N q c N = p 1 ... p N . 3) Apply the transformation P T to Y (n) : Y (n+1) = P T Y (n) = M (n+1) X. 4) If ind(P ) < ε stop iterations.

B. Rewriting according to sources

To study properties of the algorithm, let's see its implications with regard to sources. We will use the same notation for X ci extracted sequence of X corresponding to Y ci (such as Y ci = M X ci ); we will note

µ ci E(X ci ) = E(X/Y i > 0), Xci X ci -µ ci and D ci E X ci X T c i = E XX T /Y i > 0 .
Obviously, it is not possible to estimate directly these quantities linked to X from available data (i.e. the sole observations) while Γ c i can be calculated. We deduce that:

R c i = M D c i M T , and Rci = M D ci -µ ci µ T ci M T . April 11, 2003 DRAFT Then Γ c i = M -T I N -D -1 c i µ c i µ T c i M T ,
which implies

I N -D -1 c i µ c i µ T c i M T P = M T P Q.
This expression makes appear M T P which from (9) equals M (2) T . If we define V = M T P (V T = M (2) ), then V and Q are eigenvectors and eigenvalues of

I N -D -1 c i µ c i µ T c i : I N -D -1 c i µ c i µ T c i V = V Q. ( 10
)
P and Q depend on Y and are then calculable, while V is not.

From [START_REF] Cavassilas | Séparation autodidacte de sources temporellement corrélées[END_REF] we can assert that the N -1 linearly independent vectors orthogonal to

µ c i (µ ⊥ ci,k such as µ T ci µ ⊥ c i ,k = 0 ∀ k = 1, ..., N -1) are also eigenvectors of I N -D -1 ci µ ci µ T ci with the same eigenvalue 1. It is easy to verify that D -1 c i µ ci is an eigenvector of I N -D -1 c i µ ci µ T c i associated with the eigenvalue 1 -µ T c i D -1 c i µ c i .
Hence, we can write:

V T = D -1 ci µ c i µ ⊥ ci,1 ... µ ⊥ ci,N -1 T = M (2) .
If we apply the modified algorithm, instead of this expression, new mixture matrix will be

M (2) = D -1 c1 µ c1 ... D -1 c N µ c N T . The i th line of M (2) is M (2) i = D -1 ci µ c i T . ( 11
)
One line of M (2) depends on corresponding line in M (1) and PDF's of sources. Theoretical expression of M (2) is known. We will see that the difficulty lies in calculation of µ ci .

From IV-B, if PDF's of sources are symmetrical then P (Y i > 0) = 1 2 and covariance matrix D ci is diagonal of terms

D c i ,jj = σ 2 j . ( 12
)
In this case

M (2) i = µ c i ,1 Dc i ,11
...

µ c i ,N D c i ,N N .
Obviously if PDF's are identical of power σ 2 , this matrix becomes scalar: D c i ,jj = σ 2 , and

M (2) i = 1 σ 2 µ ci,1 ... µ ci,N .
remark: limitations of the study of convergence We readily see that M

(2) i

(and therefore M (2) ) closely depends on the PDF's of sources. This is the main limitation of the study of convergence of this algorithm: theoretical calculation of µ ci and D ci may hardly be done in a general way. Unfortunately proof of efficiency of the algorithm must be provided for each specific sort of PDF. This observation will prompt us to focus our attention on some noteworthy cases of PDF's, even if general considerations can be brought.

April 11, 2003 DRAFT

C. Remark about a whitening stage of observations

Before running the algorithm, it is always possible to whiten observations. Then at each step, the transformation matrix P must comply with these constraint: P must be a unitary matrix. Hence M (2) is also unitary. This procedure comes down to make orthogonal lines of M (2) calculated by the algorithm. Practically expressions R -1 ci q ci are calculated and normalized to 1, then P is calculated by

P = R -1 c 1 q c 1 ... R -1 c N q c N
and forced to be orthogonal by columns. This is only a guarantee not to converge twice or more to the same source for different initial conditions.

Practically the matrix P obtained is indeed nearly but not exactly orthogonal. However, because of this very small distance to a unitary matrix, simulations show that, without imposing a constraint of orthogonality to P at each step, the algorithm derives slowly to a situation where mixture matrix becomes non-invertible, and certain sources are found twice or more while some others are lost. So it becomes necessary to impose such a constraint to P (in practice, vectors (columns) P i of P are normalized by a Gram-Schmidt algorithm, and its normalization is never a problem because sources are retrieved except for their power). By constraining P to be orthogonal, the correction brought is not as much important than could appear. This constraint just imposes that two lines of new mixture matrix are orthogonal.

However it is not necessary to pre-whiten observations. When observations Y are whitened, that means that a linear transformation W = √ D -1 U T is applied to Y where U and D are respectively the eigenvectors and eigenvalues matrices of covariance matrix E Y Y T . If observations are whitened, transformation matrix P must comply with some constraints; let's note Φ = P T W -1 . A procedure is used to make Φ orthogonal. Transformation applied to Y will be ΦW in lieu of P T . This is a manner of avoiding pre-whitening stage of observations, but keeping the structure of invertible mixture matrix.

VI. Analysis of convergence according to PDF

In the light of previous sections, throughout this section we state that sources have symmetric PDF's.

A. Assumptions and notations

Let us consider one observation

Y i = N k=1 α ik X k , i.e. one line of the mixture matrix M M (1) i = α i1 ... α iN ; M (1) i can be considered as a point in R N .
After one step of our algorithm, the new observation Y

(2) i corresponds to the i th line of the new mixture matrix, which from (11) verifies:

M (2) i = m i1 m i2 ... m iN = D -1 ci µ c i ,
where µ ci and D ci are given by ( 3) and (4); their components are respectively calculated by ( 5) and [START_REF] Cardoso | Independent component analysis, a survey of some algebraic methods[END_REF].

m ij = µ c i ,j D c i ,jj = Σi x j p 1 (x 1 )...p N (x N )dx 1 ...dx N Σi x 2 j p 1 (x 1 )...p N (x N )dx 1 ...dx N j = 1, ...N. ( 13
)
April 11, 2003 DRAFT As PDF's of sources are symmetrical,

m ij = µ ci,j σ 2 j = 2 Σ i x j p 1 (x 1 )...p N (x N )dx 1 ...dx N σ 2 j . ( 14
)
Study consists in comparing the i th line of M (1) and the i th line of M (2) . This comparison is made trough the index performance naturally adapted to one line. A constraint exists nevertheless: M is invertible; at each step this constraint must be verified to ensure coherence of columns of new mixture matrix. If a pre-whitening stage is processed, this constraint is simply: 'P is unitary'.

B. Miscellaneous considerations

From the point of view of the algorithm there is no restriction to normalize M

(1)

i , e.g. by its largest value: if

|α i1 | > |α ij | ∀j, this line is equivalent to 1 α i2 αi1 ... α iN αi1 (the hyperplane of equation N k=1 α ik X k = 0 is the same than N k=1 α ik α i1 X k = 0).

B.1 Fixed-points of the algorithm

We define a fixed-point of the algorithm as a point such as, after one step, the new observation Y

(2) i corresponds to the following (normalized) line of new mixture matrix:

1 γ i2 ... γ iN , verifying ∀j = 2, ...N, γ ij = α ij
αi1 . We want to study the behavior of our algorithm around its fixed-points. For each fixed-point we want to determine if it is attractive or repulsive.

B.2 Absence of one source in observation

Let us consider the case where one source X i is not present in observation Y 1 (i.e. α 1i = 0). The initial mixture line is α 11 ... α 1,i-1 0 α 1,i+1 ... α 1N . Numerator of m 1i in equation ( 13) can than be written:

• • • N P k=1 k =i α k x k >0 p 1 (x 1 )...p i-1 (x i-1 )p i+1 (x i+1 )...p N (x N )dx 1 ...dx i-1 dx i+1 ...dx N × +∞ -∞ x i p i (x i )dx i .
As sources are zero-mean, this term is null: m 1i = 0. Conclusion As expected a source cannot appear in a line if it is not present from the start. The final line will be necessarily like m 11 ... m 1,i-1 0 m 1,i+1 ... m 1N .

B.3 Observations of the form

Y i = α q k=1 X k with q ≤ N
Let us find conditions for points like α ... α 0 ... 0 with q 'α' and p '0' (p + q = N ) to be fixed-points of the algorithm. [START_REF] Pham | Blind separation of instantaneous mixture of sources via an independent component analysis[END_REF] DRAFT

The p last sources can take every value in their domain of definition. Numerator of m ij in equation ( 13) can be written:

N ij = +∞ -∞ ... +∞ -∞    • • • x 1 +...+x q ≥0
x j p 1 (x 1 )...p q (x q )dx 1 ...

dx q    p q+1 (x q+1 )...p N (x N )dx q+1 ...dx N From VI-B.2 if j ∈ {q + 1, ...N }, doubtless m ij = 0. Otherwise N ij = • • • x 1 +...+x q ≥0
x j p 1 (x 1 )...p q (x q )dx 1 ...

dx q × N k=q+1 +∞ -∞ p k (x k )dx k = • • • x 1 +...+x q ≥0
x j p 1 (x 1 )...p q (x q )dx 1 ...dx q .

For obvious reasons of symmetries

N i1 = N i2 = ... = N iq provided that p 1 (x) = p 2 (x) = ... = p q (x) ∀x.
As in this case we have seen that D ci,jj = σ 2 ∀j, we conclude from ( 14) that m i1 = m i2 = ... = m iq .

Conclusion

B.4 Separation point

From previous paragraph, it comes that 0 ... 0 1 0 ... 0 (there is only one 1 at the j th position) is a fixed-point of the algorithm. This particular point is called a separation point of the algorithm in the sense that Y i = X j ; only one source appears in observation.

B.5 Conclusion

From now we will consider that all the sources have the same symmetric PDF p(x). In fact the method proposed does not necessitate this hypothesis and a mixture of sources with different symmetric PDF's can be successfully processed (see section VII).

Power of sources will be noted σ 2 and from equation ( 12), conditional covariance matrix of sources is scalar: D ci = σ 2 I N whatever the condition taken into account.

C. Analysis of convergence according to PDF

This section is devoted to an analysis of convergence of the algorithm around fixed-points in the particular case where PDF's of sources are symmetrical and identical. As proved above, results closely depend on PDF's of sources; then three cases of PDF's will be detailed:

uniformly distributed sources (continuous or discrete, i.e. taking a limited number of values): such sources are of particular interest for digital communications (16-QAM, 4-PAM,...).

gaussian sources: we will confirm that, according to known results, it is not possible to restore gaussian sources from a linear mixture of them.

April 11, 2003 DRAFT sources with a Laplace's PDF (p(x) = e -2|x| ).

Those three examples have been carefully chosen because they have different behaviors particularly interesting: convergence to a separation point for the first one (as in contrast function [START_REF] Moreau | A generalization of joint-diagonalization criteria for source separation[END_REF], [START_REF] Pham | Contrastes pour la séparation aveugle de sources[END_REF]), no convergence for the second, and a convergence to an a priori known mixture point for the last.

To illustrate these behaviors, we start our study with the particular case N = 2 where calculation can be done completely. Following, we develop results for any value of N .

C.1 Study in the case

N = 2 Starting from   Y 1 Y 2   =   α 11 α 12 α 21 α 22     X 1 X 2 
 , let's consider in a first time the first observation as the linear condition. First line of new mixture matrix is given by equations ( 11), ( 5) and ( 6).

Supposing that α 12 = 0 (else separation is already reached

), Σ 1 = {(x 1 , x 2 )/α 11 x 1 + α 12 x 2 > 0}.
For the sake of simplicity, in the whole cases, we suppose that observations have been pre-whitened and that 0 < α 11 < α 12 (other cases are deduced from this one with no difficulty); we will note 0 < γ

α 11 α 12 < 1.
The case of whitened observations can be boiled down to

Y w W Y = M w X with M w =   α 11 α 12 -α 12 α 11   ,
an orthogonal matrix. This particular case is quite interesting because ind(M w ) is readily calculable:

ind(M w ) = min(|α 11 | , |α 12 |) max(|α 11 | , |α 12 |) 2 ,
and can be immediately compared to the performance index of new mixture matrix (which will be also orthogonal). Then γ is nothing else than the square root of initial performance index ind M (1) .

Let us study the three cases of interest described above:

C.1.a Uniformly distributed sources. Let us state that X 1 and X 2 are random signals uniformly distributed within [-1; +1], their PDF is then p(x) = 1 2 1 [-1;+1] (x) and their power σ 2 = 1 3 . Calculations can be lead completely. From equation [START_REF] Cardoso | Jacobi angles for simultaneous diagonalization[END_REF]

calculations give    µ c1,1 = 1 3 γ µ c 1 ,2 = 1 2 -1 6 γ 2 , which implies    m 11 = γ m 12 = 3 2 -1 2 γ 2 Initial performance index is γ 2 . As 0 < γ < 1, 0 < m 11 < m 12 and ind M (2) = m11 m12 2 = 2γ 3-γ 2 2 .
Then 0 < m 11 m 12 < γ < 1 and ind M (2) < ind M (1) .

April 11, 2003 DRAFT If X 1 and X 2 are gaussian random signals of power σ 2 , and in the case of whitened observations, expressions of m 11 and m 12 can be easily found in the same way of previous section. From ( 5)

   µ c1,1 = 1 √ 2π σ 1 1+γ 2 µ c1,2 = 1 √ 2π σγ 1 1+γ 2
As from equation ( 12)

D ci = σ 2 I 2 , ratio m11 m12 = µ c 1 ,1 µc 1 ,2 = |γ|.
Consequently, and as we expected to, no improvement occurs after the first step of the algorithm and M (2) is no more discriminating than M (1) w.r.t. sources. This proves that, as it is well known, it is not possible to restore two gaussian sources. C.1.c Sources with Laplace's distribution. We consider here that X 1 and X 2 are random signals with PDF: p(x) = e -2|x| (its domain of definition is infinite unlike uniform PDF). Expressions of µ c1,j are:

   µ c1,1 = 1 4 2γ+1 (γ+1) 2 µ c1,2 = 1 4 γ(γ+2) (γ+1) 2 April 11, 2003 DRAFT New ratio is m 11 m 12 = µc 1 ,1 µ c 1 ,2 = 2γ+1 γ(γ+2) . For 0 < γ < 1, m 11 < m 12
, ind(M (2) ) > ind(M (1) ). The new performance index vs. the old one is shown on figure 2. There is only one attracting fixed-point (1,1) and a repelling fixed-point (0,0). This behavior is particularly interesting since the algorithm does not converge to a separation state, but to an a priori known mixture state. It will be experimentally confirmed with more than two sources with Laplace's PDF (cf. section VII).

C.1.d Conclusion.

This section shows clearly how closely the behavior of our algorithm depends on PDF's of sources. Three different behaviors have been observed:

• convergence to a separation state,

• convergence to an a priori known mixture state,

• no convergence.

Now, let's study convergence for any value of N .

C.2 Study for any number of sources

This section is devoted to a local study of convergence of the algorithm around points of interest, say fixed-points.

Starting from one line of M (1) : M From ( 14), ∀ j = 1, ..., N

m j = 2 σ 2 +∞ -∞ ... +∞ -∞       +∞ - N P i=2 αixi x j p(x 1 )dx 1       p(x 2 )...p(x N )dx 2 ...dx N . ( 15
)
As we want to study the behavior of the algorithm close to fixed-points, let's see what happens if the

initial line is 1 α 2 ... α N instead of 1 α 2 ... α N : then the final line is m 1 ... m N instead of m 1 ... m N .
Then (see appendix A)

       m j ∼ m j + 2 σ 2 +∞ -∞ ... +∞ -∞ (∆f x (α)p(f x (α))) x j p(x 2 )...p(x N )dx 2 ...dx N , ∀j = 1 m 1 ∼ m 1 + 2 σ 2 +∞ -∞ ... +∞ -∞ (∆f x (α)f x (α)p(f x (α))) p(x 2 )...p(x N )dx 2 ...dx N . ( 16
)
with

f x (α) = - N i=2 α i x i and ∆f x (α) = f x (α) -f x (α ).
Let us study this expression in the both following cases:

ratio ε j becomes m j m 1 ∼ 2p(0) m 1 ε j .
Evolution to separation point depends on the term 2p(0) m 1 according to wether it is greater, equal or smaller than 1; these quantity only depends on PDF's of sources.

A direct application of this study around other fixed-points is given in appendix C: three particular cases (uniform, gaussian and Laplace's PDF's) are studied.

C.2.b Study of the algorithm around other fixed-points.

Let us consider one observation Y 1 , i.e. one line of the mixture matrix M , which is a linear mixture of q sources with the same coefficient:

Y 1 = q k=1 X k .
1 ... 1 1 0 ... 0 (with q 1 and (N -q) 0) is a fixed-point of the algorithm. Then we wonder what will happen if we start from 1 ... 1 1 -ε 0 ... 0 . Starting from 1 ... 1 0 ... 0 leads to m 1 ... m q-1 m q 0 ... 0 with obviously m 1 = ... = m q-1 = m q . Starting from 1 ... 1 1 -ε 0 ... 0 , one step of the algorithm leads to m 1 ... m q-1 m q 0 ... 0 with m 1 = ... = m q-1 .

April 11, 2003 DRAFT

According to paragraph VI-B.2, the presence of (N -q) 0 do not change future calculation of the following proof. There are always (N -q) 0 (at same position) in resulting mixture line. Study is the same than those of a q-length line 1 ... 1 1 -ε .

Noting m = m q m 1
, three cases are possible:

• m < 1 -ε: then algorithm tends to go away from the fixed-point; this is a repulsive fixed-point.

• m > 1 -ε: then algorithm tends to come closer to the fixed-point; this is an attractive fixed-point.

• m = 1 -ε: this is a neutral fixed-point.

In this case f x (α) = -q i=2

x i and ∆f x (α) = εx q , and from equation ( 16),

       m 1 ∼ m 1 -2 σ 2 ε +∞ -∞ ... +∞ -∞ x q q k=2 x k × p - q k=2 x k p(x 2 )...p(x q )dx 2 ...dx q , m j ∼ m 1 + 2 σ 2 ε +∞ -∞ ... +∞ -∞ x q x j × p - q k=2
x k p(x 2 )...p(x q )dx 2 ...dx q , ∀j ∈ {2, ..., q}.

(

) 17
Noting ∆m i = m i -m1 ε , we observe that ∆m 1 = - q k=2
∆m k and q k=1 m k = qm 1 . Furthermore, for obvious reasons of symmetry, necessarily ∆m 1 = ... = ∆m q-1 , implying ∆m q = -(q -1)∆m 1 .

Let us find the conditions to be verified for m to be less than 1 -ε:

m = m q m 1 < 1 -ε ⇐⇒ m 1 + ε∆m q m 1 + ε∆m 1 < 1 -ε ⇐⇒ ∆m q -∆m 1 m 1 < -1 ⇐⇒ ∆m 1 m 1 > 1 q .
Evolution to separation point depends on wether the term ∆m 1 m 1 is greater, equal or smaller than 1 q ; these quantity depends on PDF's of sources. We note that this quantity also depends on q. A study of evolution of this ratio according to q shows (see appendix D) that, according to PDF's of sources, if a fixed-point is repulsive (resp. attractive) for a value of q, it remains repulsive (resp. attractive) for all value larger than q.

A direct application of this study around other fixed-points is given in appendix E: two particular cases (uniform and Laplace's PDF's) are studied.

Conclusion

Fixed-points of the algorithm can be divided in two classes, attractive and repulsive ones. For uniform PDF, the attractive fixed-point is [1 0 ... 1], i.e. separation state. For Laplace's PDF, this is [1 ... 1]. The use of a perturbation method which consists in applying a known transformation (a rotation by example) to observations ensures that even if initial point is one of the repulsive fixed-points, convergence will lead to the attractive fixed-point. This method must be applied if after one step of algorithm, nothing has changed.

C.2.c Conclusion. This section of analysis of convergence of the algorithm stresses several important points:

April 11, 2003 DRAFT

• processing is done by line (nevertheless coherence of columns is ensured by orthogonalizing P at each step because M is invertible by assumption).

• its behavior closely depends on PDF's of sources: this limits naturally the range of proofs. Then study has been done close to fixed-points: separation points and others.

• for uniform PDF, separation state is the sole attractive fixed-point: other fixed-points are repulsive.

• for Laplace's PDF, the sole attractive fixed-point is [1 ... 1].

• a perturbation method can be used to avoid freezing on a repulsive fixed-point.

In next section, computer simulation with digital signals will illustrate the power of this method.

VII. Simulations

Algorithm has been tested with digital signals. Throughout this section all signals will be understood to have the same PSD (i.i.d.), unless the contrary is specifically stated. Theoretically, sources and mixture matrix are unknown; however, as performed usually, we will use this information to quantify the relevance of results.

Sources and observations can be represented by N -by-L matrices:

X = X 1 X 2 ... X N T and Y = Y 1 Y 2 ... Y N T where X i is an L-dimensional vector x i (1) x i (2) ... x i (L) T .

A. Practical considerations

Simulations are completed with matlab . Matlab functions like rand and randn are used respectively to create time samples of uniformly distributed (as well as binary or four states) or gaussian sources. In appendix F are given the processes to generate sources with a Laplace's PDF and to generate sources with same PDF but different PSD.

Practically, to reach the condition Y i > 0 (i.e. N k=1 α ik X k > 0), we process an extraction of a sequence from observation Y i ; for example let's use the first one Y 1 : we select indices k=1,...L for which y 1 (k) is positive, and create an L c1 -dimensional vector c 1 whose entries are these indices. Then, for j = 1, ...N , we are able to create the j th extracted observations Y c 1 ,j (which is an L c 1 -length vector) defined by

Y c 1 ,j = y j (c 1 (1)) ... y j (c 1 (L c1 )) T , and conditional observations Y c 1 (which is a N × L c 1 matrix): Y c 1 = Y c 1 ,1 • • • Y c 1 ,N T . example : Y =   1 4 -1 3 -1 -5 -4 5 1 -5 3 2   ; then c 1 = 1 2 4 and Y c 1 =   1 4 3 -4 5 -5   .

B. Results and comments

For a chosen value of N and given types of PDF (specified in each case), iterations are stopped when the value of ind(P ) is less than a small value ε chosen beforehand. To quantify results, we choose the April 11, 2003 DRAFT following criterions: the value of ind(M (∞) ) where M (∞) is the mixture matrix at the end of convergence, and the mean square error (mse) between real and restored sources. Values presented in the table below are average values obtained from 500 experiences. The mean number of iterations necessary to reach the stopping condition is also given.

For each realization, mixture matrix (invertible) and sources are chosen at random.

Sources with a uniform PDF

Results are summarized in Table I. 1.9 10 -4 ± 3.7 10 -5 10 1.9 10 -4 52 ± 15 1.9 10 -5 1. Let's note that simulations with sources with a limited number of states (digital communications) have been realized. Results are very good and convergence is all the better when the number of states decreases.

In such cases, the mean square error need not to be very small, and a reasonable small value can easily be neglected a posteriori.

remark The mean number of iterations necessary to reach the stopping condition is given with a standard deviation. Nevertheless, this information is incomplete as it does not permit to realize how performance index converges. Hereafter are given two examples of this evolution versus the iteration number in cases of 2 and 10 uniformly distributed sources (figures 3 and 4).

April 11, 2003 DRAFT II.

C. Comparison with a classical HOS method

Let's compare numerically our algorithm's behavior to a classical HOS one, say JADE algorithm. This comparison consists in launching the both algorithms for identical data (sources and mixture matrix).

The criterion to evaluate these methods is the value of ϑ = ind M (∞) , performance index of the final mixture matrix.

Tables above give results obtained with the both methods for several kinds of sources' PDF's. When an acceptable separation state is reached, CPU times are presented to give the reader a rough idea of time convergence: its significance does not lie in its absolute values but in a comparison between magnitudes for the both methods.

It sometimes happens that JADE algorithm converges to a bad solution (e.g. in Table IV for N = 40 and L = 3072 where ind M (∞) = 0.2) ; this case is indicated by an asterisk. Let's note that a value of ind M (∞) greater than 0.01 is a bad one (separation is far from being reached).

For a given number of sources N , let's see what happens as the number of time samples L decreases.

Uniformly distributed sources

Let us consider 40 uniformly distributed sources. Results are summarized in Table IV.

Conclusion

Some conclusions can be deduced from the observation of these few numerical experimental results.

For a given number of sources N , both methods work all the better that the number of time samples L of each source increases: the price to pay is increasing CPU time. In fact, when N increases, L must increase significantly for results to stay as satisfactory.

For a large value of L (relatively to N ), JADE algorithm gives better results: in such a case, JADE algorithm is very efficient and it is difficult to compete with it.

But, as said before, the advantage of our method appears when N increases while L does not (or equivalently when L decreases while N is fixed). A decreasing value of L involves a slow degradation of results; under a threshold of L, different according to the algorithm, the algorithm fails. The threshold for JADE algorithm is higher than for ours, thus there exists an area of values of L for which JADE algorithm totally fails while ours continues to give acceptable results. In this area, our algorithm brings a substantial improvement.

Naturally when L is too small, both methods fail.

VIII. Conclusions

In this paper we present a new method for blind source separation problem when PSD's of sources are the same; in such cases it is well known that second-order methods are ineffective. In fact, the method proposed does not depend on PSD's of sources and also gives good results in the case of different PSD's.

Thus this algorithm can be used for both cases.

The method proposed presents several advantages:

• Contrary to geometrical methods [START_REF] Barrere | Chabriel The Archimedes principle applied to separation of uniformly distributed sources, SEE PSIP[END_REF] [25] [START_REF] Puntonet | A competitive neural network for blind separation of sources based on geometry properties, biological and artificial computation: from neuroscience to technology[END_REF] or neural network methods [START_REF] Basak | Blind separation of a mixture of uniformly distributed source signals: a novel approach[END_REF][15], it can process sources whose PDF has an unbounded domain of definition (e.g. Laplace's PDF).

April 11, 2003 DRAFT

• Calculations are limited to only first and second-order moments.

• It is robust w.r.t. the number of sources (relatively to the number of time samples of each source).

• Despite a presentation which could seem complicated, the algorithm is very easy to program and necessitates few calculation.

For reasons explained in the paper we have limited our study to symmetrical PDF's. Two classes of signals have been more precisely detailed: uniformly distributed sources (digital communications) and sources with a Laplace's PDF.

The state reached at the end of convergence depends only on PDF's of sources: this state can be a separation state (e.g. uniformly distributed sources) or an a priori known mixture state (e.g. sources with a Laplace's PDF).

Speed of convergence is different according to PDF (in fact it depends on the distance to a gaussian PDF) and number of sources; for example, for sources with a triangular PDF, convergence is slower than with uniform ones.

Practically convergence is always reached. However when the number of sources grows, mean-square error is larger; as our objective is not to find exactly sources, but an approximation as satisfying as possible, this is quite acceptable for numerous engineering applications.

For uniformly distributed sources (with any number of states) results are very convincing: the algorithm converges directly to a separation state. In this case, it will be doubtless possible to establish a link with contrast functions [START_REF] Pham | Contrastes pour la séparation aveugle de sources[END_REF].

Numerical comparisons with a commonly used HOS method (JADE) are provided. They show that for a given number of sources, our method continues to give acceptable results while JADE algorithm when the number of time samples available decreases.

On the other hand, for sources with Laplace's PDF, as at the end of convergence we obtain an a priori known mixture state, a final operation is necessary to reach the separation state: in this case it is hardly possible to talk of contrast function. However a generalization of contrast functions could be glimpsed in the light of this particular example. As a prospect, this question needs to be examined in more details: much work remains to be done concerning a link between our method and contrast functions.

Another prospect is the study of convergence of the algorithm when PDF's of sources are no more symmetric. Also a theoretical study of links with HOS methods must be performed in the future. Obviously, application of this algorithm (with adaptations) to convolutive mixture problem with sensor noise must be the subject of a detailed study in the near future.

Notation

p i (x) PDF of source X i σ 2 i power of source X i D i domain of definition of p i (x) 1 D (x) 1 if x ∈ D, 0 otherwise. Appendix A

Study for any number of sources

Expression of m j is

m j = 2 σ 2 +∞ -∞ ... +∞ -∞       +∞ - N P i=2 α i x i x j p(x 1 )dx 1       p(x 2 )...p(x N )dx 2 ...dx N ,
which can be written making appear those of m j :

m j = m j + 2 σ 2 +∞ -∞ ... +∞ -∞       - N P i=2 α i x i - N P i=2 α i xi x j p(x 1 )dx 1       p(x 2 )...p(x N )dx 2 ...dx N , or, noting f x (α) = - N i=2 α i x i , m j = m j + 2 σ 2 +∞ -∞ ... +∞ -∞    f x (α) fx(α ) x j p(x 1 )dx 1    p(x 2 )...p(x N )dx 2 ...dx N .
April 11, 2003 DRAFT

As f x (α ) and f x (α) are close enough (∆f x (α) = f x (α) -f x (α )), m j can be approximated by

m j ∼ m j + 2 σ 2 +∞ -∞ ... +∞ -∞ (∆f x (α)x j p(x 1 )) |x 1 =f x (α) p(x 2 )...p(x N )dx 2 ...dx N . Then        m j ∼ m j + 2 σ 2 +∞ -∞ ... +∞ -∞ (∆f x (α)p(f x (α))) x j p(x 2 )...p(x N )dx 2 ...dx N , ∀j = 1 m 1 ∼ m 1 + 2 σ 2 +∞ -∞ ... +∞ -∞ (∆f x (α)f x (α)p(f x (α))) p(x 2 )...p(x N )dx 2 ...dx N . 2 Appendix B

Study close to a point of separation

Obviously, ∀j = 1, m j = 0.

It comes immediately that f x (α) = 0 and ∆f

x (α) = N i=2 ε i x i . Then      m j ∼ 2 σ 2 +∞ -∞ ... +∞ -∞ N i=2 ε i x i p(0) x j p(x 2 )...p(x N )dx 2 ...dx N , ∀j = 1 m 1 ∼ m 1 .
Then

m j ∼ 2 σ 2 p(0) N i=2 ε i +∞ -∞ ... +∞ -∞ x i x j p(x 2 )...p(x N )dx 2 ...dx N , ∀j = 1, or m j ∼ 2 σ 2 p(0) N i=2 ε i E(X i X j ), ∀j = 1. As E(X i X j ) = 0 if i = j, m j ∼ p(0) 2 σ 2 E X 2 j ε j .
The initial ratio ε j becomes

ε j = m j m 1 ∼ 2p(0) m 1 ε j . 2 Appendix C
Application of the study of the algorithm around a point of separation to particular cases example 1: uniform PDF

If p(x) = 1 2a 1 [-a;+a] (x) then p(0) = 1 2a , σ 2 = a 2 3 and m 1 = a 2σ 2 . Then 2p(0) m 1 = 2 3
. [1 0 ... 0] is an attractive point of the algorithm. April 11, 2003 DRAFT In fact, in this case it is possible to find exact expression of m 1 and m j :

     m 1 = 2 σ 2 2 N -1 3 3 2 -1 2 N i=2 ε 2 i , m j = 2 σ 2 2 N -1 3 ε j , ∀ j = 1.
Obviously

Study of the algorithm around other fixed-points

We will note m (q) j the value of m j given by [START_REF] Cruces-Alvarez | Castedo-Ribas An iterative inversion approach to blind source separation[END_REF] with α i = 1 ∀i = 1, ...q. m (q)

j = 2 σ 2 +∞ -∞ ... +∞ -∞       +∞ - q P i=2 xi x j p(x 1 )dx 1      
p(x 2 )...p(x q )dx 2 ...dx q . m (q)

j = 2 σ 2 +∞ -∞ ... +∞ -∞       +∞ - q-1 P i=2 x i x j p(x 1 )dx 1       p(x 2 )...p(x q )dx 2 ...dx q + 2 σ 2 +∞ -∞ ... +∞ -∞       - q-1 P i=2 x i - q P i=2 xi x j p(x 1 )dx 1       p(x 2
)...p(x q )dx 2 ...dx q .

First term is m (q-1) j and second one is approximately (if q is large):

+∞ -∞
...

+∞ -∞

x q x j p -q i=2

x i p(x 2 )...p(x q )dx 2 ...dx q , April 11, 2003 DRAFT i.e. -∆m j . Then m (q) j ∼ m (q-1) j -∆m j . By hypothesis, ∆m j > 0 and m (q) j > 0. Then ∆m (q) j m (q) j ∼ ∆m (q) j m (q-1) j -∆m (q) j .

Considering that ∆m (q) j ∼ ∆m (q-1) j

, it comes ∆m (q) j m (q) j ∼ ∆m (q-1) j m (q-1) j -∆m (q-1) j = ϕ 1 -ϕ , with ϕ = ∆m (q-1) j m (q-1) j < 1. If ϕ > 1 q-1 (respectively < 1 q-1 ), then ∆m (q) j m (q) j > 1 q (respectively < 1 q ). In conclusion, according to PDF's of sources, if a fixed-point is repulsive (resp. attractive) for a value of q, it remains repulsive (resp. attractive) for all value larger than q. 2 Appendix E Application of the study of the algorithm around other fixed-points to particular cases Case of Laplace's PDF:

For q = 2, from [START_REF] Jutten | Blind separation of sources, Part I : An adaptative algorithm based on neuromimatic architecture[END_REF] Calculations give m 1 = 3 8 and ∆m 1 = 1 8 , then ∆m 1 m 1 = 1 3 < 1 2 = 1 q : [1 1] is an attractive point of the algorithm. Calculations give ∆m1 m 1 = 1 4 < 1 3 = 1 q ; [1 1 1] is an attractive point of the algorithm. 2 Case of uniform PDF:

For q = 2, calculations give ∆m 1 m 1 = 1 > 1 2 ; [1 1] is a repulsive point of the algorithm. In the case of uniformly distributed sources, a systematic analysis of bounds of integration allows to calculate coefficients of new observation versus those of initial one. This approach can readily be systematically applied to find the integrals needed to calculate conditional moments. Unfortunately number of integrals inflates quickly with q (3 for q = 4, 10 for q = 5 or 6, 35 for q = 7 or 8, 126 for q = 9 or 10,...).

As only ratios of m j 's are important to determine the behavior of the algorithm, we will ignore the constant term 2 σ 2 in the expressions hereafter. Example for q = 2 The first constraint is: ∀i, -1 ≤ x i ≤ 1. The equation of Π is "x 1 +αx 2 = 0" with α = 1 -ε 1 (ε 0). x 1 + αx 2 ≥ 0 ⇐⇒ x 1 ≥ -αx 2 . Then it is necessary to have -1 ≤ -αx 2 ≤ x 1 ≤ 1 and as α < 1, -1 ≤ -αx 2 ≤ 1 is always true. Consequently, m j = 1 ...dx 1 dx 2 dx 3 .

Results obtained are the following:

• m 1 = 1 8 13 12
• m 1 = m 1 + 1

  Starting from the point α ... α 0 ... 0 ∼ 1 ... 1 0 ... 0 one step of the algorithm leads to m 11 ... m 1q 0 ... 0 . Provided that PDF's of sources are all identical, it is equivalent to 1 ... 1 0 ... 0 . In this case, this is a fixed-point of the algorithm.

Fig. 1 .

 1 Fig. 1. Uniform PDF: new performance index vs. old performance index

Fig. 2 .

 2 Fig. 2. Laplace's PDF: new performance index vs. old performance index

  ... α iN , leads to the line of M(2) : M (2) i = m i1 ... m iN . For the sake of simplicity, we will note M (1) i = α 1 ... α N , and M (2) i = m 1 ... m N . As said before, it is always possible to normalize each equation of hyperplane: we will suppose that |α 1 | > |α i | ∀i: M (1) i becomes 1 α 2 ... α N . April 11, 2003 DRAFT

  around a separation point, around other fixed-points. C.2.a Study close to a point of separation. In this case 1 α 2 ... α N = 1 0 ... 0 and 1 α 2 ... α N = 1 ε 2 ... ε N with |ε i | 1: the observation y i = x 1 + N k=2 ε k x k is almost one of sources (here x 1 ) with residual contribution of others. It comes (see appendix B) that the initial

3 10 - 3 ± 8 .6 10 - 5 †

 38105 This case illustrates that the algorithm works as well with sources which have different PSD's. remark concerning the robustness w.r.t. number of sources Results obtained for N = 10 are convincing; the algorithm is robust w.r.t. the number of sources, contrary to H.O.S. methods which fail quickly when N increases. Examples with 40 uniformly distributed sources have also been processed successfully.

Fig. 3 .Fig. 4 .

 34 Fig. 3. N=2 : example of evolution of performance index vs. number of iterations

1 √ 2 √ 2πσ 2 .

 122 |m 1 | |m j | ∀j = 2, ...N , and ratio . Then, starting from 1 ε 2 ... ε N leads approximately to 1 2 3 ε 2 ... 2 3 ε N . There is an improvement in the sense that relative distance (between starting point and separation point) has decreased. Progression of decrease is geometric of ratio 2 3 ; there is convergence to [1 0... 0] which is consequently an attractive fixed point of the algorithm. 2 example 2: gaussian PDF of power σ 2 In this case, p(0) = 2πσ 2 and m 1 = Then 2p(0) m1 = 1 and as expected to, [1 0... 0] is a neutral point of the algorithm (in fact like any point). 2 example 3: Laplace PDF In this case, p(x) = a 2 e -a|x| then p(0) = a 2 . σ 2 = 6 a 2 and m 1 = 2 a ; then 2p(0) m 1 = 3. [1 0... 0] is a repulsive point of the algorithm. 2 Appendix D

2 x 1 p 2 (x 2 )

 2122 the expressions of m 1 and ∆m 1 are m 1 = 1 p(x 1 )dx 1 p(x 2 )dx 2 and ∆m dx 2 .

For q = 3 ,x 1 x 3

 313 p(x 1 )dx 1 p(x 2 )p(x 3 )dx 2 dx 3 and ∆m 1 = (x 2 +x 3 )p(x 2 +x 3 )p(x 2 )p(x 3 )dx 2 dx 3 = x 2 + x 3 )p(x 2 )p(x 3 )dx 2 dx 3 .

2 x 1 - 1 -αx 3 1 -x 2 -αx 3 ...dx 1 dx 2 dx 3 + 1 01-αx 3 - 1 1-x 2 -αx 3 ...dx 1 dx 2 dx 3

 2111313133 j dx 1 dx 2 , for j = 1, 2.April 11, 2003 DRAFTExample for q = 3 Each conditional moment will be calculated with three triple integrals:

1 = 1 1 15 • 1 = 1 1 = 1

 111151111 -15 13 ε + O ε 2 < 1 -ε: [1 1 1] is a repulsive fixed-point. This is confirmed by the calculation of ∆m 1 m 1 : ∆m 1 m for q = 4Results obtained are the following:• m 1 = m 2 = m 3 = m 4 = 1 16 28 m 1 = m 2 = m 3 = m 1 + -8 7 ε + O ε 2 < 1 -ε: [1 1 1 1] is a repulsive fixed-point. This is confirmed by the calculation of ∆m 1 m1 : -1315 1199 ε + O ε 2 < 1 -ε: [1 1 1 1 1] is a repulsive fixed-point.This is confirmed by the calculation of ∆m1 m 1 : ∆m1 m 1 ∼ 263 360 / 1199 360 = 263 1199 > 1 5 . 2

  Result presented here aims at illustrating that the method proposed can as well process mixtures of sources with different PDF's (provided that at the most one gaussian source is present). Results are summarized in TableIII.

			TABLE II
			Sources with a Laplace's PDF
	N ind(M (∞) ) Number of	ind(P )	mse (on one source)
			iterations
	2	8.8 10 -5	23 ± 8	1.9 10 -5	7.4 10 -7 ± 5.3 10 -7
	3	1.2 10 -3	50 ± 16	1.0 10 -5	1.0 10 -5 ± 1.9 10 -5
	Sources with different PDF's
	April 11, 2003				DRAFT

  ± 6.2 10 -4 † first source has a uniform PDF, second one a Laplace's PDF. ‡ sources have respectively a uniform PDF, a uniform PDF with a different PSD, a Laplace's PDF and a gaussian PDF.

			TABLE III
			Sources with different PDF
	N	ind(M (∞) ) Number of	ind(P )	mse
			iterations
	2 † 4.1 10 -5	17 ± 5	4.4 10 -6	1.2 10 -5 ± 7.3 10 -6
	4 ‡ 6.3 10 -3	32 ± 14	4.1 10 -6	7.3 10 -4

  We can observe in this table that, as expected, for a given value of N , the quality of results decreases with L. When L is too small, both methods fail.

			TABLE IV
	Comparison with HOS method for 40 uniformly distributed sources
		JADE JADE our	our
				algo.	algo.
	L	ϑ	CPUt ϑ	CPUt
	40960 0.000011 1884	0.00009 920
	20480 0.000024 1050	0.00013 677
	10240 0.00005 503	0.00019 392
	6144	0.0001 381	0.00034 710
	3072	0.2	*	0.0007 321
	2048	0.2	*	0.06	226
	1024	0.2	*	0.12	*
	Binary sources
	Let us consider 40 binary sources, i.e. each time sample is equal to -1 or +1. Results are summarized
	in Table V.

TABLE VI

 VI Comparison with HOS method for 30 Four states sources

		JADE JADE our	our
				algo.	algo.
	L	ϑ	CPUt ϑ	CPUt
	4096	0.00007 83	0.00014 248
	2048	0.00015 63	0.00028 141
	1536	0.00023 101	0.00036 106
	1200	0.203	*	0.00047 85
	1024	0.236	*	0.00056 79
	512	0.24	*	0.0011 46
	256	0.234	*	0.154	32

[START_REF] Pham | Blind separation of instantaneous mixture of sources via an independent component analysis[END_REF] DRAFT

X, so that the PDF is unaltered. Keeping the first time sample of X, the second one swaps with another one among the p following samples: we choose the sample whose value is the closest to the first sample.

And so on, from the second time sample. p must be chosen judiciously, neither too small nor too large.

Generation of sources with a Laplace's PDF

To generate n time samples of a source Z with a Laplace's PDF, i.e. p(z) = 1 2 ae -a|z| , let's start from time samples generated by the matlab function rand.

-X=rand (1,n);

Generation of a source with same PDF and different PSD Starting from a source X with a given PDF, it is possible to generate another one X with the same PDF but a different PSD, by applying a simple process. This process consists in swapping time samples of [START_REF] Pham | Blind separation of instantaneous mixture of sources via an independent component analysis[END_REF] DRAFT