NONNEGATIVE MATRIX FACTORIZATION AND SPATIAL COVARIANCE MODEL FOR UNDER-DETERMINED REVERBERANT AUDIO SOURCE SEPARATION Simon Arberet1 , Alexey Ozerov2 , Ngoc Q.K. Duong2 , Emmanuel Vincent2 , R´emi Gribonval2 , Fr´ed´eric Bimbot3 , Pierre Vandergheynst1∗ 1

´ Signal Processing Institute, Ecole Polytechnique F´ed´erale de Lausanne (EPFL), Switzerland {simon.arberet, pierre.vandergheynst}@epfl.ch 2

INRIA, Centre de Rennes - Bretagne Atlantique, France

{alexey.ozerov, Quang-Khanh-Ngoc.Duong, emmanuel.vincent, remi.gribonval}@inria.fr 3

CNRS, IRISA - UMR 6074, France [email protected]

ABSTRACT We address the problem of blind audio source separation in the under-determined and convolutive case. The contribution of each source to the mixture channels in the time-frequency domain is modeled by a zero-mean Gaussian random vector with a full rank covariance matrix composed of two terms: a variance which represents the spectral properties of the source and which is modeled by a nonnegative matrix factorization (NMF) model and another full rank covariance matrix which encodes the spatial properties of the source contribution in the mixture. We address the estimation of these parameters by maximizing the likelihood of the mixture using an expectation-maximization (EM) algorithm. Theoretical propositions are corroborated by experimental studies on stereo reverberant music mixtures. 1. INTRODUCTION In blind source separation (BSS), we observe a multichannel signal x(t) ∈ RM which is a mixture of N source signals sn (t) ∈ R, 1 ≤ n ≤ N . In the convolutive BSS case, each source sn (t) is convolved with M filters hn (l) ∈ RM which model in the audio context the acoustic paths from source n to the M microphones. The mixture process can be expressed as: x(t) =

N X

yn (t) + n(t)

(1)

n=1

where n(t) ∈ RM is an additive noise and yn (t) ∈ RM is the spatial image of source n which is expressed as: yn (t) =

L−1 X

hn (l)sn (t − l)

(2)

l=0

where L is the filter length. The BSS problem consists in recovering either the source signals sn (t) or their spatial images yn (t) given the mixture signal x(t). In this paper we consider the later BSS problem formulation. When the number N of sources is larger than the number M of mixture channels, the mixture is said under-determined. The BSS problem is often addressed in the time-frequency (TF) domain via the short time Fourier transform (STFT), and ∗ This work was supported in part by the SMALL project and the Quaero Programme, funded by OSEO.

the convolutive mixing process is approximated by a complexvalued instantaneous mixing in each frequency bin. In other words each source spatial image in the STFT domain Yn (t, f ) is approximated by the following complex-valued multiplication: b n (f )Sn (t, f ) Yn (t, f ) ≈ h

(3) b where hn (f ) is the Fourier transform of the mixing filters hn (t) and Sn (t, f ) and Yn (t, f ) are respectively the STFT of sn (t) and yn (t). Thus, according to the model (3), if Sn (t, f ) is a zero-mean random variable with variance vn (t, f ), the covariance of Yn (t, f ) is given by: (4) Ryn (t, f ) = vn (t, f )Rn (f ) H b b where Rn (f ) = hn (f )hn (f ) (H denotes the matrix conjugate transposition) is a rank-1 matrix. This rank-1 model holds only when the filter length L is short compared to the STFT window size [1]. A particular case is when the mixture is instantaneous, i.e. the filters hn (l) have length L = 1, then approximation (3) becomes an equality. However in an environment with realistic reverberation time, the filter length L is usually longer than the STFT window size. Assuming the rank-1 model (3), BSS can be achieved by estimating a mixing matrix [2, 3, 4] in each frequency bin f (whose b n (f ), 1 ≤ n ≤ N ) and then recovercolumns are the vectors h ing the source coefficients Sn (t, f ) assuming independence of the sources and some sparse prior distributions [5]. However, if the mixing matrices are estimated independently in each freb n (f ) of these mixing matrices are quency bin, the columns h arbitrary permuted in each frequency, leading to the well-known permutation problem. Recently, some spectral approaches using the rank-1 model (3) and modeling the structure of the source variances vn (t, f ) in the TF plane, with a Gaussian mixture model (GMM) [6] or nonnegative matrix factorisation (NMF) [7, 8] have been proposed. These spectral approaches have shown to provide better performance [9] than classical sparse approaches like binary masking [2], l1 -norm [10] or lp -norm [11] minimization. In the case of the NMF approach of [8], as there is a coupling of the frequency bins due to the structure of vn (t, f ), and b n (f ) are jointly estimated, we are able to avoid as vn (t, f ) and h the permutation problem. To model reverberation efficiently, Duong et al.[12] proposed recently to consider Rn (f ) as a full-rank (unconstrainted) matrix. They showed that this model led to better results than the rank-1 model on reverberant mixtures in oracle context where

