Fangchen Feng

Matthieu Kowalski

Underdetermined Reverberant Blind Source Separation: Sparse Approaches for Multiplicative and Convolutive Narrowband Approximation

We consider the problem of blind source separation for underdetermined convolutive mixtures. Based on the multiplicative narrowband approximation in the time-frequency domain with the help of Short-Time-Fourier-Transform (STFT) and the sparse representation of the source signals, we formulate the separation problem into an optimization framework. This framework is then generalized based on the recently investigated convolutive narrowband approximation and the statistics of the room impulse response. Algorithms with convergence proof are then employed to solve the proposed optimization problems. The evaluation of the proposed frameworks and algorithms for synthesized and live recorded mixtures are illustrated. The proposed approaches are also tested for mixtures with input noise. Numerical evaluations show the advantages of the proposed methods.

I. INTRODUCTION A. Time model

Blind source separation (BSS) recovers source signals from a number of observed mixtures without knowing the mixing system. Separation of the mixed sounds has several applications in the analysis, editing, and manipulation of audio data [START_REF] Comon | Handbook of Blind Source Separation: Independent component analysis and applications[END_REF]. In the real-world scenario, convolutive mixture model is considered to take the room echo and the reverberation effect into account:

x m (t) = N n=1 a mn (t) * s n (t) + n m (t), (1)
where s n is the n-th source and x m is the m-th mixture. N and M are the number of sources and microphones respectively. a mn (t) is the room impulse response (RIR) from the nth source to the m-th microphone. n m (t) is the additive white Gaussian noise at the m-th microphone. We denote also s img mn (t) = a mn (t) * s n (t) the image of the n-th source at the m-th microphone.

B. Multiplicative narrowband approximation

The source separation for convolutive mixtures is usually tackled in the time-frequency domain with the help of STFT (Short-Time-Fourier-Transform) [START_REF] Winter | MAP-based underdetermined blind source separation of convolutive mixtures by hierarchical clustering and 1 -norm minimization[END_REF], [START_REF] Sawada | Underdetermined convolutive blind source separation via frequency bin-wise clustering and permutation alignment[END_REF], [START_REF] Duong | Under-determined reverberant audio source separation using a full-rank spatial covariance model[END_REF]. With the narrowband assumption, the separation can be performed in each Fangchen Feng is with Laboratoire Astroparticule et Cosmologie, Université Paris Diderot, CNRS/IN2P3, Sorbonne Paris Cité, 75205, Paris, France (email: fangchen.feng@apc.in2p3.fr) Matthieu Kowalski is with Laboratoire des signaux et systèmes, CNRS, Centralesupélec, Université Paris-Sud, Université Paris-Saclay, 91192, Gifsur-Yvette, France (email: matthieu.kowalski@l2s.centralesupelec.fr) frequency band [START_REF] Kellermann | Wideband algorithms versus narrowband algorithms for adaptive filtering in the DFT domain[END_REF]. Because of the permutation ambiguity in each frequency band, the separation is then followed by a permutation alignment step to regroup the estimated frequency components that belong to the same source [START_REF] Sawada | Measuring dependence of binwise separated signals for permutation alignment in frequency-domain bss[END_REF]. In this paper, we concentrate on the separation step.

The multiplicative narrowband approximation [START_REF] Winter | MAP-based underdetermined blind source separation of convolutive mixtures by hierarchical clustering and 1 -norm minimization[END_REF], [START_REF] Sawada | Underdetermined convolutive blind source separation via frequency bin-wise clustering and permutation alignment[END_REF] deals with the convolutive mixtures in each frequency using the complex-valued multiplication in the following vector form:

x(f, τ ) = N n=1 ãn (f )s n (f, τ ) + ñ(f, τ ), (2)
where x(f, τ ) = [x 1 (f, τ ), . . . , xM (f, τ )] T and sn (f, τ ) are respectively the analysis STFT coeffcients of the observations and the n-th source signal. ãn (f ) = [ã 1n (f ), . . . , ãMn (f )]

T is a vector that contains the Fourier transform of the RIR associated with the n-th source. ñ(f, τ ) = [ñ 1 (f, τ ), . . . , ñM (f, τ )] T consistes not only the analysis STFT coefficients of the noise, but also the error term due to the approximation. The formulation [START_REF] Winter | MAP-based underdetermined blind source separation of convolutive mixtures by hierarchical clustering and 1 -norm minimization[END_REF] approximates the convolutive mixtures by using instantaneous mixture in each frequency. This approximation therefore largely reduces the complexity of the problem and is valid when the RIR length is less than the STFT window length. The sparsity assumption is largely utilized for source separation problem [START_REF] Winter | MAP-based underdetermined blind source separation of convolutive mixtures by hierarchical clustering and 1 -norm minimization[END_REF], [START_REF] Sawada | Underdetermined convolutive blind source separation via frequency bin-wise clustering and permutation alignment[END_REF], [START_REF] Zibulevsky | Blind source separation by sparse decomposition in a signal dictionary[END_REF], [START_REF] Yilmaz | Blind separation of speech mixtures via time-frequency masking[END_REF]. Based on the model [START_REF] Winter | MAP-based underdetermined blind source separation of convolutive mixtures by hierarchical clustering and 1 -norm minimization[END_REF] and by supposing that only one source is active or dominant in each time-frequency bin (f, τ ), the authors of [START_REF] Winter | MAP-based underdetermined blind source separation of convolutive mixtures by hierarchical clustering and 1 -norm minimization[END_REF] proposed to estimate the mixing matrix in each frequency by clustering, and then estimate the source in a maximum a posteriori (MAP) sense. This method is further improved by [START_REF] Sawada | Underdetermined convolutive blind source separation via frequency bin-wise clustering and permutation alignment[END_REF] where the authors proposed to use a soft masking technique to perform the separation. The idea is to classify each time-frequency bin of the observation x(f, τ ) into N class, where N is the number of sources. Based on a complex-valued Gaussian generative model for source signals, they inferred a bin-wise a posteriori probability P (C n |x(f, τ )) which represents the probability that the vector x(f, τ ) belongs to the n-th class C n . This method obtains good separation results for speech signals, however only in low reverberation scenario [START_REF] Sawada | Underdetermined convolutive blind source separation via frequency bin-wise clustering and permutation alignment[END_REF]. The performance of these methods is limited by the multiplicative approximation whose approximation error increases rapidely as the reverberation time becomes long [START_REF] Kowalski | Beyond the narrowband approximation: Wideband convex methods for under-determined reverberant audio source separation[END_REF]. Moreover, the disjointness of the soures in the time-frequency domain is not realistic [START_REF] Yilmaz | Blind separation of speech mixtures via time-frequency masking[END_REF].

