Rémi Dubroca

Member IEEE Christophe De Luigi

Member IEEE Marc Castella

Senior Member IEEE Eric Moreau

A General Algebraic Algorithm for Blind Extraction of one Source in a MIMO Convolutive Mixture

Keywords: Contrast Functions, Blind Source Extraction, Higher Order Statistics, Tensor Decomposition, Independent Component Analysis

The paper deals with the problem of blind source extraction from a MIMO convolutive mixture. We define a new criterion for source extraction which uses higher-order contrast functions based on so called reference signals. It generalizes existing reference-based contrasts. In order to optimize the new criterion, we propose a general algebraic algorithm based on best rank-1 tensor approximation. Computer simulations illustrate the good behavior and the interest of our algorithm in comparison with other approaches.

I. INTRODUCTION

T HE problem of blind extraction of one source signal appears in a wide range of signal processing applications such as bio-medical and telecommunications. The aim is to restore one source signal from the observation of a set of mixed sources. Among existing mixing models, we consider a Multiple-Input / Multiple-Output (MIMO) context where non observable source signals are mixed through a multidimensional convolutive channel. The separation is said to be performed blindly when the mixing system is unknown and cannot be identified.

Such mixtures are considered in telecommunication applications like cellular communications. The mobile user signals are mixed through the air channel and observed on the different antennas of the base station. In this context, when the extraction of one mobile user signal has to be performed, the problem cannot be treated as Single-Input / Single-Output blind equalization because the other mobile user signals cannot be considered as noise: they are interfering signals.

The mobile user of interest is extracted from the mixture of all users. The choice of a blind method for extraction of the source potentially avoids the necessity of a training sequence and therefore lowers bandwidth loss.

The blind source extraction problem may be considered as part of the general blind source separation (BSS) problem. Iterative methods in the context of BSS, extract the sources one by one and therefore rely on the ability to extract one source. These so-called deflation methods, presented in [START_REF] Loubaton | Blind Deconvolution of Multivariate Signals: A Deflation Approach[END_REF], developed in [START_REF] Inouye | Iterative Algorithms Based on Multistage Criteria for Multichannel Blind Deconvolution[END_REF], [START_REF] Tugnait | Identification and Deconvolution of Multichannel Linear Non-Gaussian Processes Using Higher Order Statistics and Inverse Filter Criteria[END_REF], [START_REF] Delfosse | Adaptive blind separation of independent sources: A deflation approach[END_REF], [START_REF] Castella | Quadratic Higher-Order Criteria for Iterative Blind Separation of a MIMO Convolutive Mixture of Sources[END_REF], decompose the separation problem in successive stages for each source. These methods can deal both with temporally independent identically distributed (i.i.d.) sources and temporally non i.i.d. sources [START_REF] Simon | Separation of a Class of Convolutive Mixtures: A Contrast Function Approach[END_REF]. Hence the blind source extraction problem is a critical part in these methods, the separation relies on the ability to extract well each source.

The former method is in contrast to a global separation approach that has been studied in a convolutive context (see [START_REF] Comon | Contrasts for Multichannel Blind Deconvolution[END_REF], [START_REF] Moreau | Generalized Contrasts for Multichannel Blind Deconvolution of Linear Systems[END_REF], [START_REF] Jelonnek | Generalized eigenvector algorithm for blind equalization[END_REF] for the i.i.d. case and [START_REF] Comon | Blind Separation of Independent Sources from Convolutive Mixtures[END_REF] for the non i.i.d. case). The separation is not performed sequentially, but all the sources are extracted simultaneously. Such an approach faces two major obstacles in practice. First, a prewhitening is required, which is a difficult task. In addition, the optimization scheme of the criteria is likely to converge to a wrong solution because of the existence of local spurious maxima [START_REF] Castella | A family of frequency-and time-domain contrasts for blind separation of convolutive mixtures of temporally dependent signals[END_REF].

Contrast functions based on higher-order statistics have been proved to be efficient for both blind source separation or extraction [START_REF] Comon | Contrasts for Multichannel Blind Deconvolution[END_REF], [START_REF] Moreau | A Generalization of Joint-Diagonalization Criteria for Source Separation[END_REF]. Here we focus on Multiple-Input / Single-Output (MISO) extraction criteria relying on higherorder statistics. One potential drawback is the high order dependence on the parameters. Recently "reference"-based contrasts have been investigated to lower this dependence or to improve the extraction performance [START_REF] Adib | Source Separation Contrasts using a Reference Signal[END_REF], [START_REF] Castella | Quadratic Higher-Order Criteria for Iterative Blind Separation of a MIMO Convolutive Mixture of Sources[END_REF], [START_REF] Kawamoto | Eigenvector Algorithms Incorporated with Reference Systems for Solving Blind Deconvolution of MIMO-IIR Linear Systems[END_REF]. These approaches impose some constraints on the reference signals.

