Statistical Model of Speech Signals Basedon Composite Autoregressive System with
Application to Blind Source Separation
Hirokazu Kameoka, Takuya Yoshioka, Mariko Hamamura,
Jonathan Le Roux, and Kunio Kashino
NTT Communication Science Laboratories, NTT Corporation,
3-1 Morinosato Wakamiya, Atsugi, Kanagawa 243-0198, Japan{kameoka,leroux,kunio}@cs.brl.ntt.co.jp, takuya@cslab.kecl.ntt.co.jp

Abstract. This paper presents a new statistical model for speech signals, which consists of a time-invariant dictionary incorporating a set ofthe power spectral densities of excitation signals and a set of all-polefilters where the gain of each pair of excitation and filter elements is allowed to vary over time. We use this model to develop a combined blindseparation and dereverberation method for speech. Reasonably good separations were obtained under a highly reverberant condition.
Keywords: Blind source separation, composite autoregressive system.

1

Introduction

When a set of observed data is considered to be a random sample drawn from a
population that follows a particular family of probability distributions, we often
refer to the distribution family as a statistical model. In this paper, we present
a new statistical model of speech signals, suitable for applying to blind signal
processing problems.
The aim of blind source separation (BSS) or blind dereverberation is to detect each source component from observed signals captured by one or several
microphones without using any information about the transfer characteristic of
the path from each source to each microphone. We thus usually make some assumption about the sources and then formulate an optimization problem basedon a criterion that measures the consistency with the assumption. One approach
to formulating the problem would be to deﬁne a likelihood function of a source
signal by employing a statistical model assumption. When choosing which statistical model to invoke, it is important to take account of whether it agrees
well with the actual behaviors of the source of interest and whether it leads toa mathematically tractable form of the optimization problem.
Let us now brieﬂy review the BSS problem. BSS algorithms derived from a
convolutive mixture model in the time domain such as [1] are ﬁne for short mixing ﬁlters, but when it comes to realistically long ﬁlters, they can be unrealizable
because of computational requirements. In the STFT domain where the frame
V. Vigneron et al. (Eds.): LVA/ICA 2010, LNCS 6365, pp. 245–253, 2010.
c Springer-Verlag Berlin Heidelberg 2010


246

H. Kameoka et al.