C. Beyond the multiplicative narrowband model

A generalization of the multiplicative approximation is proposed in [START_REF] Duong | Under-determined reverberant audio source separation using a full-rank spatial covariance model[END_REF] by considering the spatial covariance matrix of the source signals. By modeling the sources STFT coefficients as a phase-invariant multivariate distribution, the authors inferred that the covariance matrix of the STFT coefficients of the n-th source images s img n = s img 1n , s img 2n , . . . , s img

M n

T can be factorized as:

R simg n (f, τ ) = v n (f, τ )R n (f ), (3)
where v n (f, τ ) are scalar time-varying variances of the n-th source at different frequencies and R n (f ) are time-invariant spatial covariance matrices encoding the source spatial position and spatial spread [START_REF] Duong | Under-determined reverberant audio source separation using a full-rank spatial covariance model[END_REF]. The multiplicative approximation forces the spatial covariance matrix to be of rank-1 and the authors of [START_REF] Duong | Under-determined reverberant audio source separation using a full-rank spatial covariance model[END_REF] exploited a generalization by assuming that the spatial covariance matrix is of full-rank and showed that the new assumption models better the mixing process because of the increased flexibility. However, as we show in this paper by experiments, the separation performance of this full-rank model is still limited in strong reverberation scenarios. Moreover, as both the bin-wise method [START_REF] Sawada | Underdetermined convolutive blind source separation via frequency bin-wise clustering and permutation alignment[END_REF] and the full-rank approach [START_REF] Duong | Under-determined reverberant audio source separation using a full-rank spatial covariance model[END_REF] do not take the error term ñ into consideration, they are therefore sensitive to additional noise.

Recently, the authors of [START_REF] Li | Audio source separation based on convolutive transfer function and frequency-domain Lasso optimization[END_REF] investigated the convolutive narrowband approximation for oracle source separation of convolutive mixtures (the mixing systems is known). They showed that the convolutive approximation suits better the original mixing process especially in strong reverberation scenarios. In this paper, we investigate the convolutive narrowband approximation as the generalization of the multiplicative approximation in the full blind setting (the mixing system if unknown).

The contribution of the paper is three-folds: first based on the multiplicative narrowband approximation, we formulate the separation in each frequency as an optimization problem with 1 norm penalty to exploit sparsity. The proposed optimization formulation is then generalized based on the statistics of the RIR [START_REF] Benichoux | Convex regularizations for the simultaneous recording of room impulse responses[END_REF] and the convolutive narrowband approximation model [START_REF] Li | Audio source separation based on convolutive transfer function and frequency-domain Lasso optimization[END_REF]. At last, we propose to solve the obtained optimizations with PALM (Proximal alternating linearized minimization) [START_REF] Bolte | Proximal alternating linearized minimization for nonconvex and nonsmooth problems[END_REF] and BC-VMFB (Block coordinatevariable metric forward backward) [START_REF] Chouzenoux | A block coordinate variable metric forward-backward algorithm[END_REF] algorithms which have convergence guarantee.

The rest of the article is organized as follows. We propose the optimization framework based on multiplicative approximation with 1 norm penalty and present the corresponding algorithm in Section II. The optimization framework is then generalized in Section III based on the statistics of the RIR and the convolutive approximation. The associated algorithm is also presented. We compare the separation performance achieved by the proposed approaches to that of the state-ofthe-art in various experimental settings in Section IV. Finally, Section V concludes the paper.

II. THE MULTIPLICATIVE NARROWBAND APPROXIMATION

We first rewrite the formulation (2) with matrix notations by concatenating the time samples and source indexes. In each frequency f , we have:

Xf = Ãf Sf + Ñf , (4)
where Xf ∈ C M ×L T is the matrix of the analysis STFT coefficients of the observations at the given frequency f . Ãf ∈ C M ×N is the mixing matrix at frequency f . Sf ∈ C N ×L T is the matrix of the analysis STFT coefficients of the sources at frequency f . Ñf ∈ C M ×L T is the noise term which also contains the approximation error. In the above notations, L T is the number of time samples in the time-frequency domain.

The target of the separation is to estimate Ãf and Sf from the observations Xf . However, according to the definition of the analysis STFT coefficients , the estimated Sf has to be in the image of the STFT operator (see in [START_REF] Balazs | Adapted and adaptive linear time-frequency representations: a synthesis point of view[END_REF] for more details). To avoid this additional constraint, we propose to replace the analysis STFT coefficients Sf by the synthesis STFT coefficients α f ∈ C N ×L T , which leads to:

Xf = Ãf α f + Ñf . (5)
In the following, we denote also α f,n the n-th source componant (row) of α f and α f,n (τ ) the scalar element at position τ in α f,n .

A. Formulation of the optimization

Based on the model ( 5), we propose to formulate the separation as an optimization problem as follow:

min Ãf ,α f 1 2 Xf -Ãf α 2 F + λ α f 1 + ı C ( Ãf ), (6)
where • F denotes the Frobenius norm and • 1 is the 1 norm of the matrix which is the sum of the absolute value of all the elements. ı C ( Ãf ) is an indicator function to avoid the trivial solution caused by the scaling ambiguity between Ãf and α f :

ı C ( Ãf ) = 0, if ãf,n = 1, n = 1, 2, . . . , N + ∞, otherwise (7)
with ãf,n the n-th column of Ãf . λ is the hyperparameter which balances between the data term 1 2 Xf -Ãf α f 2 F and the penalty term α f 1 .