In order to extract one source, we propose in this paper a family of contrast functions based on higher-order cumulant and using reference signals. These contrast functions impose no constraint on the reference signals. We show the usefulness of tensorial algorithms applied to independent component analysis (ICA) [START_REF] Lathauwer | Higherorder power method -application in independent component analysis[END_REF]. The optimization problem shows an equivalence to the problem of best rank-1 approximation of tensors [START_REF]On the Best Rank-1 and Rank-(R 1 ,R 2 , . . . ,R N ) Approximation of Higher-Order Tensors[END_REF], [START_REF] Lathauwer | A Multilinear Singular Value Decomposition[END_REF], [START_REF] Dubroca | Cubic Higher-Order Criterion and Algorithm for Blind Extraction of a Source Signal[END_REF]. This optimization scheme is also applied to the classical kurtosis contrast function [START_REF] Simon | Separation of a Class of Convolutive Mixtures: A Contrast Function Approach[END_REF]. The computation of these best rank-1 approximations is done through the use of the Higher-Order Power Method (HOPM) algorithm (higher-order generalization of the power method). The novel results in this paper are the following ones: we introduce a new family of reference based contrasts, we interprete tensor approximation in BSS and we propose new optimization method adapted to the proposed criteria.

The problem is formulated in Section II and Section III introduces the new contrast function. Some recalls on Tensor algebra and decomposition from [START_REF]On the Best Rank-1 and Rank-(R 1 ,R 2 , . . . ,R N ) Approximation of Higher-Order Tensors[END_REF] are done in Section IV and an optimization algorithm for our new criterion is presented in Section V. Section VI illustrates through computer simulations the good behavior of the proposed approach. Finally we conclude this paper in Section VII.

Notations. Tensors are denoted by bold-face calligraphic letters (A, B, . . .), matrices by bold-face capital letters (A, B, . . .), vectors by bold-face lower-case letters (a, b, . . .), and scalars by normal letters (a, b, . . . ; α, β, . . .). In the whole paper, n stands for a generic integer (n ∈ Z).

All quantities throughout the paper may be either real or complex-valued. (.) * stands for the complex conjugate, (.) T for matrix transpose, (.) H for matrix conjugate transpose. Finally Cum{.} stands for the the cumulant of any set of random variables, and δ i,j stands for the Kronecker symbol, i.e. δ i,j = 1 if i = j and 0 otherwise.

II. PROBLEM FORMULATION

We consider a linear time-invariant (LTI) MIMO mixing model with K inputs and N outputs

x(n) = k∈Z M(k)s(n -k), (1)
where x(n) is the (N, 1) observation vector (N 2), s(n) is the (K, 1) source vector (K 2), and M(n) is the (N, K) matrix corresponding to the impulse response of the convolutive mixing system. The transfer function of the LTI mixing system is denoted by

M(z) = n∈Z M(n)z -n , (2)
and we assume N K.

To achieve the extraction of one source blindly, that is with no access to the source signals nor to the mixing system filter, we have to make some assumptions on the LTI mixing system and the sources:

A1. The LTI mixing system is stable, left invertible. Therefore it exists a system W(n) such that the global LTI system with impulse response

G(n) k∈Z W(k)M(n -k) (3)
corresponds to the identity system.

A2. The source signals s i (n), i ∈ {1, . . . , K} are zeromean, unit-variance, random stationary signals. Moreover they are statistically mutually independent (at least up to the order of the considered cumulants). At least one of the sources has non-zero 4-th order auto-cumulant, i.e. there exists i ∈ {1, . . . , N } such that

Cum{s i (n), s i (n), s i (n), s i (n)} = 0.
In the MISO context, the aim is to estimate a (1, N ) row vector filter w(n) such that the scalar signal

y(n) = k∈Z w(k)x(n -k) (4)
restores one of the sources s i (n), i ∈ {1, . . . , K}, up to a non-zero scalar filter. It is useful to define the global vector filter g(n) by

g(n) k∈Z w(k)M(n -k), (5)
thus

y(n) = k∈Z g(n -k)s(k) (g ⋆ s)(n). (6)
As the source signals are unobservable, there exists some inherent undetermined factors in their estimation. They can be recovered only up to a permutation and to a scalar filtering. So, the extraction of one source is said to be achieved when there exists an index i 0 ∈ {1, . . . , K} such that the filter components in g(n) read ∀i ∈ {1, . . . , K} g

i (n) (g(n)) i = g(n)δ i,i0 . (7)
The above relation is called the "extraction condition" and expresses the fact that y(n) is equal to the source signal, s i0 (n) up to a filtering by the scalar filter with impulse response g(n).

The source signal correlation sequences are denoted by

γ i (k) E{s i (n)s * i (n -k)},
where k ∈ Z and i ∈ {1, . . . , K}. For any index j ∈ {1, . . . , K} and any scalar filter with impulse response h(n), we define its j-norm by:

h 2 j (k1,k2)∈Z 2 h(k 1 )h * (k 2 )γ j (k 2 -k 1 ). (8)
h 2 j is the variance of the signal obtained when filtering the source s j (n) by the filter with impulse response h(n). We define the (weighted) l 2 -norm of the (1, K) global filter g(n) by:

g 2 K i=1 g i 2 i . (9)
Since the sources are unit power, one can restrict to the unit-norm global filter case in [START_REF] Lathauwer | A Multilinear Singular Value Decomposition[END_REF] by imposing the constraint E{|y(n)| 2 } = 1. It is straithforward to show that this constraint is equivalent to g 2 = 1. In the particular case of i.i.d. sources, this constraint is equivalent to:

K i=1 k∈Z |g i (k)| 2 = 1. (10)
Moreover, in the i.i.d. case, the contrast functions (defined in the next section) allow one to equalize the sources. Consequently we define the stronger separation condition:

∃ l ∈ Z and ∃ i 0 ∈ {1, . . . , K} such that g i (n) = αδ n,l δ i,i0 . (11)
Here the scalar filter g(n) reduces to a delay l and a scaling factor α ∈ C * and ( 11) is actually a condition for separation and equalization.

The set of source signals satisfying assumption A2. will be denoted by S. The set of unit norm vector filters will be denoted by G 1 . Finally, we denote by Y the set of output of the MISO extractor y(n) when the input source signals belong to S and the global system belongs to G 1 .

III. EXTRACTION CRITERIA

A. Recalls on contrast function

A contrast is an objective function depending on the statistics of the outputs of the extracting system. Its maximization leads to the extraction of one source. In consequence, the source extraction problem becomes an optimization problem.

We can notice that a contrast is a function from Y to R but it may also be seen as a function from G 1 to R. In the i.i.d. case, the following definition is considered:

Definition III.1. Contrast function for i.i.d. sources

Let C be a function from Y to R. It is called a contrast function when it satisfies the following two properties:

p1. ∃ (i 0 , k 0 ) ∈ {1, . . . , K} × Z such that: ∀ y(n) ∈ Y C{y(n)} C{s i0 (n -k 0 )} . ( 12
)
p2. The equality holds in [START_REF] Gorokhov | Subspace-Based Techniques for Blind Separation of Convolutive Mixtures with Temporally Correlated Sources[END_REF] if and only if y(n) equals a multiple of s i (n -k) for some i ∈ {1, . . . , K} and some k ∈ Z (i.e. the global filter is extracting and statisfies [START_REF] Ghennioui | A Nonunitary Joint Block Diagonalization Algorithm for Blind Separation of Convolutive Mixtures of Sources[END_REF]).

In the case of i.i.d. source signals, the maximization of the contrast function restores one source up to a delay and a scaling factor.

However this definition cannot be used for the non i.i.d. case because the extraction of a source can be only guaranteed up to a scalar filtering. Hence we give a weaker definition for the non i.i.d. case:

Definition III.2. Contrast function for non i.i.d. sources

Let C be a function from Y to R. It is called a contrast function when it satisfies the following two properties:

p1. ∀ y(n) ∈ Y: C{y(n)} K max i=1 sup g:||gi||i=1 ||gj ||j=0 ∀j =i C{(g ⋆ s)(n)}. ( 13
)
p2. The equality holds in [START_REF] Inouye | Iterative Algorithms Based on Multistage Criteria for Multichannel Blind Deconvolution[END_REF] if and only if y(n) is a filtered version of s i (n) for some i ∈ {1, . . . , K} (i.e. the global filter statisfies [START_REF] Lathauwer | A Multilinear Singular Value Decomposition[END_REF]) for some i 0 ).

A classical contrast function defined in [START_REF] Simon | Separation of a Class of Convolutive Mixtures: A Contrast Function Approach[END_REF] is the absolute value of the kurtosis:

C 4 {y(n)} |Cum{y(n), y * (n), y * (n), y(n)}|. (14)
A contrast function based on referenced signal is defined in [START_REF] Castella | Quadratic Higher-Order Criteria for Iterative Blind Separation of a MIMO Convolutive Mixture of Sources[END_REF] by the absolute value of the fourth order cross-cumulant between the output of the extracting system and a reference signal r(n):

C 2,4,r {y(n)} |Cum{y(n), y * (n), r(n), r * (n)}|. (15)

B. New reference based contrast function

The main goal of the paper is to propose a family of contrast functions based on reference signals generalizing the contrast in [START_REF] Castella | Quadratic Higher-Order Criteria for Iterative Blind Separation of a MIMO Convolutive Mixture of Sources[END_REF]. For that, we consider the following R-th order (R 3) cross-cumulant where (.) ( * ) k , k ∈ {1, . . . , S} are optional complex conjugate:

κ S,R,r {y(n)} Cum{y ( * )1 (n), . . . , y ( * )S (n) S , r 1 (n), . . . , r R-S (n)} . (16)
The r i (n), i ∈ {1, . . . , R -S}, are given reference signals (2 < S ≤ R) which are assumed to satisfied the following assumption A3. The reference signals r i (n), i ∈ {1, . . . , R -S} are obtained by a stable MIMO filtering of the sources or of their complex-conjugates.

Note that the case S = 2 is treated in [START_REF] Castella | Quadratic Higher-Order Criteria for Iterative Blind Separation of a MIMO Convolutive Mixture of Sources[END_REF] and the case S = R is treated in [START_REF] Moreau | A Generalization of Joint-Diagonalization Criteria for Source Separation[END_REF]. In the following we consider the general case 2 < S ≤ R.

We now define the following function:

C S,R,r {y(n)} |κ S,R,r {y(n)}|. ( 17
)

1) Case of i.i.d. source signals:

The following supremum is assumed to satisfy the assumption A4. ∃(j 0 , k 0 ) such that:

N max j=1 sup k∈Z |κ S,R,r {s j (n-k)}| = |κ S,R,r {s j0 (n-k 0 )}| < ∞. ( 18
)
We can now state the following proposition (for 2 < S R):

Proposition III.3. In the case of i.i.d. source signals, under assumptions A1.-A4., the function C S,R,r is a contrast over the set Y.

The proof is given in Appendix A.

2) Case of non i.i.d. source signals:

We define the following supremum,

M max S,R K max i=1 M S,R,i (19)
where

M S,R,i sup g:||g||i=1 C S,R,r {(g ⋆ s)(n)}. (20)
We make the following assumption:

A5. For all i such that M S,R,i = M max S,R , the extremum M S,R,i is reached by at least one filter of unit i-norm defined in [START_REF]On the Best Rank-1 and Rank-(R 1 ,R 2 , . . . ,R N ) Approximation of Higher-Order Tensors[END_REF].

We can state (with 2 < S R):

Proposition III.4. In the case of non i.i.d. source signals, under assumptions A1.-A3. and A5., the function C S,R,r is a contrast over the set Y.

The proof is given in Appendix B.

We can notice that the observation signals satisfy the assumption A3. and can serve as reference signals. Note also that r 1 (n) = . . . = r R-S (n) is a possible choice. As we have seen in the previous propositions, we remark that the new contrast requires no condition on the reference signals except the assumption A3..

IV. RECALLS ON TENSOR DECOMPOSITION

Tensorial algorithms applied to ICA [START_REF] Lathauwer | Higherorder power method -application in independent component analysis[END_REF] are of interest and the optimization problem has an equivalence to the problem of best rank-1 approximation of tensors [START_REF]On the Best Rank-1 and Rank-(R 1 ,R 2 , . . . ,R N ) Approximation of Higher-Order Tensors[END_REF], [START_REF] Lathauwer | A Multilinear Singular Value Decomposition[END_REF]. We recall notions on tensor decomposition which are necessary in Section V. See [START_REF]On the Best Rank-1 and Rank-(R 1 ,R 2 , . . . ,R N ) Approximation of Higher-Order Tensors[END_REF] for more details in tensor algebra. The n-mode product of a tensor A ∈ C I1×...×IN by a matrix U ∈ C Jn×In , denoted by A• n U, is an (I 1 ×. . .×J n ×. . . I N )tensor whose entries are given by

(A • n U) i1i2...jn...iN in A i1i2...in...iN U jnin . ( 21
)
The outer product of N vectors u (1) , . . . , u (N ) of size I 1 , . . . , I N respectively is a N -order tensor in C I1×...×IN with elements defined by A i1...iN = u

(1)

i1 . . . u (N )
iN for all values of the indices. It is denoted by u (1) • . . . • u (N ) . A N -order tensor in C I1×...×IN is rank-1 if it can be written as the outer product of N vectors u (1) , . . . , u (N ) of size I 1 , . . . , I N respectively. The best rank-1 approximation of a N -order tensor can be described as follows. Given A ∈ C I1×...×IN , find a scalar λ and unit-norm u (1) , . . . , u (N ) such that the rank-1 tensor à ∈ C I1×...×IN defined by

à λ u (1) • . . . • u (N ) (22)
minimizes the least-squares cost function

f ( Ã) = A -Ã 2 . ( 23
)
The Frobenius norm of A is defined by A A, A with the scalar product A, B of two tensors A, B ∈ C I1×...×IN defined by

A, B i1 . . . iN A i1...iN B * i1...iN . ( 24
)
The minimization of the cost function f ( Ã) is equivalent to the maximization, over the unit-norm vectors u (1) , . . . , u (N ) , of the function:

g(u (1) , . . . , u (N ) ) = |A • 1 u (1) H • 2 . . . • N u (N ) H | 2 . ( 25
)
The scalar λ is reached by:

λ = A • 1 u (1) H • 2 . . . • N u (N ) H . (26)
A solution can be found by an ALS method which is a higherorder extension of the power method (HOPM) for matrices [START_REF]On the Best Rank-1 and Rank-(R 1 ,R 2 , . . . ,R N ) Approximation of Higher-Order Tensors[END_REF].

HOPM algorithm

Input: A ∈ C I1×...×IN Output: Ã ∈ C I1×...×IN 1) Initial values: u (n) 0 , n ∈ [1, . . . , N ]. 2) Repeat until convergence: • ũ(1) k+1 = A • 2 u (2) H k • 3 u (3) H k • 4 . . . • N u (N ) H k λ (1) k+1 = ũ (1)
k+1 and u

(1)

k+1 = ũ(1) k+1 /λ (1) k+1 • ũ(2) k+1 = A • 1 u (1) H k+1 • 3 u (3) H k • 4 . . . • N u (N ) H k λ (2) k+1 = ũ(2)
k+1 and u

(2)

k+1 = ũ(2) k+1 /λ (2) k+1 • . . . • ũ(N) k+1 = A • 1 u (1) H k+1 • 2 u (2) H k+1 • 3 . . . • N -1 u (N -1) H k+1 λ (N ) k+1 = ũ(N) k+1
and u

(N ) k+1 = ũ(N) k+1 /λ (N ) k+1
Solutions : u (1) , . . . , u (N ) , λ.

3) Ã = λ u (1) • . . . • u (N ) .
This iterative algorithm has to be initialized, which can be done using the Higher Order Singular Value Decomposition (HOSVD), as proposed in [START_REF]On the Best Rank-1 and Rank-(R 1 ,R 2 , . . . ,R N ) Approximation of Higher-Order Tensors[END_REF]. The stopping criterion of the algorithm is also studied in this paper.

V. PROPOSED ALGORITHM

We describe a general algorithm for the optimization of the proposed family of contrast functions.

