INFINITE KERNEL LINEAR PREDICTION FOR JOINT ESTIMATIONOF SPECTRAL ENVELOPE AND FUNDAMENTAL FREQUENCY
Kazuyoshi Yoshii

Masataka Goto

National Institute of Advanced Industrial Science and Technology (AIST){k.yoshii, m.goto}@aist.go.jp
Periodic signal

ABSTRACTThis paper presents a new probabilistic formulation of linear prediction (LP) for jointly estimating the spectral envelope and fundamental frequency (F0) of a speech signal. A main problem of classical LP
is that the peaks of the estimated envelope are highly biased toward
the harmonic partials of a speech spectrum. To solve this problem,
we propose a nonparametric Bayesian model called infinite kernel
linear prediction (IKLP) based on a Gaussian process with multiple
kernel learning. Our model can represent the periodicity of a speech
signal by using a weighted sum of infinitely many periodic kernels
that correspond to different F0s. We put a gamma process prior on
the positive weights of those kernels and perform sparse learning todetermine a predominant kernel indicating the F0 at the same time ofspectral envelope estimation. The experimental results showed that
our model can estimate spectral envelopes and F0s of speech andsinging signals while identifying pitched segments.

Time

Frequency

SourceConvoluted by
impulse response

Time

FilterFrequency

Observed signal

Time

Multiplied byfrequency response

Time

Frequency
Frequency

Fig. 1. A source-filter model of the speech production mechanism.
that imposes smoothness constraints on all-pole filters. Villavicencio et al. [5] applied iterative cepstral smoothing to all-pole filters.
In a wider context of spectral envelope estimation not limited to LP,
Kawahara et al. [6] proposed a speech analysis-and-synthesis system called STRAIGHT that can decompose a speech spectrum into
the harmonic and non-harmonic spectra and the spectral envelope.
Nakano and Goto [7] aimed to avoid the bias of fundamental frequencies (F0s) by averaging the spectral envelopes of neighboring
frames. Although these methods are considered to yield good estimates, they all require that the F0s be known in advance.
A probabilistic approach provides a solid mathematical foundation for modeling F0s (sources) and envelopes (filters). For example,
Sasou and Tanaka [8] proposed the periodicity of the source signal
could be represented by using an AR hidden Markov model (HMM)with circularly-connected states. Toda and Tokuda [9] used a trajectory HMM for capturing the temporal dynamics of harmonic components derived from the source signal, although F0 information isrequired. Kameoka et al. [10, 11] pioneered probabilistic models forjoint F0 and spectral-envelope estimation. In [10] they formulated
multiple kernel LP (MKLP) based on a Gaussian process (GP) thatis specified by a fixed number of periodic kernels corresponding todifferent F0s. To determine a predominant kernel, the maximum-aposteriori (MAP) estimate of kernel weights is computed (F0 estimation) at the same time the all-pole filter is estimated. Our study issituated as the state of the art in this thread of research.
In this paper we propose a more sophisticated model called infinite kernel LP (IKLP) that combines the strengths of Bayesian nonparametrics and kernel methods in a principled manner. We take thelimit of the MKLP model [10] as the number of kernels diverges
to infinity, and put a gamma process (GaP) prior on the weights of
those kernels. As a result of variational Bayesian (VB) inference, we
can obtain a sparse estimate of infinitely many weights. We believethat our study could contribute not only to the field of signal processing by revealing the underlying probabilistic assumptions of LP,
but also to the field of machine learning by providing a new efficientand convergence-guaranteed algorithm as a general solution to the
problem of multiple kernel learning (MKL) [12].