For instantaneous mixtures, the formulation (6) has been firstly proposed in [START_REF] Zibulevsky | Blind source separation by sparse decomposition in a signal dictionary[END_REF] and recently investigated in [START_REF] Feng | A unified approach for blind source separation using sparsity and decorrelation[END_REF]. Compared to the masking technique of separation [START_REF] Sawada | Underdetermined convolutive blind source separation via frequency bin-wise clustering and permutation alignment[END_REF], the 1 norm term exploits only sparsity which is more realistic than the disjointness assumption for speech signals. Moreover, the Lagrangian form with the data term 1 2 Xf -Ãf α f 2 F allows us to take the noise/approximation error into consideration.

B. Algorithm: N-Regu

The optimization problem ( 6) is non-convex with a nondifferentiable term. In this paper, we propose to solve the problem by applying the BC-VMFB (block coordinate variable metric forward-backward) [START_REF] Chouzenoux | Variable metric forwardbackward algorithm for minimizing the sum of a differentiable function and a convex function[END_REF] algorithm. This algorithm relies on the proximal operator [START_REF] Combettes | Proximal splitting methods in signal processing[END_REF] given in the next definition.

Definition 1. Let ψ be a proper lower semicontinuous function, the proximal operator associated with ψ is defined as:

prox ψ := argmin y ψ(y) + 1 2 y -x 2 F . (8)
When the function ψ(y) = λ y 1 , the proximal operator becomes the entry-wise soft thresholding presented in the next proposition.

Proposition 2. Let α ∈ C N ×L T . Then, α = prox λ • 1 (α) := S λ (α) is given entrywise by soft-thresholding:

αi = α i |α i | (|α i | -λ) + , (9)
where (α) + = max(0, α).

When the function ψ in Definition 1 is the indicator function ı C , the proximal operator reduces to the projection operator presented in Proposition 3.

Proposition 3. Let à ∈ C M ×N . Then  = prox ı C ( Ã) := P C ( Ã)
is given by the column-wise normalization projection:

ân = ãn ãn , n = 1, 2, . . . , N (10)
With the above proximal operators, we present the algorithm derived from BC-VMFB in Algorithm 1. We denote the data term by

Q(α f , A f ) = 1 2 Xf -Ãf α f 2 F . L (j) = Ã(j+1) H f Ã(j+1) f 2 is the Lipschitz constant of ∇ α f Q(α (j) f , Ã(j) f ) with • 2
denoting the spectral norm of matrix. Details of the derivation of this algorithm and the convergence study are given in Appendix VI-A. In the following, this algorithm is referred as N-Regu (Narrowband optimization with regularization).

Algorithm 1: N-Regu Initialisation : α (1) f ∈ C N ×L T , Ã(1) f ∈ C M ×N , L (1) = Ã(1) H f Ã(1) f 2 , j = 1; repeat ∇ α f Q α (j) f , Ã(j) f = - Ã(j) H f Xf - Ã(j) f α (j) f ; α (j+1) f = S λ/L (j) (α (j) f -1 L (j) ∇ α f Q(α (j) f , Ã(j) f ); Ã(j+1) f = P C ( Xf α (j+1) H f ); L (j+1) = Ã(j+1) H f Ã(j+1) f 2 ; j = j + 1; until convergence;

III. THE CONVOLUTIVE NARROWBAND APPROXIMATION

A. Convolutive approximation

Theoretically, the multiplicative narrowband approximation ( 2) is valid only when the RIR length is less than the STFT window length. In practice, this condition is rarely varified because the STFT window length is limited to ensure the local stationarity of audio signals [START_REF] Li | Audio source separation based on convolutive transfer function and frequency-domain Lasso optimization[END_REF]. To avoid such limitation, the convolutive narrowband approximation was proposed in [START_REF] Avargel | System identification in the short-time fourier transform domain with crossband filtering[END_REF], [START_REF] Talmon | Relative transfer function identification using convolutive transfer function approximation[END_REF]:

x(f, τ ) = N n=1 L l=1 hn (f, l)s n (f, τ -l), (11)
where hn = h1n , . . . , hMn T is the vector that contains the impulse responses in the time-frequency domain associated with the n-th source. L is the length of the convolution kernel in the time-frequency domain.

The convolutive approximation ( 11) is a generalization of the multiplicative approximation (2) as it considers the information diffusion along the time index. When the kernel length L = 1, it reduces to the multiplicative approximation. The convolution kernel in the time-frequency domain hmn (f, τ ) is linked to the RIR in the time domain a mn (t) by [START_REF] Li | Audio source separation based on convolutive transfer function and frequency-domain Lasso optimization[END_REF]:

hmn (f, τ ) = [a mn (t) * ζ f (t)] | t=τ k0 , (12)
which represents the convolution with respect to the time index t evaluated with a resolution of the STFT frame step k 0 with:

ζ f (t) = e 2πif t/L F j ϕ(j) φ(t + j), (13)
where L F is the number of frequency bands. ϕ(j) et φ(j) denote respectively the analysis and synthesis STFT window. With matrix notations, for each frequency f , the convolutive approximation ( 11) can be written as:

Xf = Hf Sf + Ñf , (14)
where Hf ∈ C M ×N ×L is the mixing system formed by concatenating the impulse responses of length L. In the following, we denote also hf,mn the vector that represents the impulse response at position (m, n) in Hf and hf,mn (τ ) the scalar element at position (m, n, τ ). The operator denotes the convolutive mixing process [START_REF] Benichoux | Convex regularizations for the simultaneous recording of room impulse responses[END_REF].

Compared to the original mixing process in the time domain (1), the convolutive approximation ( 14) largely reduces the length of the convolution kernel, thus makes the estimation of both the mixing system and the source signals practically possible.

B. Proposed optimization approach a) Basic extension of the multiplicative model: Once again, to circumvent the additional constraint brought by the analysis coefficients of the sources, we replace the analysis STFT coefficients Sf by the synthesis coefficients α f , which leads to:

Xf = Hf α f + Ñf . ( 15
)
Based on (15), we generalize ( 6) by replacing the multplicative operator by the convolutive mixing operator:

min Hf ,α f 1 2 Xf -Hf α f 2 F + λ α f 1 + ı Conv C ( Hf ), ( 16
)
where ı Conv C ( Hf ) is the normalisation constraint to avoid trivial solutions:

ı Conv C ( Hf ) =      0, if m,τ | hf,mn (τ )| 2 = 1, n = 1, . . . , N + ∞, otherwise.
(17) b) Regularization for the convolution kernel: In [START_REF] Benichoux | Convex regularizations for the simultaneous recording of room impulse responses[END_REF], the authors consider the problem of estimating the RIR supposing that the mixtures and the sources are known. They formulated the estimation problem as an optimization problem and proposed a differentiable penalty for the mixing system in the time domain:

m,n,t |a mn (t)| 2 2ρ 2 (t) , (18)
where ρ(t) denotes the amplitude envelope of RIR which depends on the reverberation time RT 60 :

ρ(t) = σ10 -3t/RT60 , (19)
with σ being a scaling factor. The penalty ( 18) is designed to force an exponential decrease of the RIR which satisfaits the acoustic statistics of the RIR [START_REF] Kuttruff | Room acoustics[END_REF].

As the convolutive kernel in the time-frequency domain is linked to the RIR in time domain by [START_REF] Bolte | Proximal alternating linearized minimization for nonconvex and nonsmooth problems[END_REF], in this paper, we consider the penalty in the time-frequency domain in the same form:

P( Hf ) = m,n,τ | hf,mn (τ )| 2 2ρ 2 (τ ) , (20)
where ρ(τ ) is the decreasing coefficients in the time-frequency domain which depends on ρ(t) and the STFT transform.

Other forms of penalty are also proposed in [START_REF] Benichoux | Convex regularizations for the simultaneous recording of room impulse responses[END_REF]. However, their adaption in the time-frequency domain is not straightforward.

c) Final optimization problem: With the above penalty term, the formulation ( 16) can be improved as:

min Hf ,α f 1 2 Xf -Hf α f 2 F + λ α f 1 + P( Hf ) + ı Conv C ( Hf ). (21)

C. Algorithm: C-PALM

We propose to use the Proximal Alternating Linearized Minimization (PALM) algorithm [START_REF] Bolte | Proximal alternating linearized minimization for nonconvex and nonsmooth problems[END_REF] to solve the problem. The derived algorithm is presented in Algorithm 2, and one can refer to Appendix VI-C for details on the derivation and the convergence study. We refer to this algorithm as C-PALM (Convolutive PALM) in the following. We denote:

Q( αf , Hf ) = 1 2 X -Hf α f 2 F + P( Hf )
and the gradient of P( Hf ) is given coordinate-wise by:

∇ Hf P( Hf ) f,mnτ = hf,mn (τ ) ρ4 (τ ) . ( 22
)
In Algorithm 2, HH f and α H f are respectively the adjoint operators of the convolutive mixtures with respect to the convolution kernel and the sources. Details of derivation of these adjoint operators are given in Appendix VI-B. L (j)

α f and L (j) Hf are respectively the Lipschitz constant of ∇ α f Q(α (j) f , H(j) f ) and ∇ Hf Q(α (j+1) f , H(j) f ). L (j)
α f can be calculated with the power iteration algorithm [START_REF] Kowalski | Beyond the narrowband approximation: Wideband convex methods for under-determined reverberant audio source separation[END_REF] shown in Algorithm 3. L (j) Hf can be approximately estimated thanks to the next proposition.

Algorithm 2: C-PALM Initialisation : α (1) f ∈ C N ×L T , H(1) f ∈ C M ×N , j = 1; repeat ∇ α f Q α (j) f , H(j) f = - H(j) H f Xf - H(j) f α (j) ; α (j+1) f = S λ/L (j) α f α (j) f -1 L (j) α f ∇ α f Q(α (j) f , H(j) f ) ; ∇ Hf Q(α (j+1) f , H(j) f ) = -( Xf - H(j) f α (j+1) f ) α (j+1) H f + ∇ Hf P( H(j) f ); H(j+1) f = P Conv ı C H(j) f -1 L (j) Hf ∇ Hf Q(α (j+1) f , H(j) f ) ; Update L (j) α f et L (j) Hf ; j = j + 1; until convergence; Algorithm 3: Power iteration for the calculation of L α f Initialisation : v f ∈ C N ×L T ; repeat W = HH f Hf v f ; L α f = W ∞ ; v f = W Lα f ; until convergence;
Proposition 4. If we suppose that the source componants α f,1 , α f,2 , . . . , α f,N are mutually independant and L L T , then L Hf , the Lipschitz constant of ∇ Hf Q(α f , Hf ) can be calculated as:

L Hf = max n (L f,n ) + max τ ( 1 ρ8 (τ ) ), (23)
where

L f,n = Γ f,n , with Γ f,n =      γ f,n (0) γ f,n (1) . . . γ f,n (L -1) γ f,n (-1) γ f,n (0) . . . γ f,n (L -2) . . . . . . . . . . . . γ f,n (1 -L) γ f,n (2 -L) . . . γ f,n (0)      , (24)
and γ f,n (τ ) is the empirical autocorrelation function of α f,n :

γ f,n = L T -1 =1 α f,n ( + τ )α * f,n ( ). (25)
Proof. The proof is postponed in Appendix VI-D.

If the independance assumption mentioned in Proposition 4 appears to be strong, it is well adapted for audio signals as it is the basic hypothesis of the FDICA (frequency domain independant component analysis) [START_REF] Sawada | A robust and precise method for solving the permutation problem of frequency-domain blind source separation[END_REF] used for source separation of determined convolutive mixtures. Although we do not have any guarantee of independence in the proposed algorithm, the experiments show that good performances are obtained.

