Expectation-Maximization for Speech Source Separation Using Convolutive Transfer Function

TL;DR

This paper introduces an EM-based method for speech source separation in reverberant environments using convolutive transfer functions in the STFT domain, improving accuracy over narrowband models.

Contribution

It proposes a novel joint estimation approach of CTF coefficients and source signals using EM, enhancing separation performance in reverberant conditions.

Findings

Effective in highly reverberant environments

Outperforms narrowband models

Provides accurate room filter estimation

Abstract

This paper addresses the problem of under-determinded speech source separation from multichannel microphone singals, i.e. the convolutive mixtures of multiple sources. The time-domain signals are first transformed to the short-time Fourier transform (STFT) domain. To represent the room filters in the STFT domain, instead of the widely-used narrowband assumption, we propose to use a more accurate model, i.e. the convolutive transfer function (CTF). At each frequency band, the CTF coefficients of the mixing filters and the STFT coefficients of the sources are jointly estimated by maximizing the likelihood of the microphone signals, which is resolved by an Expectation-Maximization (EM) algorithm. Experiments show that the proposed method provides very satisfactory performance under highly reverberant environments.

Tables

Table 1. TABLE I: Real-time factor of the four methods.

BM	FR-SCM	CTF-Lasso	Proposed
0.01	81.7	54.2	630.0

Equations

$$y(n)=a(n)\star s(n),$$

$$y(p,k)=a(p,k)\star s(p,k)=\sum _{p^{\prime }=0}^{Q}a(p^{\prime },k)s(p-p^{\prime },k),$$

$$a(p,k)=(a(n)\star {\zeta }_k(n)){|}_{n=pL}.$$

$${\zeta }_k(n)=e^{j\frac{2\pi }{N}kn}\sum _{m=-\infty }^{+\infty }{\tilde{\omega }}_a(m){\tilde{\omega }}_s(n+m),$$

$$x_i(p,k)=\sum _{j=1}^J{a}_{ij}(p,k)\star s_j(p,k)+e_i(p,k),$$

$$\displaystyle\mathcal{N}_{c}(s_{j}(p,k);0,v_{j}(p,k))=\frac{1}{\pi v_{j}(p,k)}\exp\big{(}\!-\!\frac{|s_{j}(p,k)|^{2}}{v_{j}(p,k)}\big{)}.$$

$$\displaystyle\mathcal{N}_{c}(\mbox{\boldmath$e$}(p,k);0,\mathbf{\Sigma_{e}}(k))=\frac{1}{\pi^{I}|\mathbf{\Sigma_{e}}(k)|}e^{-\mbox{\boldmath$e$}(p,k)^{H}\mathbf{\Sigma_{e}}(k)^{-1}\mbox{\boldmath$e$}(p,k)},$$

$$\displaystyle\mbox{\boldmath$s$}_{j}(p)$$

$$\displaystyle\mbox{\boldmath$s$}(p)$$

$${\mathbf{a}}_{ij}$$

$${\mathbf{A}}_j$$

$$\mathbf{A}$$

$$\displaystyle\mbox{\boldmath$x$}(p)=\sum_{j=1}^{J}\mathbf{A}_{j}\mbox{\boldmath$s$}_{j}(p)+\mbox{\boldmath$e$}(p)=\mathbf{A}\mbox{\boldmath$s$}(p)+\mbox{\boldmath$e$}(p).$$

$$\displaystyle\tilde{\mbox{\boldmath$s$}}_{j}$$

$$\displaystyle\tilde{\mbox{\boldmath$e$}}_{i}$$

$$\displaystyle\tilde{\mbox{\boldmath$x$}}_{i}$$

$$\displaystyle\tilde{\mbox{\boldmath$s$}}$$

$$\displaystyle\tilde{\mbox{\boldmath$e$}}$$

$$\displaystyle\tilde{\mbox{\boldmath$x$}}$$

$${\mathcal{A}}_{ij}=\left[\begin{array}{cccccc}{a}_{ij}(0)& 0& \cdots & \cdots & \cdots & 0\\ ⋮& \ddots & \ddots & \ddots & \ddots & ⋮\\ a_{ij}(Q)& \ddots & a_{ij}(0)& 0& \ddots & 0\\ ⋮& \ddots & \ddots & \ddots & \ddots & ⋮\\ 0& \cdots & 0& a_{ij}(Q)& \cdots & a_{ij}(0)\end{array}\right],$$

$${\mathcal{A}}_i$$

$$\mathcal{A}$$

$$\displaystyle\tilde{\mbox{\boldmath$x$}}_{i}$$