Index Terms— Linear prediction, source-filter model, Bayesiannonparametrics, kernel methods, Gaussian and gamma processes.
1. INTRODUCTION
Spectral envelope estimation forms the basis for speech and singing
analysis. The speech production mechanism is considered to be well
explained by the source-filter theory, which assumes that an excitation signal generated by the vocal cords is modified acoustically by
the vocal tract, i.e., the source signal is convoluted by an impulse
response of the filter. In the frequency domain, the fine structure andspectral envelope of a speech spectrum reflect the frequency characteristics of the source and filter, respectively (Fig. 1).
Linear prediction (LP) is a popular parametric approach to spectral envelope estimation [1]. In general, we assume speech signals to
be autoregressive (AR), i.e., the filter response can be represented byan all-pole transfer function. This assumption is widely accepted asreasonable because most phonemes have no anti-resonance and thehuman auditory system is sensitive to spectral peaks (formants) corresponding to poles. If the source signal is a white Gaussian noise,
we can correctly estimate the filter coefficients by maximum likelihood estimation in a probabilistic framework [1]. When we analyze
a pitched signal, however (i.e., when the source signal is periodic),the estimated filter (spectral envelope) unnecessarily has very sharp
peaks at the harmonic partials of the spectrum.
A lot of effort has been devoted to solving this problem. For example, El-Jaroudi and Makhoul [2] proposed a basic approach thatfits the frequency response of an all-pole filter to a discrete set of harmonic partials. This discrete all-pole (DAP) modeling was extended
by Badeau and David [3] for AR and moving-average (ARMA) modeling. Alternatively, Oudot et al. [4] took a regularization approachThis study was supported in part by the JSPS KAKENHI Grant Number
23700184 and the JST OngaCREST project.

978-1-4799-0356-6/13/$31.00 ©2013 IEEE

White noise

Autoregressive system

463

ICASSP 2013

we instead discuss a linear regression problem that aims to modelan excitation signal (continuous function) (t) over time t. Let t =M{tm }M
m=1 be a set of times at which {m }m=1 are sampled from
(t). Our goal is to approximate (t) by the weighted sum of J basis
functions {φj (t)}Jj=1 as follows:(6)
(t) = Jj=1 wj φj (t) + η(t) = φ(t)T w + η(t),2. LINEAR PREDICTION
This section introduces relevant probabilistic models for estimating
spectral envelopes (coefficients of AR filters) of speech and singing
signals. First we revisit a classical AR model based on a strong
assumption that excitation signals are white Gaussian noise. We thenexplain kernelized models that can perform robust estimation evenwhen excitation signals have a periodic nature.

where η(t) is an error function and φ(t) = (φ1 (t), · · · , φJ (t))T .
Let η = (η1 , · · · , ηM )T be a set of values sampled from η(t) at t
and let Φ = (φ(t1 ), · · · , φ(tM ))T be an M -by-J design matrix.
Then, we can write the “marginal” of Eq. (6) on t as follows:

2.1. Probabilistic Formulation
We review a standard formulation that assumes a target audio signal
to locally follow an AR process. Let x = (x1 , x2 , · · · , xM )T be M
consecutive samples contained in a short segment (shifting window).
If the local signal x follows a P -order AR process, we can write(1)
xm = P
p=1 ap xm−p + m ,

 = Φw + η.
w ∼ N (0, νw I), η ∼ N (0, νe I),

where a = (a1 , · · · , aP ) is a set of P coefficients of the AR filter
called predictor coefficients, and  = (1 , · · · , M )T is a noise term.
In the source-filter modeling of speech signals,  represents an excitation signal generated by the vocal cords and a represents the resonance characteristics of the vocal tract. This AR model can be regarded as a linear system that takes m as input and then outputs xm
according to an all-pole transfer function given by A(z) = 1/(1 −a1 z −1 · · · − aP z −P ), i.e., X(z) = E(z)A(z), where X(z) andE(z) are z transforms of x and . The frequency response of the allpole filter (spectral envelope) is therefore given by |A(e2πim/M )|2 ,where m here is the index of a frequency bin (Fig. 1).
Given the observed audio signal x, our objective is to estimate
the coefficients a in a probabilistic framework. To do this, we needto specify the characteristics of the excitation signal . A standard
assumption is that  is a white Gaussian noise given by ∼ N (0, νw ΦΦT + νe I).

Km,m = φ(tm )T φ(tm ).

x ∼ N (0, Ψ−1 (νw K + νe I)Ψ−T ).

2.3. Multiple Kernel Learning
We discuss how to design a kernel matrix K that reflects the characteristics of the excitation signal . If the observation x is a pitchedsignal having the F0, a natural choice is using a periodic kernel. Forexample, k(t, t ) = exp(−2 sin2 (π t−t)/l2 ) is a well-known kerTnel having the period T . This means that all basis functions are
implicitly assumed to have the period T (the F0 is given by 1/T ).
Alternatively, Kameoka et al. [10] designed a basis function having
H harmonic components with equal power as follows:
t−cj(cj is a phase)(12)
φj (t) = H
h=1 sin 2πh T.
.
.
.
xM −1 · · · xM −P

