Composite Autoregressive System for Sparse
Source-Filter Representation of Speech
Hirokazu Kameoka, and Kunio Kashino
NTT Communication Science Laboratories, NTT Corporation
3-1 Morinosato Wakamiya, Atsugi, Kanagawa 243-0198, Japan
Abstract— This paper presents a new generative model for
speech signals called a “composite autoregressive system”. Thismodel consists of a composite dictionary incorporating a setof the power spectral densities (PSDs) of excitation sources
and a set of all-pole filters where the gain of each pair of
excitation and filter elements is allowed to vary over time. We
use this model to develop a computationally efficient scheme forgenerating a sparse mixture representation of speech based on the
Expectation-Maximization algorithm. The algorithm iteratively
updates the excitation PSDs and the gains through the update
formulae, which reduce under a particular condition to the multiplicative update rule for non-negative matrix factorization with
the Itakura-Saito distance criterion, and the all-pole parameters
using the Levinson-Durbin algorithm.

I. I NTRODUCTION
We live in a world containing a huge variety of sounds
emanating from various combinations of sound sources. In
recent years there has been a growing interest in the searchfor sparse representations of acoustic signals that allow us
to analyze complex acoustic scenes in a data-driven manner.
The assumption behind this approach is that the underlying
model generating a time series comprising a wide varietyof observations is simply a linear combination of a compact
set of static dictionary vectors each of which is associatedwith a time-varying coefficient. Given a set of observed
vectors (which can be either signals or spectra), the goal
is to find a set of ‘atoms’, each defined by a dictionary
vector, such that any observed vector can be parsimoniously
represented using only a small subset of active atoms. As aconsequence of this sparsity requirement, each of the atoms iscapable of representing an independently recurrent pattern that
underlies the observations. Sparse coding [1], non-negative
matrix factorization (NMF) [2], shift-invariant sparse coding[3], shifted NMF [4], non-negative matrix factor deconvolution(NMFD) [5, 6], NMF2D [7], non-negative tensor factorization(NTF) [8], semi-NMF [9], shift-invariant semi-NMF [10], andcomplex NMF [11] can be designed to produce various typesof sparse representations of acoustic signals, many of whichhave been applied extensively with notable success. The choiceof which method to invoke depends on what kind of structural
regularity or property we would like to exploit to obtain thedesired form of signal decomposition.
This paper describes a signal model that allows us to obtain
yet another type of sparse representation suited to analyzing
speech. For this purpose we first consider factorizing speech
into the product of prosodic and phonemic factors on thebasis of the speech production mechanism. It has been shownthat the speech production system can be well modeled bya linear system comprising a glottal excitation input and a

978-1-4244-3828-0/09/$25.00 ©2009 IEEE

vocal-tract resonance filter that respectively determine thedegree of periodicity (pitch) and the phoneme (timbre) of
the voice. Although the characteristic of a vocal-tract filter
and the periodicity of a glottal excitation vary over time
during an entire speech utterance, the phoneme number andperiodicity range are both usually limited. It is therefore
reasonable to assume that observed speech signals (of one
or more concurrent voices) have been generated from a composite linear system that is composed of the direct productof a limited set of excitation sources and a limited set of
vocal-tract filters where each excitation and filter element
pair is associated with a time-varying activation coefficientthat determines the instantaneous loudness of each of the
concurrent voice components. We then assume that each vocaltract filter is restricted to an all-pole filter. Since the output ofan all-pole filter with a Gaussian white noise input is usuallyreferred to as an autoregressive (AR) model, which is usedin diverse applications of speech analysis, we will term thepresent generative model a “composite autoregressive system”.
In the long term, we are planning to apply the present model
to various blind signal processing problems including speechdenoising, speech dereverberation and speech separation forwhich a time-varying all-pole filter model has proved to bevery useful [14]. As a preliminary to these applications, thispaper presents the basic formulation of the model, an efficientalgorithm for estimating the model parameters, and someexperimental results to show how well the present model
represents speech signals.
II. G ENERATIVE M ODEL
A. Autoregressive system with colored noise input
We assume that a speech signal in a short-term segment,{xt }K
t=1 , is a sampled sequence drawn from the P -order
autoregressive process such that
xt =

P


ap xt−p + t ,

(1)

p=1