$$\displaystyle\text{or}\qquad\tilde{\mbox{\boldmath$x$}}$$

$${\mathbf{\Psi }}_{\mathbf{s}}=\text{Diag}\{[v_1(1),\dots ,v_1(P),\dots ,v_J(1),\dots ,v_J(P){]}^{\top }\}$$

$${\mathbf{\Psi }}_{\mathbf{e}}((i_1-1)P+p_1,(i_2-1)P+p_2)=\{\begin{array}{ll}{\mathbf{\Sigma }}_{\mathbf{e}}(i_1,i_2),& \text{mbox}if{p}_1=p_2,\\ 0,& \text{mbox}otherwise.\end{array}$$

$$\displaystyle-(\tilde{\mbox{\boldmath$x$}}-\mathcal{A}\tilde{\mathbf{s}})^{H}\mathbf{\Psi_{e}}^{-1}(\tilde{\mbox{\boldmath$x$}}-\mathcal{A}\tilde{\mathbf{s}})-\tilde{\mbox{\boldmath$s$}}^{H}\mathbf{\Psi_{s}}^{-1}\tilde{\mbox{\boldmath$s$}},$$

$${\hat{\mathbf{\Sigma }}}_{\mathbf{s}}$$

$$\displaystyle\widehat{\tilde{\mbox{\boldmath$s$}}}$$

$$\displaystyle\mathbf{\widehat{\tilde{R}}_{s}}=\widehat{\tilde{\mbox{\boldmath$s$}}}\widehat{\tilde{\mbox{\boldmath$s$}}}^{H}+\mathbf{\widehat{\Sigma}_{s}}.$$

$$p(\text{mbox}\mathbf{X},$$

$$\displaystyle=\prod_{p=1}^{P}\mathcal{N}_{c}(\mbox{\boldmath$x$}(p);\mathbf{A}\mbox{\boldmath$s$}(p),\mathbf{\Sigma_{e}})\prod_{j=1}^{J}\prod_{p=1}^{P}\mathcal{N}_{c}(s_{j}(p);0,v_{j}(p)).$$

$$\displaystyle\log(p(\mbox{\bf X},\mbox{\bf S}|\Theta))=-\sum_{p=1}^{P}\Big{(}\text{log}(|\mathbf{\Sigma_{e}}|)+(\mbox{\boldmath$x$}(p)-\mathbf{A}\mbox{\boldmath$s$}(p))^{H}\mathbf{\Sigma_{e}}^{-1}$$

$$\displaystyle(\mbox{\boldmath$x$}(p)-\mathbf{A}\mbox{\boldmath$s$}(p))\Big{)}-\sum_{j=1}^{J}\sum_{p=1}^{P}\Big{(}\log(v_{j}(p))+\frac{|s_{j}(p)|^{2}}{v_{j}(p)}\Big{)}+\text{const}.$$

$$Q(\Theta ,{\Theta }^{\text{old}})=$$

$$\displaystyle\qquad\qquad+\mbox{\boldmath$x$}(p)\widehat{\mbox{\boldmath$s$}}(p)^{H}\mathbf{A}^{H}-\mathbf{A}{\mathbf{\widehat{R}_{s}}(p)}\mathbf{A}^{H})\}\Big{)}$$

$$\displaystyle-\sum_{j=1}^{J}\sum_{p=1}^{P}\Big{(}\log(v_{j}(p))+\frac{\widehat{v}_{j}(p)}{v_{j}(p)}\Big{)}+\text{const},$$

Full text

I Introduction

Most of speech source separation techniques are designed in the short time Fourier transform (STFT) domain where the narrowband assumption is generally used, e.g. [1-4]. In the narrowband assumption, at each frequency band, the time-domain filter is represented by the acoustic transfer function (ATF), and the time-domain convolutive process is transformed to a product between the ATF and the STFT coefficients of the source signal. This assumption is also referred to as the multiplicative transfer function (MTF) approximation [5]. Based on the ATF or its variant, e.g. relative transfer functions [6, 7], beamforming techniques are widely used for multichannel speech source separation and speech enhancement. Popular beamformers for multisource separation include linearly constrained minimum variance/power (LCMV/LCMP) [6, 8]. Furtherly, because of the spectral sparsity of speech, the microphone signal can be assumed to be dominated by only one speech source in each time-frequency (TF) bin. This is refered to as the W-disjoint orthogonality (WDO) assumption [1]. The binary masking method [1, 2] and $\ell_{1}$-norm minimization method [3] exploit such WDO assumption. More examples of narrowband assumption-based techniques can be found in [9] and references therein.

