Speech Enhancement Using an Iterative Posterior NMF

Abstract

Over the years, miscellaneous methods for speech enhancement have been proposed, such as spectral subtraction (SS) and minimum mean square error (MMSE) estimators. These methods do not require any prior knowledge about the speech and noise signals nor any training stage beforehand, so they are highly flexible and allow implementation in various situations. However, these algorithms usually assume that the noise is stationary and are thus not good at dealing with nonstationary noise types, especially under low signal-to-noise (SNR) conditions. To overcome the drawbacks of the above methods, nonnegative matrix factorization (NMF) is introduced. NMF approach is more robust to nonstationary noise. In this chapter, we are actually interested in the application of speech enhancement using NMF approach. A speech enhancement method based on regularized nonnegative matrix factorization (NMF) for nonstationary Gaussian noise is proposed. The spectral components of speech and noise are modeled as Gamma and Rayleigh, respectively. We propose to adaptively estimate the sufficient statistics of these distributions to obtain a natural regularization of the NMF criterion.

1. Introduction

Over the past several decades, there has been a large research interest in the problem of single-channel sound source separation. Such work focuses on the task of separating a single mixture recording into its respective sources and is motivated by the fact that real-world sounds are inherently constructed by many individual sounds (e.g., human speakers, musical instruments, background noise, etc.). While source separation is difficult, the topic is highly motivated by many outstanding problems in audio signal processing and machine learning, including the following:

1. Speech denoising and enhancement—the task of removing background noise (e.g., wind, babble, etc.) from recorded speech and improving speech intelligibility for human listeners and/or automatic speech recognizers

2. Content-based analysis and processing—the task of extracting and/or processing audio based on semantic properties of the recording such as tempo, rhythm, and/or pitch

3. Music transcription—the task of notating an audio recording into a musical representation such as a musical score, guitar tablature, or other symbolic notations

4. Audio-based forensics—the task of examining, comparing, and evaluating audio recordings for scientific and/or legal matters

5. Audio restoration—the task of removing imperfections such as noise, hiss, pops, and crackles from (typically old) audio recordings

6. Music remixing and content creation—the task of creating a new musical work by manipulating the content of one or more previously existing recordings

2. Nonnegative matrix factorization

2.1 NMF model

Nonnegative matrix factorization is a process that approximates a single nonnegative matrix as the product of two nonnegative matrices. It is defined by

$$V\approx \mathit{WH}$$ E1

$V\in R_{+}^{N_f\times N_t}$ is a nonnegative input matrix. $W\in R_{+}^{N_f\times N_z}$ is a matrix of basis vectors, basis functions, or dictionary elements; $H\in R_{+}^{N_z\times N_t}$ is a matrix of corresponding activations, weights, or gains; and $N_f$ is the number of rows of the input matrix. $N_t$ is the number of columns of the input matrix; $N_z$ is the number of basis vectors [1].

$V\in R_{+}^{N_f\times N_t}$ —original nonnegative input data matrix

• Each column is an $N_f$ -dimensional data sample.

• Each row represents a data feature.

$W\in R_{+}^{N_f\times N_z}$ —matrix of basis vectors, basis functions, or dictionary elements.

• A column represents a basis vector, basis function, or dictionary element.

• Each column is not orthonormal, but commonly normalized to one.

$H\in R_{+}^{N_z\times N_t}$ —matrix of activations, weights, encodings, or gains.

• A row represents the gain of a corresponding basis vector.

• Each row is not orthonormal, but sometimes normalized to one.

When used for audio applications, NMF is typically used to model spectrogram data or the magnitude of STFT data [2]. That is, we take a single-channel recording, transform it into the time-frequency domain using the STFT, take the magnitude or power V, and then approximate the result as $V\approx \mathit{WH}$ . In doing so, NMF approximates spectrogram data as a linear combination of prototypical frequencies or spectra (i.e., basis vectors) over time.

This process can be seen in Figure 1 [3], where a two-measure piano passage of “Mary Had a Little Lamb” is shown alongside a spectrogram and an NMF factorization. Notice how W captures the harmonic content of the three pitches of the passage and H captures the time onsets and gains of the individual notes. Also note that $N_z$ is typically chosen manually or using a model selection procedure such as cross-validation and $N_f$ and $N_t$ are a function of the overall recording length and STFT parameters (transform length, zero-padding size, and hop size).

Figure 1.

NMF of a piano performing “Mary had a little lamb” for two measures with $N_z$ = 3. Notice how matrix W captures the harmonic content of the three pitches of the passage and matrix H captures the time onsets and gains of the individual notes [3].

This leads to two related interpretations of how NMF models spectrogram data. The first interpretation is that the columns of V (i.e., short-time segments of the mixture signal) are approximated as a weighted sum of basis vectors as shown in Figure 2 and Eq. (2):

$$V\approx \left[\begin{array}{ccccc}\left|\right)& \left|\right)& \left|\right)& & \left|\right)\\ V1& V2& V_3& \dots \dots & V_{N_t}\\ \left|\right)& \left|\right)& \left|\right)& & \left|\right)\end{array}\right]\approx \left[{\sum }_{j=1}^{N_z}{H}_{j1}{W}_j{\sum }_{j=1}^K{H}_{j2}{W}_j{\sum }_{j=1}^K{H}_{j{N}_t}{W}_j\right]$$ E2