where t is an excitation source signal that is assumed to bea zero-mean stationary Gaussian noise and its autocorrelation
function is ft . It is important to note that t does not have to bea white noise. Let x = (x1 , · · · , xK )T ∈ RK and its discrete
Fourier transform (DFT) be X = Dx = (X1 , · · · , XK )T ∈CK where D is a DFT matrix. Then, according to Eq. (1),the linearity of DFT, and the stationarity and the Gaussianityof t , X follows a multivariate complex Gaussian distribution

2477

X ∼ NC (0, Λ),

(2)



with a diagonal covariance matrix Λ = diag λ1 , · · · , λK
whose diagonal elements are equal toFk,|A(ej2πk/K )|2
A(z) = 1 − a1 z −1 − · · · − aP z −P ,
λk =

(3)
(4)

where k is the frequency index. Fk is the DFT of ft andrepresents the power spectral density (PSD) of the excitation
source signal t (namely, the spectral fine structure), that canhave any shape and is not necessarily flat. 1/|A(ej2πk/K )|2 , on
the other hand, corresponds to the spectral envelope expressedas the PSD of the all-pole transfer function.
B. Composite autoregressive system
The composite autoregressive system is assumed to consistof a dictionary of I excitation source PSDs and a dictionaryof J all-pole filters. Subsequently, we use superscripts iand j to denote the indices of the excitation PSDs and the
all-pole filters, respectively, and let the ith excitation PSD
and the jth all-pole transfer function be denoted by Fki
and 1/Aj (ej2πk/K ). The system is able to generate I × J
different voice components each of which is characterized by
combining elements drawn from the respective dictionaries.
The gain of each voice component, i.e. the voice activity, isallowed to vary individually over time. We thus use Uni,j todenote the voice activity of the {i, j}th voice component at
the nth frame (or segment). Following the discussion in II-A,
the DFT of the {i, j}th voice component at the nth frame,
i,j
i,j T
X i,j
n = (X1,n , · · · , XK,n ) , follows a multivariate complex
Gaussian distribution
NC (0, Λi,jn ) with a diagonal covariance
 i,j
i,j
matrix Λn = diag λ1,n , · · · , λi,jK,n whose elements are
Fki Uni,j,|Aj (ej2πk/K )|2
Aj (z) = 1 − aj1 z −1 − · · · − ajP z −P .
λi,j
k,n =

(5)
(6)

I,J
Now, if we assume that X 1,1
n , · · · , Xn

are mutually independent and i.i.d., it follows that




Yn=
X i,j
Λi,j
n ∼ NC 0, Φn , Φn =n , (7)ijijIII. R EGULARIZED I NVERSE P ROBLEM
A. Maximum A-Posteriori (MAP) estimation
Given a set of observed STFT components Y , we would
like to find the estimate of θ that maximizes the posterior
density P (θ|Y ) ∝ P (Y |θ)P (θ), or equivalently,
log P (Y |θ) + log P (θ).

As seen from an inspection of Eq. (8), log P (Y |θ) is equal up
to constant terms to the goodness of fit between |Yk,n |2 and
φk,n defined by the Itakura-Saito distance [12, 13]. We arethus led to obtain a PSD model with as small a reconstruction
error as possible. On the other hand, as with the sparse coding
concept [1, 15], we would like to keep the voice activities assparse as possible. For this purpose, the term log P (θ) canbe designed to have a sparsifying effect. In the subsequent
analysis, for convenience we will use an inverse Gamma priordefined over Uni,j > 0


β
1
(12)
exp − i,j ,P (θ) ∝
i,j α+1Un
i,j,n Un
which promotes sparsity when α is large. The reason for
its convenience will become apparent later. Maximizing Eq.
(11) therefore combines the goals of a small reconstruction
error and sparseness. As a consequence, the more frequently
a certain spectral fine/envelope structure emerges in the observed spectra, the more likely it is to be captured in the
excitation/filter dictionary.
B. Expectation-Maximization (EM) algorithm
Although it is difficult to obtain a closed-form solution for
the MAP estimate of θ, we can develop a computationallyefficient scheme for its estimation based on the ExpectationMaximization (EM) algorithm [16]. In this subsection, we
briefly review the general principle of the EM algorithm.
We use y to denote the observed “incomplete” data, that
follows P (y|θ) indexed by the parameter θ, and let ŷ denotethe “complete” data, related to y by a certain noninvertible
transform H : ŷ → y. The complete data pdf is written as
P (ŷ|θ) = P (ŷ|y, θ)P (y|θ),(13)

