A unified approach for underdetermined blind signal separation and source
activity detection by multichannel factorial hidden Markov models
Takuya Higuchi1 , Hirofumi Takeda1 , Tomohiko Nakamura1 and Hirokazu Kameoka1,2
1

Graduate School of Information Science and Technology, The University of Tokyo
2
NTT Communication Science Laboratories, NTT Corporation{higuchi, h-takeda, nakamura, kameoka}@hil.t.u-tokyo.ac.jp

Abstract
This paper proposes to introduce a new model called “the multichannel factorial hidden Markov Model (MFHMM)” for underdetermined blind signal separation (BSS). For monaural source
separation, one successful approach involves applying nonnegative matrix factorization (NMF) to the magnitude spectrogram of a mixture signal, interpreted as a non-negative matrix.
Up to now, multichannel extensions of NMF, which allow for
the use of spatial information as an additional clue for source
separation, have been proposed by several authors and provento be an effective approach for underdetermined BSS. This approach is based on the assumption that an observed signal is a
mixture of a limited number of source signals each of which
has a static power spectral density scaled by a time-varying amplitude. However, many source signals in real world are nonstationary in nature and the variations of the spectral densities
are much richer in time. Moreover, many sources including
speech tend to stay inactive for some while until they switch toan active mode, implying that the total power of a source may
depend on its underlying state. To reasonably characterize sucha non-stationary nature of source signals, this paper proposes toextend the multichannel NMF model by modeling the transitionof the set consisting of the spectral densities and the total power of each source using a hidden Markov model (HMM). By
letting each HMM contain states corresponding to active and inactive modes, we will show that voice activity detection andsource separation can be solved simultaneously through parameter inference of the present model. The experiment showedthat the proposed algorithm provided a 7.65 dB improvementcompared with the conventional multichannel NMF in terms of
the signal-to-distortion ratio.
Index Terms: blind signal separation, source activity detection,a hidden Markov model, non-negative matrix factorization

1. Introduction
Blind source separation (BSS) refers to a technique for separating out individual source signals from microphone array inputs when the transfer characteristics between the sources andmicrophones are unknown. The best-known commercial application of BSS techniques is their use in teleconferencing systems. To solve BSS problems, it is generally necessary to makesome assumptions about the sources, and formulate an appropriate optimization problem based on criteria designed according to those assumptions. For example, if the observed signals
outnumber the sources, we can employ independent component
analysis (ICA) [1] by assuming that the sources are statistically independent of each other. However, in an underdeterminedcase, the independence assumption is too weak to allow us todetermine a unique solution and so directly applying ICA will
not work well.
For monaural source separation, one successful approach
involves applying non-negative matrix factorization (NMF) tothe magnitude spectrogram of a mixture signal, interpreted asa non-negative matrix [2, 3]. With this approach, the spectrogram of a mixture signal is factorized into the product of a basis matrix consisting of basis spectra and an activation matrix
consisting of time-varying amplitudes associated with the basis