Rn (f ) and vn (t, f ) are known, in semi-blind context where Rn (f ) is known but vn (t, f ) is a free variance estimated from the mixture. They also formulated an expectation-maximization (EM) algorithm [13] to blindly estimate Rn (f ) and vn (t, f ) in each frequency bin. However as the parameters Rn (f ) and vn (t, f ) are estimated independently in each frequency bin, the permutation problem has to be solved a posteriori. Duong et al. [13] applied a DOA-based algorithm to solve the permutation problem. However, in order to deploy this DOA-based algorithm, it is imperative to know the inter-microphone distance beforehand. Motivated by the effectiveness of the full rank spatial model of Duong et al.[13] and the NMF spectral model [7, 8], we investigate in this paper the modeling of each spatial source image with a combination of these two models. We describe the proposed source spatial image model in Section 2. Section 3 addresses the proposed inference method which consists of maximizing the likelihood of the mixture data using an EM algorithm [14]. In Section 4, we compare the source separation performance achieved by our full-rank NMF method with the rank-1 NMF method [8] and with other state-of-the-art algorithms over stereo music data. Finally, we conclude in Section 5.

Thus the noise can be seen as a particular source with only one component which is time-invariant. In other words relation (4) applied to the noise “source” becomes Rb (t, f ) = Rb (f ). 2.3. Mixture Since the STFT is a linear transform, the mixing process (1) can be rewritten as: X X(t, f ) = Yn (t, f ) + N(t, f ). n

The sources and noise are assumed to be independent of each other. Thus, the model of the mixture STFT X(t, f ) is a zero-mean Gaussian vector with covariance matrix: X (9) Rx (t, f ) = Ryn (t, f ) + Rb (f ). n

3. ESTIMATION OF THE MODEL PARAMETERS We wish to estimate in the maximum likelihood (ML) sense, the mixing parameters Rn (f ) of each source, the source variances vn (t, f ) under the constraint given by (6) and the noise covariance Rb (f ).

2. MODEL 3.1. Criterion and indeterminacies 2.1. Source spatial image We assume that each spatial source image Yn (t, f ) at TF point (t, f ) is a zero-mean complex random vector: Yn (t, f ) ∼ Nc (0, Ryn (t, f )),

(5)

where Nc (µ, Σ) is a proper complex distribution with probability density function (pdf): h i 4 Nc (Y; µ, Σ) = |π Σ|−1 exp − (Y − µ)H Σ−1 (Y − µ) , where |A| denotes the determinant of a square matrix A. µ and Σ are respectively, the M -dimensional mean vector and the M × M covariance matrix of Y. Covariance matrix Ryn (t, f ) is given by (4), where Rn (f ) is a full-rank unconstrained timeinvariant covariance matrix which encodes the spatial properties of the source [13] and vn (t, f ) is a time-varying source variance which is an assumed sum of K components: vn (t, f ) =

K X

n wf,k hn k,t

(6)

k=1 n + where wf,k , hn k,t ∈ R . Thus, the power spectrum Vn = [vn (t, f )]f,t of each source n is structured as a product of two n nonnegative matrices Wn = [wf,k ]f,k and Hn = [hn k,t ]k,t :

Vn = Wn Hn . According to (4), (5), (6), each source spatial image Yn (t, f ) can be seen P as a sum of K independant zero-mean Gaussians Yn (t, f ) = K k=1 Yn,k (t, f ) with the respective covariances: n Ryn,k (t, f ) = wf,k hn k,t Rn (f ).

Let Θ = {Rn (f ), Wn , Hn , Rb (f ), ∀f, n} be the set of all the parameters we wish to estimate. As P (X(t, f )|Θ) is a zeromean Gaussian P according to section 2.3, maximizing the loglikelihood t,f log P (X(t, f )|Θ) is equivalent to minimizing the cost: X C(Θ) = X(t, f )H R−1 x (t, f )X(t, f ) + log |Rx (t, f )| t,f

