W Souidène
email: souidene@tsi.enst.fr

A Aïssa-El-Bey

K Abed-Meraim

A Beghdadi
email: beghdadi@l2ti.univ-paris13.fr

BLIND IMAGE SEPARATION USING SPARSE REPRESENTATION

Keywords: Separation, image restoration, sparse matrices

This paper focuses on the blind image separation using their sparse representation in an appropriate transform domain. A new separation method is proposed that proceeds in two steps: (i) an image pretreatment step to transform the original sources into sparse images and to reduce the mixture matrix to an orthogonal transform (ii) and a separation step that exploits the transformed image sparsity via an p-norm based contrast function. A simple and efficient natural gradient technique is used for the optimization of the contrast function. The resulting algorithm is shown to outperform existing techniques in terms of separation quality and computational cost.

INTRODUCTION

Blind source separation (BSS) is an important research field in signal and image processing. In particular, separating linear mixtures of several 'independent' images has application in biomedical imaging [1-3], in cosmology and multispectral imaging [START_REF] Naceur | The contribution of the sources separation method in the decomposition of mixed pixels[END_REF][START_REF] Bijaoui | Blind source separation of multispectral astronomical images[END_REF], in polarimetric imaging [START_REF] Bronstein | Separation of reflections via sparse ICA[END_REF], etc. Recently, an important research activity has been observed for solving the BSS problem using a sparse representation of the source signals. Solution for the blind separation of image sources using sparsity include the wavelet-transform domain method in [START_REF] Bronstein | Separation of reflections via sparse ICA[END_REF] and the method in [START_REF] Zibulevsky | Blind source separation by sparse decomposition in signal dictionary[END_REF] using projection onto sparse dictionaries. In this work, we propose a new solution based on the transformed image sparsity. The new BSS algorithm is shown to be more efficient that other existing techniques in the literature and leads to improved separation quality with lower computational cost.

NOTATIONS AND DATA MODEL

We assume that N images f1, • • • , f N each of size (m f , n f ) are merged and M linear mixtures of these original images are observed. The latter mixtures can be modeled by the following linear system:

g(m, n) = Af (m, n) + w(m, n) (1) where, f (m, n) = [f 1 (m, n), • • • , f N (m, n)] T is a N × 1 im- age source vector consisting of the stack of corresponding pixels of source images, w(m, n) = [w1(m, n), • • • , wM (m, n)] T is the M × 1 gaussian complex noise vector which affects each image mixture pixel, A is the M × N full column rank mixing matrix (i.e., M ≥ N ), g(m, n) = [g 1 (m, n), • • • , g M (m, n)]
T is an M × 1 vector of mixture image pixels and the superscript T denotes the transpose operator. The purpose of blind image separation is to find a separating matrix, i.e. a N × M matrix B such that

b f (m, n) = Bg(m, n
) is an estimate of original images. In practice, the separating matrix estimation is performed up to a permutation and a certain fixed scalar, i.e. B is a separating matrix iff:

BA = P Λ ( 2
)
where P is a permutation matrix and Λ a non-singular diagonal matrix that represent the inherent ambiguities of the BSS problem.

SEPARATION ALGORITHM

As shown in [START_REF] Bronstein | Separation of reflections via sparse ICA[END_REF][START_REF] Zibulevsky | Blind source separation by sparse decomposition in signal dictionary[END_REF], exploiting the sparsity of some representations of the original images afford us to achieve the BSS problem. Indeed, the mixture destroys or 'reduces' the sparsity of the considered signals that is restored after source separation. Reversely, it is shown in [START_REF] Bronstein | Separation of reflections via sparse ICA[END_REF][START_REF] Zibulevsky | Blind source separation by sparse decomposition in signal dictionary[END_REF] that restoring (maximizing) the sparsity leads to the desired source separation. Based on this, we propose in the sequel a two-step BSS solution consisting in a linear pre-treatment that transforms the original sources into sparse signals followed by a BSS algorithm that minimizes the p norm of the transformed image mixtures using natural gradient technique.

Image pre-treatment