In real scenarios, the time-domain filter, i.e. the room impulse response (RIR), is normally much longer than the STFT window (frame) length, since the latter should be set to be sufficiently short to account for the local stationarity of speech. For this case, the narrowband assumption is no longer valid, and thus leads to unsatisfied speech source separation performance. In the literature, only a few studies had questioned the validity of the narrowband assumption and attempted to tackle this problem. Based on the narrowband assumption, the theoretical covariance matrix of one source image is a rank-one matrix [4]. To mitigate the invalidity of this matrix in practice, a full-rank spatial covariance matrix was adopted in [10], even if the narrowband assumption is used. To circumvent the inaccurancy of the narrowband assumption, the wideband time-domain convolution model was used in [11, 12], where the source STFT coefficients are recovered by minimizing the fit cost between the time-domain mixtures and the time-domain convolution model. Meanwhile, based on the Lasso technique, the $\ell_{1}$-norm of the source STFT coefficients is minimized as well to impose the sparsity of speech spectra. This method achieves good performance, but its computational complexity is very large due to the high cost of time-domain convolution operation. In [13, 14], based on the criterion of likelihood maximization, a variational EM algorithm was proposed also using the time-domain convolution model and the STFT domain signal model.

In [15], the time-domain convolution can be ideally represented in the STFT-domain by the cross-band filters. More precisely, the STFT coefficients of the source image can be computed by summing multiple convolutions (over frequencies) between the STFT coefficients of the source signal and the STTF-domain filter. Note that the convolution is conducted along the frame axis. To simplify the analysis, for each frequency, the band-to-band filter, i.e. the CTF model [16], is used, while with the cross-band information omitted. Compared to the narrowband assumption that uses a frequency-wise scalar product to approximate the time-domain convolution, the CTF model uses a frequency-wise convolution and thus is more accurate. Following the principle of the wide-band Lasso [11], based on the CTF model, a subband Lasso technique was proposed in [17], which largely reduces the complexity relative to the wide-band Lasso technique. In [18], two CTF inverse filtering methods were proposed based on the multiple-input/output inverse theorem (MINT). In [19], the CTF was integrated into the generalized sidelobe canceler beamformer. A CTF-based EM algorithm was proposed in [20] for single-source dereverberation, in which the Kalman filter was exploited to achieve online EM update. The cross-band filters were adopted in [21], combined with a non-negative matrix factorization model for the source signal. To estimate the source signals, the likelihood of the microphone signals is maximized via a variational EM algorithm. In [22], also based on likelihood maximization and EM, a STFT-domain convolutive model was used for source separation, combined with an HMM model for source activity. Even though this STFT-domain convolutive model was not named as CTF in [22], it actually plays the same role as CTF.

Due to the high model complexity, the above mentioned source separation techniques that are beyond the narrowband assumption actually can not be performed in a blind manner, in other words, some prior knowledge are required. For example, the RIRs for the wide-band Lasso techniques, and the CTFs for the CTF-Lasso and MINT techniques, are required to be knwon or well estimated. Both the mixing filters and source parameters are required as a good initialization for the (variational) EM techniques [13, 14, 21, 22]. A blind multichannel CTF identification method was proposed in [23] and the identified CTF can be fed in the semi-blind methods. However, this CTF identification method was only suitable for the single-source case.

In the present paper, based on the CTF model, we propose a likelihood maximization method for speech source separation. First, the CTF model is presented in a source mixture probabilistic framework. Then, an EM algorithm is proposed for resovling the likelihood maximization problem. The STFT coefficients of the source signals are taken as hidden variables, and are estimated in the expectation step (E-step). The CTF coefficients and source parameters are estimated in the maximization step (M-step). Experiments show that the proposed method performs better than the narrowband assumpation based methods [1, 10] and the CTF-lasso method [17] within a semi-blind setup where the mixing filters are initialized with a perturbed version of the ground-truth CTF.

The rest of this paper is organized as follows. Section II presents the CTF formulation, which is plugged in a probabilistic framework in Section III. The proposed EM is given in Section IV. Experiments are presented in Section V. Section VI concludes the paper. This paper is an extension of a conference paper [24]. The main improvements over [24] consists of i) we present the methodology in more detail, such as the two vector/matrix formulations in Section III, the detailed derivation of the EM algorithm in Section IV and the execution process of EM in Algorithm 1; ii) in Section V, we add the experiments with CTF perturbations, and analyze the computational complexity of the proposed method.

II CTF Formulation