where Rx (t, f ) is defined in (9). Thus, the ML criterion suffers from scaling indeterminacies because for any α ∈ R+,∗ : `1 ´ Ryn (t, f ) = (αvn (t, f )) α Rn (f ) , and also for any αk ∈ “ ” P n ) α1k hn R+,∗ : vn (t, f ) = k (αk wf,k k,t . In order to remove these scaling indeterminacies, we normalize Rn (f ) according to the Frobenius norm kRn (f )kF =P 1 (and scale vn (t, f )) acn = 1 (and scaling cordingly) and impose the condition, f wf,k accordingly) as in [8]. hn k,t 3.2. Algorithm We derive an EM algorithm [14] based on the complete data {Yn,k (t, f ), N(t, f ) ∀t, f, n, k}, that is the set of the STFT coefficients of all the source spatial image components and the noise. Each iteration of the EM algorithm is composed of two steps: the E-step and the M-step. The E-step consists of computb yn (t, f ), R b y (t, f ), ing the expectation of the natural statistics R n,k b Rb (t, f ), that is, the covariances of Yn (t, f ), Yn,k (t, f ) and N(t, f ), conditionally on the mixture data and the current parameter estimates Θ. The M-step consists in re-estimating the parameters Θ using the updated natural statistics. 3.2.1. E-step: Conditional expectation of natural statistics

(7) b yn (t, f ) = Y b n (t, f )Y b nH (t, f ) + (I − Gn (t, f )) Ryn (t, f ) R (10)

2.2. Noise Let N(t, f ) be the STFT of n(t). We assume that the noise is a stationary zero-mean Gaussian process with covariance Rb (f ): N(t, f ) ∼ Nc (0, Rb (f )).

(8)

H b y (t, f ) = Y b n,k (t, f )Y b n,k R (t, f ) n,k

+ (I − Gn,k (t, f )) Ryn,k (t, f ) bH

b b (t, f ) = N(t, b f )N (t, f ) + (I − Gb (t, f )) Rb (f ) R

(11)

4. EXPERIMENTAL RESULTS

where b n (t, f ) = Gn (t, f )X(t, f ) Y

(12)

b n,k (t, f ) = Gn,k (t, f )X(t, f ) Y

4.1.1. Datasets

b f ) = Gb (t, f )X(t, f ) N(t, Gn (t, f ) = Ryn (t, f ) (Rx (t, f ))−1 Gn,k (t, f ) = Ryn,k (t, f ) (Rx (t, f ))−1 Gb (t, f ) = Rb (f ) (Rx (t, f ))−1 . 3.2.2. M-step: Update of the parameters The re-estimation of Rn (f ) in the ML sense is equivalent to minimizing the sum over all the ˛ t of the ”Kullback“ time frames b yn (t, f )˛˛ Ryn (t, f ) , with reLeibler (KL) divergence DKL R spect to (w.r.t.) Rn (f ), between two zero-mean Gaussian disb yn (t, f ) and Ryn (t, f ) tributions with covariances matrices R defined in (10) and (4) respectively and with: DKL ( R1 | R2 ) =

4.1. Setting

´ ` ´ ´ 1` ` tr R1 R−1 − log det R1 R−1 −M . 2 2 2

Given Wn and Hn , this minimization has a closed-form representation which is [13]

We evaluated our NMF full rank algorithm with the NMF rank-1 algorithm of Ozerov and F´evotte [8] and the rank-1 and full rank EM algorithms of Duong et al [13], over stereo noiseless music mixtures under various mixing conditions. For each experiment 10 mixtures were generated by convolving different 10 seconds source signals sampled at 16 kHz with room impulse responses obtained with the Roomsim toolbox [15]. The microphones are omnidirectional and the room dimensions are 4.45 m x 3.55 m x 2.5 m. The number of sources is set to either 3 or 4, the reverberation time (RT) to either 130 ms or 250 ms and the distance between the two microphones to either 5 cm or 1 m, resulting in 8 mixing conditions overall. 4.1.2. Algorithms setting and evaluation criterion The STFT was computed with a sine window of length 1024. The number of components per source of the NMF models was set to K = 5, the number of iterations for each EM algorithm was 50. Separation performance was evaluated using the signalto-distortion ratio (SDR) criterion averaged over all sources. 4.1.3. Initialization

T 1 X 1 b y (t, f ) Rn (f ) = R T t=1 vn (t, f ) n