Finally, we must stress that the BC-VMFB algorithm is not suitable for [START_REF] Sawada | A robust and precise method for solving the permutation problem of frequency-domain blind source separation[END_REF] as it relies on the second derivative of Q( αf , Hf ) w.r.t Hf , which does not necessarily simplify the algorithm.

IV. EXPERIMENTS A. Permutation alignment methods

For the proposed approaches, we use the existing permutation alignment methods. For N-Regu, we compare the approach based on TDOA (Time Difference Of Arrival) used in Full-rank method [START_REF] Duong | Under-determined reverberant audio source separation using a full-rank spatial covariance model[END_REF] and the approach based on interfrequency correlation used in the Bin-wise approach [START_REF] Sawada | Underdetermined convolutive blind source separation via frequency bin-wise clustering and permutation alignment[END_REF]. For the inter-frequency correlation permutation, we use the power ratio [START_REF] Sawada | Measuring dependence of binwise separated signals for permutation alignment in frequency-domain bss[END_REF] of the estimated source to present the source activity. For C-PALM, as the TDOA permutation is not adapted, we use only the correlation permutation.

For the proposed approaches (N-Regu and C-PALM) and the reference algorithms (Bin-wise and Full-rank), we also designed an oracle permuation alignment method. In each frequency, we look for the permutation that maximizes the correlation between the estimated and the original sources. Such permutation alignment is designed to show the best permutation possible in order to have a fair comparison of the separation approaches instead of the choice made for solving the permutation problem.

B. Experimental setting

We first evaluated the proposed approaches with 10 sets of synthesized stereo mixtures (M = 2) containing three speech sources (N = 3) of male/female with different nationalities. The mixtures are sampled at 11 kHz and truncated to 6 s. The room impulse response were simulated via the toolbox [START_REF] Lehmann | Prediction of energy decay in room impulse responses simulated with an image-source model[END_REF]. The distance between the two microphone is 4 cm. The reverberation time is defined as 50 ms, 130 ms, 250 ms and 400 ms. The Fig. 1 illustrates the room configuration. For each mixing situation, the mean values of the evaluation results over the 10 sets of mixtures are shown.

We then evaluated the algorithm C-PALM with the live recorded speech mixtures from the dataset SiSEC2011 [START_REF] Araki | The 2011 signal separation evaluation campaign (sisec2011):-audio source separation[END_REF]. Music mixtures are avoided because the instrumental sources are often synchronized to each other and this situation is difficult for the permutation alignment based on inter-frequency correlation [START_REF] Sawada | Underdetermined convolutive blind source separation via frequency bin-wise clustering and permutation alignment[END_REF]. An effective alternative way is to employ nonnegative matrix factorization [START_REF] Feng | Sparsity and low-rank amplitude based blind source separation[END_REF]. The parameters of STFT for the synthesized and live recorded mixtures are summarized in Table I. The STFT window length (and window shift) for synthesized mixtures are chosen to preserve local stationarity of audio sources without bringing too much computational costs. The parameters for the live recorded mixtures are the same as the reported reference algorithm Bin-wise [START_REF] Sawada | Underdetermined convolutive blind source separation via frequency bin-wise clustering and permutation alignment[END_REF]. The separation performance is evaluated with the signal to distortion ratio (SDR), signal to interference ratio (SIR), source image to spatial distortion ratio (ISR) and signal to artifact ratio (SAR) [START_REF] Vincent | First stereo audio source separation evaluation campaign: data, algorithms and results[END_REF]. The SDR reveals the overall quality of each estimated source. SIR indicates the crosstalk from other sources. ISR measures the amount of spatial distortion and SAR is related to the amount of musical noise.

N-Regu is initialized with Gaussian random signals. C-PALM is initialized with the results of N-Regu with 1000 iterations. This choice of initialization for C-PALM compensates the flexibility of the convolutive model (and then the number of local minima in ( 21)) without bringing too much computational cost. We use the stopping criteria α (j+1) f α j f F < 10 -4 for both algorithms.

C. Tuning the parameters

For the proposed methods, we chose several pre-defined hyperparameter λ and select the λ which corresponds to the best SDR. Even though such a way of choosing this hyper-parameter is not possible for real applications, such evaluation offers a "fair" comparison with the state-of-the-art approaches and gives some empirical suggestions of choosing this parameter in practice. We implement the continuation trick, also known as warm start or fixed point continuation [START_REF] Hale | Fixed-point continuation for 1minimization: Methodology and convergence[END_REF] for a fixed value of λ: we start the algorithms with a large value of λ and iteratively decrease λ to the desired value.

It is also important to mention that the hyperparameter λ should be theoretically different for each frequency since the sparsity level of the source signals in each frequency can be very different (for speech signals, the high frequency componants are usually sparser than the low frequency componants). Therefore, different λ should be determined for each frequency. However, in this paper, we used a single λ for all the frequencies and the experiments show that this simplified choice can achieve acceptable results if we perform a whitening pre-processing for each frequency.

For C-PALM, as the reverberation time is unknown in the blind setting, we pre-define the length of the convolution kernel in the time-frequency domain L = 3 as well as the penalty parameter ρ(τ ) = [1.75, 1.73, 1.72]

T . Although these parameters should vary with the reverberation time, we show in the following that the proposed pre-defined parameters work well in different strong reverberation conditions.

D. Synthesized mixtures without noise

We first evaluate the algorithms with synthesized mixtures in the noiseless case as a function of the reverberation time RT 60 . The results are shown in Fig. 2.

For RT 60 = 50 ms, it is clear that the Full-rank method performs the best in terms of all four indicators. Its good performance is due to the fact that the full-rank spatial covariance model suits better the convolutive mixtures than the multiplicative approximation and the fact that the TDOA permutation alignement has relatively good performance in low reverberation scenario. N-Regu outperforms Bin-wise only in SDR and SAR. It is because that N-Regu has better data fit than the masking-based Bin-wise method while Binwise obtains time-frequency domain disjoint sources which have lower inter-source interference. C-PALM is dominated by other methods in SDR, SIR and ISR. We believe it is because that the pre-defined penalty parameter ρ(τ ) does not fit the low reverberation scenario. The advantages of C-PALM can be seen in relatively stronger reverberation scenarios (especially RT 60 = 130, 250 ms) where C-PALM outperforms other methods in SDR and SIR. For RT 60 = 400 ms, all the presented algorithms have similar performance while C-PALM performs slightly better in SIR. To compare the two permutation methods used for N-Regu, TDOA permutation performs better than inter-frequency correlation permutation in SDR, SIR and SAR.