In an indoor (room) environment, a speech source signal propagates to the receivers (microphones) through the room effect. In the time domain, the received source image ${y}(n)$ is given by the linear convolution between the speech source signal ${s}(n)$ and the RIR ${a}(n)$:

[TABLE]

where $\star$ denotes convolution. Applying STFT, the narrowband approximation usually approximates the time domain convolution (1) as the product $y(p,k)=a(k)s(p,k)$, where $y(p,k)$ and $s(p,k)$ are the STFT coefficients of the corresponding signals, and $a(k)$ is the ATF. Let $N$ denote the STFT window (frame) length, then $k\in[0,N-1]$ is the frequency bin index, and $p\in[1,P]$ is the frame index. This narrowband approximation is no longer valid when ${a}(n)$ is long relative to the STFT window.

We use the CTF model to circumvent the inaccurancy of the narrowband assumption, then $y(p,k)$ can be presented as [16]:

[TABLE]

where $a(p,k)$ denotes the CTF, which can be derived from the time-domain filter ${a}(n)$ by:

[TABLE]

This equation is the convolution with respect to the time index $n$ evaluated at multiples of the frame step $L$, with

[TABLE]

where $\tilde{\omega}_{a}(n)$ and $\tilde{\omega}_{s}(n)$ are respectively the STFT analysis and synthesis windows. The length of CTF, i.e. $Q+1$, approximately equals the length of RIR divided by $L$.

III Mixture Model Formulations

III-A Basic Formuation for Mixture Model

We consider a source separation problem with $J$ sources and $I$ sensors, which could be either underdetermined ($I<J$) or (over)determined ($I\geq J$). Using the CTF formulation (2), in the STFT domain, the microphone signal $x_{i}(p,k)$ is:

[TABLE]

where $a_{ij}(p,k)$ is the CTF from source $j,\ j=1,\dots,J$ to sensor $i,\ i=1,\dots,I$, and $e_{i}(p,k)$ denotes the noise signal.

III-B Probabilistic Model

In the literature of source separation, each source signal $s_{j}(p,k)$ is normally assumed to be independent to other sources, and is also independent across STFT frames and frequencies. Each STFT coefficient $s_{j}(p,k)$ is assumed to follow a complex Gaussian distribution with a zero mean and variance $v_{j}(p,k)$ [4, 10], i.e. its probability density function (pdf) is:

[TABLE]

The noise signal is assumed to be stationary, temporally uncorrelated, and independent to the speech source signals. We define the noise vector accross microphones as \mbox{\boldmathe}(p,k)=[e_{1}(p,k),\dots,e_{I}(p,k)]^{\top}\in\mathbb{C}^{I\times 1}. This vector is assumed to follow a zero-mean complex Gaussian, with a non-diagonal covariance matrix denoted as $\mathbf{\Sigma_{e}}(k)$. This convariance matrix encodes the spatial correlation of the noise signals. The pdf is:

[TABLE]

where H denotes complex transpose, $|\cdot|$ the determinant of matrix.

Since the proposed source separation method is carried out independently at each frequency, hereafter, we omit the frequency index $k$ for notational simplicity.

III-C Vector/matrix Formulation 1

To formulate the mixture model (4) more compactly, we have several different choices to organize the signals and the convolution operation in vector/matrix forms. To facilitate the derivation of the following EM algorithm, we use two different vector/matrix formulations. In this section, Formulation 1 will be presented to enable us to easily derive the M-step, and Formulatioin 2 will be presented used for the derivation of the E-step in the next section. This two formulations are different just in the organization of the variables and parameters, thence transforming from the E-step to the M-step, and vice-versa, will only necessitate reorganizing the vector/matrix elements.

In Formulation 1, we define the source signals in vector form as, for $p\in[1,P]$:

[TABLE]

where ⊤ denotes vector/matrix transpose. If $p\leq q$, we set $s_{j}(p-q)=0$. Define the CTF in vector/matrix form as:

[TABLE]

We already defined \mbox{\boldmathe}(p)\in\mathbb{C}^{I\times 1}. Similarly, the microphone signal is \mbox{\boldmathx}(p,k)=[x_{1}(p,k),\dots,x_{I}(p,k)]^{\top}\in\mathbb{C}^{I\times 1}. Finally, we can rewrite (4) as:

[TABLE]