and, taking the logarithm on both sides, we obtain

where Y n ∈ CK denotes the DFT of the observed signal at
the nth frame. Assuming the observations to be generated independently from the generative model, the probability density
K×N
function (pdf) of Y = (Yk,n )K,N = (Y 
1, · · · , Y N ) ∈ C
can be written concisely as P (Y |θ) = k,n P (Yk,n |θ) with
1
|Yk,n |2,
(8)
P (Yk,n |θ) =
exp −
πφk,n
φk,n
where θ contains all the unknown parameters of the system:Fki , ajp , Uni,j .
(9)θ=log P (y|θ) = log P (ŷ|θ) − log P (ŷ|y, θ).

The diagonal element of Φn , i.e. φk,n , corresponds to the PSDof the output signal produced by the system such that  i,jλk,n .
(10)
φk,n =(14)

Now, taking the conditional expectation of both sides given yand the value θ = θ , it follows thatlog P (y|θ) = E log P (ŷ|θ)|y, θ = θ
− E log P (ŷ|y, θ)|y; θ = θ .

(15)

Adding log P (θ) to both sides and putting L(θ) ≡ log P (y|θ)+ log P (θ), Q(θ, θ ) ≡ E[log P (ŷ|θ)|y, θ = θ ] + log P (θ)
and V (θ, θ ) ≡ E[log P (ŷ|y, θ)|y, θ = θ ], it reads
L(θ) = Q(θ, θ ) − V (θ, θ ).

i,j,k,n

i

(11)

(16)

We notice here that L(θ) has the same form as Eq. (11).
By invoking the Jensen’s inequality, V (θ, θ ) ≤ V (θ , θ ).
Therefore, if we obtain θ such that Q(θ, θ ) > Q(θ , θ ) then

j

2478

L(θ) > L(θ  ).

(17)

Based on this result, we can construct an iterative algorithm
that ensures the monotonic increase of L(θ) as follows:
Set at 1 and θ() at a certain value.
Compute Q(θ, θ() ).
Update θ by θ(+1) ← argmax Q(θ, θ() ).

Initial step:E step:M step:θ

Set ← + 1 and return to E step.
C. Definition of incomplete data
When applying the EM algorithm to the current MAP
estimation problem, the first step is to define the “completedata”. As the observed STFT component Yk,n is assumed
to contain I × J concurrent voice components, a natural
choice for the complete data Ŷ k,n is the corresponding hidden
1,1
I,J T
, · · · , Xk,n) withcomponents, that is, Ŷ k,n = (Xk,ni,jXk,n∼ NC (0, λi,jk,n ).

D. M-step update formulae
We can derive closed-form expressions of the M-step update
formulae for all the model parameters thanks to the definitionsof the complete data and the prior. Differentiating Q(θ, θ )partially with respect to Fki and Uni,j and setting them at zero,
we obtain the following update formulaei,jj j2πk/K 2)|
1   Ψk,n |A (e,i,j
NJ n jUn
j j2πk/K 2
 Ψi,j)|1