Fig. 3 compares the presented algorithms with oracle permutation alignment. For RT 60 = 50 ms, once again, Full-rank has the best performance in all four indicators. This confirms the advantages of the full rank spatial covariance model. In high reverberation conditions, C-PALM performs better than others in SDR and SIR. In particular, for RT 60 = 130, 250 ms, C-PALM outperforms Full-rank by more than 1 dB in SDR and outperforms Bin-wise by about 1.2 dB in SIR. N-Regu performs slightly better than Bin-wise in SDR for all reverberation conditions.

The above observations show the better data fit brought by the optimization framework used in N-Regu (and C-PALM) and confirm the advantages of convolutive narrowband approximation used in C-PALM for high reverberation conditions (especially RT 60 = 130, 250 ms). Fig. 4 illustrates the performance of the presented algorithm as a function of the sparsity level1 of the estimated synthesis coeffcients of the sources for RT 60 = 130 ms. As the sparsity level is directly linked to the hyperparameter λ in the proposed algorithms, this comparison offers some suggestions of choosing this hyperparameter. Full-rank method does not exploits sparsity, thus has 0% as sparsity level. As the number of sources N = 3, the sparsity level of the masking-based Bin-wise method is 66.6%.

C-PALM performs better than N-Regu in terms of SDR, SIR and SAR when the sparsity level is less than 60% and its best performance is achieved when the sparsity level is around 40%. For N-Regu, in terms of SDR, SAR and ISR, the best performance is achieved with the least sparse result.

E. Synthesized mixtures with noise

In this subsection, we evaluate the proposed methods with synthesized mixtures with additive white Gaussian noise. The noise of different energy is added which leads to different input SNR. Fig. 5 reports the separation performance as a function of input SNR with the reverberation time fixed to RT 60 = 130 ms.

It is clear that N-Regu with TDOA permutation outperforms other methods in terms of SDR and SIR. In particular, it performs better than others by about 1 dB in SIR for all the input SNR tested. C-PALM outperforms the state-of-the-art approaches only in SDR. We believe that it is due to the fact that the freedom degree of the convolutive narrowband approximation used in C-PALM could be sensitive to input noise. Another reason is that the inter-frequency correlation based permutation could be sensitive to input noise. The latter conjecture is supported by the observation that, in terms of SDR and SIR, the gap between N-Regu with TDOA permutation and with correlation permutation increases as the input noise becomes stronger. Further evidence can be found by the comparisons between the presented algorithm with oracle permutation alignment in Fig. 6.

In Fig. 6, in terms of SDR and SIR, it is clear that the gap between N-Regu and C-PALM decreases as the input noise gets stronger. This remark confirms that the separation step of C-PALM is sensitive to input noise. Moreover, in terms of SIR, C-PALM with oracle permutation performs consistently better than N-Regu with oracle permutation, while C-PALM with correlation permutation is dominated by N-Regu with TDOA permutation by about 1 dB (Fig. 5). This observation shows that the performance of C-PALM can be largely improved for noisy mixtures if better permutation alignment method is developped. Fig. 7 reports the separation performance as a function of the sparsity level of the estimated synthesis coefficients of the sources. RT 60 = 130 ms and the input SNR is 15 dB. The results of Full-rank and Bin-wise method are also shown.

In terms of SDR and SIR, N-Regu with TDOA permutation consistently outperforms the other methods and achieves its best performance when the sparsity level is about 78%. Compared to Bin-wise method, this observation coincides with the intuition that, for noisy mixtures, the coefficients of the noise in the observations should be discarded to achieve better separation. C-PALM achieves its best performance in terms of SDR and SIR when the sparsity level is about 75%.

Fig. 8 illustrates the results of separation as a function of the reverberation time for a fixed input SNR (SNR=15 dB).

We can see that N-Regu with TDOA permutation has the best performance in terms of SDR.

F. Synthesized mixtures with different sources positions

In this subsection, we tested the robusteness of the proposed algorithms w.r.t the sources positions. The same room setting as shown in Fig. 1 is used. Fig. 9 illustrates the four tested In these experiments, the reverberation time is fixed to RT 60 = 130 ms and no noise is added to the mixtures. Fig. 10 shows the separation performance.

It is clear that in terms of SDR, SIR and ISR, all the presented algorithms have the worst performance in setting 3.

This remark shows that having two sources close to each other and one source relatively far (setting 3) could be a more difficult situation for blind source separation than having three sources close to each other (setting 4). For C-PALM, it has the best performance in terms of SDR, SIR and ISR for all the settings. This observation shows that C-PALM (and the pre-defined penalty parameter) is robust to sources positions,

G. Live recorded mixtures without noise

This subsection reports the separation results of C-PALM for publicly avaiable benchmark data in SiSEC2011 [START_REF] Araki | The 2011 signal separation evaluation campaign (sisec2011):-audio source separation[END_REF]. We used the speech signals (male3, female3, male4 and female4) from the first development data (dev1.zip) in "Under-determined speech and music mixtures" data sets. Table II shows the separation results. For C-PALM, we chose the hyperparameter λ such that the sparsity level of the estimated coefficients of the sources is about 20%, 60% for RT 60 = 130, 250 ms respectively. Compared to the performances reported in [START_REF] Araki | The 2011 signal separation evaluation campaign (sisec2011):-audio source separation[END_REF], C-PALM obtains relatively good separation results epsecially when the number of sources N = 3.

H. Computational time

We terminate the expriment section by presenting the computational time of the presented algorithm for the synthesized mixtures in Table III. C-PALM is of relative big computational cost mainly because of the convolution operator in each iteration of the algorithm.