In Formulation 1, the source vector \mbox{\boldmaths}(p) follows a zero-mean complex Gaussian distribution with $J(Q+1)\times J(Q+1)$ diagonal covariance matrix $\mathbf{R_{s}}(p)$ where the first $Q+1$ diagonal entries (for the first source) are $v_{1}(p),\dots,v_{1}(P-Q)$, the next $Q+1$ diagonal entries (for the second source) are $v_{2}(1),\dots,v_{2}(P)$, and so on. The pdf of the mixture given the sources is \mathcal{N}_{c}(\mbox{\boldmathx}(p);\mathbf{A}\mbox{\boldmaths}(p),\mathbf{\Sigma_{e}}).

III-D Vector/matrix Formulation 2

In Formulation 2, let

[TABLE]

denote the $j$-th source vector and the $i$-th noise and microphone vectors, all involving all $P$ frames. Concatenate them along source/microphone as:

[TABLE]

Define the CTF convolution matrix in $\mathbb{C}^{P\times P}$:

[TABLE]

where the CTF $\{a_{ij}(p)\}$ is first flipped and then duplicated as the row vectors, with one element shift per row. Concatenate it along source and microphone as:

[TABLE]

Then we can rewrite (4) as:

[TABLE]

In Formulation 2, the pdf of source vector \tilde{\mbox{\boldmaths}} is a zero-mean complex Gaussian distribution with $JP\times JP$ diagonal covariance matrix

[TABLE]

where $\text{Diag}\{\cdot\}$ denotes diagonal matrix of a vector. The noise vector \tilde{\mbox{\boldmathe}} follows a zero-mean complex Gaussian distribution with $IP\times IP$ covariance matrix $\mathbf{\Psi_{e}}$. The entries of $\mathbf{\Psi_{e}}$ are

[TABLE]

where the arguments in parentheses denotes the row and column indices; $i_{1},i_{2}\in[1,I],\ p_{1},p_{2}\in[1,P]$. The pdf of the mixture given the sources is \mathcal{N}_{c}(\tilde{\mbox{\boldmathx}};\mathbf{\mathcal{A}}\tilde{\mbox{\boldmaths}},\mathbf{\Psi_{e}}).

IV Expectation-Maximization Algorithm

Collect the source variances for all sources and frames, we have the source variance set $\mathbf{V}=\{v_{j}(p)\}_{j\in[1,J],p\in[1,P]}$. The parameters set in the present problem is $\Theta=\{\mathbf{V},\mbox{\bf A},\mathbf{\Sigma_{e}}\}$ in Formulation 1 or $\Theta=\{\mathbf{\Psi_{s}},\mathbf{\mathcal{A}},\mathbf{\Psi_{e}}\}$ in Formulation 2. The likelihood of the mixture will be maximized by an EM algorithm, in which the parameters $\Theta$ are the optimization variables. Meanwhile, the the STFT coefficients of the source signals, i.e. $\{s_{j}(p)\}_{j,p}$, are taken as hidden variables, whose posterior statistics will be inferred, and the posterior mean is taken as the estimation of the source signals. The proposed EM algorithm is summarized in Algorithm 1.

IV-A E-step

The E-step will be derived based on Formulation 2. Using the parameters estimates $\Theta$, given in the preceding M-step by (IV-B) in Formulation 1, we construct the CTF convolution matrix $\mathbf{\mathcal{A}}$, the source covariance matrix $\mathbf{\Psi_{s}}$ and the noise covariance matrix $\mathbf{\Psi_{e}}$ following (III-D), (8) and (9), respectively.

In Formulation 2, the posterior distribution of the source signals is p(\tilde{\mbox{\boldmaths}}|\tilde{\mbox{\boldmathx}},\Theta)\propto p(\tilde{\mbox{\boldmathx}}|\tilde{\mbox{\boldmaths}},\Theta)p(\tilde{\mbox{\boldmaths}}|\Theta). Since both p(\tilde{\mbox{\boldmathx}}|\tilde{\mbox{\boldmaths}},\Theta) and p(\tilde{\mbox{\boldmaths}}|\Theta) are Gaussian, p(\tilde{\mbox{\boldmaths}}|\tilde{\mbox{\boldmathx}},\Theta) is also Gaussian. Let \mathbb{E}_{\tilde{\mbox{\boldmaths}}|\tilde{\mbox{\boldmathx}},\Theta}[\ .\ ] denotes the expectation in the sense of the posterior distribution p(\tilde{\mbox{\boldmaths}}|\tilde{\mbox{\boldmathx}},\Theta). From the exponent of p(\tilde{\mbox{\boldmaths}}|\tilde{\mbox{\boldmathx}},\Theta), i.e.

[TABLE]