Figure 2.

NMF interpretation I. the columns of V (i.e., short-time segments of the mixture signal) are approximated as a weighted sum or mixture of basis vectors W [3].

The second interpretation is that the entire matrix V is approximated as a sum of matrix “layers,” as shown in Figure 3 and Eq. (3).

$$V\approx \left[\begin{array}{ccccc}\left|\right)& \left|\right)& \left|\right)& & \left|\right)\\ V1& V2& V_3& \dots \dots & V_{N_t}\\ \left|\right)& \left|\right)& \left|\right)& & \left|\right)\end{array}\right]\approx \left[\begin{array}{ccccc}\left|\right)& \left|\right)& \left|\right)& & \left|\right)\\ W1& W2& W_3& \dots \dots & W_{N_z}\\ \left|\right)& \left|\right)& \left|\right)& & \left|\right)\end{array}\right]\left[\begin{array}{c}{h}_1^T\\ h_2^T\\ h_3^T\\ .\\ h_{N_z}^T\end{array}\right]$$

$$V\approx W_1{h}_1^T+W_2{h}_2^T+W_3{h}_3^T+..\dots +W_{N_z}{h}_{N_z}^T$$ E3

Figure 3.

NMF interpretation II. The matrix V (i.e., the mixture signal) is approximated as a sum of matrix “layers” [3].

The application of NMF on noisy speech can be seen in Figure 4.

Figure 4.

Applying NMF on noisy speech.

2.2 Optimization formulation

To estimate the basis matrix W and the activation matrix H for a given input data matrix V, NMF algorithm is formulated as an optimization problem. This is written as:

$$\underset{\mathrm{W},\mathrm{H}}{\mathit{\text{argmin}}}D\left(V\left|\mathrm{WH}\right)\right)$$

$$W\ge 0,H\ge 0$$ E4

where $D\left(V\left|\mathrm{WH}\right)\right)$ is an appropriately defined cost function between V and W H and the inequalities $\ge$ are element-wise. It is also common to add additional equality constraints to require the columns of W to sum to one, which we enforce. When $D\left(V\left|\mathrm{WH}\right)\right)$ is additively separable, the cost function can be reduced to

$$D\left(V\left|\mathrm{WH}\right)\right)=\sum _{f=1}^{N_f}\sum _{t=1}^{N_t}d\left(V_{\mathit{\text{ft}}}\left|{\left[\mathrm{WH}\right]}_{\mathit{\text{ft}}}\right)\right)$$ E5

where ${\left[\right]}_{\mathit{\text{ft}}}$ indicates its argument at row f and column t and $D\left(V\left|\mathrm{WH}\right)\right)$ is a scalar cost function measured between V and WH.

Popular cost functions include the Euclidean distance metric, Kullback-Liebler (KL) divergence, and Itakura-Saito (IS) divergence. Both the KL and IS divergences have been found to be well suited for audio purposes. In this work, we focus on the case where $d\left(q\left|p\right)\right)$ is generalized (non-normalized) KL divergence:

$$d_{\mathit{KL}}\left(q\left|p\right)\right)=\mathit{q} \mathit{ln}\frac{q}{p}-q+p$$ E6

where ${\left[\right]}_{\mathit{\text{ft}}}$ indicates its argument at row f and column t and $d\left(q\left|p\right)\right)$ is a scalar cost function measured between q and p.

This results in the following optimization formulation:

$$\underset{\mathrm{W},\mathrm{H}}{\mathit{\text{argmin}}}\sum _{f=1}^{N_f}\sum _{t=1}^{N_t}-V_{\mathit{\text{ft}}} \mathit{ln}\left|{\left[\mathrm{WH}\right]}_{\mathit{\text{ft}}}\right)+\left|{\left[\mathrm{WH}\right]}_{\mathit{\text{ft}}}\right)+\mathit{\text{const}}$$

Subject to

$$W\ge 0,H\ge 0$$ E7

Given this formulation, we notice that the problem is not convex in W and H, limiting our ability to find a globally optimal solution to Eq. (7). It is, however, biconvex or independently convex in W for a fixed value of H and convex in H for a fixed value of W, motivating the use of iterative numerical methods to estimate locally optimal values of W and H.

2.3 Parameter estimation

To solve Eq. (7), we must use an iterative numerical optimization technique and hope to find a locally optimal solution. Gradient descent methods are the most common and straightforward for this purpose but typically are slow to converge. Other methods such as Newton’s method, interior-point methods, conjugate gradient methods, and similar [4] can converge faster but are typically much more complicated to implement, motivating alternative approaches.

The most popular alternative that has been proposed is by Lee and Seung [1, 5] and consists of a fast, simple, and efficient multiplicative gradient descent-based optimization procedure. The method works by breaking down the larger optimization problem into two subproblems and iteratively optimizes over W and then H, back and forth, given an initial feasible solution. The approach monotonically decreases the optimization objective for both the KL divergence and Euclidean cost functions and converges to a local stationary point.