We can compactly rewrite Eq. (1) as  = Ψx, which means(4)

Using Eqs. (2) and (4) gives a likelihood function of x as follows:x ∼ N (0, νΨ−1 Ψ−T ).

(11)This is a GP regression model [10], and it includes Eq. (5) as a special case. In fact, when ΦΦT = I, i.e., when J basis functions are
independent of each other like a series of Dirac delta functions, we
can recover Eq. (5) by assuming ν = νw + νe .

(2)

x = Ψ−1 .

(10)Instead, any positive semidefinite matrix can be used as a kernel matrix. This enables the kernel trick; we can directly calculate K without explicitly specifying basis functions as Km,m = k(tm , tm ),where k(tm , tm ) is a kernel function. From this kernelization andEq. (4), we derive a likelihood function of x as follows:0−aP · · · −a1 1(9)This is indeed a GP [13] because any marginal  sampled from (t) is
Gaussian distributed with a positive semidefinite covariance matrix
specified by a term K = ΦΦT called a kernel matrix. Each elementof K is given by the inner product of basis functions as follows:

where ν is a noise variance and I an identity matrix. This means that
M elements of  are independent and identically distributed according to a Gaussian N (0, ν). We let Ψ be an M -by-M approximate
circulant matrix and X be an M -by-P matrix as follows:⎡ 0 ··· 0 ⎤
⎤
⎡ 1
.
..
.
⎢ x
⎥
⎥
⎢ −a . . .
.
.
⎢ 1
⎥
⎥
⎢ 1
⎢
⎥
.
⎥
⎢ . .
⎢ .. . . . 0 ⎥⎥⎢ . .. ...
⎥. (3)
Ψ=⎢ .
⎥, X = ⎢
⎢ .
⎥
⎥
⎢ −aP
⎢ ..
⎥
⎥
⎢
x1 ⎥
⎢
..
.. ..
⎦
⎣
⎣
⎦
.
.
.
.
.

0

(8)

where νw and νe are scales of Gaussian variance. Then, we getT ∼ N (0, νI),

(7)

We assume that both ω and η are Gaussian distributed as follows:

and used a kernel matrix K calculated using Eq. (10). A problem isthat the true period T of the excitation signal  is unknown. Given
the observation x, one must therefore estimate K itself.
Multiple kernel learning (MKL) [12] is a powerful solution tothis problem. More specifically, the kernel matrix K is defined as
the weighted sum of I kernel matrices as follows:(13)
K = Ii=1 θi K i ,(5)This is a standard probabilistic model of LP. In maximum-likelihood(ML) estimation, the optimal values of a are obtained by solving anormal equation given by X T Xa = X T x.
This model works well when the excitation signal  is well approximated by an isotropic Gaussian defined as Eq. (2). Speech andsinging signals, however, exhibit a clear periodicity that is derived
from the periodic vibration of the vocal cords (source).

where θ = {θ1 , · · · , θI } is a set of kernel weights to be estimatedfrom the observed signal x and K i is a periodic kernel having theperiod Ti . The value of θi indicates a degree of the predominanceof period Ti in x (the F0 is given by 1/Ti ). Plugging Eq. (13) into
Eq. (11) gives a likelihood function of x as follows:(14)x ∼ N 0, Ψ−1 νw Ii=1 θi K i + νe I Ψ−T .

2.2. Kernelization based on a Gaussian Process
Kameoka et al. [10] generalized the probabilistic model specified as
Eq. (5) in terms of Gaussian process (GP) regression. In Section 2.1,
we assumed that  is a white Gaussian noise given by Eq. (2). Here

464

The parameters a and θ can be estimated by using the expectationmaximization (EM) algorithm. In [10], several hundreds of kernels
having different periods T = {T1 , · · · , TI } are prepared, and MAP
estimation is then performed by putting a generalized Gaussian distribution on each θi as a prior distribution. However, the values of θ
do not become truly sparse in MAP estimation.