the posterior mean \widehat{\tilde{\mbox{\boldmaths}}}=\mathbb{E}_{\tilde{\mbox{\boldmaths}}|\tilde{\mbox{\boldmathx}},\Theta}[\tilde{\mbox{\boldmaths}}] and covariance matrix \mathbf{\widehat{\Sigma}_{s}}=\mathbb{E}_{\tilde{\mbox{\boldmaths}}|\tilde{\mbox{\boldmathx}},\Theta}[(\tilde{\mbox{\boldmaths}}-\widehat{\tilde{\mbox{\boldmaths}}})(\tilde{\mbox{\boldmaths}}-\widehat{\tilde{\mbox{\boldmaths}}})^{H}] can be derived by reorganizing the quadratic and linear forms in \tilde{\mbox{\boldmaths}} . We obtain:

[TABLE]

and then the posterior second-order moment matrix \mathbf{\widehat{\tilde{R}}_{s}}=\mathbb{E}_{\tilde{\mbox{\boldmaths}}|\tilde{\mbox{\boldmathx}},\Theta}[\tilde{\mbox{\boldmaths}}\tilde{\mbox{\boldmaths}}^{H}] can be computed as

[TABLE]

As derived based on the narrowband assumption in [4], Eq. (IV-A) is also in the form of classical Wiener filtering. The difference is, here the interframe elements in the posterior covariance matrix $\mathbf{\widehat{\Sigma}_{s}}$ are nonzero, which means the correlation between frames due to the convolution is encoded. As a result, the posterior mean of source, i.e. \widehat{\tilde{\mbox{\boldmaths}}}, is recovered by deconvoluting the mixture.

IV-B M-step

The M-step will be derived based on Formulation 2. Collect the multichannel mixture vectors and source signal vectors along frames, we have the observation set \mbox{\bf X}=\{\mbox{\boldmathx}(p)\}_{p\in[1,P]} and the source signal set \mbox{\bf S}=\{\mbox{\boldmaths}(p)\}_{p\in[1,P]}. The complete-data (including observations and hidden variables) likelihood function is:

[TABLE]

Considering only the terms related to the parameters and hidden variables, the corresponding loglikelihood writes:

[TABLE]

Denote the auxiliary function for likelihood maximization as $Q(\Theta,\Theta^{\text{old}})=\mathbb{E}_{\mathbf{S}|\mathbf{X},\Theta^{\text{old}}}[\log(p(\mbox{\bf X},\mbox{\bf S}|\Theta))]$, where $\Theta^{\text{old}}$ denotes the parameters estimated at the previous iteration. From the loglikelihood (IV-B), the auxiliary function can be derived as:

[TABLE]

where $\text{Trace}\{\cdot\}$ denotes matrix trace, and

[TABLE]

are the posterior statistics of the source signal, namely the posterior mean, the posterior second-order moment matrix and the element-wise posterior second-order moment, respectively. Actually, $\widehat{v}_{j}(p)$ is the $((j-1)(Q+1)+1)$-th diagonal entry of $\mathbf{\widehat{R}_{s}}(p)$.

Reformulation: \widehat{\mbox{\boldmaths}}(p) and $\mathbf{\widehat{R}_{s}}(p)$ can be obtained by reformulating \widehat{\tilde{\mbox{\boldmaths}}} and $\mathbf{\widehat{\tilde{R}}_{s}}$ derived in the preceding E-step. The reformulation is mainly to find the elements with the same source and frame indices:

[TABLE]

where the subscript {j,p} denote “the $j$-th source at $p$-th frame”. In a vector, or in a row/column of a matrix, {j,p} represents i) the $\{(j-1)P+p\}$-th element in \widehat{\tilde{\mbox{\boldmaths}}} and $\mathbf{\widehat{\tilde{R}}_{s}}$, based on Formulation 2, and ii) the $\{(j-1)(Q+1)+q+1\}$-th ($q\in[0,Q]$) element in \widehat{\mbox{\boldmaths}}(p) and $\mathbf{\widehat{R}_{s}}(p)$, based on Formulation 1.

With respect to $\mathbf{A}^{*}$ (∗ denotes conjugate), $v_{j}(p)$ and $\mathbf{\Sigma_{e}}$, the (complex) derivative of $Q(\Theta,\Theta^{\text{old}})$ are respectively:

[TABLE]

To maximize $Q(\Theta,\Theta^{\text{old}})$, the three derivatives are set equal to zero, then $\mathbf{A}$, $v_{j}(p)$ and $\mathbf{\Sigma_{e}}$ can be estimated as, respectively:

[TABLE]

V Experiments

V-A Experimental Configuration

The binaural (two-channel) simulated signals were used to evaluate the proposed EM algorithm. The experiments were conducted under various acoustic conditions, in terms of reverberation time, number of sources and intensity of room filter perturbation.