(13)

P n n where vn (t, f ) = K k=1 wf,k hk,t as defined in (6). n n The re-estimation of wf,k and “ hk,t in the˛ML sense is equiv” P b y (t, f )˛˛ Ry (t, f ) w.r.t. alent to minimizing DKL R t,f

n,k

n,k

n b wf,k and hn k,t , where Ryn,k (t, f ) and Ryn,k (t, f ) are defined in (11) and (7) respectively. Given Rn (f ), these minimizations have closed-form representations which are:

n wf,k =

T 1 X vˆn,k (t, f ) , T t=1 hn k,t

hn k,t =

F 1 X vˆn,k (t, f ) n F wf,k

As the EM algorithm is very sensitive to the initialization and to be sure to have a “good initialization”, we provide it with perturbed oracle initializations, where the parameters Rn (f ) and vn (t, f ) are estimated form the original source spatial images as in [12] and then perturbed with a high level additive noise (SNR n , hn of 3 dB) as in [8]. Parameters wf,k k,t of the NMF approaches are then computed with NMF decomposition using multiplicative update (MU) rules and KL divergence as in [8]. For the rank-1 methods (i.e. binary masking, Duong et al. rank-1 and b n (f ) by calculating the first princiNMF rank-1) we compute h pal component of Rn (f ) using the principal component analysis (PCA).

(14)

f =1

4.1.4. Noise

Assuming that the noise is spatially uncorrelated, we set the off-diagonal coefficients of Rb (f ) to zero.

When the noise tends to zero, the estimation of the mixture parameters using the NMF rank-1 algorithm gets stuck [8] and when the noise is small, the convergence of this EM algorithm is very slow. Thus, the authors of [8] proposed a strategy called noise annealing with noise injection where the noise covariance N(t, f ) = σb2 (f )I is initialized with a large value of σb2 (f ) and instead of being re-estimated at each iteration, is gradually decreased through iterations to a small value. Noise injection means that a random noise with a covariance N(t, f ) is added to X(t, f ) at each EM iteration. This technique accelerates the overall global convergence [8]. Although this stuck problem doesn’t hold in our full rank NMF algorithm, we used this noise annealing with noise injection scheme for both the NMF rank-1 algorithm and our NMF full rank algorithm1 .

3.3. Estimation of the sources

4.2. Results

After the convergence of the EM algorithm, the source spatial images are estimated in the TF domain with the Wiener estimator as in (12). The estimation of source spatial images in the time domain are then obtained via inversion of the STFT map using the overlappadd technique.

The results corresponding to reverberation times of 130 ms and 250 ms are respectively shown in Table 1 and Table 2.

with:

” 1 “ −1 b y (t, f ) vˆn,k (t, f ) = tr Rn (f )R (15) n,k M where tr(.) denotes the trace of a square matrix. To remove the scaling ambiguity between the components n wf,k and hn k,t , we normalize them as explained in section 3.1. The process of re-estimating of the noise covariance involves a similar set of steps as in (13) and is given by: Rb (f ) =

T 1 Xb Rb (t, f ). T t=1

1 We also noticed that there is a marginal increase of the performance (between 0 and 0.5 dB of the SDR) of NMF full rank when using this noise annealing with noise injection scheme.

Table 1. Source separation performance, RT = 130 ms Reverberation Time Microphone distance Number of sources

130 ms 5 cm 3

1m

4

Appoaches

3

4

SDR in dB

6. REFERENCES

Binary masking

0.4

0.0

3.6

0.9

Duong et al. rank-1

1.2

0.2

2.2

0.8

Duong et al. full rank

9.5

8.3

9.8

8.3

NMF rank-1

8.7

7.1

7.3

4.6

NMF full rank

9.1

7.5

10.2

8.5

Table 2. Source separation performance, RT = 250 ms Reverberation Time

250 ms

Microphone distance Number of sources

5 cm 3

Appoaches

4

1m 3

GMM) to the full rank model and validation of the proposed method over real-word recordings with more than four sources. As the EM algorithm is sensitive to parameter initialization, it is important to investigate blind initialization procedures of the model parameters and particularly, the initialization of the spatial covariance matrices.

4

SDR in dB

Binary masking

0.8

-0.3

2.9

0.7

Duong et al. rank-1

0.1

-0.1

1.5

0.5

Duong et al. full rank

8.1

7.2

8.1

7.1

NMF rank-1

7.7

6.4