3.2.1. Deriving Matrix-variate Inequalities
To derive the tractable lower bound L , we need to use some inequalities. We start with the following definitions:Definition 1 (Positive semidefinite matrix). We say a symmetric matrix A is positive semidefinite iff z T Az ≥ 0 for any vector z or iff
A = Z T Z for some real matrix Z.
Definition 2 (Convex/concave matrix functions). We say a matrixvariate function f (·) is convex iff λf (A)+(1−λ)f (B) ≥ f (λA+(1 − λ)B) for any number 0 ≤ λ ≤ 1. Conversely, we say f (·) is
concave iff λf (A) + (1 − λ)f (B) ≤ f (λA + (1 − λ)B).
Let V be a positive semidefinite (PSD) matrix and z be any vector.
Using these two definitions, we can derive the following lemmas (weomit the proofs because of space limitations):
Lemma 1. A matrix-variate function f (V ) = log |V | is concave.

3.1. Nonparametric Bayesian Formulation(15)

We here put a gamma process (GaP) prior on the infinite-dimensional
vector θ. More specifically, we introduce independent gamma priors
on elements of θ as follows:(16)

Lemma 2. A matrix-variate function g(V ) = z T V −1 z is convex.
These lemmas lead to useful inequalities of matrix-variate functions.
First, we can apply first-order Taylor expansion to f (V ) around an
arbitrary PSD matrix Ω as follows:

As the truncation level I diverges to infinity, the vector θ approximates an infinite sequence drawn from a GaP with shape parameter
α. It is proven that the effective number of elements, I + , such that
θi >  for some number  > 0 is almost surely finite. If I is sufficiently larger than α, we can expect that only a few of the I elementsof θ will be substantially greater than zero. This property gives atheoretical basis to sparse learning in an infinite space.
To complete the Bayesian formulation, we put gamma priors on
the positive weights νw and νe as follows:νw ∼ Gamma(aw , bw ), νe ∼ Gamma(ae , be ),log |V | ≤ log |Ω| + tr(Ω−1 V ) − M,(21)

where M is the size of V . Second, we can apply a matrix inequality
proposed by Sawada et al. [14] to g(V ) as follows:
−1

I
zT
z ≤ Ii=1 z T ΥTi V −1
(22)i Υi z,i=1 V i(17)where {V i }Ii=1 is a set of PSD matrices and {Υi }Ii=1 is a set of
auxiliary matrices that sum to the identity matrix.

where a∗ and b∗ are hyperparameters indicating the shape and rate
parameters of the gamma distribution. We also put a Gaussian prior
on a thusly: a ∼ N (0, λI), where λ is a hyperparameter.
The F0 is given by identifying a predominant kernel K i with thehighest value of E[θi ] (the F0 is 1/Ti ). The spectral envelope is specified by the filter coefficients E[a]. In addition, our model can distinguish whether the target signal x is pitched or unpitched according
to a degree of periodicity E[νw ]/E[νw + νe ] as in MKLP [10].

3.2.2. Deriving Evidence Lower Bound and Updating Formulas
We derive updating formulasby computing the tractable ELBO L .
Let K be K = νw i θi K i + νe I. Using Eqs. (21) and (22)enables us to compute the lower bound of E[log p(x|θ, a, νw , νe )](the first term of L given by Eq. (19)) as follows:11
M
log(2π) − E[log |K|] − E[xT ΨT K −1 Ψx]
2
2
2
1
1
1
E[νw θi ] tr(Ω−1 K i ) − E[νe ] tr(Ω−1 ) + const.
≥ − log |Ω| −
2
2 i
2
1  
1   1  T T T −1
x Ψ Υi K i Υi Ψx − E ν1 xT ΨT ΥTE−
0 Υ0 Ψx,e
2 i νw θ i2(23)

E[log p(x|·)] = −3.2. Variational Bayesian Inference
Given an observed signal x, our goal is to compute a posterior distribution over random variables p(θ, a, νw , νe |x) by using the Bayes
rule p(θ, a, νw , νe |x) = p(x, θ, a, νw , νe )/p(x). Since analytical
calculation of p(x) is infeasible, we use a variational Bayesian (VB)
method for approximating p(θ, a, νw , νe |x) by a variational distribution that can be factorized into four posteriors as follows:q(θ, a, νw , νe ) = q(a)q(νw )q(νe )(20)

where q(x|θ, a, νw , νe ) is a lower bound of p(x|θ, a, νw , νe ) (seeEq. (23)). For tractability we assume q(a) = δa∗ (a), where δa∗ is
a Dirac delta function taking infinity at a MAP point estimate a∗ .