V-A1 Simulation set-up

We use a KEMAR dummy head [25] with one microphone embedded in each ear as the recording system. The head related impulse responses (HRIRs) for a large grid of directions were measured in advance. The ROOMSIM simulator [26] simulates the binaural room impulse responses using these HRIRs as both the direct-path wave and the reflections. Four reverberant conditions were simulated, with the reverberation time $T_{60}$ as [math] s (the anechoic case), $0.22$ s, $0.5$ s and $0.79$ s, respectively. The TIMIT [27] speech signals were used as the speech source signals, and were convolved with the simulated BRIRs to generate microphone signals (mixtures). The sampling rates of source signals and microphone signals are both 16-kHz, and the length of source signals are about 3 s. The speech sources were set to locate at different directions in front of the dummy head. The noisy microphone signals is generated by adding a spatially uncorrelated stationary speech-like noise to the noise-free signals. One SNR (signal-to-noise ratio) condition, i.e. 20 dB, is tested. The STFT uses Hamming window with a length of $1,\!024$ samples ($64$ ms) and frame step of $256$ samples ($16$-ms). In this experiment, the ground-truth noise covariance matrix $\mathbf{\Sigma_{e}}(k)$ is used, and fixed during EM iterations.

V-A2 EM initialization

Depends on what types of prior knowledges are available, the EM algorithm can be initialized either from the E-step or from the M-step. For both choices, the accuracy of the initialization is crucial for the EM iterations to converge to a good solution. In this experiment, we initialize EM from the M-step. Due to the difficulty of the blind initialization, we consider a semi-blind initialization scheme. The time-domain filters are assumed to be known, from which the CTFs are computed by (3). To fit the realistic situation that the time-domain filters (or CTFs) should actually be blindly estimated and suffer from some estimation error, a proportional Gaussian random noise is added to the time-domain filters to generate the perturbed filters. The normalized projection misalignment (NPM) [28] in decibels (dB) is used to measure the intensity of the perturbation. The lower NPM is, the intenser the perturbation will be. To have a good initialization for the source variance, the CTF-Lasso method proposed in [17] is first applied, which solves the problem:

[TABLE]

where $\|\cdot\|_{1}$ denotes $\ell_{1}$-norm, and is used to impose the spectral spasity of speech sources. For more details, please refer to [17]. Then the source variance is initialized as the magnitude square of each source coefficient estimate. It is found that, for most of acoustic conditions and frequency bins, the number of EM iterations required for convergency is less than 10, thence in this experiment we use a constant number, i.e. 7, of EM iterations.

V-A3 Baseline methods

Three baseline methods were used for comparison: i) The CTF-Lasso method used for initialization; ii) The binary masking (BM) method [1], which is based on the narrowband approximation. To make a fair comparison, the narrowband mixing filters are also computed using the knwon perturbed time-domain filters. However, to compute the one-order mixing filters, for the high reverberation time cases, the time-domain filters should be first truncated to have a length being equal to (or less than) the STFT window length. Based on some pilot experiments, we use the HRIRs (without reverberation) as the truncated filters, which achieves the best results. For source separation, each TF bin is assigned to one of the sources based on the mixing filters; iii) the full-rank spatial covariance matrix (FR-SCM) method [10]. The full-rank spatial covariance matrix for each source was separately estimated using the corresponding source image, and kept fixed during the EM, following the line of the semi-oracle experiments in [10].

V-A4 Performance metrics

The signal-to-distortion ratio (SDR), signal-to-interferences ratio (SIR) and signal-to-artifacts ratio (SAR) [29], all in decibels (dB), are used as the separation performance metric. In the following, three sets of experiments will be conducted: i) for various reverberation time, ii) for various numbers of sources, i.e. with $2$, $3$, $4$ and $5$ sources, iii) with various NPM settings. For each condition, the metric scores are averaged over 20 mixtures.

The computation complexity of each method is measured with the real-time factor, which is the processing time of one method divided by the length of the processed signal. Note that all the methods are implemented in MATLAB.

V-B Results as a Function of Reverberation Time