spectra. An important feature of this approach is that it is capable of finding a finite set of basis spectra that are considered
to be the dominant elements constituting the observed spectrogram, in an unsupervised manner. Up to now, several attempts
have been made to extend this approach to a multichannel casein order to allow for the use of spatial information as an additional clue for separation, which have opened a door to a newpromising approach for underdetermined BSS [4, 5]. This approach is based on the assumption that an observed signal is a
mixture of a limited number of source signals each of which has
a static power spectral density (i.e., the basis spectrum) scaled
by a time-varying amplitude. However, many source signals
in real world are non-stationary in nature and the variations of
the spectral densities are much richer in time. For example, thesound of a piano note would be more accurately characterized
by a succession of several basis spectra corresponding to “attack,” “decay,” “sustain” and “release” segments than only bya single basis spectrum. Another important fact is that manysources including speech tend to stay inactive for some while
until they switch to an active mode. This implies that the total power of a source may depend on its underlying state. To
reasonably characterize such a non-stationary nature of source
signals, this paper proposes to extend the multichannel NMF
model by modeling the transition of the set consisting of thespectral densities and the total power of each source using ahidden Markov model (HMM). With this model, we would be
able to flexibly reflect the time-varying nature of each source
by appropriately specifying or training the transition probabilities prior to analysis. We formulate an entire generative modelof a multichannel mixture signal by incorporating a generativemodel of the spectral densities of each source, described using
an HMM. Each hidden state of the HMM represents the state of
the corresponding source. We call this model “the multichannel
factorial hidden Markov model (MFHMM)”.
In general, simultaneous estimation is preferable when several estimation problems are interdependent. If we knew wheneach source is active and inactive, source separation would become a relatively simple matter. On the other hand, if all thesources were already separated, voice activity detection would
become a relatively simple matter. This simply implies that the
problems of source separation and voice/source activity detection are interdependent of each other. It is important to note thatthrough parameter inference of the present model, we would
be able to simultaneously perform voice activity detection andsource separation based on a unified maximum likelihood criterion.
The remainder of this paper is organized as follows. Sec. 2
formulates a generative model of a multichannel observed signal and source signals based on the factorial HMM, Sec. 3 describes a parameter inference algorithm for the present modeland Sec. 4 presents some experimental results.

2. Multichannel factorial HMM
2.1. Mixing modelFirst we consider a situation where I source signals are recorded
by M microphones. Here, let ym (ωk , tl ) ∈ C be the shorttime Fourier transform (STFT) component observed at the m-th
microphone, and si (ωk , tl ) ∈ C be the STFT component of thei-th source. 1 ≤ k ≤ K and 1 ≤ l ≤ L are the frequency andtime indices, respectively. If we assume that the length of theimpulse response from a source to microphones is sufficiently
shorter than the frame length of the STFT, the observed signal
can be approximated fairly well by an instantaneous mixture in
the frequency domain:y(ωk , tl ) =I∑ai (ωk )si (ωk , tl ),

(1)

i=1