The algorithm proposed in this article is efficient for separating sparse sources. For some signals, one can assume that the spatial or temporal representation is naturally sparse, whereas for natural scenes, this assumptions falls down. We propose to make the image sparse by simply taking into account its Laplacian transform:

F = ∇f = ∂ 2 f ∂x 2 + ∂ 2 f ∂y 2 , ( 3
)
or, in discrete form

F (m, n) = f (m + 1, n) + f (m -1, n) + f (m, n + 1) +f (m, n -1) -4f (m, n) .
Our motivation for choosing this transformation is two fold. First the Laplacian transform is a sparse representation of the image since it acts as an edge detector which provides a two-level image, the edges and the homogeneous background. Second, the Laplacian is a linear transformation. This latter property is 'critical' since the separation matrix estimated to separate the image mixtures is the same to separate the mixture of Laplacian images:

G = ∂ 2 Af ∂x 2 + ∂ 2 Af ∂y 2 = AF ( 4
)
where G is the Laplacian transform of the mixtures. In the literature, some other linear transformations were proposed in order to make the image sparse, including the projection into a sparse dictionary [START_REF] Zibulevsky | Blind source separation by sparse decomposition in signal dictionary[END_REF]. In Figure 1, the original cameraman image is displayed as well as its Laplacian transform and their respective histograms that clearly show the sparsity of the latter. In the pre-treatment phase, we also propose an optional whitening step which aims to set the mixtures to the same energy level. Furthermore, this procedures reduces the number of parameters to be estimated. More precisely, the whitening step is applied to the Laplacian image mixtures before using our separation algorithm. The whitening is achieved by applying a N × M matrix Q to the Laplacian image mixtures in such a way Cov(QG) = I in the noiseless case, where Cov(•) stands for the covariance operator. As shown in [START_REF] Belouchrani | A blind source separation technique using second-order statistics[END_REF], Q can be computed as the inverse square root of the noiseless covariance matrix of the Laplacian image mixtures (see [START_REF] Belouchrani | A blind source separation technique using second-order statistics[END_REF] for more details). In the following, we apply our separation algorithm on the whitened data:

G w (m, n) = QG(m, n).

Sparsity-based BSS algorithm