This section proposes a nonparametric Bayesian model, called infinite kernel linear prediction (IKLP), for jointly estimating the spectral envelope and F0 of a speech signal in a principled manner. First
we take the limit of Eq. (14) as the number of kernels goes to infinity,i.e., I → ∞, in the framework of Bayesian nonparametrics. Thenwe put appropriate priors on unknown variables for a full Bayesiantreatment. To approximate a posterior distribution of those variables,
we derive an efficient and convergence-guaranteed algorithm.

θi ∼ Gamma(α/I, α).

q(θ) ∝ p(θ) exp(Eq(a,νw ,νe ) [log q(x|θ, a, νw , νe )]),q(νw ) ∝ p(νw ) exp(Eq(θ,a,νe ) [log q(x|θ, a, νw , νe )]),q(νe ) ∝ p(νe ) exp(Eq(θ,a,νw ) [log q(x|θ, a, νw , νe )]),

3. INFINITE KERNEL LINEAR PREDICTION

A likelihood function of the observed signal x is given by−T.
x ∼ N 0, Ψ−1 νw I→∞
i=1 θi K i + νe I Ψ

However, the first term is is still intractable. We therefore take a further lower bound L such that L ≥ L . Note that L can be indirectly
maximized by maximizing L . The updating formulas areiq(θi ).

where Ω is a PSD matrix and Υ = {Υi }I→∞
is a set of auxiliaryi=0matrices that sum to unity. By setting the partial derivatives equal tozero, we obtain the following optimal values of Ω and Υ:(18)Ω = E[νw ] i E[θi ]K i + E[νe ]I,−1 −1 −1
K i S −1 , Υ0 = E ν1S ,Υi = E ν 1θEach posterior is then iteratively updated such that an evidence lower
bound (ELBO) L is monotonically increased, where L is given bylog p(x) ≥ E[log p(x|θ, a, νw , νe )]

(19)



w i

e

−1
1
K i + E[ ν1e ]−1 I.
νw θ i(24)(25)

where S =Note that Eq. (23)iE
involves the expectations both of the parameters and of their reciprocals, i.e., the sufficient statistics are x and 1/x. Since the sufficient+ E[log p(θ)] + E[log p(a)] + E[log p(νw )] + E[log p(νe )]− E[log q(θ)] − E[log q(a)] − E[log q(νw )] − E[log q(νe )] ≡ L.

465

8kHz

8kHz

8kHz

400Hz

1.0

Observations

8kHz

Reconstructions

8kHz

Estimated spectral envelopes

8kHz

Estimated F0s (predominant kernels)

Original F0s

Detected pitched frames (indicated by

)

1s

2s

3s

4s

5s

400Hz

1.0
6s

Observations

Reconstructions

Estimated spectral envelopes

Estimated F0s (predominant kernels)

Detected pitched frames (indicated by
1s

2s

)
3s

4s

Fig. 3. Estimation results for a female-speech signal.

Fig. 2. Estimation results for a male-singing signal.
statistics of the gamma priors are log(x) and x, the generalized inverse Gaussian (GIG) distribution is a convenient choice as the functional form of posteriors [15]. We thus assume

4.2. Experimental Results
The experimental results (Fig. 2 and Fig. 3) showed the potential of
IKLP not only for joint spectral-envelope and F0 estimation but also
for detection of pitched segments. We conjecture that speech signals
can be analyzed more correctly than singing signals. In Fig. 2, IKLP
failed to detect pitched segments around 4.2s–4.6s because of relatively weak harmonic partials. The F0s were estimated accurately inthe pitched segments of both signals, even in vibrato around 6.0s in
Fig. 2. The estimated filter spectra (envelopes) were little skewed byharmonic partials because the source spectra were explicitly modeled, as shown in the second and third rows of Fig. 2 and Fig. 3.
Although the IKLP model itself is theoretically sound, we foundthat half pitch errors tend to occur because the optimization method
easily gets stuck in bad local optima when using periodic kernels.
The likelihood function given by Eq. (15) imposes a larger penalty
when the model variance Ψ−1 KΨ−T underestimates the observed
variance xxT , resulting in false (overestimated) harmonic partials
with a smaller penalty. This implies a deep connection of Eq. (15)to the Itakura-Saito divergence (ISD), which in fact acts as a costfunction when we restrict K to an identity matrix (IKLP reduces to
ISD-based AR modeling [1]). It is known that the ISD is harder to
optimize because of its nonconvexity than the Kullback-Leibler divergence (KLD) [18]. This problem could be solved by representing
the temporal dynamics of F0s and/or learning multiple models withdifferent initializations based on rough F0 estimates.