where y(ωk , tl ) = (y1 (ωk , tl ), . . . , yM (ωk , tl ))T ∈ CM .
ai (ωk ) denotes the frequency array response for source i at
frequency ωk . For convenience of notation, we hereafter use
subscripts k and l to indicate ωk and tl respectively.
2.2. Generative process of observed signals
Here we describe the generative process of an observed signal
based on Eq. (1). If we assume that each source signal follows
a piecewise stationary Gaussian process, then si,k,l follows a2complex normal distribution with mean 0 and covariance σi,k,l,2
2),
si,k,l |σi,k,l∼ NC (si,k,l ; 0, σi,k,l(2)2where σi,k,ldenotes the power spectral density of i-th source at
frequency k and time l. From Eq. (1) and Eq. (2), y k,l is also
normally distributed such that∑2
y k,l |a1:I,k , σ1:I,k,l ∼ NC (y k,l ; 0,C i,k σi,k,l), (3)i
2
conditioned on a1:I,k = {a1,k , . . . , aI,k }, σ1:I,k,l=22{σ1,k,l , . . . , σI,k,l } where C i,k=ai,k aHandi,k
NC (x; µ, Σ) ∝ exp(−(x − µ)H Σ−1 (x − µ)). C i,krepresents a spatial correlation matrix of the i-th source at
frequency k.

2.3. Generative model for multichannel NMF [4, 5]
In the regular NMF model (see [6]), the power spectra of asource signal is assumed to be static up to a scale factor. We
can incorporate this assumption into the above model by setting
2
σi,k,l= wi,k hi,l ,(4)2
where σi,k,l
is assumed to be factorized into the product of thestatic power spectrum wi,k and the time-varying amplitude hi,l .
The generative model of si,k,l is thus rewritten as

si,k,l |wi,k , hi,l ∼ NC (si,k,l ; 0, wi,k hi,l ),(5)conditioned on wi,k and hi,l . Since the generative model of
Y = {y k,l }k,l under the assumption Eq. (5) can be viewed
as a natural extension of NMF to a multichannel case, a BSS
approach based on this model is called the multichannel NMF[4, 5].
2.4. Generative modeling of source signals using HMMs
As described above, the multichannel NMF model roughly assumes that the power spectra of each sound source is static up
to a scale factor. However, many sound sources exhibit different spectra according to underlying “states” of the sources.
For example, the spectra of the sound of a piano note would
be different in “attack,” “decay,” “sustain” and “release” segments. Another important fact is that many sources including
speech tend to stay inactive for some while until they switch toan active mode. This implies that the total power of a source
also depends on its underlying state. To reasonably characterizesuch a non-stationary nature of source signals, here we propose
to model the sequence of the power spectra and the total powersof each source using an HMM.

Figure 1: Illustration of the concept of the FHMM.
Now we introduce a latent variable zi,l ∈ {1, . . . , D} todenote a state of the i-th source at time l. The state sequencezi,1 , . . . , zi,L is assumed to follow a Markov chain:zi,l |zi,l−1 ∼ Categorical(zi,l ; ρzi,l−1 ),(6)where Categorical(x; y) = yx , ρd = (ρd,1 , . . . , ρd,D ) denotes the transition probability of state d to each state 1, . . . , D,and ρ = (ρd,d′ )D×D denotes the transition matrix. Here, we
assume that hi,l follows a gamma distribution with hyperparameters determined according to zi,l ,
hi,l |zi,l ∼ Gamma(hi,l ; αzi,l βzi,l ),(7)where α1:D and β1:D are the shape and scale parameters of aα−1 −x/β
e
gamma distribution, and Gamma(x; α, β) = x Γ(α)β. As
αwe want hi,l to take a small value when zi,l is the “inactive”(i.e., silent) state, we set the hyperparameters of the gamma distribution of that state so that it becomes a sparsity-inducing distribution. As regards the gamma distributions of the remaining
states, we consider setting the hyperparameters so that they become uniform distributions. In particular, we consider setting
α and β at 1 and 10−2 respectively for the inactive state, andat 1 and 1020 respectively for the active states. We expect that
this setting allows us to solve source separation and voice activity detection in a cooperative manner. Let us use wi,k,d todenote the power spectrum of the i-th source at state d. The
power spectrum of the i-th source at time l is also assumed to
be determined according to zi,l . Thus, the generative model ofsi,k,l is eventually written as
si,k,l |wi,k,1:D , hi,l , zi,l∼ NC (si,k,l ; 0, wi,k,zi,l hi,l ).

(8)

Since the generative of y k,l contains multiple HMMs associatedwith the underlying sources, the overall model can be viewed asa Factorial HMM. Fig. 1 shows an illustration of the proposed
FHMM. Note that we can integrate an assumption that the signal tend to stay non-active for some time by setting the transition probabilities of self-transitions at reasonably large values
in a specific state of the HMM. If we assume the emission probabilities of the HMMs as uniform distributions, our overall generative model is given by Eqs. (6) andy k,l |a1:I,k , w1:I,k,1:D , h1:I,l , z1:I,l∑∼ NC (y k,l ; 0,C i,k wi,k,zi,l hi,l ),(9)iconditioned on a1:I,k ,w1:I,k,1:D ,h1:I,l and z1:I,l .

3. Algorithm for parameter estimation
3.1. Objective function
In this section, we describe a parameter estimation algorithm for our generative model based on an auxiliary function
method. The random variables of interest in our model are
W = w1:I,1:K,1:D , H = h1:I,1:L ,C = C 1:I,1:K and Z =z1:I,1:L . We denote the entire set of the above parameters as Θ.
In the following, ρ is constants that is determined experimentally. Our goal is to compute the posteriorp(Θ|Y ) =p(Y , Θ)
,
p(Y )

(10)

where Y = y 1:K,1:L is a set consisting of the time-frequency
components of observed multichannel signals. By using theconditional distributions defined in Sec. 2, we can write thejoint distribution p(Y , Θ) asp(Y , Θ) ∝p(Y |Θ)p(H|Z)p(Z).

(11)

The objective function is defined as L(Θ) = log p(Θ|Y ). Our
goal is to obtain Θ̂ such thatΘ̂ = argmax log p(Θ|Y ).

Figure 2: A spectrogram of the mixed signal.
set of the auxiliary variables as Λ. tr(·) is the trace of a matrix.
The equality L(Θ) = L+ (Θ, Λ) is satisfied when

(12)

Θ

−1

By using Eqs. (10), (11) and (12), the current optimization
problem can be rewritten as(Θ̂ = argmax log p(Y |Θ)Θ)+ log p(H|Z) + log p(Z) .
(13)
According to the generative model defined in Sec.
log p(Y |Θ) is written as

2,

log p(Y |Θ)
∑ M
−1
1
1
=
log 2π − log |X̂ k,l | − y k,l H X̂ k,l y k,l ),
(−
2
2
2
k,l

(14)
where X̂ k,l =∑i C i,k wi,k,zi,l hi,l .

3.2. Optimization algorithm based on an auxiliary function
methodThe optimization problem of maximizing L(Θ) with respect to Θ is difficult to solve analytically. However, we can invoke an auxiliary function approach to derive an iterative algorithm that searches for the estimate of Θ, as with [5]. To
apply an auxiliary function approach to the current optimization problem, the first step is to construct an auxiliary function
L+ (Θ, Λ) satisfying L(Θ) = maxΛ L+ (Θ, Λ). We refer to Λas an auxiliary variable. It can then be shown that L(Θ) is nonincreasing under the updates Θ ← argmaxΘ L+ (Θ, Λ) andΛ ← argmaxΛ L+ (Θ, Λ). The proof of this shall be omitted
owing to space limitations. Thus, L+ (Θ, Λ) should be designed
as a function that can be maximized analytically with respect toΘ and Λ. Such a function can be constructed as follows.
L(Θ)≥L+ (Θ, Λ)( (−1
H1 ∑ ∑ tr(y k,l y k,l Ri,k,l C i,k Ri,k,l )=−2
wi,k,zi,l hi,l
i
k,l
)
)
h
+
log
|U
|
−
M
+ tr(U −1
C
)w
k,l
i,k
i,k,zi,l i,l
k,l
)
∑(
+
(αzi,l − 1) log hi,l − hi,l /βzi,l − αzi,l log βzi,li,l+ log p(Z),(15)

where Ri,k,l and U k,l are auxiliary∑ variables that satisfy Hermitian positive definiteness and i Ri,k,l = I. We denote theRi,k,l = C i,k wi,k,zi,l hi,l X̂ k,l ,(16)

U k,l = X̂ k,l .

(17)

Therefore, we can monotonically increase L by repeating
the following two steps.
1. Maximizing L+ with respect to R and U .
2. Maximizing L+ with respect to W , H, C and Z.
Step 1 consists in updating R and U using Eqs. (16) and(17). In step 2, we can obtain update rules of W , H, C by
setting the partial derivative of L+ with respect to each of theparameters at zero. The partial derivatives of L+ with respect
to W and H are given by∑( tr(y k,l H y k,l Ri,k,l C −1
∂L+i,k Ri,k,l )=2∂wi,k,zi,lwi,k,zi,l hi,ll)− tr(U −1
C
)h
(18)
i,k
i,l ,
k,l
∑( tr(y k,l H y k,l Ri,k,l C −1
∂L+i,k Ri,k,l )
=
2
∂hi,l
wi,k,zi,l hi,lk)
− tr(U −1C)wi,ki,k,zi,lk,l+ (αzi,l − 1)/hi,l − 1/βzi,l ,(19)respectively. By setting them at zero, we obtain the following
update rules:v
u ∑ tr(y H y R−1uk,lk,l i,k,l C i,k Ri,k,l )
u l
h
i,l
wi,k,zi,l ← t
,
(20)
∑
−1
l tr(U k,l C i,k )hi,l√(αzi,l − 1) + (αzi,l − 1) + 4µi,l νi,l
, (21)
hi,l ←
2µi,l
whereµi,l =∑ tr(y k,l H y k,l Ri,k,l C −1
i,k Ri,k,l )kνi,l =∑,(22)

tr(U −1
k,l C i,k )wi,k,zi,l + 1/βzi,l .

(23)

wi,k,zi,lk

The partial derivatives of L+ with respect to C is given by−1
H∑( C −1∂L+i,k Ri,k,l y k,l y k,l Ri,k,l C i,k=∂C i,kwi,k,zi,l hi,llFigure 3: Examples of spectrogram of signal of a bell (top),that of separated signal (middle), and its activity detection result(bottom). Black indicates the state is assigned at the time.
Table 1: The average output SDRs and standard deviations of
the three sources by the conventional and proposed methods.
SDR(±SD) [dB]
Proposed
Conventionalbell
19.92(±11.50)
13.33(±8.22)

whistle31.21(±6.69)
23.37(±4.57)

stapler
16.81(±10.09)
8.28(±5.98))− U −1w
hk,l i,k,zi,l i,l .

(24)

By setting this at zero, we obtain an algebraic Riccati equation:
C i,k Ai,k C i,k = B i,k ,where
Ai,k =∑−1

wi,k,zi,l hi,l X̂ k,l , B i,k(25)∑= C i,k ( wi,k,zi,l hi,ll
−1
−1X̂ k,l y k,l H y k,l X̂ k,l )C i,k .

l

(26)

We can solve this equation by using a method in [5]. We perform an eigenvalue decomposition of a 2M × 2M matrix
[
]
0
−Ai,k
,
(27)
−B
0
i,k

and let e1,i,k . . . eM,i,k be eigenvectors with negative eigenvalues. It is guaranteed that there are exactly M negative eigenvalues. Then, let us decompose the 2M -dimensional eigenvectors
as[]f
em,i,k = g m,i,k,
(28)
m,i,k

for m = 1 . . . M where f m,i,k and g m,i,k are M -dimensional
vectors. We obtain the update rule for C i,k as
C i,k ← Gi,k F −1
i,k ,

(29)

where F i,k=[f 1,i,k , . . . , f M,i,k ] and Gi,k=[g 1,i,k , . . . , g M,i,k ].
L+ is equal up to constant terms to the sum of the log posteriors of HMMs, when viewed as a function of Z. Thus, we
can invoke the Viterbi algorithm to search for the optimal pathzi,1 , . . . , zi,L for each i individually. Note that updating W ,H, C and Z corresponds to solving the problems of blind signal separation and source activity detection based on a unified
objective function.

4. Experimental evaluation
We evaluated the performance of the present method in termsof the ability to separate sources and detect their activity. We
used a mixed stereo signal as the experimental data, each ofwhich we obtained by mixing the non-speech signals (soundsof a whistle, a bell and a stapler) from the RWCP database[8]and was convolved with the measured room impulse response
from the RWCP database [8] (in which the distance between
the microphones was 5.85 cm and the reverberation time was 0ms). Thus, the three signals were artificially located 20◦ , 60◦
and 100◦ from the microphones respectively. The sampling rate
was 32 kHz. To compute the STFT components of the observed
signal, the STFT frame length was set at 16 ms and a Hamming
window was used with an overlap length of 8 ms. Fig. 2 showsa spectrogram of the mixed signal. We set the number of statesof HMMs D as 3. We expected that d = 1 is an inactive state and d = 2 and 3 are active states, by setting α1 and β1
as 1 and 10−2 respectively, and α2:3 and β2:3 1 and 1020 respectively. the transition probability ρ as ρ1 = (0.9, 0.1, 0),
ρ2 = (0, 0.5, 0.5) and ρ3 = (0.5, 0, 0.5) based on the assumption that a sound source tends to stay inactive.
The diagonal ele√ments of C were initially all set to 1/ M , and the off-diagonal
elements were initially set to zero. We first performed the single channel NMF based on IS-divergence starting from different random initial matrices W and H, then we set the results of
W , H as the initial parameters of W and H. The parameter
estimation algorithm was run for 100 iterations. We chose themethod proposed in [5] as a comparison. The separated signalŷ i,k,l was obtained by Wiener filtering−1ŷ i,k,l = wi,k,zi,l hi,l C i,k X̂ k,l y k,l .

(30)

As an evaluation measure, we used the signal-to-distortion ratio(SDR) [9]. The SDR is expressed in decibels (dB), and a higher
SDR indicates superior quality. The input SDRs [dB] of thebell, the whistle and the stapler were -8.36, 6.64 and -12.78[dB], respectively.
Table 1 shows the average SDRs and standard deviations for the ten trials obtained by the conventional and proposed
methods. The average SDRs of the proposed method were superior to those of the conventional method for each signal. The
total average of the SDRs obtained with the proposed method
was 7.65 dB more than that obtained with the conventional approach. These results show the effectiveness of the proposed
method for BSS. Examples of separated signals are available
at http://www.hil.t.u-tokyo.ac.jp/˜higuchi/
demo/Examples_interspeech.htm. Fig. 3 shows examples of spectrogram of signal of a bell (top), that of separated
signal (middle), and its activity detection result (bottom). The
activities of the bell were almost all estimated correctly. However, our results also showed that abilities of signal separationand source activity detection by our method depended on thefirst condition of W and H to some extent. As future work, weplan to solve this kind of problem by learning W and ρ froma clean signal and apply our proposed method to acoustic event
detection.

5. Conclusion
In this paper we extend our previous approach to a multichannel observed signal and proposes a comprehensive approach todeal with an underdetermined BSS problem and source activity detection. Specifically, we describe a generative model of
mixture signals by incorporating a generative model of spectra for each source, using MFHMMs. Each hidden state of thepresent HMM represents states of each source such as its activities. Through the estimation of the parameters of the overall
generative model, we can simultaneously the underdeterminedblind signal separation and source activity detection. The experiment showed that the proposed algorithm provided a 7.65 dB
improvement compared with the conventional method in termsof the signal-to-distortion ratio.

6. Acknowledgement
This work was supported by JSPS KAKENHI Grant Number
26730100.

7. References[1] A. Hyvärinen, J. Karhunen, and E. Oja, Independent Component
Analysis, John Wiley & Sons, 2001.
[2] D. D. Lee, and H. S. Seung, “Learning the parts of objects with
nonnegative matrix factorization,” Nature, vol. 401, pp.788–791,
1999.
[3] P. Smaragdis, and J. C. Brown, “Non-negative matrix factorization for polyphonic music transcription,” Proc. WASPAA 2003,
Oct. 2003, pp. 177–180.
[4] A. Ozerov, and C. Févotte, “Multichannel nonnegative matrix
factorization in convolutive mixtures for audio source separation,”
IEEE Trans. Audio, Speech and Language Processing, vol. 18, no.
3, pp. 550-563, Mar.2010.
[5] H. Sawada, H. Kameoka, S. Araki and N. Ueda, “Efficient algorithms for multichannel extensions of Itakura-Saito nonnegative
matrix factorization,” IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP) 2012, pp. 261-264,
2012.
[6] C. Févotte, N. Bertin, and J.L. Durrieu, “Nonnegative matrix factorization with the Itakura-Saito divergence: With application tomusic analysis,” Neural computation, vol. 21, no. 3, pp. 793–830,
2009.
[7] Masahiro Nakano, Le Roux Jonathan, Hirokazu Kameoka, Yu Kitano, Nobutaka Ono, and Shigeki Sagayama, “Nonnegative matrix factorization with Markov-chained bases for modeling timevarying patterns in music spectrograms,” Latent Variable Analysis
and Signal Separation, vol. 6365, pp. 149–156, 2010.
[8] S. Nakamura, K. Hiyane, F. Asano, T. Nishiura, and T. Yamada,“Acoustical sound database in real environments for sound scene
understanding and hands-free speech recognition,” in Proc. 2nd
International Conference on Language Resources & Evaluation(LREC 2000), pp. 965–968, 2000.
[9] E. Vincent, R. Gribonval, and C. Févotte, “Performance measurement in blind audio source separation,” IEEE Transactions on
Audio, Speech, and Language Processing, pp. 1462–1469, 2006.