The approach is justified using the machinery of majorization-minimization (MM) algorithms [6]. MM algorithms are closely related to expectation maximization (EM) algorithms. In general, MM algorithms operate by approximating an optimization objective with a lower bound auxiliary function. The lower bound is then maximized instead of the original function, which is usually more difficult to optimize.

Algorithm 1 shows the complete iterative numerical optimization procedure applied to Eq. (7) with the KL divergence, where the division is element-wise, $\otimes$ is an element-wise multiplication, and 1 is a vector or matrix of ones with appropriately defined dimensions [5].

Algorithm 1 KL-NMF parameter estimation

Procedure KL-NMF ( $V\in R_{+}^{N_f\times N_t}$ //input data matrix.

$N_z$ //number of basic vectors.)

Initialize: $W\in R_{+}^{N_f\times N_z}$ , $H\in R_{+}^{N_z\times N_{t.}}$

repeat

Optimize over W

$$W\leftarrow W\otimes \frac{\left(\frac{V}{\mathrm{WH}}\right)H^T}{1H^T}$$ E8

Optimize over H

$$H\leftarrow H\otimes \frac{W^T\left(\frac{V}{\mathrm{WH}}\right)}{W^{T}1}$$ E9

until convergence

return: W and H

NMF is an optimization technique using EM algorithm in terms of matrix, whereas probabilistic latent component analysis (PLCA) is also an optimization technique using EM algorithm in terms of probability. In PLCA, we are going to incorporate probabilities of time and frequency. In the next section, the development of PLCA-based algorithm to incorporate time-frequency constraints is discussed.

3. A probabilistic latent variable model with time-frequency constraints

Considering this approach, we now develop a new PLCA-based algorithm to incorporate the time-frequency user-annotations. For clarity, we restate the form of the symmetric two-dimensional PLCA model we use:

$$p\left(f,t\right)={\sum }_{z}p\left(z\right)p\left(f\left|z\right)\right)p\left(t\left|z\right)\right)$$ E10

Compared to a modified NMF formulation, incorporating optimization constraints as a function of time, frequency, and sound source into the factorized PLCA model is particularly interesting and motivating to our focus.

Incorporating prior information into this model, and PLCA in general, can be done in several ways. The most commonly used methods are by direct observations (i.e., setting probabilities to zero, one, etc.) or by incorporating Bayesian prior probabilities on model parameters. Direct observations do not give us enough control, so we consider incorporating Bayesian prior probabilities. For the case of Eq. (10), this would result in independently modifying the factor terms $p\left(f\left|z\right)\right)$ , $p\left(t\left|z\right)\right)$ , or $p\left(z\right)$ . Common prior probability distributions used for this purpose include Dirichlet priors [7], gamma priors [8], and others.

Given that we would like to incorporate the user-annotations as a function of time, frequency, and sound source, however, we notice that this is not easily accomplished using standard priors. This is because the model is factorized, and each factor is only a function of one variable and (possibly) conditioned by another, making it difficult to construct a set of prior probabilities that, when jointly applied to $p\left(f\left|z\right)\right)$ , $p\left(t\left|z\right)\right)$ , and/or $p\left(z\right)$ , would encourage or discourage one source or another to explain a given time-frequency point. We can see this more clearly when we consider PLCA to be the following simplified estimation problem:

$$X\left(f,t\right)\approx \phi \left(z\right)\phi \left(f,z\right)\phi \left(t,z\right)$$ E11

where $X\left(f,t\right)$ is the observed data that we model as the product of three distinct functions or factors $\phi \left(z\right)$ , $\phi \left(f,z\right)$ , and $\phi \left(t,z\right)$ . Note, each factor has different input arguments and each factor has different parameters that we wish to estimate via EM. Also, forget for the moment that the factors must be normalized probabilities.

Given this model, if we wish to incorporate additional information, we could independently modify:

• $\phi \left(z\right)$ to incorporate past knowledge of the variable z

• $\phi \left(f,z\right)$ to incorporate past knowledge of the variable f and z

• $\phi \left(t,z\right)$ to incorporate past knowledge of the variable t and z

This way of manipulation allows us to maintain our factorized form and can be thought of as prior-based regularization. If we would like to incorporate additional information/regularization that is a function of all three variables z, f, and t, then we must do something else. The first option would be to try to simultaneously modify all factors together to impose regularization that is a function of all three variables. This is unfortunately very difficult—both conceptually difficult to construct and practically difficult to algorithmically solve.

This motivates the use of posterior regularization (PR). PR provides us with an algorithmic mechanism via EM to incorporate constraints that are complementary to prior-based regularization. Instead of modifying the individual factors of our model as we saw before, we directly modify the posterior distribution of our model. The posterior distribution of our model, very loosely speaking, is a function of all random variables of our model. It is natively computed within each E step of EM and is required to iteratively improve the estimates of our model parameters. In this example, the posterior distribution would be akin to $\phi \left(z,f,t\right)$ , which is a function of t, f, and z, as required. We now formally discuss PR below, beginning with a general discussion and concluding with the specific form of PR we employ within our approach.