q(θi ) = GIG(θi |γi , ρi , τi ),q(νw ) = GIG(νw |γw , ρw , τw ), q(νe ) = GIG(νe |γe , ρe , τe ), (26)γ/2(ρ/τ )√
xγ−1 e−(ρx+τ /x)/2 . The updatwhere GIG(x|γ, ρ, τ ) = 2K
γ ( ρτ )ing formulas are given by

γi = α/I, ρi = 2α + E[νw ] tr(Ω−1 K i ), 
τi = E ν1w xT ΨT ΥTi K −1i Υi Ψxγw = aw , ρw = 2bw + i E[θi ] tr(Ω−1 K i ),  
τw = i E θ1i xT ΨT ΥTi K −1i Υi Ψxγe = ae , ρe = 2be + tr(Ω−1 ), τe = xT ΨT ΥT0 Υ0 Ψx.(27)The MAP estimate of a is determined such that the partial derivativeof L is equal to zero. More specifically, a∗ is a solution of theregularized normal(X T Σ−1 X + λI)a = X T Σ−1 x, equation−1T
1
1
where Σ = i E νw θi ΥTi K −1i Υi + E νe Υ0 Υ0 .
4. EXPERIMENTS
This section reports experiments that were conducted to evaluate thebasic performance of the proposed model called IKLP.

5. CONCLUSION

4.1. Experimental Conditions

We presented a nonparametric Bayesian model for joint spectralenvelope and F0 estimation in the solid framework of multiple kernel
learning. The experimental results showed that our model is robust to
high-pitched signals and can detect unpitched segments. We plan to
conduct more comprehensive experiments for comparison. There areseveral interesting directions of this research. One possibility would
be to use more precise models of the source signal  [19–21] when
designing basis functions Φ. The kernel parameters T could be optimized with type-II ML estimation (empirical Bayes). In addition,
the deep connection of IKLP to classical ISD-based LP [1] opens up
a door to the fundamental generalization of ISD-based nonnegative
matrix factorization [15, 22] (known as useful for music signal separation) in that the covariance structure within the elements of eachbasis vector can be considered by using various kernels as in [23].

We used two audio signals sampled at 16 [kHz]; One is an unaccompanied solo of a male singer (RWC-MDB-G-2001 No.91) excerpted
from the RWC Music Database [16]. We analyzed and resynthesized
the signal with STRAIGHT [6] by using ground-truth F0s. The otheris a female speech (FSUSA101) excerpted from the ATR Japanese
speech database [17]. Since these signals have high F0s over 300Hz,
higher-order LP is more difficult. We applied IKLP on a frame-byframe basis with a window size of 2048 samples (M = 2048) anda shifting interval of 160 samples. The hyperparameters were set as
α = 1.0, aw = bw = ae = be = 1.0, P = 30, and λ = 0.1. We
prepared sufficiently many periodic kernels {K i }Ii=1 that represent
different F0s ranging from 100 [Hz] to 400 [Hz] at 6-cent intervals,i.e., the truncation level was set as I = 400α.

466

[15] M. Hoffman, D. Blei, and P. Cook, “Bayesian nonparametric matrix factorization for recorded music,” in International
Conference on Machine Learning (ICML), 2010, pp. 439–446.

6. REFERENCES[1] F. Itakura and S. Saito, “Analysis synthesis telephony based on
the maximum likelihood method,” in International Congresson Acoustics (ICA), 1968, pp. C17–C20.

[16] M. Goto, H. Hashiguchi, T. Nishimura, and R. Oka, “RWC
music database: Popular, classical, and jazz music database,”in International Conference on Music Information Retrieval(ISMIR), 2002, pp. 287–288.

[2] A. El-Jaroudi and J. Makhoul, “Discrete all-pole modeling,”
IEEE Transactions on Signal Processing, vol. 39, no. 2, pp.
411–423, 1991.

[17] A. Kuramatsu, K. Takeda, Y. Sagisaka, S. Katagiri,
H. Kawabara, and K. Shikano, “ATR Japanese speech databaseas a tool of speech recognition and synthesis,” Speech Communication, vol. 9, no. 4, pp. 357–363, 1990.