A. Optimization method

We assume that the mixing filter is a FIR filter with impulse response of length L and that M(z) is irreducible. Then it admits a MIMO-FIR left inverse filter of length D (see [START_REF] Gorokhov | Subspace-Based Techniques for Blind Separation of Convolutive Mixtures with Temporally Correlated Sources[END_REF]), which is assumed causal with no loss of generality. The row vectors which define the impulse response can be stacked in the following (1, N D) row vector:

w (w(0) . . . w(D -1)). (27)
We also define the (N D, 1) column vector

x(n) (x(n) T x(n -1) T . . . x(n -D + 1) T ) T . (28)
It is then easily seen that

y(n) = w x(n). (29)
Let the covariance matrix be R = E{x(n)x(n) H }, so we have E{|y(n)| 2 } = wRw H . Now using the assumption A3., the multilinearity property of cumulants and (29), we have

κ S,R,r {y(n)} = i1,...,iS w ( * )1 i1 . . . w ( * )S iS Cum{x ( * )1 i1 (n), . . . , x ( * )S iS (n), r 1 (n), . . . , r R-S (n)}, (30)
where w

( * ) k i k , k ∈ {1, .
. . , S}, stands for w i k conjugated in the same way as y ( * ) k in [START_REF] Loubaton | Blind Deconvolution of Multivariate Signals: A Deflation Approach[END_REF]. Thus this relation can be written as a S-order tensor decomposition ( [START_REF]On the Best Rank-1 and Rank-(R 1 ,R 2 , . . . ,R N ) Approximation of Higher-Order Tensors[END_REF])

κ S,R,r {y(n)} = C S • 1 w ( * )1 • 2 . . . • S w ( * )S ( 31
)
where the tensor C S is defined component wise as The covariance matrix can have some of its eigenvalues equal to zero. It is important to notice that for any row vector such that w H 0 ∈ ker R we have w 0 Rw H 0 = 0. The solution of the above maximization problem (33) is not unique, it follows that we may impose in addition w H ∈ (ker R) ⊥ . Thus, w is projected onto the signal subspace. So, a SVD, we decompose

R = UDU H , (34)
and we define

P = UD 1/2 and Q = D -1/2 U H (35)
in order to project onto the signal subspace:

w = wP and CS = C S • 1 Q ( * )1 • 2 . . . • S Q ( * )S . (36)
We obtain CS and w and finally the problem in (33) reduces to the following one:

max| CS • 1 w( * )1 • 2 . . . • S w( * )S | with w wH = 1. ( 37
)
This projection onto the signal subspace is equivalent to a prewhitening step on the observations.

In Section IV, we have seen that the above maximization is equivalent to the minimization of the cost function for the best rank-1 approximation of tensor CS [START_REF]On the Best Rank-1 and Rank-(R 1 ,R 2 , . . . ,R N ) Approximation of Higher-Order Tensors[END_REF]. So we compute here the best rank-1 approximation of tensor CS , using the HOPM algorithm (it can be seen as a higher order extension of the power method for best rank-1 approximation of matrix), to maximize the proposed contrast.

B. A "fixed-point" like method

It has been proposed in a previous work [START_REF] Castella | Quadratic Higher-Order Criteria for Iterative Blind Separation of a MIMO Convolutive Mixture of Sources[END_REF] an iterative method to improve the performances of the extraction algorithm for the case S = 2. This "fixed-point" like method can be also generalized to the considered family of contrasts. Here, we recall briefly the method.

We assume, in this paragraph, that all references are the same. That is for all i ∈ {1, . . . , R-S}, we have r i (n) = r(n) (or possibly r i (n) = r * (n) in the complex case). The output obtained after maximization of C S,R,r should be closer to the extracted source than the reference signal. Basically, the output which has been previously obtained by the maximization of C S,R,r serves as a new reference signal. The fixed point principle is illustrated for three iterations in the figure 1 and explained below: "fixed-point" like algorithm (for l max iterations) Input: x,w 0 . Output: w lmax .

1) for l ∈ {1, . . . , l max }:

• fix the reference to r(n) = w l-1 x(n),

• w l = arg max w C S,R,r {y(n)} where y(n) = w x(n).

2) The separating filter is given by the coefficients w lmax .

VI. SIMULATIONS

We now propose computer simulations in order to illustrate the usefulness of our criterion w.r.t. other approaches.

We focus on the fourth order cumulant based contrasts in order to compare the generalized approach to the contrast functions defined in ( 14) and [START_REF] Kawamoto | Eigenvector Algorithms Incorporated with Reference Systems for Solving Blind Deconvolution of MIMO-IIR Linear Systems[END_REF]. In this context, we define the following contrast function for R = 4 and S = 3:

C 3,4,r {y(n)} |Cum{y(n), y * (n), y * (n), r(n)}|. ( 38
)
We compare three different algorithms:

1) The algorithm E 2 maximizes the quadratic criterion defined in [START_REF] Kawamoto | Eigenvector Algorithms Incorporated with Reference Systems for Solving Blind Deconvolution of MIMO-IIR Linear Systems[END_REF]. The solution is given by the eigenvector of C2 associated to the eigenvalue with largest modulus.

In practice it is obtained by the SVD of the considered matrix (second order tensor).