3.1 Posterior regularization

The framework of posterior regularization, first introduced by Graca, Ganchev, and Taskar [9, 10], is a relatively new mechanism for injecting rich, typically data-dependent constraints into latent variable models using the EM algorithm. In contrast to standard Bayesian prior-based regularization, which applies constraints to the model parameters of a latent variable model in the maximization step of EM, posterior regularization applies constraints to the posterior distribution (distribution over the latent variables, conditioned on everything else) computation in the expectation step of EM. The method has found success in many natural language processing tasks, such as statistical word alignment, part-of-speech tagging, and similar tasks that involve latent variable models.

In this case, what we do is constrain the distribution q in some way when we maximize the auxiliary bound $F\left(q,\Theta \right)$ with respect to q in the expectation step of an EM algorithm, resulting in

$$q^{n+1}=\underset{q}{\mathit{\text{argmin}}}\mathit{KL}\left(q\left.‖p\right)\right)+\Omega \left(q\right)$$ E12

where $\Omega \left(q\right)$ constrains the possible space of q.

Note, the only difference between Eq. (12) and our past discussion on EM is the added term $\Omega \left(q\right)$ . If $\Omega \left(q\right)$ is set to zero, we get back the original formulation and easily solve the optimization by setting q = p without any computation (except computing the posterior p). Also note to denote the use of constraints in this context, the term “weakly supervised” was introduced by Graca [11] and is similarly adopted here.

This method of regularization is in contrast to prior-based regularization, where the modified maximization step would be

$${\Theta }^{n+1}=\underset{\Theta }{\mathit{\text{argmax}}}F\left(q^{n+1},\Theta \right)+\Omega \left(\Theta \right)$$ E13

where $\Omega \left(\Theta \right)$ constrains the model parameter $\Theta$ .

3.2 Linear grouping expectation constraints

Given the general framework of posterior regularization, we need to define a meaningful penalty $\Omega \left(q\right)$ for which we map our annotations. We do this by mapping the annotation matrices to linear grouping constraints on the latent variable z. To do so, we first notice that Eq. (12) decouples for each time-frequency point for our specific model. Because of this, we can independently solve Eq. (12) for each time-frequency point, making the optimization much simpler. When we rewrite our E step optimization using vector notation, we get

$$\underset{q}{\mathit{\text{argmin}}}-q_{\mathit{\text{ft}}}^T\mathit{ln}{p}_{\mathit{\text{ft}}}+q_{\mathit{\text{ft}}}^T\mathit{ln}{q}_{\mathit{\text{ft}}}$$

subject to $$q_{\mathit{\text{ft}}}^{T}1=1,q_{\mathit{\text{ft}}}\ge 0$$ E14

where q and $p\left(z\left|f,t\right)\right)$ for a given value of f and t is written as $q_{\mathit{\text{ft}}}$ and $p_{\mathit{\text{ft}}}$ without any modification; we note q is optimal when equal to $p\left(z\left|f,t\right)\right)$ as before.

We then apply our linear grouping constraints independently for each time-frequency point:

$$\underset{q}{\mathit{\text{argmin}}}-q_{\mathit{\text{ft}}}^T\mathit{ln}{p}_{\mathit{\text{ft}}}+q_{\mathit{\text{ft}}}^T\mathit{ln}{q}_{\mathit{\text{ft}}}+q_{\mathit{\text{ft}}}^T{\lambda }_{\mathit{\text{ft}}}$$

Subject to $$q_{\mathit{\text{ft}}}^{T}1=1,q_{\mathit{\text{ft}}}\ge 0,$$ E15

where we define ${\lambda }_{\mathit{\text{ft}}}={\left[{\Lambda }_{\mathit{\text{ft}}1}..\dots \dots {\Lambda }_{\mathit{\text{ft}}1}{\Lambda }_{\mathit{\text{ft}}2}\dots \dots \dots {\Lambda }_{\mathit{\text{ft}}2}\right]}^T\in R^{N_z}$ as the vector of user-defined penalty weights, $T$ is a matrix transpose, $\ge$ is element-wise greater than or equal to, and 1 is a column vector of ones. In this case, positive-valued penalties are used to decrease the probability of a given source, while negative-valued coefficients are used to increase the probability of a given source. Note the penalty weights imposed on the group of values of z that correspond to a given source s are identical, linear with respect to the z variables, and applied in the E step of EM, hence the name “linear grouping expectation constraints.”

To solve the above optimization problem for a given time-frequency point, we form the Lagrangian

$$L\left(q_{\mathit{\text{ft}}},\gamma \right)=-q_{\mathit{\text{ft}}}^T\mathit{ln}{p}_{\mathit{\text{ft}}}+q_{\mathit{\text{ft}}}^T\mathit{ln}{q}_{\mathit{\text{ft}}}+q_{\mathit{\text{ft}}}^T{\lambda }_{\mathit{\text{ft}}}+\gamma \left(1-q_{\mathit{\text{ft}}}^{T}1\right)$$ E16