[3] R. Badeau and B. David, “Weighted maximum likelihood autoregressive and moving average spectrum modeling,” in International Conference on Acoustics, Speech, and Signal Processing (ICASSP), 2008, pp. 3761–3764.
[4] M. Oudot, O. Cappé, and E. Moulines, “Estimation of the spectral envelope of voiced sounds using a penalized likelihood approach,” IEEE Transactions on Speech and Audio Processing,vol. 9, no. 5, pp. 469–481, 2001.

[18] N. Bertin, C. Févotte, and R. Badeau, “A tempering approachfor Itakura-Saito non-negative matrix factorization. with application to music transcription,” in International Conference on
Acoustics, Speech, and Signal Processing (ICASSP), 2009, pp.
1545–1548.

[5] F. Villavicencio, A. Röbel, and X. Rodet, “Improving LPC
spectral envelope extraction of voiced speech by true-envelope
estimation,” in International Conference on Acoustics, Speech,
and Signal Processing (ICASSP), 2006, pp. 869–872.

[19] M. Rothenberg, “A new inverse-filtering technique for deriving
the glottal air flow waveform during voicing,” Journal of the
Acoustical Society of America, vol. 53, no. 6, pp. 1632–1645,
1973.

[6] H. Kawahara, I. Masuda-Katsuse, and A. de Cheveigne, “Restructuring speech representations using a pitch-adaptive timefrequency smoothing and an instantaneous-frequency-based
F0 extraction: Possible role of a repetitive structure in sounds,”
Speech Communication, vol. 27, no. 3–4, pp. 187–207, 1999.

[20] D. H. Klatt and L. C. Klatt, “Analysis, synthesis and perceptionof voice quality variations among female and male talkers,”
Journal of the Acoustical Society of America, vol. 87, no. 2,
pp. 820–857, 1990.
[21] G. Fant, J. Liljencrants, and Q. Lin, “A four-parameter modelof glottal flow,” STL-QPSR, vol. 26, no. 4, pp. 1–13, 1985.

[7] T. Nakano and M. Goto, “A spectral envelope estimation
method based on F0,” in Statistical and Perceptual Audition (SAPA) - Speech Communication with Adaptive Learning(SCALE) Conference, 2012, pp. 11–16.

[22] C. Févotte, N. Bertin, and J.-L. Durrieu, “Nonnegative matrix
factorization with the Itakura-Saito divergence: With application to music analysis,” Neural Computation, vol. 21, no. 3,
pp. 793–830, 2009.

[8] A. Sasou and K. Tanaka, “Robust LP analysis using glottal
source HMM with application to high-pitched and noise corrupted speech,” in European Conference on Speech Communication and Technology (Eurospeech), 2001, vol. 4, pp. 2443–
2446.

[23] A. Liutkus, R. Badeau, and G. Richard, “Gaussian processes
for underdetermined source separation,” IEEE Transactions
on Audio, Speech, and Language Processing, vol. 59, no. 7,
pp. 3155–3167, 2011.

[9] T. Toda and K. Tokuda, “Statistical approach to vocal tract
transfer function estimation based on factor analyzed trajectory
HMM,” in International Conference on Acoustics, Speech, andSignal Processing (ICASSP), 2008, pp. 3925–3928.
[10] H. Kameoka, Y. Ohishi, D. Mochihashi, and J. Le Roux,“Speech analysis with multi-kernel linear prediction,” in
Acoustical Society of Japan (ASJ) Spring Meeting, 2-Q-24 (inJapanese), 2010.
[11] H. Kameoka, N. Ono, and S. Sagayama, “Speech spectrum
modeling for joint estimation of spectral envelope and fundamental frequency,” IEEE Transactions on Audio, Speech, andLanguage Processing, vol. 18, no. 6, pp. 1507–1516, 2010.
[12] G. Lanckriet, N. Cristianini, P. Bartlett, L. El Ghaoui, and M. I.
Jordan, “Learning the kernel matrix with semidefinite programming,” Journal of Machine Learning Research, vol. 5,
pp. 27–72, 2004.
[13] C. E. Rasmussen and C. K. I. Williams, Eds., Gaussian Processes for Machine Learning, MIT Press, 2006.
[14] H. Sawada, H. Kameoka, S. Araki, and N. Ueda, “Efficientalgorithms for multichannel extensions of Itakura-Saito nonnegative matrix factorization,” in International Conference on
Acoustics, Speech, and Signal Processing (ICASSP), 2012, pp.
261–264.

467