2) The algorithm E 3 maximizes the cubic criterion defined in (38). The solution is given by the best rank-one approximation of the tensor C3 . In practice it is obtained by a HOPM applied to the considered tensor (third order tensor).

3) The algorithm E 4 maximizes the classical quartic criterion with no reference defined in [START_REF] Jelonnek | Generalized eigenvector algorithm for blind equalization[END_REF]. The solution is given by the best rank-one approximation of the tensor C4 . In practice it is obtained by a HOPM applied to the considered tensor (fourth order tensor). The quality of extraction is measured thanks to an error index suggested in [START_REF] Ghennioui | A Nonunitary Joint Block Diagonalization Algorithm for Blind Separation of Convolutive Mixtures of Sources[END_REF] and defined by:

ind(g) 1 - max i∈{1,...,N } k∈Z |g i (k)| 2 i∈{1,...,N } k∈Z |g i (k)| 2 .
(39)

(n) = βs 1 (n) + (1 -β)s 2 (n)).
This error index is non negative and takes its values between 0 and 1. The zero value represents a perfect extraction. In all simulations the mean of the error index over 200 Monte-Carlo realizations is considered. At each run, the mixing system and the sources have been drawn randomly. The reference signal is chosen as the first observation from the mixing system (except in section VI-A).

In the following simulations, we have choosen to fix in the "fixed-point" method the number of iterations to two for i.i.d. sources and five for non i.i.d. sources to compare the referenced contrasts to the classical one (except in section VI-A).

A. Experiment 1 -Influence of the reference signal

We consider the case of K = 2 i.i.d. real-valued binary source signals (they take their values in {-1, 1} with equal probabilities), N = 5 observation signals and a real-valued mixing filter of length L = 3. Its coefficients are drawn randomly according to a zero-mean and unit variance normal distribution. To evaluate the influence of the reference signal, we choose it as r(n

) = βs 1 (n) + (1 -β)s 2 (n) with β ∈ [0, 1].
In Figure 2, for N e = 5000, the mean error index is plotted versus β for the E 2 and the E 3 algorithms (there is no reference in the E 4 algorithm). Better results are obtained for values of β near 0 or 1 that is to say when the reference signal is close to one source. The performance of the E 2 algorithm fails for β = 0.5 whereas the E 3 algorithm still works. This shows that our approach imposes nearly no constraint on the reference signal. E 3 is less sensitive to the chosen reference signal as E 2 .

B. Experiment 2 -fixed point method

The mixing system is complex-valued and its coefficients are drawn randomly according to a zero-mean and unit variance normal distribution. We consider here a mixture of K = 3 i.i.d. source signals. The length of the mixing filter is L = 3 and the number of observations is N = 6. In figure 3 the error index is drawn versus the number Fig. 3. Fixed point method for E 3 algorithm, complex i.i.d. case: Performance versus number of samples without noise. of sources samples N e of i.i.d. complex-valued Quadrature Amplitude Modulation sources with 16 states (16-QAM). In this figure, we give the performance of the E 3 algorithm for different number of iterations in the "fixed-point" method. Performing more than two "fixed-point" iterations does not improve the performance sharply. Therefore, in the following experiments for i.i.d. sources, we will fix, for the E 2 and E 3 algorithms, the number of iterations of the "fixed-point" method to two. Now we consider non i.i.d. complex-valued Continuous Phase Modulated sources. Pseudo-symbols of CPM signals are unit modulus signals defined by the relation