With $\gamma$ being a Lagrange multiplier, take the gradient with respect to q and $\gamma$ :

$${\nabla }_{q_{\mathit{\text{ft}}}}L\left(q_{\mathit{\text{ft}}},\gamma \right)=-\mathit{ln}{p}_{\mathit{\text{ft}}}+1+\mathit{ln}{q}_{\mathit{\text{ft}}}+{\lambda }_{\mathit{\text{ft}}}-\gamma 1=0$$ E17

$${\nabla }_{a}L\left(q_{\mathit{\text{ft}}},\gamma \right)=\left(1-q_{\mathit{\text{ft}}}^{T}1\right)=0$$ E18

set Eqs. (17) and (18) equal to zero, and solve for $q_{\mathit{\text{ft}}}$ , resulting in

$$q_{\mathit{\text{ft}}}=\frac{P_{\mathit{\text{ft}}}\otimes \mathit{exp}\left(-{\lambda }_{\mathit{\text{ft}}}\right)}{P_{\mathit{\text{ft}}}^T\mathit{exp}\left(-{\lambda }_{\mathit{\text{ft}}}\right)}$$ E19

where exp{} is an element-wise exponential function.

Notice the result is computed in closed form and does not require any iterative optimization scheme as may be required in the general posterior regularization framework [9], minimizing the computational cost when incorporating the constraints. Also note, however, that this optimization must be solved for each time-frequency point of our spectrogram data for each E step iteration of our final EM parameter estimation algorithm.

3.3 Parameter estimation

Now knowing the posterior-regularized expectation step optimization, we can derive a complete EM algorithm for a posterior-regularized two-dimensional PLCA model (PR-PLCA):

$$p\left(z\left|f,t\right)\right)\leftarrow \frac{p\left(z\right)p\left(f\left|z\right)\right)p\left(t\left|z\right)\right){\overline{\Lambda }}_{\mathit{ftz}}}{{\sum }_{z\prime }p\left(z\prime \right)p\left(f\left|z\right)\prime \right)p\left(t\left|z\right)\prime \right){\overline{\Lambda }}_{\mathit{ftz}\prime }}$$ E20

where $\overline{\Lambda }=\text{exp}\left\{-\Lambda \right\}$ . The entire algorithm is outlined in Algorithm 2. Notice we continue to maintain closed-form E and M steps, allowing us to optimize further and draw connections to multiplicative nonnegative matrix factorization algorithms.

Algorithm 2 PR-PLCA with linear grouping expectation constraints

Procedure PLCA (

$V\in R_{+}^{N_f\times N_t}$ //observed normalized data

$N_z$ //number of basis vectors

$N_s$ //number of sources

$\Lambda \in R^{N_f\times N_t\times N_z}$ //penalties

)

Initialize: feasible $p\left(z\right)$ , $p\left(f\left|z\right)\right)$ and $p\left(t\left|z\right)\right)$

$$\mathbf{\text{Precompute}}\mathbf{:}\overline{\Lambda }\leftarrow \mathit{exp}\left(-\Lambda \right)$$ E21

repeat

Expectation step

for all z, f, t do

$$p\left(z\left|f,t\right)\right)\leftarrow \frac{p\left(z\right)p\left(f\left|z\right)\right)p\left(t\left|z\right)\right){\overline{\Lambda }}_{\mathit{ftz}}}{{\sum }_{z\prime }p\left(z\prime \right)p\left(f\left|z\right)\prime \right)p\left(t\left|z\right)\prime \right){\overline{\Lambda }}_{\mathit{ftz}\prime }}$$ E22

end for

Maximization step

for all z, f, t do

$$p\left(f\left|z\right)\right)=\frac{{\sum }_t{V}_{\mathit{\text{ft}}}p\left(z\left|f,t\right)\right)}{{\sum }_{f\prime }{\sum }_{t\prime }{V}_{f\prime t\prime }p\left(z\left|f\prime ,t\prime \right)\right)}$$ E23

$$p\left(t\left|z\right)\right)=\frac{{\sum }_f{V}_{\mathit{\text{ft}}}p\left(z\left|f,t\right)\right)}{{\sum }_{f\prime }{\sum }_{t\prime }{V}_{f\prime t\prime }p\left(z\left|f\prime ,t\prime \right)\right)}$$ E24

$$p\left(z\right)=\frac{{\sum }_f{\sum }_t{V}_{\mathit{\text{ft}}}p\left(z\left|f,t\right)\right)}{{\sum }_{z\prime }{\sum }_{f\prime }{\sum }_{t\prime }{V}_{f\prime t\prime }p\left(z\left|f\prime ,t\prime \right)\right)}$$ E25

end for

until convergence

return: $p\left(f\left|z\right)\right)$ , $p\left(t\left|z\right)\right)$ , $p\left(z\right)$ and $p\left(z\left|f,t\right)\right)$

• Multiplicative Update Equations

We can rearrange the expressions in Algorithm 2 and convert to a multiplicative form following similar methodology to Smaragdis and Raj [12].

Rearranging the expectation and maximization steps, in conjunction with Bayes’ rule, and