size is assumed to be suﬃciently larger than the ﬁlter length, a convolutive mixture signal can be approximated by an instantaneous mixture in the frequency
domain, thus allowing for an eﬃcient implementation of BSS algorithms [2]. However, when reverberation comes into play, the ﬁlter length can often be larger
than the frame size and so the above approximation becomes relatively less accurate. According to several studies such as [3], reverberation can be modeled
fairly well as a convolution for each frequency-band in the STFT domain. Therefore, under highly reverberant conditions, we can expect a convolutive mixture
model in the time-frequency domain [4] to be a better approximation.
We consider the situation where we observe signals emanating from M sourceslbe the STFT of the signal observedand captured by M microphones. Let Yk,n1
M T
) be a set of observed data,at the l-th microphone and Y k,n = (Yk,n , · · · , Yk,nwhere k and n are the frequency and time indices, respectively. We use S k,n =
1
M T
(Sk,n
, · · · , Sk,n
) to denote a set of M source components. In this paper, we use
a separation system that has a multichannel ﬁnite-impulse-response form in thetime-frequency domain such that:
S k,n =nk

W k,τ Y k,n−τ ,

(1)

τ =0

where W k,τ , 0 ≤ τ ≤ nk are matrix coeﬃcients of size M ×M . Let us assume that

m
m
Sk,n
is a random variable and that Sk,nand Skm ,n are statistically independentwhen (k, n, m) = (k  , n , m ). The deﬁnition of the probability density functionmmmm |θ m (s(pdf) for Sk,n, denoted by fSk,nk,n |θ ), is the main subject of the following
section. This pdf corresponds to a statistical model of the m-th source signal
and θm is the set of parameters characterizing its distribution. Once we deﬁne
its speciﬁc form, the joint pdf of Y := {Y k,n }k,n can be written in terms ofmm |θ m ( · |θfSk,n
) explicitly as
fY |Θ (Y |Θ) =k| det Wk,0 |
t,m

 m m
m |θ m S̄
,
fSk,n
k,n|θ

(2)


m
where S̄k,nis the m-th element of the separated signal vector, given by τ W k,τ
Y k,n−τ . The joint pdf fY |Θ (Y |Θ) is the likelihood of the unknown variables
Θ := {{θm }m , {W k,τ }k,τ } given observation Y , which is an objective function
for achieving separation and dereverberation in a joint manner, as with [4].
It has been shown that the speech production system can be well modeled
on a frame-by-frame basis by a linear system comprising a glottal excitation
input and a vocal-tract resonance ﬁlter that respectively determine the degreeof periodicity and the phoneme of the voice. One of the most frequently usedmodels for short-term speech signals is the autoregressive (AR) model, whichmodels the signal as the output of an all-pole system. [5] was among the ﬁrst
to propose a BSS system in which a Gaussian AR source model is incorporatedmm
m |θ (sin fSk,nk,n |θ ), where a complex spectrum at each frame is assumed to be a
set of random samples drawn from a diﬀerent AR system with a Gaussian white
noise input. This type of statistical model for STFT spectrograms of speech has

Statistical Model of Speech Signals

247

later been shown to work successfully for BSS in highly reverberant environments[4]. The objective of this paper is to investigate the possibility to improve theperformance of this state-of-the-art BSS system [4] by replacing the Gaussian
AR source model with the speech model we proposed previously [6], which will
be reviewed in the next section.

2

Proposed Statistical Model of Speech

The white noise assumption as regards the excitation inputs underlying the standard AR model is known to not hold especially for voices with high fundamental
frequencies (F0 s). This is because when F0 increases the spacing between theharmonics of the excitation spectrum increases correspondingly, thus departingfrom a white (ﬂat) spectrum. Rather than ﬁxing the characteristics of the excitation input, we would therefore, if possible, like to estimate them in the same way
as the vocal-tract characteristics. However, if we simply treat both of the characteristics as separate variables for each frame, it would no longer be possible
to determine these characteristics uniquely since there are an inﬁnite number of
combinations giving the same ﬁlter output. Some additional constraint is needed
to avoid this indeterminacy.
There are a wide variety of regularities in speech that can be exploited toconstrain speech models. For example, it may be reasonable to assume that every short-term signal of a speech utterance can be represented by a combinationof elements drawn from two dictionaries, one consisting of a small number of
vocal-tract characteristics, and the other of a small number of excitation characteristics. This is because the phoneme number and periodicity range of speechduring an entire utterance are both usually limited. We thus assume here thatspeech signals have been generated by a compound linear system composed ofthe direct product of a limited set of time-invariant excitation characteristics anda limited set of time-invariant vocal-tract characteristics where each output associated with an excitation and ﬁlter element pair is activated by a time-varying
gain. A signal at a particular frame is thus assumed to be characterized by thevolume levels of all the ﬁlter outputs. Note that since the dictionary elements
are in general unknown, they need to be estimated in a data-driven manner.
In this section, we focus on a particular source and so the index m will beomitted for simplicity of notation. We start by modeling an output signal characterized by a particular excitation and ﬁlter element pair. Let us assume thata signal in the n-th frame, {xn [t]}K
t=1 , is a sampled sequence drawn from the
P -order AR process with a set of AR parameters common over n such thatxn [t] =P

ap xn [t − p] + n [t],

(3)

p=1

where n [t] is an excitation input signal that is assumed to be a zero-mean
stationary Gaussian noise. Its autocorrelation function, νn [t], is constant up to
the gain over all the frames such that νn [t] = Un h[t]. Un is assumed to be the248

H. Kameoka et al.

energy of n [t] within the frame. We note here that n [t] is not restricted tobeing white noise. Let xn = (xn [1], · · · , xn [K])T ∈ RK and its discrete Fouriertransform (DFT) be X n = (X1,n , · · · , XK,n )T ∈ CK . Then, according to Eq. (3),distribution X n ∼ NC (0, Λn ), with
X n follows a multivariate complex Gaussiana diagonal covariance matrix Λn = diag λ1,n , · · · , λK,n whose elements are
Hk U n,|A(ej2πk/K )|2(4)

A(z) = 1 − a1 z −1 − · · · − aP z −P ,(5)λk,n =

where k is the frequency index and j is the imaginary unit. {Hk }K
k=1 is the DFTof {h[k]}K,
which
represents
the
power
spectraldensity(PSD)
ofthe excitation
k=1
source signal n [t] (namely, the spectral ﬁne structure), which can have any shapeand is not necessarily ﬂat. On the other hand, 1/|A(ej2πk/K )|2 corresponds to aspectral envelope expressed as the PSD of the all-pole transfer function.
We can now construct a “composite” autoregressive model by extending theabove model. The composite autoregressive system is assumed to consist of adictionary of I excitation PSDs and a dictionary of J all-pole ﬁlters. Subsequently, we use superscripts i and j to denote the indices of the excitation PSDs
and the all-pole ﬁlters, respectively, and we denote the ith excitation PSD andthe jth all-pole transfer function by Hki and 1/Aj (ej2πk/K ). The system is able
to generate I × J diﬀerent voice components each of which is characterized by
combining elements drawn from the respective dictionaries. If we assume that
only one of the I × J voice components is active at each frame, the source
pdf can be deﬁned using a Gaussian scaled mixture model [7]. By contrast, we
assume that all the voice components are always active with diﬀerent volumelevels. Let Uni,j denote the volume level of the {i, j}-th voice component at then-th frame. By following the discussion above, the DFT of the {i, j}-th voice
i,j
i,j T
component at the nth frame, X i,jn = (X1,n , · · · , XK,n ) , follows a multivariatei,jcomplex Gaussian
distribution
NC (0, Λn ) with a diagonal covariance matrix
 i,j
i,j 
Λi,j
n = diag λ1,n , · · · , λK,n whose elements are
Hki Uni,j,|Aj (ej2πk/K )|2(6)

Aj (z) = 1 − aj1 z −1 − · · · − ajP z −P .

(7)

λi,j
k,n =

I,J
If we now assume that X 1,1
n , · · · , X n are mutually independent, it follows that




X i,j
Λi,j
(8)
Sn =
n ∼ NC 0, Φn , Φn =
n ,
i

j

i

j

where S n ∈ CK denotes the DFT of the speech signal at the n-th frame. The
statistical model, fSk,n |θ (Sk,n |θ), is thus given concisely as1|Sk,n |2
fSk,n |θ (Sk,n |θ) =
exp −
,
(9)
πφk,n
φk,n

Statistical Model of Speech Signals

249

where θ contains all the unknown parameters of the present system: θ := Hki , ajp ,Uni,j i,j,p,k,n . The diagonal element of Φn , i.e. φk,n , corresponds to the PSD of
the output signal produced by the present system such that
φk,n =ij

Hki Uni,j.
|Aj (ej2πk/K )|2(10)

It is important to note that in the special case where J =
1 and P = 0, the PSD
φk,n reduces to the form of matrix multiplication φk,n = i Hki Uni,1 and Eq. (9)thus reduces to the likelihood function under the statistical interpretation of
non-negative matrix factorization (NMF) given in [8]. This fact suggests thatthe present entire BSS system is structurally related to multichannel NMF [9].

3

Optimization Process

Given a set of observed STFT components Y , we want to ﬁnd the estimate ofΘ = {{θm }m , {W k,τ }k,τ } that maximizes the posterior density fΘ|Y (Θ; Y ) ∝fY |Θ (Y ; Θ) fΘ (Θ), or equivalently, the log posterior density
L(θ) := log fY |Θ (Y |Θ) + log fΘ (Θ).

(11)

Eq. (11) can be iteratively increased by using a coordinate descent method inwhich each iteration comprises the following three maximization steps: (S1)θm ← argmaxθm L(θ) for all m, (S2) Wk,0 ← argmaxWk,0 L(θ) for all k,k← argmax{Wk,τ }nk L(θ) for all k. If we were able to oband (S3) {Wk,τ }nτ =1
τ =1tain an estimate of the PSD of each source, namely φmk,n , we could invoke [4]to perform (S2) and (S3). Therefore, obtaining the update formula of (S1) will
suﬃce to complete the derivation of the entire optimization process. It should
be noted that (S1) amounts to maximizingmmmm |θ m (S̄log fSk,n(12)k,n |θ ) + log fθ m (θ )k,nmis the m-th vector element of τ Wk,τ Y k,n−τ . Aswith respect to θm where S̄k,nthis maximization is carried out for each m separately, the index m is omitted
again in the following.
log fSk,n |θ (S̄k,n |θ) is equal up to constant terms to the goodness of ﬁt between|S̄k,n |2 and φk,n deﬁned by the Itakura-Saito divergence. We are thus led toobtain a PSD model with as small a reconstruction error as possible. On theother hand, as with the sparse coding concept [10], we would like to keep thevoice components as sparse as possible. The prior term log fθ (θ) can be used topromote the sparseness of Uni,j . In the subsequent analysis, for convenience we
use an exponential prior (a folded Laplacian prior) deﬁned over Uni,j ≥ 0fθ (θ) =α exp(−αUni,j ),
(13)
i,j,n

250

H. Kameoka et al.

which promotes sparsity when α is large. Maximizing Eq. (12) thus combines
the goals of a small reconstruction error and sparseness. As a consequence, themore frequently a certain spectral ﬁne/envelope structure emerges in |S̄k,n |2 , themore likely it is to be captured in the excitation/ﬁlter dictionary.
Although it is diﬃcult to obtain a closed-form solution for maximizing Eq. (12),we can develop a computationally eﬃcient scheme for its estimation based on
the Generalized Expectation-Maximization (GEM) algorithm. When applying
the GEM algorithm to the current “partial” MAP estimation problem (that is,(S1)), the ﬁrst step is to deﬁne the “complete data”. As the “observed” sourcecomponent S̄k,n is assumed to contain I × J concurrent voice components, anatural choice for the complete data is the corresponding hidden components,
1,1
I,J T
i,j
, · · · , Xk,n) with Xk,n∼ NC (0, λi,jthat is, X k,n = (Xk,nk,n ). From Eq. (8),the many-to-one relationship between X k,n and S̄k,n is described as S̄k,n =i,j
1T X k,n with 1 = (1, · · · , 1)T . In Sec. 2, we have already assumed that Xk,n



is independent of all other Xki ,j,n , so the log-likelihood of the complete data
X := {X k,n }k,n is
log fX|θ (X|θ) = −





H
log det πΛk,n + tr Λ−1
k,n X k,n X k,n ,

(14)

k,n


I,J 
where Λk,n = diag λ1,1k,n , · · · , λk,n . Taking the conditional expectation of Eq. (14)given S̄k,n and θ = θ and then adding log fθ (θ) to both sides, we obtain
Q(θ, θ ) = log fθ (θ) −k,nlog det πΛk,nH, (15)+ tr Λ−1
k,n E X k,n X k,n |Sk,n = S̄k,n , θ = θ T −1 T 
where E X k,n X H
1 Λk,n +k,n |Sk,n = S̄k,n , θ = θ = Λk,n − Λk,n 1 1 Λk,n 1
 T 
−1  T 
−1 T 2 |S̄k,n | Λk,n 1 1 Λk,n 1
1 Λk,n 1
1 Λk,n . Writing it in an element-wise expression, we obtain


i,j


|Aj (ej2πk/K )|2
Ψk,n
 c
i i,j
Q(θ, θ ) = −Uni,j , (16)log Hk Un +−αi U i,jHnkn i,jk,n i,ji,jrepresents the PSD estimate of the {i, j}th voice components,
where Ψk,n
i,j
=
Ψk,n
c

i,j
i,j

λ k,n  
λ k,n i,j−λ+|S̄k,n |2 .
φ
k,n
k,n


φ k,n
φ k,n

(17)

The notation = denotes equality up to constant terms. By setting the partial
derivatives of Q(θ, θ ) with respect to Hki and Uni,j at zero, we obtain the following
update formulae for Hki and Uni,jStatistical Model of Speech Signals

1   i,j j j2πk/K 2 i,jΨ |A (e)| /Un ,N J n j k,n i,jαUni,j 2 + KUni,j −Ψk,n |Aj (ej2πk/K )|2 Hki = 0, Uni,j ≥ 0

Hki =

251

(18)
(19)

k

By setting the partial derivatives of Q(θ, θ ) with respect to aj1 , · · · , ajP at zero,
we obtain the Yule-Walker equations
rpj =

P


j
ajq rp−q(p = 1, · · · , P ),

(20)

q=1

where rpj is deﬁned by the inverse DFT of the average spectral envelope over all
the voice components with index j such that
rpj =


k

i,j
Ψk,n

H  ik U  i,jnn,iepj2πk/K .

(21)

The update formula for the autoregressive parameters of the jth all-pole ﬁlter
can therefore be calculated using the well-known Levinson-Durbin algorithm.
Here it is important to note that when a sparse constraint comes into play,there is a need for some constraint on the scales of the factorized elements inorder to avoid an indeterminacy in the scaling. We thus adopt a simple procedurethat consists of calculating
Eq. (18) and then projecting it onto the unit normspace: Hki ← Hki / k Hki  .

4Experimental Results

We present here the separation results we obtained with the present BSS algorithm. All of the examples use the same two-input four-output impulse response,
which was measured in a varechoic chamber where the reverberation time was
0.5 sec. With this impulse response, we mixed two speech signals into four mixtures. The two speech signals of female speakers, taken from the ATR speech
database, were sampled at 16 kHz and band limited to the 50 Hz to 7 kHz
frequency range. The input Signal-to-Interference ratios (SIRs) are shown in
Tab. 1. Time-frequency representations were obtained using the polyphase ﬁlterbank analysis with a frame length of 32 ms and a hop size of 8 ms. The ﬁlter
length nk was set as follows: nk = 25 for Fk < 0.8; nk = 20 for 0.8 ≤ Fk < 1.5;
nk = 15 for 1.5 ≤ Fk < 3; nk = 10 for Fk ≥ 3, where Fk is the frequency in kHzof the k-th frequency bin. The AR order P was set at 12, and α at 10K. The
iterative algorithm comprising (S1)–(S3) was run for 3 iterations. For each stepof (S1), the GEM algorithm was run for 100 iterations.
We tested the performance of the present method with diﬀerent I and J
settings. Tab. 2 lists the SIRs obtained with the proposed method for various
settings of I and J, along with those obtained with the baseline methods, where

252

H. Kameoka et al.

Table 1. Input SIRs (dB)

Table 2. Output SIRs (dB) for various settings

ChannelSource #1 #2 #3 #4#1 -0.6 -0.3 -0.1 0.6#2
0.6 0.3 0.1 -0.6

I=2
I=5
I = 10 Base- BaseSource J = 12 J = 15 J = 8 J = 10 J = 8 J = 12 line1 line2#1
18.0 16.8 19.9 18.5 17.7 17.6 11.6 17.2
#2
11.7 10.9 14.1 13.9 13.6 14.5 11.0 13.9

Baseline1 and Baseline2 refer to Sawada’s method [11] and Yoshioka’s method[4], respectively. For Baseline1, we performed time-frequency analysis with aframe length of 256 ms and a hop size of 64 ms. For Baseline2, the frame lengthand the hop size were set at 16 ms and 8 ms. The best SIR result by the present
method was obtained when I and J were set at 5 and 8, which signiﬁcantly
outperforms both of the baseline methods. The results are very preliminary andthey need to be conﬁrmed by a more thorough analysis in the future.

5

Conclusion

This paper described a statistical model called the “composite autoregressive
system”, which consists of a time-invariant dictionary incorporating a set of
PSDs of excitation signals and a set of all-pole ﬁlters. Under this model, speech
signals are assumed to be characterized by the volume levels of the excitationand ﬁlter element pairs, that vary over time. We proposed to use this model
to develop a combined blind separation and dereverberation method for speech
and reasonably good separations were obtained for four mixtures of two speech
signals under a reverberant condition.

References
1. Douglas, S., Sawada, H., Makino, S.: Natural gradient multichannel blind deconvolution and speech separation using causal FIR filters. IEEE Trans. Speech, Audio
Process. 13(1), 92–104 (2005)
2. Smaragdis, P.: Blind separation of convolved mixtures in the frequency domain.
Neur. Comp. 22, 21–34 (1998)
3. Nakatani, T., Yoshioka, T., Kinoshita, K., Miyoshi, M., Juang, B.-H.: Blind speech
dereverberation with multi-channel linear prediction based on short time Fouriertransform representation. In: Proc. Int’l. Conf. Acoust., Speech, Signal Process.,
pp. 85–88 (2008)
4. Yoshioka, T., Nakatani, T., Miyoshi, M., Okuno, H.G.: Blind separation and dereverberation of speech mixtures by joint optimization. IEEE Trans. Audio, Speech,
Language Process (2010) (accepted for publication)
5. Dégerine, S., Zaı̈di, A.: Separation of an instantaneous mixture of Gaussian autoregressive sources by the exact maximum likelihood approach. IEEE Trans. Signal
Processing 52(6), 1499–1512 (2004)
6. Kameoka, H., Kashino, K.: Composite Autoregressive System for Sparse SourceFilter Representation of Speech. In: Proc. 2009 IEEE International Symposium on
Circuits and Systems (ISCAS 2009), pp. 2477–2480 (2009)

Statistical Model of Speech Signals

253

7. Benaroya, L., Bimbot, F., Gribonval, R.: Audio source separation with a single
sensor. IEEE Trans. Audio Speech Language Processing 14(1), 191–199 (2006)
8. Févotte, C., Bertin, N., Durrieu, J.-L.: Nonnegative matrix factorization,with
the Itakura-Saito divergence. With application to music analysis. Neural Comput. 21(3), 793–830 (2009)
9. Ozerov, A., Févotte, C.: Multichannel nonnegative matrix factorization in convolutive mixtures for audio source separation. IEEE Trans. Audio, Speech, Language
Process. 18(3), 550–563 (2010)
10. Olshausen, B.A., Field, D.J.: Emergence of simple-cell receptive field properties bylearning a sparse code for natural images. Nature 381, 607–609 (1996)
11. Sawada, H., Araki, S., Makino, S.: Measuring dependence of binwise separated
signals for permutation alignment in frequency-domain BSS. In: Proc. Int’l. Symp.
Circ., Syst., pp. 3247–3250 (2007)