s k (n + 1) = exp(ıπh k a k n )
where h k ∈]0, 1[ is a modulation index and a k n is the i.i.d. binary symbol sequence transmitted by th k-th user. We fix the modulation index h to 0.2. In figure 4 the error index is drawn versus the number of sources samples N e . In this figure, we give the performance of the E 3 algorithm for different number of iterations in the "fixed-point" method. After five iterations, the performance does not improve sharply, then in the following experiments for non i.i.d. sources, we will fix, for the E 2 and E 3 algorithms, the number of iterations of the "fixed-point" method to five.

C. Experiment 3 -i.i.d. real-valued source signals

We consider real-valued binary sources. The mixing system is also real-valued and its coefficients are drawn randomly according to a zero-mean and unit variance normal distribution. We consider here a mixture of K = 3 i.i.d. source signals.

The length of the mixing filter is L = 3 and the number of observations is N = 6. In figure 5 the error index is drawn versus the number of sources samples N e . In this figure, we compare the performance of the E 2 , E 3 and E 4 algorithms. Our approach performs better than E 4 for small amount of samples and similarly for large amount of samples. Here, even if E 4 is much more complex than E 3 (fourth order tensor versus third order tensor), the latter performs slighty better. The extraction performance of the E 2 algorithm is slightly Fig. 4.

Fixed point method for E 3 algorithm, complex non i.i.d. case: Performance versus number of samples without noise. lower in this configuration but it is important to note that its performances can be improved by increasing the number of iteration in the fixed point method.

D. Experiment 4 -i.i.d. complex-valued source signals

We consider complex-valued 16-QAM sources. The mixing system is complex-valued and its coefficients are drawn randomly according to a zero-mean and unit variance normal distribution. We consider here a mixture of K = 3 i.i.d. source signals. The length of the mixing filter is L = 3 and the number of observations is N = 6. In figure 6 the error index is drawn versus the number of sources samples N e . In this figure, we compare the performance of the E 2 , E 3 and E 4 algorithms. Our approach performs again better than E 4 for small amount of samples. E 2 is performing better for number of samples below 2000. With larger amount of samples, the E 4 and the E 3 algorithms are performing very similarly. The extraction performance of the E 2 algorithm is again slightly lower in this configuration.

In order to compare different optimization schemes, we put in the table I the performances of the cubic criterion and the contrast defined in ( 14) optimized using a gradient ascent method with an adaptatively adjusted step-size (initially set to one and divided by two each time the contrast does not increase), called Kurtosis. The Kurtosis and E 3 algorithms perform very similarly. The results in table II show up that 2 fixed point optimizations of the cubic criterion can be performed about 6 times faster than a single optimization of the classical kurtosis contrast (the execution times have been estimated on a dual processor Intel Xeon Dual Core 5150 running at 2.667GHz clock frequency and with 16GB RAM). Even if we have higher value of the mean error index than in the previous experiment for the three algorithms, the extraction is well performed. We illustrate this fact in the figure 7, with the constellation of a extracted 16-QAM source of 5000 samples. The algorithm used for this extraction is E 3 with N e = 5000 and the error index is equal to 0.002.

E. Experiment 5 -non i.i.d. complex-valued source signals

Now we consider non i.i.d. complex-valued CPM sources. We still consider a mixture of K = 3 source signals, a mixing filter length of L = 3 and N = 6 observations. We fix the modulation index h to 0.2. In figure 8 the error index is drawn versus the number of sources samples N e and with no noise. In figure 8, we compare the performance of the E 2 , E 3 and E 4 algorithms. For non i.i.d. complex-valued sources, E 4 and E 3 perform equally well and better than E 2 .

F. Experiment 6 -additive noise signal

We consider here an additive noise signal in the model (1) that is:

x

(n) = k∈Z M(k)s(n -k) + n(n), (40)
where n(n) is the (N, 1) noise vector. The noise signals n i (n), i ∈ {1, . . . , N } are additive, complex-valued and follow the same normal law with zero-mean and N 0 variance. We define the SNR in the worst case i.e. through the power ratio of the contributions of all the sources to the additive noise.

SNR dB = ( M 2 ) dB -(N 0 ) dB , (41)
where M(n) is the convolutive mixing system matrix.

M 2 = k∈Z N i=1 K j=1 |M i,j | 2 (k), (42)
We consider here a mixture of K = 3 source signals in two situations: with i.i.d. and non i.i.d. source signals.

The length of the mixing filter is L = 3 and the number of observations is N = 6. In the figures 9 and 10 the error index is drawn versus the SNR with N e = 5000 and we compare the performance of the E 2 , E 3 and E 4 algorithms.

1) i.i.d. case:

In the figure 9, we consider a mixture of 16-QAM source signals. Our approach performs better than E 4 for low SNR values, under 10 dB. After that point, we get back the same results as in the figure 6 for N e = 5000: E 3 and E 4 are performing similarly and better than E 2 .

2) non i.i.d. case: In the figure 10, we consider a mixture of CPM source signals. We fix the modulation index h to 0.2. In this case all the algorithms are performing well for SNR values above 15 dB. After 50 dB, our algorithm and E 4 are performing similarly better than the E 2 . The E 3 algorithm is less computationally intensive than E 4 (fourth order tensor versus third order tensor), and the extraction thanks to E 3 is done in allmost all the situation presented here.

VII. CONCLUSION

We have proposed generalized contrast functions within the framework of blind extraction of one source from a MIMO convolutive mixture. This family of contrasts holds for both i.i.d. and non i.i.d. sources. In order to optimize this criterion, we have developped a new generalized algebraic algorithm based on best rank-1 tensor decomposition.

Unlike other referenced approaches, our new contrasts require no constraint on the reference signals. The new approach realizes a good compromise between generality and implementation simplicity. Computer simulations illustrate interesting features and good performance in comparison with quadratic and quartic criteria.

(47) corresponds to the equalization condition.

APPENDIX B PROOF OF PROPOSITION III.4

Proof: We want to prove that at most one component of the norm of the global filter is non-zero. For any i ∈ {1, . . . , K}, we write the i-th component of the unit norm vector filter g: g i = g i i gi where gi reads g i / g i i if g i i = 0 and gi is zero if g i i = 0. Defining ỹi (n) (g i ⋆ s i ) (n), from the assumption A3., the cumulant multilinearity and the mutual independance of the signals ỹi (n), i ∈ {1, . . . , K}, yield:

κ S,R,r {y(n)} = K i=1 g i S i κ S,R,r {ỹ i (n)}. (48)
Then we have the following inequality:

C S,R,r {(g ⋆ s)(n)} K i=1 g i S i
κ S,R,r {ỹ i (n)} .

Using [START_REF] Simon | Separation of a Class of Convolutive Mixtures: A Contrast Function Approach[END_REF], gi i = 1 and g = 1, we get:

C S,R,r {(g ⋆ s)(n)} K i=1 g i S i M S,R,i M max S,R K i=1 g i S i M max S,R K i=1 g i 2 i M max S,R (50)
The last inequality holds because, as S > 2,

K i=1 g i S i K i=1 g i 2 i = 1 .
The previous inequality proves that the property p1. of the definition (III.2) is satisfied. Moreover, if the equality is satisfied, we can write