$$Z\left(f,t\right)={\sum }_{z}p\left(z\right)p\left(f\left|z\right)\right)p\left(t\left|z\right)\right){\overline{\Lambda }}_{\mathit{ftz}},$$

we get

$$p\left(z\left|f,t\right)\right)=\frac{p\left(f\left|z\right)\right)p\left(t,z\right){\overline{\Lambda }}_{\mathit{ftz}}}{Z\left(f,t\right)}$$ E26

$$p\left(t,z\right)={\sum }_f{V}_{\mathit{\text{ft}}}q\left(z\left|f,t\right)\right)$$ E27

$$p\left(f\left|z\right)\right)=\frac{{\sum }_t{V}_{\mathit{\text{ft}}}q\left(z\left|f,t\right)\right)}{{\sum }_{t}p\left(t,z\right)}$$ E28

$$p\left(z\right)={\sum }_{t}p\left(t,z\right)$$ E29

Rearranging further, we get

$$p\left(f\left|z\right)\right)=\frac{p\left(f\left|z\right)\right){\sum }_t\frac{V_{\mathit{\text{ft}}}{\overline{\Lambda }}_{\mathit{ftz}}}{Z\left(f,t\right)}p\left(t,z\right)}{{\sum }_{t}p\left(t,z\right)}$$ E30

$$p\left(t,z\right)=p\left(t,z\right){\sum }_{f}p\left(f\left|z\right)\right)\frac{V_{\mathit{\text{ft}}}{\overline{\Lambda }}_{\mathit{ftz}}}{Z\left(f,t\right)}$$ E31

which fully specifies the iterative updates. By putting Eqs. (30) and (31) in matrix notation, we specify the multiplicative form of the proposed method in Algorithm 3.

Algorithm 3. PR-PLCA with linear grouping expectation constraints in matrix notation

Procedure PLCA (

$V\in R_{+}^{N_f\times N_t}$ //observed normalized data

$N_z$ //number of basis vectors

$N_s$ //number of sources

${\Lambda }_s\in R^{N_f\times N_t},\forall \in \left\{1,.\dots N_s\right\}$ //penalties

)

Initialize: $W\in R_{+}^{N_f\times N_z}$ , $H\in R_{+}^{N_z\times N_t}$

Precompute:

For all s do

$${\overline{\Lambda }}_s\leftarrow \mathit{exp}\left\{-{\Lambda }_s\right\}$$ E32

$$X_s\leftarrow V\otimes {\overline{\Lambda }}_s$$ E33

End for

Repeat

$$\Gamma \leftarrow {\sum }_s\left(W_s{H}_s\right)\otimes {\overline{\Lambda }}_s$$ E34

For all s do

$$Z_s\leftarrow \frac{X_s}{\Gamma }$$ E35

$$W_{\left(s\right)}\leftarrow W_s\otimes \frac{Z_s{H}_s^T}{1H_s^T}$$ E36

$$H_{\left(s\right)}\leftarrow H_s\otimes \left(W_s^T{Z}_s\right)$$ E37

End for

until convergence

return: W and H

4. An iterative posterior NMF method for speech enhancement in the presence of additive Gaussian noise (proposed algorithm)

Over the past several years, research has been carried out in single-channel sound source separation methods. This problem is motivated by speech denoising, speech enhancement [13], music transcription [14], audio-based forensic, and music remixing. One of the most effective approach is nonnegative matrix factorization (NMF) [5]. The user-annotations can be used to obtain the PR terms [15]. If the number of sources is more, then it is difficult to identify sources in the spectrogram. In such cases, the user interaction-based constraint approaches are inefficient.

In order to avoid the previous problem, in the proposed method, an automatic iterative procedure is introduced. The spectral components of speech and noise are modeled as Gamma and Rayleigh, respectively [16].

4.1 Notation and basic concepts

Let noisy speech signal x[n] be the sum of clean speech s[n] and noise d[n] and their corresponding magnitude spectrogram be represented as

$$\left|\mathit{X}\left(\mathit{f},\mathit{t}\right)\right|=\left|\mathit{S}\left(\mathit{f},\mathit{t}\right)\right|+\left|\mathit{D}\left(\mathit{f},\mathit{t}\right)\right|$$ E38

where f represents the frequency bin and t the frame number. The observed magnitudes in time-frequency are arranged in a matrix X $\in R_{+}^{f\times t}$ of nonnegative elements. The source separation algorithms based on NMF pursue the factorization of X as a product of two nonnegative matrices, $\mathit{W}=\left[{\mathit{w}}_1,{\mathit{w}}_2,\dots ,{\mathit{w}}_R\right]\in R_{+}^{f\times R}$ in which the columns collect the basis vectors and $\mathit{H}={\left[h_1^T,h_2^T,.\dots ,h_R^T\right]}^T\in R_{+}^{R\times t}$ that collects their respective weights, i.e.,

$$X=\mathrm{WH}=\sum _{z=1}^R{W}_z{H}_z$$ E39

where R denotes the number of latent components.

4.2 Proposed regularization