V. CONCLUSION

In this paper, we have developped several approaches for blind source separation with underdetermined convolutive mixtures. Based on the sparsity assumption for the source signals and the statistics of the room impulse response, we developed the N-Regu with multiplicative narrowband approximation and C-PALM with convolutive narrowband approximation. The numerical evaluations show the advantages of C-PALM for noiseless mixtures in strong reverberation scenarios. The experiments also show the good performance of N-Regu for noisy mixtures. The penalty parameter ρ(τ ) in C-PALM has to be predefined, which makes C-PALM not suitable for low reverberation condition. Future work will concentrate on the estimation of ρ(τ ). In this paper, we used inter-frequency correlation permutation alignment for C-PALM. It would be interesting to exploit TDOA based permutaiton method for convolutive narrowband approximation to improve C-PALM.

VI. APPENDIX

A. Derivation of N-Regu

We consider the following optimization problem: F is a constant and does not change the minimizer. The reason of adding this term is purely algorithmic. We then solve the optimization [START_REF] Hale | Fixed-point continuation for 1minimization: Methodology and convergence[END_REF] with BC-VMFB [START_REF] Chouzenoux | A block coordinate variable metric forward-backward algorithm[END_REF].

min Ãf ,α f 1 2 Xf -Ãf α f 2 F + µ 2 Ãf 2 F + λ α f 1 + ı C ( Ãf ). ( 26
Let the general optimization min

x,y

F (x) + Q(x, y) + G(y) , (27)
where F (x) and G(y) are lower semicontinuous functions, Q(x, y) is a smooth function with Lipschitz gradient on any bounded set. BC-VMFB uses the following update rules to solve (27):

x (j+1) = argmin x F (x) + x -x (j) , ∇ x Q(x (j) , y (j) ) + t 1,(j) 2 x -x (j) 2 U 2,(j) , (28)
y (j+1) = argmin y G(y) + y -y (j) , ∇ y Q(x (j+1) , y (j) ) + t 2,(j) 2 y -y (j) 2 U 2,j , (29)
where U 1,(j) and U 2,(j) are positive definite matrices. x 2 U denotes the variable metric norm:

x 2 U = x, Ux . (30)
With the variable metric norm, the proximal operator (8) can be generalized as:

prox U,ψ := argmin y ψ(y) + 1 2 y -x 2 U . (31)
Then ( 28) and (29) can be rewritten as follow:

x (j+1) =prox U 1,(j) ,F/t 1,(j) x (j) - 1 t 1,(j) U 1,(j) -1 ∇ x Q(x (j) , y (j) ) , (32)
y (j+1) =prox U 2,(j) ,G/t 2,(j) y (j) - 1 t 2,(j) U 2,(j) -1 ∇ y Q(x (j+1) , y (j) ) . (33)
It is shown in [START_REF] Chouzenoux | A block coordinate variable metric forward-backward algorithm[END_REF] that the sequence generated by the above update rules converges to a critical point of the problem (27).

For the problem (26), we make the following substitutions:

F (α f ) = λ α f 1 , Q(α f , Ãf ) = 1 2 Xf -Ãf α f 2 F + µ 2 Ãf 2 F , G( Ãf ) = ı C ( Ãf ), (34)
Denoting by L

(j) the Lipschitz constant of ∇ α f Q(α (j) f , Ã (j)
f ), we have chosen:

U 1,(j) = L (j) I, U 2,(j) = ∂Q( Ãf , α (j+1) f ) 2 ∂ 2 Ãf = α (j+1) f α (j+1) H f
+ µI, t 1,(j) = t 2,(j) = 1.

(35)

The update step of the mixing matrix can be written as: (36)

Ã(j+1/2) f = Xf α (j+1) H f (α (j+1) f α (j+1) H f + µI) -1 , Ã (j+1)
As the choice of the parameter µ does not change the minimizer of [START_REF] Hale | Fixed-point continuation for 1minimization: Methodology and convergence[END_REF], by choosing µ sufficiently large, the update step of Ãf becomes:

Ã(j+1/2) f = P C Xf α (j+1) H f . ( 37
)
We obtain the N-Regu as shown in Algorithm 1.

B. Convolutive mixing operator and its adjoint operators

Given a signal s ∈ C T , and a convolution kernel h ∈ C L , the convolution can be written under the matrix form:

x = Hs = Sh , (38)
H ∈ C T ×T and S ∈ C T ×L being the corresponding circulant matrices.

The convolutive mixing operator can then be represented by

where s 1 , s 2 , . . . , s N ∈ C T are N source signals and x 1 , x 2 , . . . , x M ∈ C T are M observations. H mn is the convolution matrix from the n-th source to the m-th microphone. Thanks to these notations, the adjoint operator of convolutive mixing with respect to the mixing system is a linear operator C M ×T → C N ×T and can be represented by the following matrix multiplication:

     s 1 s 2 . . . s N      =    
          x 1 x 2 . . . x M      . (40)
In order to coincide with the previous notations in [START_REF] Balazs | Adapted and adaptive linear time-frequency representations: a synthesis point of view[END_REF], we denote the above formulation as:

S = H H X. (41)
The adjoint operator of the convolutive mixture with respect to the sources can then be written as:

H = X S H , (42)
with

h mn = S H n x m . (43)
C. Derivation of C-PALM

The PALM algorithm [START_REF] Bolte | Proximal alternating linearized minimization for nonconvex and nonsmooth problems[END_REF] is designed to solve the nonconvex optimization problem in the general form (27) by the following update rules:

x (j+1) = argmin x F (x) + xx (j) , ∇ x Q(x (j) , y (j) ) + t 1,(j) 2

xx (j) 2 2 , (44) y (j+1) = argmin y G(y) + yy (j) , ∇ y Q(x (j+1) , y (j) ) + t 2,(j) 2 yy (j) 2 2 , (45

)
where j is the iteration index and t 1,(j) et t 2,(j) are two step parameters.

It is shown in [START_REF] Bolte | Proximal alternating linearized minimization for nonconvex and nonsmooth problems[END_REF] that the sequence generated by the above update rules converges to a critical point of the problem (27).

From the general optimization (27), we do the following substitutions:

F (α f ) = λ α f 1 , Q(α f , Hf ) = 1 2 Xf -Hf α f 2 F + P( Hf ), G( Hf ) = ı Conv C ( Hf ), (46)
and the particular choices: t 1,(j) = L 1,(j) , t 2,(j) = L 2,(j) , where L 1,(j) and L 2,(j) are respectively the Lipschitz constant of ∇ α f Q(α (47)

Let Ψ n denotes the circulant matrix associated with α n . If the synthesis coefficients of different sources are independent, we have E [Ψ i Ψ j ] = 0, for i = j.

Then, using similar notations as in Appendix VI-B, one can write Ĥ as:

ĥmn = Ψ H n Ψ n hm , (48)
Finally, Proposition 4 comes from the definition of the Lipschitz constant.

Fig. 1 .

 1 Fig. 1. Room configuration for synthesized mixtures

Fig. 3 .

 3 Fig. 3. Separation performance of different algorithms with oracle permutation alignment as a function of the reverberation time RT 60 in noiseless case

Fig. 4 .Fig. 5 .

 45 Fig. 4. performance of different algorithms as a function of the sparsity level in noiseless case. RT 60 = 130 ms

Fig. 6 .Fig. 7 .

 67 Fig. 6. Separation performance of different algorithms with oracle permutation alignment as a function of the input SNR for RT 60 = 130 ms

Fig. 8 .

 8 Fig. 8. Separation performance of different algorithms as a function of the reverberation time RT 60 with input SNR=15 dB

4 Fig. 9 .

 49 Fig. 9. Different settings of source positions for synthesized mixtures without input noise

Fig. 10 .

 10 Fig. 10. Separation performance of different algorithms for different sources positions in noiseless case. RT 60 = 130 ms.

  ) This optimization is equivalent to the problem (6): the indicator function ı C ( Ãf ) forces the normalization on each column of Ãf , therefore the term µ 2 Ãf 2

f∈

  prox U 2,(j) ,ı C ( Ã(j+1/2) f ).

  We obtain the C-PALM algorithm presented in Algorithm 2.D. Calculation of the Lipschitz constant in C-PALMWe present the calculation of the Lipschitz constant of the functionI( Hf ) := Hf α f α H f + ∇ Hf P( Hf ) = Ĥf + ∇ Hf P( Hf ).

TABLE I

 I

	EXPERIMENTAL CONDITIONS
		synthesized	live recorded
	Number of microphones	M = 2	M = 2
	Number of sources	N = 3	N = 3, 4
	Duration of signals	6 s	10 s
	Reverberation time (RT 60 )	50, 130, 250, 400 ms	130, 250 ms
	Sample rate	11 KHz	16 kHz
	Microphone distance	4 cm	5 cm, 1 m
	STFT window type	Hann	Hann
	STFT window length	512 (46.5 ms)	2048 (128 ms)
	STFT window shift	256 (23.3 ms)	512 (32 ms)

  Fig. 2. Separation performance as a function of the reverberation time RT 60 in noiseless case

		10									15
	SDR (dB)	8 4 6				N-Regu with TDOA permutation N-Regu with correlation permutation C-PALM with correlation permutation Bin-wise Full-rank	SIR (dB)	10 5
		2									0
		50 0	100	150	200	250	300	350	400		50 -5	100	150	200	250	300	350	400
					RT 60 (ms)								RT 60 (ms)
		15									14
											12
	ISR (dB)	5 10								SAR (dB)	6 8 10
											4
		50 0	100	150	200	250	300	350	400		50 2	100	150	200	250	300	350	400
					RT 60 (ms)								RT 60 (ms)
		10									15
		8				N-Regu with oracle permutation C-PALM with oracle permutation
	SDR (dB)	4 6				Bin-wise with oracle permutation Full-rank with oracle permutation	SIR (dB)	5 10
		2
		50 0	100	150	200	250	300	350	400		50 0	100	150	200	250	300	350	400
					RT 60 (ms)								RT 60 (ms)
		16									14
		14									12
	ISR (dB)	8 10 12								SAR (dB)	6 8 10
		6									4
		50 4	100	150	200	250	300	350	400		50 2	100	150	200	250	300	350	400
					RT 60 (ms)								RT 60 (ms)

TABLE II SEPARATION

 II RESULTS OF C-PALM FOR LIVE RECORDED MIXTURES FROM SISEC2011 (SDR/SIR/ISR/SAR IN DB)

		RT 60 = 130 ms	RT 60 = 250 ms
	microphone space	5 cm	1 m	5 cm	1 m
	male3	7.65 / 11.38 / 12.10 / 10.65 7.53 / 11.27 / 11.77 / 10.58	5.20 / 7.67 / 9.01 / 8.62	4.98 / 10.62 / 6.67 / 7.04
	female3	6.69 / 9.81 / 10.90 / 10.52	9.77 / 14.49 / 14.13 / 13.02	5.29 / 9.16 / 7.77 / 8.75 7.34 / 11.22 / 10.97 / 11.02
	male4	3.25 / 4.65 / 6.09 / 6.01	2.34 / 2.15 / 5.16 / 5.47	2.10 / 1.79 / 4.63 / 5.49	3.08 / 4.22 / 6.00 / 6.11
	female4	2.36 / 2.05 / 5.37 / 6.53	3.66/ 6.05 / 6.80 / 7.15	2.39 / 2.20 / 5.27 / 6.51	3.12 / 4.51 / 6.07 / 6.84
	x-axis (m)

  H 11 H 12 . . . H 1N H 21 H 22 . . . H 2N

		x 1										s 1	
	   	x 2 . . .	   	=	   	. . .	. . .	. . .	. . .	   	   	s 2 . . .	   	,
	x M										s N

H M 1 H M 2 . . . H M N

In this paper, the sparsity level is the percentage of zero elements in a vector or matrix. A higher sparsity level means a sparser vector or matrix.