5.5

3.8

NMF full rank

8.8

7.5

9.6

8.0

Unsurprisingly, when the number of sources increases as well as when the reverberation time increases, the performance of all the tested algorithms degrades. NMF full rank outperforms NMF rank-1 and NMF rank-1 outperforms Duong et al. rank-1 by between 3 dB and 7 dB. In the “low” reverberant setting (RT = 130 ms), Duong et al. full rank performs better than NMF full rank when the microphones distance is 5 cm, but less than NMF full rank when the microphone distance is 1 m. When the reverberation time is longer (RT = 250 ms), the NMF full rank outperforms all the other tested methods. Thus it shows that combining the full rank spatial covariance model with the NMF spectral model improves the separation in realistic reverberant environment. 5. CONCLUSION In this paper we have introduced a new model for convolutive blind source separation that combines the advantages of the two existing models. The source spectrum is modeled via nonnegative matrix factorization (NMF) and the convolutive mixing process is modeled using a full rank spatial covariance instead of a rank-1. We addressed the estimation of the model parameters by maximizing the likelihood of the observed mixture using an EM algorithm. Experimental results over music data in different settings (number of sources, microphone distance, reverberation time) validate that our model outperfoms the NMF rank-1 approach [8], the full rank method of Duong et al. [13] and binary masking, when the reverberation time is realistic (RT of 250 ms). Future works include the extension of other spectral models (e.g.

[1] L. Parra and C. Spence, “Convolutive blind separation of non-stationary sources,” IEEE Transactions on Speech and Audio Processing, vol. 8, no. 3, pp. 320–327, 2000. [2] O. Yilmaz and S. Rickard, “Blind separation of speech mixtures via time-frequency masking,” IEEE Transactions on Signal Processing, vol. 52, no. 7, pp. 1830–1847, July 2004. [3] S. Arberet, R. Gribonval, and F. Bimbot, “A robust method to count and locate audio sources in a multichannel underdetermined mixture,” IEEE Transactions on Signal Processing, vol. 58, no. 1, January 2010. [4] S. Winter, W. Kellermann, H. Sawada, and S. Makino, “Map-based underdetermined blind source separation of convolutive mixtures by hierarchical clustering and l 1norm minimization,” EURASIP Journal on Advances in Signal Processing, vol. 2007, pp. 1–12, 2007. [5] B. Pearlmutter P. O’Grady and S. Rickard, “Survey of sparse and non-sparse methods in source separation,” IJIST, March 2005. [6] S. Arberet, A. Ozerov, R. Gribonval, and F. Bimbot, “Blind spectral-GMM estimation for underdetermined instantaneous audio source separation,” in ICA, 2009. [7] C. F´evotte, N. Bertin, and J.L. Durrieu, “Nonnegative matrix factorization with the itakura-saito divergence: With application to music analysis,” Neural Computation, vol. 21, no. 3, pp. 793–830, 2009. [8] A. Ozerov and C. F´evotte, “Multichannel nonnegative matrix factorization in convolutive mixtures for audio source separation,” IEEE Transactions on Audio, Speech and Language Processing, vol. 18, no. 3, pp. 550–563, 2010. [9] E. Vincent, S. Araki, and P. Bofill, “The 2008 signal separation evaluation campaign: A community-based approach to large-scale evaluation,” in ICA, 2009, pp. 734–741. [10] P. Bofill and M. Zibulevsky, “Underdetermined blind source separation using sparse representations,” Signal Processing, vol. 81, no. 11, 2001. [11] E. Vincent, “Complex nonconvex lp norm minimisation for underdetermined source separation,” in ICA, 2007. [12] N.Q.K. Duong, E. Vincent, and R. Gribonval, “Spatial covariance models for under-determined reverberant audio source separation,” in WASPAA, 2009. [13] N.Q.K. Duong, E. Vincent, and R. Gribonval, “Underdetermined convolutive blind source separation using spatial covariance models,” in ICASSP, 2010. [14] A.P. Dempster, N.M. Laird, D.B. Rubin, et al., “Maximum likelihood from incomplete data via the EM algorithm,” Journal of the Royal Statistical Society. Series B (Methodological), vol. 39, no. 1, pp. 1–38, 1977. [15] D. Campbell, K. Palom¨aki, and G. Brown, “A matlab simulation of ’shoebox’ room acoustics for use in research and teaching,” Computing and Information Systems Journal, vol. 9, no. 3, pp. 48–51, October 2005.