There are several ways to incorporate the user-annotations into latent variable models, for instance, by using the suitable regularization functions. For expectation maximization (EM) algorithms, posterior regularization was introduced by [9, 11]. This method is data dependent. This method gives richness and also gives the constraints on the posterior distributions of latent variable models. The applications of this method is used in many natural language processing tasks like statistical word alignment, part-of-speech tagging. The main idea is to constrain on the distribution of posterior, when computing expectation step in EM algorithm.

The prior distributions for the magnitude of the noise spectral components are modeled as Rayleigh probability density function (PDF) with scale parameter $\sigma$ , which is fitted to the observations by a maximum likelihood procedure [16, 17], i.e.,

$$f\left(x;\sigma \right)=\frac{x}{{\sigma }^2}{e}^{-x^2/2{\sigma }^2}\text{for}x\ge 0\text{with}{\sigma }^2=\frac{1}{2N}\sum _{i=1}^N{x}_i^2$$ E40

The above equation can be written as

$$f\left(x;\sigma \right)=e^{\mathit{log}\left(\frac{x}{{\sigma }^2}\right)}{e}^{-x^2/2{\sigma }^2}=e^{\mathit{log}\left(\frac{x}{{\sigma }^2}\right)-\frac{x^2}{2{\sigma }^2}}$$ E41

By applying negative logarithm on both sides of (41), we will get

$$-\mathit{log}\left(f\left(x;\sigma \right)\right)=-\mathit{log}\left(e^{\mathit{log}\left(\frac{x}{{\sigma }^2}\right)-\frac{x^2}{2{\sigma }^2}}\right)=\frac{x^2}{2{\sigma }^2}-\mathit{log}\left(\frac{x}{{\sigma }^2}\right)$$ E42

Then, the regularization term for the noise is defined as

$${\mathrm{\Lambda }}_N\equiv {\mathrm{\Lambda }}_{S_1}=-\text{log}f\left(x;\sigma \right)=\frac{x^2}{2{\sigma }^2}-\text{log}\frac{x}{{\sigma }^2}.$$ E43

The spectral components of speech modeled as Gamma probability density function [16, 18]

$$f\left(x;k,\theta \right)=\frac{x^{k-1}{e}^{\frac{-x}{\theta }}}{{\theta }^k\Gamma \left(k\right)}$$ E44

with shape parameter $k>0$ and scale parameter $\theta >0$ :

$$\theta =\frac{1}{\mathit{kN}}\sum _{i=1}^N{x}_i \text{and} k\approx \frac{3-s+\sqrt{{\left(s-3\right)}^2+24s}}{12s}$$ E45

where the auxiliary variable s is defined as $s=\mathit{ln}\left(\frac{1}{N}{\sum }_{i=1}^N{x}_i\right)-\frac{1}{N}{\sum }_{i=1}^N\mathit{ln}\left(x_i\right)$ .

The regularization term for the speech samples is defined as (by applying negative logarithm in both sides of (44))

$${\mathrm{\Lambda }}_{\mathit{SP}}\equiv {\mathrm{\Lambda }}_{S_2}=-\text{log}f\left(x;k,\theta \right)=\frac{x}{\theta }-\text{log}\left(\frac{x^{k-1}}{{\theta }^k\mathrm{\Gamma }\left(k\right)}\right),x\ge 0$$ E46

Special case: When we fix k = 1, the Gamma density simplifies to the exponential density and

$$f\left(x;1,\theta \right)=\frac{1}{\theta }{e}^{\frac{-x}{\theta }},{\mathrm{\Lambda }}_{\mathit{SP}}\equiv {\mathrm{\Lambda }}_{S_2}=\frac{x}{\theta },x\ge 0$$ E47

The proposed multiplicative nonnegative matrix factorization method is summarized in Algorithm 4 [16]. In general, like in the specific case of Algorithm 4, one can only guarantee the monotonous descent of the iteration through a majorization-minimization approach [19] or the convergence to a stationary point [20].

The subscript(s) with parenthesis represents corresponding columns or rows of the matrix assigned to a given source. 1 is an approximately sized matrices of ones, and $\otimes$ represents element-wise multiplication.

Algorithm 4: Gamma-Rayleigh regularized PLCA method (GR-NMF)

Procedure

$X\in R_{+}^{f\times t}$ % Observed normalized data

${\Lambda }_S\in R_{+}^{f\times t}$ , $s\in \left\{1,..\dots ,N_S\right\}$ % ${\Lambda }_S$ -Penalties, $N_S$ -Number of sources

$${\Lambda }_{s\left(\mathit{NEW}\right)}=0$$

$${\Lambda }_{S_1}={\Lambda }_{N\left(\mathit{OLD}\right)}=\frac{X^2}{2{\sigma }^2}-\mathit{log}\frac{X}{{\sigma }^2}\text{and}{\Lambda }_{S_2}={\Lambda }_{\mathit{SP}\left(\mathit{OLD}\right)}=\frac{X}{\theta }-\mathit{log}\left(\frac{X^{k-1}}{{\theta }^k\Gamma \left(k\right)}\right)$$ E48