Fig. 1 plots the performance measures obtained for the $4$ reverberation times. The number of sources is $3$. NPM is set to -35 dB, for which the filter pertrubation is light, and the CTFs used for CTF-Lasso and the proposed method are very accurate. Thence, for this experiment, the CTFs are fixed during EM iterations. For the anechoic case, all the four methods largely improve the SDR and SIR scores compared to the unprocessed signals, however slightly reduces the SAR score. This indicates that the three methods can efficiently separate the multiple sources, but introduces some artifacts. Especially, BM suffers from more artifacts than the other methods, since the hard assignment of the time-frequency bins to the dominant source largely distorts the less dominant sources. As $T_{60}$ increases, the performance measures of BM and FR-SCM dramatically decrease, since the length of RIRs is much larger than the STFT window and the narrowband approximation is not valid any more for the high reverberation cases. FR-SCM outperforms BM due to the use of the full-rank convariance matrix, which however is only suitable for the low reverberation cases. This can be testified by the fact that the most prominent advantage of FR-SCM over BM presents when $T_{60}$ is 0.22 s. In contrast to BM and FR-SCM, CTF-Lasso and the proposed EM method achieve good performances. It is a bit surprising that the performance measures actually increases with the increase of reverberation time. The possible reason is that the long filters involve more information to differentiate and separate the multiple sources. Of course, to satisfy this assertion the filters should be accurate enough. Compared to CTF-Lasso whose outputs are taken as the initial point for the EM algorithm, the EM algorithm improves the SDR by about 1.5 dB for every reverberation times, which indicates that the EM iterations are able to refine the quality of the source estimate.

V-C Results as a Function of Number of Sources

Fig. 2 plots the results for various number of sources, for $T_{60}=0.5$ s. NPM is set to -35 dB, and the CTFs are fixed during EM iterations. As expected, the performance measures of all methods degrade when the number of sources increases. For BM, the WDO assumption for speech sources becomes less valid when more sources present. For FR-SCM, CTF-Lasso and the proposed EM method, the mutual confusion of sources gets larger with increasing number of sources. The performance degradation rate of the four methods are similar to the one for the unprocessed signals. Overall, when the CTFs are properly initialized, CTF-Lasso and the proposed EM method achieve good source separtion performance even for the mixtures with five sources using only two microphones.

V-D Results as a Function of NPM

To evaluate the proposed method under the conditions that the initialized CTFs suffer from large estimation error, we conducted the experiments with various NPM settings, and the results are shown in Fig. 3. Since the initialization is not accurate, in this experiment, the CTFs are updated during EM interations to refine the CTFs. To demonstrate the efficiency of CTF updation, the results with the fixed CTFs are also given. As NPM increases, the performance of CTF-Lasso and the proposed EM method largely degrade. When NPM is larger than -20 dB, the proposed method with fixed CTFs does not improve the performance of CTF-Lasso, which means the source estimation can not be refined based the inaccurate CTFs. In contrast, the proposed method with updating CTFs improves the performance of CTF-Lasso, due to the refining of CTFs in the M-step. The performance measures of CTF-Lasso and the proposed method become close to the ones of BM and FR-SCM, which indicates that the CTF-based methods are more sensitive to the filter perturbations than the narrowband assumption based methods.

V-E Computational Complexity Analysis

Table I shows the real-time factor of the four methods for the case with three sources and $T_{60}$ of 0.5 s. BM is the fastest, since it is an one-step method. The other three are all iterative methods, whose computational complexity is proportional to the number of iterations. FR-SCM and the proposed method are similar in the sense that they both are based on EM iteration: estimating the source statistics using a Wiener-like filter in the E-step, and estimating the mixing filters in the M-step. The main difference between them is that the proposed method uses the CTF with a length of $Q$, while FR-SCM uses the mixing filter with a length of $1$. As a result, the proposed method has a much larger complexity than FR-SCM. CTF-Lasso also uses the CTF. However, unlike FR-SCM and the proposed method that several matrix inverse operations are performed (as shown in (IV-A) and (IV-B)), Lasso optimization only executes the first-order convolution operation. Thence, the real-time factor of CTF-Lasso is lower than the factor of FR-SCM and the proposed method.

VI Conclusion

In this work, an EM algorithm has been proposed for speech separation. The subband convolutive model, i.e. CTF model, was adopt. To concisely derive the M-step and E-step, two convolution vector/matrix formulations were used. The CTF model based methods, i.e. CTF-Lasso and the proposed EM method, outperform the narrowband assumption based method, i.e. binary masking and FR-SCM, for the reverberant case. The proposed EM algorithm is capable to refine the CTFs and source estimates starting with the output of CTF-Lasso, and thus improves the source separation performance. Only the semi-blind experiments were conducted in this work due to the difficulty of EM initialization. In the future, a blind CTF identification method could be developed to enable the blind initialization of EM. To this aim, the CTF identification methods proposed in [23] and [30] could be combined, which are in the contexts of single-source dereverberation and multisource localization, respectively.