K i=1 g i S i = K i=1 g i 2 i = 1. (51)
It comes that for all i ∈ {1, . . . , K}, at most one g i i is non zero. Additionally this value of i is such that M S,R,i = M max S,R . Thus, the property p2. of (III.2) is satisfied. The converse is easily deduced, and the proposition (III.2) is satisfied.

(

  C S ) i1,...,iS = Cum{x ( * )1 i1 (n), . . . , x ( * )S iS (n), r 1 (n), . . . , r R-S (n)} . (32) Hence the optimization of the contrast function in (17) under the unit power constraint reads max |C S • 1 w ( * )1 • 2 . . .• S w ( * )S | with wRw H = 1 . (33)

Fig. 1 .

 1 Fig. 1. Fixed point like method for 3 iterations.

Fig. 2 .

 2 Fig. 2. Performance versus coefficient β for two sources (where the reference signal is r(n) = βs 1 (n) + (1 -β)s 2 (n)).

Fig. 5 .

 5 Fig. 5. Real-valued i.i.d. case: Performance versus number of samples without noise.

Fig. 6 .

 6 Fig. 6. Complex-valued i.i.d. case: Performance versus number of samples without noise.

Fig. 7 .

 7 Fig. 7. Complex-valued i.i.d. case: Extracted QAM-16 constellation, error index is equal to 0.002.

Fig. 9 .

 9 Fig. 9. Complex-valued i.i.d. case: Performance versus SNR for Ne = 5000.

Fig. 10 .

 10 Fig. 10. Complex-valued non i.i.d. case: Performance versus SNR for Ne = 5000.

1 . 2 .

 12 APPENDIX A PROOF OF PROPOSITION III.3Proof: Using the assumption A3., the multilinearity property of cumulants and the i.i.d. assumption of source signals, we haveκ S,R,r {y(n)} = κ S,R,r {s j (n -k)} . (43) Hence C S,R,r {y(n)} = κ S,R,r {s j (n -k)} , (44) and thus C S,R,r {y(n)} j,k |g j (k)| S |κ S,R,r {s j (n -k)}| . (45)As y(n) is unit power, we havej,k |g j (k)| 2 =Hence |g j (k)| 2 1 and as S > 2 we have for all j and k|g j (k)| S |g j (k)|Using this result and the assumption A4., we haveC S,R,r {y(n)} j,k |g j (k)| 2 |κ S,R,r {s j (n -k)}| N max j=1 sup k∈Z |κ S,R,r {s j (n -k)}| R,r {s j (n -k)}| .(46)Finally considering the equality, we have |g j (k)| S = |g j (k)| 2 only if it exists a given couple (j 0 , k 0 ) such that g j0 (k 0 ) = 1 g j (k) = 0 ∀(j, k) = (j 0 , k 0 ) .

TABLE I COMPARISON

 I OF THE AVERAGE ERROR INDEX FOR 1) ONE OPTIMIZATION OF THE KURTOSIS CONTRAST 2) 2 FIXED POINT ITERATIONS OPTIMIZATION

				OF CUBIC CRITERION
		Number of samples	5000	6000	7000	8000	9000	10000
		Kurtosis	12.38	17.38	20.55	24.41	28.53	32.71
		Cubic Criterion	2.18	2.55	2.99	3.52	4.13	4.80
				TABLE II
	COMPARISON OF THE AVERAGE EXECUTION TIME (IN S) FOR 1) ONE OPTIMIZATION OF THE KURTOSIS CONTRAST 2) 2 FIXED POINT ITERATIONS
				OPTIMIZATION OF CUBIC CRITERION
	Fig. 8.	Complex-valued non i.i.d. case: Performance versus number of
	samples without noise.

the field of signal processing from the University of Toulon, France. From 2000 to 2003, he was in postdoctoral position with the University of Toulon (2 years) and with STMicroelectronics (1 year) in the field of wireless telecommunications. Since 2003, he has been an assistant professor within ISITV's Telecommunications Department. His main research interests are in statistical signal processing using higher-order statistics.

Marc Castella was born in 1976 in Courbevoie,

France. In 2000, he received the "Agrégation" degree in the field of applied physics and in 2001 he received the M.Sc. degree in electrical engineering, from both the " École Normale Supérieure de Cachan" and the "Université Paris-Sud, Orsay" (France). He obtained his Ph.D. degree from "Université de Marne-la-Vallée" (France) in 2004. Since then, he has been "Maître de Conférence" in Evry (France) at "Télécom & Management SudParis" (formerly "Institut National des Télécommunications (INT)"). He is also a member of the UMR-CNRS 5157 research team SAMOVAR. His research activities have focused on the problem blind source separation in general.

Eric Moreau (M'96-SM'08) was born in Lille, France. He graduated from the "Ecole Nationale Supérieure des Arts et Métiers" (ENSAM), Paris, France, in 1989 and received the "Agrégation de Physique" degree from the "Ecole Normale Supérieure de Cachan" in 1990. He received the DEA degree in 1991 and the Ph.D. degree in 1995, both in the field of signal processing and from the University of Paris-Sud, France. From 1995 to 2001, he was assistant professor within the Telecommunications Department of the Engineering School "Institut des Sciences de l'Ingénieur de Toulon et du Var" (ISITV), La Valette, France. He is currently a Professor with the University of Toulon, France. His main research interests are in statistical signal processing using high-order statistics.