$${\stackrel{\sim }{\Lambda }}_{s\left(\mathit{OLD}\right)}\leftarrow \mathit{exp}\left\{-{\Lambda }_s\right\}$$ E49

Repeat

For all s do

$${\stackrel{\sim }{\Lambda }}_s=\left(1-\mu \right){\Lambda }_{s\left(\mathit{OLD}\right)}+\mu {\Lambda }_{s\left(\mathit{NEW}\right)}\%\text{Update penalties using}\mathrm{LMS}$$ E50

$${\Lambda }_{s\left(\mathit{OLD}\right)}={\Lambda }_{s\left(\mathit{NEW}\right)}$$ E51

$$X_s\leftarrow X\otimes {\stackrel{\sim }{\Lambda }}_S$$ E52

End for

$$\Gamma \leftarrow {\sum }_s\left(W_{\left(s\right)}{H}_{\left(s\right)}\right)\otimes {\stackrel{\sim }{\Lambda }}_S$$ E53

$$Z_s\leftarrow \frac{X_s}{\Gamma }$$ E54

For all s do

$$W_{\left(s\right)}\leftarrow W_{\left(s\right)}\otimes \frac{Z_s{H}_{\left(s\right)}^T}{1{\mathrm{H}}_{\left(s\right)}^T}$$ E55

$$H_{\left(s\right)}\leftarrow H_{\left(s\right)}\otimes \left(W_{\left(s\right)}^T{Z}_s\right)$$ E56

End for

Reconstruction

For all s do

$$M_s\leftarrow \frac{W_{\left(s\right)}{H}_{\left(s\right)}}{\mathrm{WH}}\%\text{Compute Filter}$$ E57

$${\widehat{X}}_s\leftarrow M_s\otimes X\%\text{Filter Mixture}$$ E58

$$x_s\leftarrow \mathit{\text{ISTFT}}\left({\widehat{X}}_s,\angle \mathrm{X},P\right)\%\text{P- STFT parameters}$$ E59

if update k % Gamma model

$$s=\mathit{ln}\left(\frac{1}{N}\sum {\widehat{X}}_s\right)-\frac{1}{N}\sum \mathit{ln}\left({\widehat{X}}_s\right),k\approx \frac{3-s+\sqrt{{\left(s-3\right)}^2+24s}}{12s}$$ E60

else % Exponential model

k = 1,

end

$$\theta =\frac{1}{\mathit{kN}}\sum {\widehat{X}}_s$$ E61

$${\Lambda }_{S_1}={\Lambda }_{N\left(\mathit{OLD}\right)}=\frac{{\widehat{X}}_{s_1}^2}{2{\sigma }^2}-\mathit{log}\frac{{\widehat{X}}_{s_1}}{{\sigma }^2}$$

$${\Lambda }_{S_2}={\Lambda }_{\mathit{SP}\left(\mathit{OLD}\right)}=\frac{{\widehat{X}}_{s_2}}{\theta }-\mathit{log}\left(\frac{{\widehat{X}}_{s_2}^{\mathrm{k}-1}}{{\theta }^k\Gamma \left(k\right)}\right)$$ E62

$${\Lambda }_{s\left(\mathit{NEW}\right)}=\mathit{exp}\left(-{\mathrm{\Lambda }}_{s\left(\mathit{OLD}\right)}\right)\%{\Lambda }_{s\left(\mathit{OLD}\right)}\text{represents both}{\Lambda }_{\mathit{SP}}\text{and}{\Lambda }_N$$ E63

End for

Until Convergence

Return: Time domain signals $x_s$

5. Experimental results

The speech and noise audio samples were taken from NOIZEUS [21]. Sampling frequency is 8 KHz. The algorithm is iterated until convergence [16]. The proposed method was compared with Euclidean NMF (EUC-NMF) [5], Itakura-Saito NMF (IS-NMF) [22], posterior regularization NMF (PR-NMF) [15], Wiener filtering [23], and constrained version of NMF (CNMF)[24]. These methods are implemented by considering nonstationary noise, babble noise and street noise. The performance of proposed method was evaluated by using perceptual evaluation of speech quality (PESQ) [25] and source-to-distortion ratio (SDR) [26]. SDR gives the average quality of separation on dB scale and considers signal distortion as well as noise distortion. For PESQ and SDR, the higher value indicates the better performance. Tables 1 and 2 show the PESQ and SDR values of different NMF algorithms evaluated. The experimental results show that proposed method performs better than other existing methods in terms of the PESQ and SDR indices.

Table 1.

PESQ and SDR for babble noise.

Table 2.

PESQ and SDR for street noise.

6. Conclusion

A novel speech enhancement method based on an iterative and regularized NMF algorithm for single-channel source separation is proposed. The clean speech and noise magnitude spectra are modeled as Gamma and Rayleigh distributions, respectively. The corresponding log-likelihood functions are used as penalties to regularize the cost function of the NMF. The estimation of basis matrices and excitation matrices are calculated by using proposed regularization of multiplicative update rules. The experiments reveal that the proposed speech enhancement method outperforms other existing benchmark methods in terms of SDR and PESQ values.