In this section, we propose an iterative algorithm for the separation of sparse signals, namely the ISBS for Iterative Sparse Blind Separation algorithm. It is well known that Laplacian image transform is characterized by its sparsity property in the spatial domain [1,[START_REF] Zibulevsky | Sparse source separation with relative Newton method[END_REF]. This property can be measured by the p norm where 0 ≤ p < 2. More specifically, one can define the following sparsity based contrast function,

G p (F ) = N X i=1 [J p (F i )] 1 p ( 5
)
where

Jp(Fi) = 1 m f n f m f X n=1 n f X m=1 |Fi(m, n)| p . ( 6
)
The algorithm finds a separating matrix B such as,

B = arg min B {Gp(B)} (7)
where Gp(B) Gp(H) [START_REF] Belouchrani | A blind source separation technique using second-order statistics[END_REF] and H(m, n) BG w (m, n) represents the estimated image sources Laplacian. The approach we choose to solve ( 7) is inspired from [START_REF] Pham | Blind separation of mixture of independent sources through a quasi-maximum likelihood approach[END_REF]. It is a block technique based on the processing of m f n f observed image pixels and consists in searching the minimum of the sample version of [START_REF] Zibulevsky | Blind source separation by sparse decomposition in signal dictionary[END_REF]. Solutions are obtained iteratively in the form:

B (k+1) = (I + (k) )B (k) ( 9
)
H (k+1) (m, n) = (I + (k) )H (k) (m, n) . ( 10
)
At iteration k, a matrix (k) is determined from a local linearization of G p (BG w ). It is an approximate Newton technique with the benefit that (k) can be very simply computed (no Hessian inversion) under the additional assumption that B (k) is close to a separating matrix. This procedure is illustrated in the following. At the (k+1) th iteration, the proposed criterion ( 6) can be developed as follows:

J p (H (k+1) i ) = = 1 m f n f m f X m=1 n f X n=1 ˛H(k) i (m, n) + N X j=1 (k) ij H (k) j (m, n) ˛p .
Under the assumption that B (k) is close to a separating matrix, we have

| (k) ij | 1
and thus, a first order approximation of J p (H (k+1) i

) is given by:

Jp(H (k+1) i ) ≈ 1 m f n f m f P m=1 n f P n=1 n |H (k) i (m, n)| p + p N P j=1 (k) ij " |H (k) i (m, n)| p-1 sgn " H (k) i (m, n) " H (k) j (m, n) "o (11)
where sgn(•) represents the sign value operator. Using equation [START_REF] Beghdadi | A new image distortion measure based on wavelet decomposition[END_REF], equation ( 5) can be rewritten in more compact form as:

G p " B (k+1) " = G p " B (k) " + T r " (k) R (k)T D (k) " ( 12
)
where T r(•) is the matrix trace operator, the ij th entry of matrix R (k) is given by:

R (k) ij = 1 m f n f m f X m=1 n f X n=1 |H (k) i (m, n)| p-1 sgn " H (k) i (m, n) " H (k) j (m, n) .
and

D (k) = h diag " R (k) 11 , . . . , R (k) N N "i 1 p -1 . ( 13
)
Using a gradient technique, (k) can be written as:

(k) = -µD (k) R (k) , ( 14
)
where µ > 0 is the descent step. Replacing ( 14) into (12) leads to,

Gp " B (k+1) " = Gp " B (k) " -µ D (k) R (k) 2 , ( 15
)
so µ controls the decrement of the criterion. Now, to avoid the algorithm's convergence to the trivial solution B = 0, one normalizes the outputs of the separating matrix to unit-power, i.e. ρ (k+1)

H i E " |H (k+1) i (m, n)| 2 " = 1, ∀ i.
Using first order approximation, this normalization leads to:

(k) ii = 1 -ρ (k) H i 2ρ (k) H i . ( 16
)
The final estimated separation matrix B = B (K) Q is applied to the image mixtures g to obtain an estimation of the original images. K denotes here the number of iterations that can be either chosen a priori or given by a stopping criterion of the form B (k+1) -B (k) < δ where δ is a small threshold value.

PERFORMANCE EVALUATION

All simulations are carried on 256 × 256 parrot and cameraman images. The number of observed mixtures is M = 2. The algorithms are developed on MATLAB environment. Monte Carlo simulations are carried over 200 random realizations of the additive gaussian noise for variable Signal to Noise Ratios (SNR). In order to objectively evaluate the performance of the proposed algorithm, we consider two different criteria, the first one is the Interference to Signal Ratio (ISR) criterion [START_REF] Belouchrani | A blind source separation technique using second-order statistics[END_REF] defined as:

ISR N X i=1 N X j =i E `|(BA) ij | 2 ´ρj E (|(BA) ii | 2 ) ρ i ( 17
)
where

ρ i = E(|f i (m, n)| 2
) is the i th source power. The second one is an objective image quality measure inspired from the Human Visual System (HVS) properties and developed in [START_REF] Beghdadi | A new image distortion measure based on wavelet decomposition[END_REF]. It is called P SN R -W AV for Peak Signal to Noise Ratio based on Wavelet decomposition. The separation result of the proposed algorithm is depicted on Fig. 2 where we represent the two original images (f 1 , f 2 ), the mixtures (g 1 , g 2 ) and the recovered ones

( b f 1 , b f 2
) by the proposed algorithm in the noiseless case. In Fig. 3, we compare the performance of the proposed algorithm to the Relative Newton algorithm developed by Zibulevsky et al. in [START_REF] Zibulevsky | Sparse source separation with relative Newton method[END_REF] where the case of sparse sources is considered. We plot the residual interference between separated images (ISR) versus the SNR. It is clearly shown that our algorithm (ISBS) performs better in terms of ISR especially for low SNRs. We plot on Fig. 4 the objective distortion measure between the original and separated images versus the SNR. One can observe that, for each image, we reach the same conclusion for the P SN R -W AV as for the ISR. In Fig. 5, we represent the evolution of the ISR as a function of the iteration number. A fast convergence rate is observed. Moreover, the complexity of the proposed algorithm is equal to 2N 2 m f n f + O(N 2 ) flops per iteration whereas the complexity of the Relative Newton algorithm in [START_REF] Zibulevsky | Sparse source separation with relative Newton method[END_REF] is

2N 4 + N 3 m f n f + N 6 /6.

CONCLUSION

This article deals with a simple and efficient two-step algorithm of blind image separation. The proposed method consists in a sparsification of the natural observed mixtures followed by a blind separation of the original images. The sparsification is simply the Laplacian transform and has a low computational cost. The separation is performed using an iterative algorithm based on the minimizing of

[1] A. Cichocki and S. Amari, Adaptive Blind Signal and Image Processing, Wiley & Sons, Ltd., UK, 2003.

[2] N. Zhang, J. LU, and T. Yahgi, "Nonlinear blind source separa-

Fig. 1 .

 1 Fig. 1. (a) Original image, (b) Laplacian transform, (c) Original image histogram, (d) Sparse Laplacian transform histogram

Fig. 2 .

 2 Fig. 2. (a)-(b) original images, (c)-(d) M = 2 observed mixtures, (e)-(f) restored images using ISBS algorithm.

Fig. 3 .

 3 Fig. 3. Interference to Signal Ratio (ISR) versus SNR

Fig. 4 .Fig. 5 .

 45 Fig. 4. Objective image quality measure (P SN R -W AV )


Abdeldjalil Aïssa-El-Bey

Karim Abed-Meraim

Yves Grenier
email: grenier@tsi.enst.fr

UNDERDETERMINED BLIND SEPARATION OF AUDIO SOURCES FROM THE TIME-FREQUENCY REPRESENTATION OF THEIR CONVOLUTIVE MIXTURES

Keywords: Separation, deconvolution, time-frequency analysis, identification

This paper considers the blind separation of nonstationary sources in the underdetermined convolutive mixture case. We introduce two methods based on the sparsity assumption of the sources in the timefrequency (TF) domain. The first one assumes that the sources are disjoint in the TF domain; i.e. there is at most one source signal present at a given point in the TF domain. In the second method, we relax this assumption by allowing the sources to be TF-nondisjoint to a certain extent. In particular, the number of sources present (active) at a TF point should be strictly less than the number of sensors. In that case, the separation can be achieved thanks to subspace projection which allows us to identify the active sources and to estimate their corresponding time-frequency distribution (TFD) values.

INTRODUCTION

The blind source separation of more sources than sensors (referred to as UBSS for underdetermined blind source separation) is still a challenging problem especially in the convolutive mixtures case. In the instantaneous mixture case, some methods exploiting the sparseness of the sources in certain transform domain have been proposed for UBSS [START_REF] Bofill | Underdetermined blind source separation using sparse representations[END_REF][START_REF] Linh-Trung | Separating more sources than sensors using time-frequency distributions[END_REF][START_REF] Yilmaz | Blind separation of speech mixtures via time-frequency masking[END_REF][START_REF] Linh-Trung | Underdetermined blind source separation of non-disjoint nonstationary sources in time-frequency domain[END_REF]. These methods proceed 'roughly' as follows: The mixtures are first transformed to an appropriate representation domain; the transformed sources are then estimated using their sparseness, and finally one recovers their time waveforms by source synthesis (for more information, see the recent survey work [START_REF] O'grady | Survey of sparse and nonsparse methods in source separation[END_REF]). UBSS methods for nonstationary sources have been proposed, given that these sources are sparse in the time-frequency (TF) domain [START_REF] Linh-Trung | Separating more sources than sensors using time-frequency distributions[END_REF][START_REF] Yilmaz | Blind separation of speech mixtures via time-frequency masking[END_REF]. The main assumption used in these methods is that the sources are TF-disjoint. In other words, there is at most one source present at any point in the TF domain. This assumption is rather restrictive, though the methods have also showed that they worked well under a quasi sparseness condition, i.e. sources are TF-almost-disjoint. In this paper we focus on the UBSS in convolutive mixture case and target the relaxation of the TF-disjointness condition by allowing the sources to be nondisjoint in the TF domain; that is, multiple sources are possibly present at any point in the TF domain. This case has been considered in [START_REF] Linh-Trung | Underdetermined blind source separation of non-disjoint nonstationary sources in time-frequency domain[END_REF] for the separation of instantaneous mixtures, in [START_REF] Rosca | Generalized sparse signal mixing model and application to noisy blind source separation[END_REF] for the deconvolution of single-path channels with non-zero delays and in [START_REF] Araki | Blind separation of more speech than sensors with less distortion by combining sparseness and ICA[END_REF] where binary TF-masking and ICA technique are jointly used. The main contribution of this paper consists in two new algorithms (TF-CUBSS for Time-frequency convolutive underdetermined blind source separation) for solving the UBSS in the TF domain; the first one uses vector clustering while the other uses subspace projection.

PROBLEM FORMULATION

Data model

Let s 1 (t), . . . , s N (t) be the desired sources to be recovered from the convolutive mixtures x 1 (t), . . . , x M (t) given by:

x(t) = K k=0 H(k)s(t -k) + w(t) (1)
where s(t) = [s1(t), . . . , sN (t)] T is the source vector, x(t) = [x1(t), . . . , xM (t)] T is the mixture vector (with M < N ), w(t) is the observation noise, and

H(k) are M × N matrices for k ∈ [0, K]
representing the impulse response coefficients of the channel that satisfies:

Assumption 1 The channel is such that each column vector of

H(z) def = K k=0 H(k)z -k def = [h1(z), . . . , hN (z)]
is irreducible, i.e. the entries of hi(z) denoted hji(z), j = 1, . . . , M , have no common zeros ∀i. Moreover, any M column vectors of H(z) form a polynomial matrix H(z) that it full rank over the unit-circle, i.e. rank( H(f )) = M ∀f .

TF conditions on the sources

In order to deal with UBSS, one often seeks for a sparse representation of the sources [START_REF] Bofill | Underdetermined blind source separation using sparse representations[END_REF]. In other words, if the sources can be sparsely represented in some domain, then their separation can be carried out in that domain by exploiting their sparseness.

TF-disjoint sources

Recently, there have been several UBSS methods, notably those in [START_REF] Linh-Trung | Separating more sources than sensors using time-frequency distributions[END_REF] and [START_REF] Yilmaz | Blind separation of speech mixtures via time-frequency masking[END_REF], in which the TF domain has been chosen to be the underlaying sparse domain. These two papers have based their solutions on the assumption that the sources are disjoint in the TF domain. Mathematically, if Ω 1 and Ω 2 are the TF supports of two sources s 1 (t) and s 2 (t) then the sources are said TF-disjoint if Ω 1 ∩ Ω 2 = ∅. However, this is a rather strict assumption. A more practical assumption is that the sources are almost-disjoint in the TF domain [START_REF] Linh-Trung | Separating more sources than sensors using time-frequency distributions[END_REF], allowing some small overlapping in the TF domain, for which the above two methods also worked.

TF-nondisjoint sources

In this paper, we want to relax the TF-disjoint condition by allowing the sources to be nondisjoint in the TF domain. Therefore, we will allow the sources to be nondisjoint in the TF domain; that is, multiple sources are allowed to be present at any point in the TF domain. However, instead of being inevitably nondisjoint, we limit ourselves by making the following constraint:

Assumption 2 The number of active sources (i.e. sources that overlap) at any TF point is strictly less than the number of sensors.

In other words, for the configuration of M sensors, there exists at most (M -1) overlapping sources at any point in the TF domain.

For the special case when M = 2, Assumption 2 reduces to the disjoint condition.

TF-CUBSS ALGORITHM

In order to solve the UBSS problem in the convolutive case, we propose to identify first the impulse response of the channels. This problem in overdetermined case is very difficult and becomes almost impossible in the underdetermined case without side information on the considered sources. In this work and similarly to [START_REF] Huang | A blind channel identificationbased two-stage approach to separation and dereverberation of speech signals in a reverberant environment[END_REF], we exploit the sparseness property of the audio sources by assuming that from time to time only one source is present. In other words, we consider the following assumption:

Assumption 3 There exists, periodically, time intervals where only one source is present in the mixture. This occurs for all source signals of the considered mixtures.

To detect these time intervals, we propose to use information criteria based testing for the estimation of the number of sources present in the signal (see Section 3.1 for more details).

Channel estimation

Based on assumption 3, we propose here to apply SIMO (Single Input Multiple Output) based techniques to blindly estimate the channel impulse response. Regarding the problem at hand, we have to solve three different problems: first, we have to select time intervals where only one source signal is effectively present; then, for each selected time interval one should apply an appropriate blind SIMO identification technique to estimate the channel parameters; finally, the way we proceed, the same channel may be estimated several times and hence one has to group together (cluster) the channel estimates into N classes corresponding to the N source channels.

Source number estimation

Let define the spatio-temporal vector:

x d (t) = [x T (t), . . . , x T (t-d+1)] T = N k=1 H k s k (t)+w d (t), (2)
where H k are block-Sylvester matrices of size dM × (d + K),

s k (t) def = [s k (t), . . . , s k (t -K -d + 1)
] T and d is a chosen processing window size. Under the data model assumption and for large window sizes (see [START_REF] Wax | Detection of signals by information theoretic criteria[END_REF] for more details), matrices H k are full column rank. Hence, in the noiseless case, the rank of the data covariance matrix

R def = E[x d (t)x H d (t)] is equal to min(p(d + K), dM )
where p is the number of sources present in the considered time interval over which the covariance matrix is estimated. In particular, for p = 1, one has the minimum rank value equal to (d + K). Therefore, our approach consists in estimating the rank of the sample averaged covariance matrix R over several time slots (intervals) and select those corresponding to the smallest rank value r = d+K. The estimation of the rank value is done here by Akaike's criterion [START_REF] Wax | Detection of signals by information theoretic criteria[END_REF] according to:

r = arg min k       -2 log      M d i=k+1 λ 1/(M d-k) i 1 M d-k M d i=k+1 λ i      (M d-k)Ts + 2k(2M d -k)       , (3)
where λ 1 ≥ . . . ≥ λ M d represent the eigenvalues of R and T s is the time slot size. Note that it is not necessary at this stage, to know exactly the channel degree K as long as d > K (i.e. an over-estimation of the channel degree is sufficient) in which case the presence of one signal source is characterized by: d < r < 2d .

Blind channel identification

To perform the blind channel identification, we have used in this paper the Cross-Relation (CR) technique described in [START_REF] Xu | A least-squares approach to blind channel identification[END_REF]. This method is used on the time slots, where only one source signal is active. The latter are selected using the previously described Akaike's criterion. Note that there exist an improved, but more expensive, version of the CR method exploiting the quasi-sparse nature of acoustic impulse response [START_REF] Ahmad | Proportionate frequency domain adaptive algorithms for blind channel identification[END_REF] which can be used as well at this stage.

Clustering of channel vector estimates

The first step of our channel estimation method consists in detecting the time slots where only one single source signal is 'effectively' present. However, the same source signal s i may be present in several time intervals leading to several estimates of the same channel vector h i def = = [h 1i (0) . . . h M i (0) . . . h 1i (K) . . . h M i (K)] T . We end up, finally, with several estimates of each source channel that we need to group together into N classes. This is done by clustering the estimated vectors using k-means algorithm [START_REF] Frank | The data analysis handbook[END_REF]. The i th channel estimate is evaluated as the centroid of the i th class.

UBSS algorithm with TF-disjoint assumption

In this section, we propose a new cluster-based TF-CUBSS algorithm using the STFT (Short Time Fourier Transform) for convolutive mixture case. After transformation into the TF domain using the STFT, the model in (1) becomes (in the noiseless case):

Sx(t, f ) = H(f )Ss(t, f ), (4)
where Sx(t, f ) is the mixture STFT vector, Ss(t, f ) is the source STFT vector and

H(f ) = [h 1 (f ) . . . h N (f )]
is the channel Fourier transform matrix. Under the assumption that all sources are disjoint in the TF domain, (4) reduces to

S x (t, f ) = h i (f )S s i (t, f ), ∀(t, f ) ∈ Ω i , ∀i ∈ N , ( 5
)
where N = {1, . . . , N }. Consequently, two TF points (t 1 , f 1 ) and (t2, f2) belonging to the same region Ωi (i.e. corresponding to the source signal si) are 'associated' with the same channel hi. This latter observation is used next to cluster together the TF points of a given source signal. More precisely the algorithm proceeds as follows: First, we compute the STFT of the mixtures according to:

S x i (t, f ) = ∞ -∞ x i (τ )w(τ -t)e -j2πf τ dτ, i = 1, . . . , M, S x (t, f ) = [S x 1 (t, f ), . . . , S x M (t, f )] T .
Then, we apply a noise thresholding procedure which mitigates the noise effect but also reduces the computational cost as only the selected TF points are further treated by our algorithm. In particular, for each frequency-slice (t, fp) of the TFD representation, we apply the following criterion for all the time points t k belonging to this frequency-slice

If Sx(t k , fp) max t { S x (t, f p ) } > , then keep (t k , f p ), (7)
where is a small threshold (typically, = 0.01). Then, the set of all selected points, Ω, is expressed by Ω = N i=1 Ω i , where Ω i is the TF support of source s i . Note that, the effects of spreading the noise energy while localizing the source energy in the time-frequency domain amounts to increasing the robustness of the proposed method with respect to noise. Hence, by equation ( 7), we would keep only time-frequency points where the signal energy is non-negligible, the other time-frequency points are rejected, i.e. not further processed, since considered to represent noise contribution only. Also, due to the noise energy spreading, the contribution of the noise in the source time-frequency points is relatively, negligible at least for moderate and high SNRs. Finally,the clustering procedure can be done as follows: For each TF point, we obtain the spatial direction vectors by:

v(t, f ) = Sx(t, f ) S x (t, f ) , (t, f ) ∈ Ω, (8)
and force them, without loss of generality, to have the first entry real and positive.These vectors are clustered into N classes {C i | i ∈ N } by minimizing the criterion:

v(t, f ) ∈ Ci ⇐⇒ i = arg min k v(t, f ) - h k (f )e -jθ k h k (f ) 2 (9)
where h k (f ) is the Fourier Transform of the k th channel vector estimate and θ k is the phase argument of h 1k (f ) (this is to force the first entry to be real positive as for v(t, f )).

The collection of all points, whose vectors belong to the class C i , form the TF support Ω i of source s i . Therefore, we can estimate the STFT of each source s i by:

S s i (t, f ) = h H i (f ) h i (f ) 2 S x (t, f ), ∀ (t, f ) ∈ Ω i , 0, otherwise, (10)
since, from (5), we have

h H i (f ) h i (f ) 2 Sx(t, f ) = h H i (f )h i (f ) h i (f ) 2 Ss i (t, f ) ≈ Ss i (t, f ), ∀ (t, f ) ∈ Ω i .

UBSS algorithm with TF-nondisjoint assumption

We have seen the cluster-based TF-CUBSS method, using the STFT. This method relies on the assumption that the sources are TF-disjoint, which led to the TF-transformed structure in [START_REF] O'grady | Survey of sparse and nonsparse methods in source separation[END_REF]. The latter is no longer valid, when the sources are nondisjoint in the TF domain.

Under the TF-nondisjointness condition, stated in Assumption 2, we propose in this section an alternative method using subspace projection. Recall that the first two steps of the cluster-based quadratic TF-CUBSS algorithm do not rely on the assumption of TF-disjoint sources. Therefore, we can reuse these steps to obtain the channel estimates and all the TF points of the sources, i.e. Ω. Under the TF-nondisjointness condition, consider a TF point (t, f ) ∈ Ω at which there are J < M sources1 sα 1 (t), . . . , sα J (t) present, where α 1 , . . . , α J ∈ N denote the indices of the active sources at (t, f ). Our goal is to identify the sources that are present at (t, f ), i.e. α 1 , . . . , α J , and to estimate the STFT of each of these contributing sources. We define the following:

s = [s α 1 (t), . . . , s α J (t)] T , (11a) Hα(f ) = [hα 1 (f ), . . . , hα J (f )]. (11b)
Then, ( 4) is reduced to the following

S x (t, f ) = Hα (f )Ss(t, f ). (12)
Let Hβ (f ) = [h β 1 (f ), . . . , h β J (f )] and Q β (f ) be the orthogonal projection matrix onto the noise subspace of Hβ (f ) expressed by:

Q β (f ) = I -Hβ (f ) HH β (f ) Hβ (f ) -1 HH β (f ). (13)
We have the following observation:

Q β (f )h i (f ) = 0, i ∈ {β 1 , . . . , β J } Q β (f )h i (f ) = 0, i ∈ N \{β 1 , . . . , β J } . ( 14
)
Consequently, as Sx(t, f ) ∈ Range{ Hα(f )}, we have

Q β (f )S x (t, f ) = 0, if {β 1 , . . . , β J } = {α 1 , . . . , α J } Q β (f )Sx(t, f ) = 0, otherwise .
(15) Since H(f ) has already been estimated by the method presented in Section 3.1, then this observation gives us the criterion to detect the indices α1, . . . , αJ and hence, the contributing sources at the considered TF point (t, f ). In practice, to take into account noise, one detects the column vectors of Hα(f ) by minimizing:

{α 1 , . . . , α J } = arg min β 1 ,...,β J { Q β (f )S x (t, f ) } . ( 16
)
Next, TFD values of the J sources at the considered TF point are estimated by:

Ss(t, f ) ≈ H# α (f )S x (t, f ), (17)
where the superscript (•) # represents the Moore-Penrose's pseudoinversion operator.

In the simulation, the optimization problem of ( 16) is solved using exhaustive search. This is computationally tractable for small vector array sizes but would be prohibitive if M is very large.

SIMULATION RESULTS

In the simulations, we have considered an array of M = 3 sensors, that receives signals from N = 4 independent speech sources. The filter coefficients are chosen randomly and the channel order is K = 6. The sample size is T = 8192 samples (corresponding approximately to 1 second recording of speech signals sampled at 8 KHz). The separation quality is measured by the normalized mean squares estimation errors (NMSE) of the sources evaluated over N r = 200 Monte-Carlo runs and defined as:

NMSE i def = 1 N r Nr r=1 min α α si,r -si 2 s i 2 (18) NMSEi = 1 Nr N r r=1 1 - s i,r s H i si,r si 2 (19) NMSE = 1 N N i=1 NMSEi . ( 20
)
where s i def = [s i (0), . . . , s i (T -1)], s i,r ( defined similarly) represents the r th estimate of source s i and α is a scalar factor that compensates for the scale indeterminacy of the BSS problem. In Fig. 1, we compare the separation performance obtained by the subspacebased algorithm with J = 2 and the cluster-based algorithm. It is observed that subspace-based algorithm provides much better separation results than those obtained by the cluster-based algorithm. This is mainly due to the high occurrence of overlapping sources in the TF domain for this type of signals so that the 'TF-disjointness' assumption used by the TF-CUBSS algorithm is poorly satisfied. The plot in Fig. 2 presents the separation performance of the sub- space method when using the exact matrix H compared to that obtained with the proposed estimate H. The observed performance loss is due to the channel estimation error which is relatively high for low SNRs and becomes negligible for high SNRs.

CONCLUSION

This paper introduces new methods for the UBSS of TF-disjoint and TF-nondisjoint nonstationary sources in the convolutive mixture case using their time-frequency representations. The first proposed method has the advantage of simplicity while the second uses a weaker assumption on the source 'sparseness', i.e. the sources are not necessarily TF-disjoint, and proposes an explicit treatment of the overlapping points using subspace projection, leading to significant performance improvements.

Fig. 1 .

 1 Fig. 1. Comparison between subspace-based and cluster-based TF-CUBSS algorithms: normalized MSE (NMSE) versus SNR for 4 speech sources and 3 sensors.

Fig. 2 .

 2 Fig. 2. Comparison, for the subspace-based TF-CUBSS algorithm, when the mixing channel H is known or unknown: NMSE of the source estimates.

To avoid the difficult problem of estimating the number of active sources at each TF point, we have chosen in this paper to set J to a fixed value in the range 1 < J < M .