k,n |A (e
i,j
β+
.
Un =
K +α+1Fki

Fki =Differentiating Q(θ, θ ) partially with respect to aj1 , · · · , ajP
and setting them at zero, we obtain the normal equations

(18)

(19)

with H = 1, · · · , 1 ∈ RIJ . In II-B we have already seen thati,jis mutually independent across i, j, k and n, so the logXk,nlikelihood of the complete data Ŷ = (Ŷ 1,1 T , · · · , Ŷ K,N T )Tis given by
log det πΛk,n + Ŷ k,n H Λ−1Ŷlog P (Ŷ |θ) = −
k,n
k,n
=−



k,n



H
, (20)
log det πΛk,n + tr Λ−1
Ŷ
Ŷ
k,n
k,n
k,n

k,n


I,J 
where Λk,n = diag λ1,1k,n , · · · , λk,n .
Taking the conditional expectation of Eq. (20) given Yk,nand θ = θ and then adding log P (θ) to both sides, we obtain
log det πΛk,nQ(θ, θ  ) = log P (θ) −


k,n


H
+ tr Λ−1
k,n E Ŷ k,n Ŷ k,n |Yk,n , θ = θ, (21)

where
E Ŷ k,n Ŷ k,n H |Yk,n , θ = θ =
(22)




T

T −1

Λk,n − Λk,n H HΛk,n H
HΛk,n +



−1
−1

T

T
HΛk,n H T|Yk,n |2 Λk,n H HΛk,n HHΛk,n .
Writing it in an elementwise expression and using c to denotethe terms that do not depend on θ, we obtain
j j2πk/K 2

)|
Ψi,j
k,n |A (ei i,jlog Fk Un +Q(θ, θ ) = −Fki Uni,jk,n i,jβ(α + 1) log Uni,j + i,j + c, (23)
−
Un
n i,j
where Ψi,jk,n represents the PSD estimate of the {i, j}th voicecomponents, which is given by


i,j 
λ k,n |Yk,n |2 − φk,n
i,j
 i,j
.
(24)
Ψk,n = λ k,n 1 +
φ2
k,n

(26)

k

rpj =

From Eq. (7), the many-to-one relationship between the complete data Ŷ k,n and the incomplete data Yk,n is described as
Yk,n = H Ŷ k,n ,

(25)

P


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

(27)

q=1

where rpj is defined by the inverse DFT of the average spectralenvelope over all the voice components with index j such that  Ψi,j
k,n
−1
j
.
(28)
rp = DFTi i,jni Fk U n
The update formula for the autoregressive parameters of thejth all-pole filter can therefore be calculated via⎡ j⎤ ⎡ j⎤−1 ⎡ j ⎤j
· · · r1−P
r0
a1
r1
⎢ .. ⎥ ⎢ ..
⎥
⎢
.
.
..
.. ⎦ ⎣ ... ⎥(29)⎣ . ⎦=⎣ .
⎦,ajP

rPj −1

···

r0j

rPj

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 in order to avoid anindeterminacy in the scaling. We thus adopt a simple procedurethat consists of calculating Eq. (25)and then projecting it ontoi2
the unit norm space: Fki ← Fki /k Fk as with [15].
E. Reduction to NMF with Itakura-Saito distance
When P = 0, J = 1, and P (θ) = const., the present
algorithm reduces to the NMF algorithm with the ItakuraSaito distance criterion discussed thoroughly by Févotte in[13]. Since P = 0 ⇒|Aj (ej2πk/K )|2 = 1, we verify that,
under this condition, Eqs. (25), (26) become1  Ψk,n
1  Ψk,n, Uni,1 =,i,1N n Un
K
Fki
i,1

Fki =

i,1

(30)

k

which amount to the update formulae developed in [13].
IV. E XPERIMENTIn this section we report results for speech data excerpted
from the ATR B-set speech database. Speech data was monaural and sampled at 16kHz. STFT was computed using a
Hanning window that was 32ms long with a 16ms overlap.
P was set at 12. The algorithm was run for 300 iterations.
Fki , ajp and Uni,j were initialized with random values.

2479

1
Sparseness score

Frequency bin #

250
200
150
100
50

0.6
0.4
0.2
0

50

100

150

200

250

300

350

400

450

50

100

150

200

Sparseness score

Observed spectrogram |Yk,n |2 .

250
200
150

350

400

450

300

350

400

450

0.6
0.4
0.2
0

50

100

150

200

250
Frame #

Fig. 4. Sparseness scores obtained with the present methods using an inverse
Gamma prior (top) and using a uniform prior (bottom)

50
50

100

150

200

250

300

350

400

450

TABLE IS IGNAL - TO - NOISE RATIO ( D B) UNDER EACH CONDITION FOR I AND J .

Frame #

Fig. 2.







Reconstructed spectrogram φk,n .



 

(TGSWGPE[KPFGZ

Fig. 3.

300

0.8

100



250
Frame #

1

Frame #

Fig. 1.

Frequency bin #

0.8









I \J

1

3

5

5
10
15

3.07
3.78
5.18

3.54
4.63
6.32

4.91
5.81
6.65

V. C ONCLUDING R EMARKS
We presented a new generative model for speech signals
called a “composite autoregressive system”. It consists of acomposite dictionary of excitation source PSDs and all-polefilters where the gain of each excitation and filter element pairis allowed to vary over time. In addition, we used this model
to develop a computationally efficient scheme for generatinga sparse mixture representation based on the EM algorithm.

(TGSWGPE[KPFGZ

Estimates of the excitation PSDs (left) and all-pole filters (right).

For the first experiment we tested the present method undervarious conditions for I and J on a single voice signal froma female speaker, whose spectrogram is shown in Fig. 1. For
each condition we calculated the signal-to-noise ratio,
2k,n |Yk,n |(31)
SNR(dB) = 10 log10   ,
|Yk,n | − φk,n 2k,nto see how well the present model represented speech. The
results can be seen in Tab. I. As examples, Fig. 2 and Fig.
3, respectively, show the reconstructed spectrogram and thedictionaries of the excitation source PSDs and all-pole filters
obtained with the present method with I = 15 and J = 5.
For the second experiment we tested the present algorithm
to ascertain whether the inverse Gamma prior had a sparsifying
effect on the atomic activation. We used

( i,j |Uni,j |d )/( i,j Uni,j 2 )d/2 − (IJ)1−d/2
sparseness[n] =1 − (IJ)1−d/2
to measure the framewise sparseness of the voice activities
Uni,j across i and j. This function evaluates to unity if andonly if {Uni,j }i,j contains only a single non-zero component,
and takes a value of zero if and only if all components areequal. Fig. 4 shows the sparseness scores with d = 1 obtainedwith the present method and with the present method run
without including the prior. The result confirmed that theinverse Gamma prior had a certain effect on the present model
as regards obtaining a parsimonious mixture representation.

R EFERENCES[1] B.A. Olshausen, and D.J. Field, “Emergence of simple-cell receptive field
properties by learning a sparse code for natural images,” Nature, 381,
pp. 607–609, 1996.
[2] D.D. Lee, and H.S. Seung, “Learning the parts of objects by non-negative
matrix factorization,” Nature, 401, pp. 788–791, 1999.
[3] M.S. Lewicki, and T.J. Sejnowski, “Coding time-varying signals using
sparse, shift-invariant representations,” Proc. NIPS’99, 11, 1999.
[4] M. Mørup, K.H. Madsen, L.K. Hansen, “Shifted non-negative matrix
factorization,” Proc. MLSP’07, pp.139–144, 2007.
[5] P. Smaragdis, “Discovering auditory objects through non-negativity constraints,” Proc. SAPA’04, in CD-ROM, 2004.
[6] T. Virtanen, “Separation of sound sources by convolutive sparse coding,”
Proc. SAPA’04, in CD-ROM, 2004.
[7] M.N. Schmidt and M. Mørup, “Nonnegative matrix factor 2-D deconvolution for blind single channel source separation,” Proc. ICA’06, 2006.
[8] A. Shashua, T. Hazan, “Non-negative tensor factorization with applications to statistics and computer vision,” Proc. ICML’05, 792–799, 2005.
[9] C. Ding, T. Li, and M.I. Jordan, “Convex and semi-nonnegative matrix
factorization for clustering and low-dimension representation,” Lawrence
Berkeley National Laboratory, Tech. Rep. LBNL-60428, 2006.
[10] J. Le Roux, A. de Cheveigné, and L.C. Parra, “Adaptive templatematching with shift-invariant semi-NMF,” Proc. NIPS’08, 2008.
[11] H. Kameoka, N. Ono, K. Kashino, S. Sagayama, “Complex NMF: A newsparse representation for acoustic signals,” Proc. ICASSP’09, to appear.
[12] F. Itakura and S. Saito, “Analysis synthesis telephony based on the
maximum likelihood method,” Proc. ICA’68, pp. C-17–C-20, 1968.
[13] C. Févotte, N. Bertin, and J.-L. Durrieu, “Nonnegative matrix factorization with the Itakura-Saito divergence with application to music analysis,”
Tech. Rep. TELECOM ParisTech 2008D006, 2008.
[14] T. Yoshioka, T. Nakatani, and M. Miyoshi, “An integrated method forblind separation and dereverberation of convolutive audio mixtures,” Proc.
EUSIPCO’08, in CD-ROM, 2008.
[15] P.O. Hoyer, “Non-negative sparse coding,” Neural Networks for Signal
Processing XII, pp. 557–565, 2002.
[16] M. Feder, E. Weinstein, “Parameter estimation of superimposed signals
using the EM algorithm,” IEEE Trans. ASSP, 36(4), pp. 477–489, 1988.

2480