Available online at www.sciencedirect.com

Speech Communication 53 (2011) 658–676
www.elsevier.com/locate/specom

Computational auditory induction as a missing-data model-ﬁtting
problem with Bregman divergence q
Jonathan Le Roux a,b,⇑, Hirokazu Kameoka b, Nobutaka Ono a, Alain de Cheveigné c,
Shigeki Sagayama aaGraduate School of Information Science and Technology, The University of Tokyo, 7-3-1, Hongo, Bunkyo-ku, Tokyo 113-8656, Japanb
NTT Communication Science Laboratories, NTT Corporation, 3-1, Morinosato Wakamiya, Atsugi, Kanagawa 243-0198, Japanc
Centre National de la Recherche Scientiﬁque, Université Paris 5, and Ecole Normale Supérieure, 29 rue d’Ulm, 75230 Paris Cedex 05, France
Available online 1 September 2010

Abstract
The human auditory system has the ability, known as auditory induction, to estimate the missing parts of a continuous auditory
stream brieﬂy covered by noise and perceptually resynthesize them. In this article, we formulate this ability as a model-based spectrogram
analysis and clustering problem with missing data, show how to solve it using an auxiliary function method, and explain how thismethod is generally related to the expectation–maximization (EM) algorithm for a certain type of divergence measures called Bregman
divergences, thus enabling the use of prior distributions on the parameters. We illustrate how our method can be used to simultaneously
analyze a scene and estimate missing information with two algorithms: the ﬁrst, based on non-negative matrix factorization (NMF), performs analysis of polyphonic multi-instrumental musical pieces. Our method allows this algorithm to cope with gaps within the audio
data, estimating the timbre of the instruments and their pitch, and reconstructing the missing parts. The second, based on a recently
introduced technique for the analysis of complex acoustical scenes called harmonic-temporal clustering (HTC), enables us to performrobust fundamental frequency estimation from incomplete speech data.
 2010 Elsevier B.V. All rights reserved.
Keywords: Auditory induction; Acoustical scene analysis; Missing data; Auxiliary function; Bregman divergence; EM algorithm; Non-negative matrix
factorization; Harmonic-temporal clustering

1. Introduction
The main goal of Computational Auditory Scene Analysis (CASA) is to enable computers to imitate human auditory segregation abilities. CASA has been an area ofqPart of this work was presented at the 2008 ISCA Workshop on
Statistical and Perceptual Audition (SAPA 2008).
⇑ Corresponding author at: NTT Communication Science Laboratories,
NTT Corporation, 3-1, Morinosato Wakamiya, Atsugi, Kanagawa 2430198, Japan. Tel.: +81 46 240 4029; fax: +81 46 240 4708.
E-mail addresses: leroux@cs.brl.ntt.co.jp, leroux@hil.t.u-tokyo.ac.jp(J. Le Roux), kameoka@cs.brl.ntt.co.jp (H. Kameoka), onono@hil.t.
u-tokyo.ac.jp (N. Ono), Alain.de.Cheveigne@ens.fr (A. de Cheveigné),sagayama@hil.t.u-tokyo.ac.jp (S. Sagayama).
0167-6393/$ - see front matter  2010 Elsevier B.V. All rights reserved.
doi:10.1016/j.specom.2010.08.009

intensive research in recent years. Particular attention has
been given to solving the so-called “cocktail party problem”, the computational counterpart of the “cocktail party
eﬀect” (Cherry, 1953; von Helmholtz, 1954), i.e., the abilityof the human auditory system to focus on a single talker
within a mixture of conversations and background noise.
This has lead to the development of methods for multipitch estimation, noise canceling, or source separation(Wang and Brown, 2006). Less emphasis has been put on
the computational realization of another remarkable ability of the human auditory system, auditory induction.
Humans are able, under certain conditions, to estimate
the missing parts of a continuous acoustic stream brieﬂy
covered by noise, to perceptually resynthesize and clearlyhear them (Bregman, 1990; Kashino, 2006; Warren, 1970,

J. Le Roux et al. / Speech Communication 53 (2011) 658–676

1982). They are thus able to simultaneously analyze an
auditory scene, as in the cocktail party eﬀect, in the presence of gaps and to perceive the underlying acoustic events
as if the information inside those gaps had not been missing(whether the incomplete stimuli are actually reconstructed
at low levels in human perception or not is a diﬀerent issue,
which we shall not address here).
An eﬀective computational counterpart to this ability
would have many important engineering applications, frompolyphonic music recording analysis and restoration tomobile communications robust to both packet-loss andbackground noise. Attempts to combine scene analysis on
incomplete data and reconstruction of the missing parts
are rare, with the notable exception of (Ellis, 1993, 1996).
There have been few attempts to address this problem
through a statistical approach.
This article aims at developing such a computationalcounterpart to auditory induction, by simultaneously performing a decomposition of the magnitude wavelet spectrogram of an acoustical scene with missing or corrupted
samples, and ﬁlling in the gaps into that spectrogram. Various approaches have emerged recently which attempt toanalyze the structure of the spectrogram of an acousticalscene (Kameoka et al., 2007; Schmidt and Mørup, 2006;
Smaragdis, 2004), while on the other side gap interpolation
techniques have been the subject of research for many years(Achan et al., 2005; Cemgil and Godsill, 2005; Clark andAtlas, 2008; Esquef and Biscainho, 2006; Godsill andRayner, 1998; Lu et al., 2003; Wolfe and Godsil, 2005).
However, only few models so far try to deal with bothissues. One example is the framework developed byReyes-Gomez et al. (2004) that relies on local regularitiesof the spectrogram. The framework that we introduce
can use both local and global regularities.
We show here how statistical models that globally model
acoustical scenes can be extended for the analysis of scenes
with incomplete data. We ﬁrst derive the method for a general class of distortion functions that measure the goodnessof ﬁt between the model and the observed data. We thenshow how, for a particular class of functions called Bregman divergences (Banerjee et al., 2005; Bregman, 1967),
the method can be interpreted in terms of the expectation–maximization (EM) algorithm, enabling the use ofprior distributions on the parameters, for example toenforce local smoothness or other regularities. To illustrate
the concept, we apply it to the non-negative matrix factor
2D deconvolution algorithm (NMF2D) (Schmidt and
Mørup, 2006), and evaluate its performance on a polyphonic multi-instrumental musical piece: the proposed
method is able to analyze the scene in spite of the presenceof gaps, i.e., it can estimate the timbre of the instruments,
their pitch and the time of their activation, and separate
their contributions from those of other instruments, while
simultaneously reconstructing the missing parts. We ﬁnallyshow how to apply this method to the harmonic-temporalclustering framework (HTC) (Kameoka et al., 2007; Le
Roux et al., 2007a), and how the obtained algorithm can659

be used to perform robust fundamental frequency (F0) estimation of speech on incomplete data.
2. Computational auditory induction
2.1. Problem setting
We consider the problem of interpolating gaps in audio
signals by ﬁlling in the gaps in their magnitude spectrogram. We will not consider here the reconstruction of thephase: if the magnitude spectrogram can be accurately
reconstructed, other methods could be used to obtain aphase consistent with it (Griﬃn and Lim, 1984; Le Rouxet al., 2008b). We are interested in using local and global
regularities in the spectrogram to simultaneously analyze
the acoustical scene and ﬁll in gaps that may have occurred
into it, or in other words to perform “audio inpainting” byreconstructing missing regions of the spectrogram in thesame spirit as what is done in image inpainting (Bertalmio
et al., 2000), where diﬀusion-based (local) and exemplarbased (global) techniques are used to restore missing partsof an image (Criminisi et al., 2004).
This is analogous to what is performed by humans inauditory induction, when for example phonemes deleted
from a speech signal and replaced by louder broadband
noise can be perceptually synthesized by the brain and subjectively heard as if they were present (Bregman, 1990;
Kashino, 2006; Warren, 1970, 1982). The auditory induction phenomenon is a striking illustration of the law of closure of Gestalt psychology, according to which the humanperception system has a tendency to close “strong” perceptual forms which are incomplete, such as a circle partiallyoccluded by an irregular form. More generally, it can beconsidered as an expression of Mach’s “economy ofthought” (Barlow, 2001), in that it is more rewarding interms of simplicity of explanation to assume that someparts are occluded but actually present and need to bereconstructed (either at the primitive grouping stage or at
higher levels) than to assume that the stimuli is actuallycomposed of several disconnected parts. The fact that auditory induction does not occur if the phonemes are replacedby silence can be linked to the point made by Bregman(1990) that our perceptual system needs to be shown thatsome evidence is missing: humans can indeed see ﬁgureswith actual gaps in them, as with no special hint or trigger
mechanism, they have no reason to believe that the “missing” parts are not missing but merely hidden. The localization of the gaps will thus be considered known in thefollowing.
Most previous methods for gap interpolation focus on a
local modeling of the signal around the gaps. An important
corpus of work is for example based on auto-regressive(AR) modeling, stemming from work by Janssen et al.
(1986) on an AR interpolator which alternately maximizes
the likelihood w.r.t. the missing data and the model parameters, and later extended by Vaseghi and Rayner (1990) to
consider samples which are a pitch period apart, by Rayner

660

J. Le Roux et al. / Speech Communication 53 (2011) 658–676

and Godsill (1991) to cope with the tendency of the AR
interpolator to lead to over-smoothed interpolants whose
amplitude decreases at the center of the gap, or by Rajan
et al. (1997) to consider time-varying AR processes, among
many others. Apart from AR, sinusoidal modeling(McAulay and Quatieri, 1986) was also used by Maher(1993) for missing-data interpolation by performing theinterpolation directly on the parameters of the sinusoidal
model. More details and references can be found in(Veldhuis, 1990; Godsill and Rayner, 1998).
While previous works on gap interpolation mainly focus
on local regularities and do not attempt to explicitly modeland exploit the underlying structure of the scene, ourapproach is more global and conceptually closer to missing-feature approaches to automatic speech recognition(Barker et al., 2005; Cooke et al., 2001; Raj et al., 2004).
Trying to understand speech in the presence of noise canbe considered as a particular type of missing-data scene
analysis, in which time-frequency regions which are dominated by noise are assumed missing. There as well, the goal
is to analyze a scene (recognize its speech content) in spiteof the presence of unreliable data, and high-order knowledge given by the acoustic and language models can be
exploited to estimate the unreliable parts. We shall referto the very good reviews by Raj and Stern (2005) and Barker (2006) for more details on these methods.
Assuming we have at hand a statistical framework
which globally models an acoustical scene, we explain inthis article how to use it on scenes with incomplete data
and reconstruct to some extent the missing information.
We present two examples of statistical frameworks which
can be used in this context, Schmidt and Mørup’s NMF2D
algorithm (Mørup and Schmidt, 2006; Schmidt and
Mørup, 2006), and the Harmonic-Temporal Clustering
framework introduced in (Kameoka et al., 2007; Le Rouxet al., 2007a).
2.2. General method and applicability
The general idea is simple: given a statistical model thatcan be matched to observed data, we show how it can beused on incomplete data by iterating between analysis steps
and reconstruction steps. Furthermore, if the model is sufﬁciently speciﬁed so as to describe the underlying data
everywhere, it can be used to reconstruct the missing parts
as well. The procedure goes as follows: during a reconstruction step, the missing data are estimated based on the current value of the model; during an analysis step, the modelis updated based on the data completed during the reconstruction step. We show in the following subsection how
this iterative algorithm can be interpreted as using an auxiliary function method to optimize the ﬁtting of the statistical model on the regions where data were actuallyobserved.
Situations in which such an incomplete data frameworkneeds to be used are quite varied. One can cite for example
situations where a portion of the power spectrogram

1. is lost, for example after a packet loss during a network
communication,
2. has been discarded, for example by a binary mask
designed to suppress noise or select a particular speakerinside an acoustical scene,
3. or is simply not observed, for example because it lies
outside the observed frequency band or the time intervalof analysis.
The method we introduce can be used in general tomatch a statistical model to incomplete data based on theﬁtting on observed regions even if the original optimization
algorithm was designed to be eﬀective only on complete
data, such as Gaussian ﬁtting for example. For reconstruction purposes, however, an important point to ensure isthat the statistical model has suﬃcient “prediction power”
to interpolate the missing parts. The capacity to reconstruct the missing-data regions will indeed depend on thedesign of the model, and especially on the constraints introduced: the only guaranty is to obtain, on the wholedomain, a complete model which ﬁts the data were they
were observed. If the model is designed in such a way thatit inherently encompasses the same regularities as the datathat it is supposed to ﬁt, then we can expect that what canbe inferred on the missing data, based on these regularities,
from the observed parts of the data will naturally be reconstructed by the model in the course of the optimization.
For example, models which enforce continuity constraints would ensure a reconstruction with smooth transitions over the missing-data regions. Decomposition
models, such as NMF2D, which use information fromthe whole domain to build a lower-dimensional representation of the acoustical scene, will ensure a reconstructionthat conforms to the underlying representation. We willshow in particular in Section 4 how incomplete polyphonicmusic scenes can be analyzed with NMF2D on the basis of
information on the spectro-temporal envelopes of the notesof each instrument gathered from the non-missing portionsof the music. Similarly, a model such as HTC, through theuse of relevant prior distributions, will lead to reconstructions that are inherently guaranteed to respect Bregman’sgrouping cues (Bregman, 1990; Kameoka, 2007). We willshow in Section 6 in particular that the continuity constraint on the pitch contour of the HTC model enables us
to perform robust F0 estimation on incomplete speech data.
2.3. Auxiliary function method
Suppose one wants to ﬁt a parametric distribution to an
observed contour which is incomplete, in the sense that its
values are only known on a subset I  D # Rn , where D isthe domain of deﬁnition of the problem of interest. Suppose also that if the data were complete, the ﬁtting could
be performed (e.g., Gaussian distribution ﬁtting, etc.).
Then we show that using an iterative procedure based on
the auxiliary function method, the ﬁtting to the incomplete
data can also be performed.

J. Le Roux et al. / Speech Communication 53 (2011) 658–676

661

b ¼ argmin Lþ ðH; ^hÞ:H

ð5ÞLet f be the observed contour, and g(; H) a model
parameterized by H such that the ﬁtting of this model toan observed contour deﬁned on the whole domain D canbe performed.
We consider a distortion function d : S  S ! ½0; þ1Þ
where S # Rn , such that dðx; yÞ P 0; 8x; y 2 S and equality
holds if and only if x = y. As this function d is not required
to respect the triangle inequality, it is not necessarily a metric. For such a distortion function, we can introduce a measure of the distance between the observed data and themodel by integrating d between f and g(; H) on the subset
I:ZLðHÞ ¼
dðf ðxÞ; gðx; HÞÞdx:ð1Þ
I

In this kind of situation, it is often preferable, instead ofdeﬁning an “incomplete model” whose estimation may becumbersome, to try to fall back on a complete data estimation problem. This is what we do here by introducing an
auxiliary function. For any function h taking values in S
and deﬁned on DnI, let us deﬁneZLþ ðH; hÞ ¼ LðHÞ þdðhðxÞ; gðx; HÞÞdx:ð2Þ
DnI

As the second term on the right-hand side is itself derived
from the distortion measure, it is non-negative, and thusLðHÞ 6 Lþ ðH; hÞ;

8h:

ð3Þ

Moreover, there is equality in the inequality for h = g(; H).
The minimization procedure can now be described asfollows. After initializing H for example by performing
the distribution ﬁtting on the observed data completed by
0 on DnI, one then iteratively performs the following
updates:

The optimization process is illustrated in Fig. 1.

3. Probabilistic interpretation for Bregman divergences
We investigate in this section the probabilistic interpretation of the auxiliary function framework introduced
above in the particular case where the distortion function
is a Bregman divergence.
3.1. Relation between Bregman divergence-based
optimization and maximum-likelihood estimation
We follow (Banerjee et al., 2005; Grünwald, 2007) togive a brief overview of the concepts of exponential family
and Bregman divergence and to present the relation
between them. As a complete presentation would take us
too far from the purpose of the present discussion, we shallrefer to them for more details and for rigorous derivations.
We tried however to keep this article as self-contained as
possible.
Exponential families form a group of probability distributions which comprise many common families of probability distributions such as the normal, gamma, Dirichlet,
binomial and Poisson distributions, among others. They
are deﬁned as follows.
Deﬁnition 1. Let K be an open convex subset of Rd and let
M ¼ fP b jb 2 Kg be a family of probability distributions ona sample space X . M is an exponential family if there exist afunction f ¼ ðf1 ; . . . ; fd Þ : X ! Rd and a non-negative function r : X ! ½0; þ1Þ such that, for all b 2 K,
P w;b ðX Þ , ehb;fðX ÞiwðbÞ rðX Þ;þ

Step 1 Estimate h such that LðHÞ ¼ L ðH; hÞ:^
h ¼ gð; HÞ:H

ð4Þ

Step 2 Update H with ^h ﬁxed:ð6Þ

where hb; f(X)i is the inner product between b and f(X),andZwðbÞ ¼ logexpðhb; fðxÞiÞrðxÞ dx < þ1:X

An exponential family deﬁned in terms of a functionf = (f1, . . . , fd) is called a regular exponential family if therepresentation (6) is minimal, i.e., there exists noa0 ; a1 ; . . . P; ad 2 Rdþ1 n f0g such that for all x withdrðxÞ > 0; j¼1 aj fj ðxÞ ¼ a0 .

Fig. 1. Optimization through the iterative procedure. During the step 1,
the auxiliary parameter h is updated to ^h so that LðHÞ ¼ Lþ ðH; ^hÞ. Then,
during the step 2, Lþ ðH; ^hÞ is optimized w.r.t. H, ensuring thatb ^b 6 Lþ ð H;
hÞ < Lþ ðH; ^hÞ ¼ LðHÞ. The minimization of LðHÞ canLð HÞthus be performed through the minimization of the auxiliary functionLþ ðH; hÞ alternately w.r.t. h and H.

As an example, we consider the family of Poisson distributions {Phjh 2 (0, +1)} on the sample space X ¼ N
deﬁned as P h ðxÞ ¼ x!1 eh hx . We see that it is an exponential
family, with b = logh, f(X) = X, w(b) = eb, and r(x) = 1/x!.
The function f is not always the identity function as in the
Poisson case, as can be seen with the family of normal
distributions ffl;r2 jðl; r2 Þ 2 R  ½0; þ1Þg with fl;r2 ðxÞ ¼ðxlÞ2

1 ﬃ  2r2
pﬃﬃﬃﬃﬃﬃe2pr2, which can be seen to be an exponential family
by setting b = (l/r2, 1/(2r2)), f(X) = (X, X2) and rðxÞ ¼
pﬃﬃﬃﬃﬃﬃ
1= 2p.

662

J. Le Roux et al. / Speech Communication 53 (2011) 658–676

It can be shown (Banerjee et al., 2005; Grünwald, 2007)
that exponential families can actually be parameterized by
the mean value lðbÞ ¼ EðfðX ÞÞ of f(X). If we let
Kmean , {lj$b 2 K such that l(b) = l}, then l() is a
1-to-1 mapping from K to Kmean. Moreover, there exists
a function / : Kmean ! R such that for all b 2 K and for
all l 2 Kmean such that l = l(b),/ðlÞ þ wðbÞ ¼ hb; li;

ð7Þ

from which one can deduce in particular that b(l) = $/(l).
Altogether, by noticing thathb; fðX Þi  wðbÞ ¼ hb; li  wðbÞ þ hb; fðX Þ  li
¼ /ðlÞ þ hr/ðlÞ; fðX Þ  li;

ð8Þ

Pw,b can be rewritten parameterized by l = l(b), leading to
the so-called mean-value parameterization of the exponential family:
P /;l ðX Þ , P w;bðlÞ ðX Þ ¼ e/ðlÞhr/ðlÞ;fðX Þli rðX Þ:

ð9Þ

We will call l the expectation parameter of the exponential
family, which will be denoted by F / .
We are now ready to introduce the concept of Bregman
divergence and to derive its relation with the exponential
families.
Deﬁnition 2. Let / : S ! R be a strictly convex function
deﬁned on an open convex set S # Rd such that / isdifferentiable
on
S.
The
Bregmandivergenced / : S  S ! ½0; þ1Þ is deﬁned as

is strict, in particular when the support of the carrier r(x) is
strictly smaller than dom(/). Indeed, for all x outside that
support, Eq. (10) is veriﬁed as both members are equal tozero, but it is not informative on the relation between P/,l(x) and d/(f(x),l) as the right-hand side member is zero
only because b/(x) is. This is what happens for example
for the I -divergence (with /(l) = l log l  l) for whichdomð/Þ ¼ Rþ (extending the deﬁnition of / for l = 0).
The corresponding exponential family is the Poisson family, for which the set of instances and the support of thecarrier are only N.
The relation (10) builds a bridge between optimization
based on Bregman divergences and Maximum-Likelihood(ML) estimation with exponential families. As distribution-ﬁtting problems usually involve only a ﬁrst-order statistic, we will focus on the case f(X) = X. Trying to ﬁt a
model g(; H), deﬁned on a domain D with parameter H,
to an observed distribution f with a measure of distance
between the two based on a Bregman divergenceR d/ thenamounts to looking for H minimizing D d / ðf ðxÞ;gðx; HÞÞdx. But according to (10), this is equivalent (up tosome precautions which may have to be taken because of
the misﬁt between the domains of deﬁnition of the
Bregman divergence and the exponential family evoked
above)
to maximizing w.r.t. H the log-likelihoodRPðf ðxÞÞdx where the observed data points f(x) at
D /;gðx;HÞ
point x are assumed to have been independently generated
from P/,g(x;H).

d / ðx; yÞ ¼ /ðxÞ  /ðyÞ  hx  y; r/ðyÞi;3.2. Relation to the EM algorithm

where $/(y) is the gradient vector of / evaluated at y.

As we showed above, optimization based on a Bregman
divergence corresponds to an ML problem in which thedata are supposed to have been generated independently
from probability distributions of an associated exponential
family with expectation parameters g(x, H). We investigate
here the relation between the application of the EM algorithm to this ML problem and the auxiliary function
framework of Section 2.3, in the particular case where
the distortion function is a Bregman divergence d/ suchthat the associated exponential family F / veriﬁes f(X) = X.
The EM algorithm is based on the derivation of a socalled Q-function, which is classically obtained by considering the expectation of the log-likelihood logP(fjH) of
the observed data f against the conditional probability of
the unobserved data h with respect to the observed data
and the model with parameter H:

Bregman divergences include a large number of usefulloss functions such as squared loss, KL-divergence, logistic
loss, Mahalanobis distance, Itakura–Saito distance, andthe I -divergence. They verify a non-negativity property:d / ðx; yÞ P 0; 8x; y 2 S, and equality holds if and only if
x = y.
Banerjee et al. (2005) showed that the following informal derivation can be rigorously justiﬁed for a wide subclass of Bregman divergences, which includes inparticular all the loss functions cited above.
If P/,l is the probability density function of the regularexponential family F / (in its mean-value parameterization)
associated to the function / deﬁning the Bregman divergence d/, from (9) we have,
P /;l ðxÞ ¼ e/ðlÞþhr/ðlÞ;fðxÞli rðxÞ ¼ ed / ðfðxÞ;lÞþ/ðfðxÞÞ rðxÞlog P ðf jHÞ ¼ Eðlog P ðf jHÞÞP ðhjf ;HÞ

and eventually
P /;l ðxÞ ¼ e¼ Eðlog P ðf ; hjHÞÞP ðhjf ;HÞ

d / ðfðxÞ;lÞ

b/ ðxÞ

/(f(x))

ð10Þ

r(x). This relation holds for all
where b/(x) = e
x 2 dom(/), which can be shown (Banerjee et al., 2005)
to include the set of all the instances that can be drawn
from the distribution P/,l. However, one must be carefulwhen using this relation in certain cases where the inclusion Eðlog P ðhjf ; HÞÞP ðhjf ;HÞ
¼ QðH; HÞ  H ðH; HÞ;ð11Þ

where the functions Q and H were deﬁned in the obvious
way from the previous line. One notices through Jensen’s
inequality that

J. Le Roux et al. / Speech Communication 53 (2011) 658–676

8H; H ðH; HÞ 6 H ðH; HÞ;such that if one can update H such that QðH; HÞ >QðH; HÞ, then log P ðf jHÞ > log P ðf jHÞ.
In the problem we consider, we can show that there isactually a correspondence between the Q-function andthe auxiliary function Lþ that we introduced in Section
2.3. The computations of Appendix A indeed lead to thefollowing relation:
QðH; HÞ ¼ Lþ ðH; gðx; HÞÞ þ Cðf ; HÞ;ð12Þ

where Cðf ; HÞ does not depend on H. Computing the Qfunction, i.e., the E-step of the EM algorithm, corresponds
to computing the auxiliary function, which is done byreplacing the unknown data by the model at the current
step. Maximizing the Q-function w.r.t. H, i.e., the M-stepof the EM algorithm, corresponds to minimizing the auxiliary function w.r.t. H. This shows how to derive the auxiliary function in an EM point of view, and enables us forexample to consider prior distributions on the parametersand perform a MAP estimation.

663

In applications to audio, the matrix V to decompose is usually the magnitude or power spectrogram of the observed
signal. The horizontal and vertical dimensions of the matrices then respectively represent time and frequency (orlog-frequency), and the non-negative factorization of V is
expected to lead to a decomposition of the spectrogram inspectral templates W and their activations H. The assumption behind this decomposition is that the spectrogram ofan acoustical scene can be modeled as the repetition throughtime of characteristic spectral templates with varying amplitudes, the shape of these templates being time-invariant, i.e.,
invariant with the time at which they appear. NMFD
extends NMF by introducing a convolution in the time
direction, and looks for a decomposition of V asX!s
W s H;
ð14Þ
V K¼
s

3.3. Remark on the limitations of this interpretation

where each Ws is a set of bases and !s denotes the rightshift operator which moves each element in a matrix s columns to the right, e.g.:
0
1
0
1
0 1 2
1 2 3
!1
B
C
B
C
A ¼ @ 4 5 6 A; A ¼ @ 0 4 5 A:

We showed that the auxiliary function method in Section 2.3 could be derived through the EM algorithm inthe special case of the function d being a Bregman divergence d/ such that the associated exponential family veriﬁes f(X) = X. We shall note however that one has to pay
attention to the support of the probability distributionsof the exponential family. Indeed, as noted earlier, it may
happen that these distributions have a smaller support thanthe original set on which the Bregman divergence isdeﬁned. This is for example the case for the I -divergence,
which is deﬁned on Rþ but is associated to the Poisson distribution, whose support is N. The formulation presented
in Section 2.3 is thus more general than its EM counterpart, although it does not justify the use of penalty functions as prior distributions on the parameters. In theparticular case of the I -divergence, it is actually possible
to justify the use of the ML interpretation with continuous
data (Le Roux et al., 2008a).

NMFD thus also enables the representation of time structure in the extracted templates W, as, for each k, the set
W k ¼ ðW sm;k Þm;s of kth columns of all Ws can be considered
as a spectro-temporal template whose activation is determined by the kth row of H. NMF2D generalizes thisapproach to the frequency direction through a 2D convolution. Using a log-frequency spectrogram and assuming thatthe spectral patterns to be modeled are roughly pitchinvariant, i.e., that the spectral patterns of similar sounds
diﬀering only by their pitch are approximately equal up
to a shift on the frequency axis, NMF2D can account for
the repetition of a similar spectro-temporal structure atvarious instants and with various frequency shifts as the
convolution of a single spectro-temporal template with
the information on the time and height of its activations.
Concretely, the NMF2D model is
X X #/ !sV K¼
W s H /;
ð15Þ

7

8 9

s

4. Missing-data non-negative matrix factor 2D deconvolution
4.1. Overview of the original algorithm
The NMF2D algorithm is an extension of Smaragdis’s
non-negative matrix factor deconvolution (NMFD) (Smaragdis, 2004), itself an extension of the original non-negative
matrix factorization (NMF) (Lee and Seung, 1999). NMF
is a general tool which attempts to decompose a non-negative matrix V 2 RP0;MN in the product of two usually
lower-rank non-negative matrices W 2 RP0;MR andH 2 RP0;RN ,V  WH :

ð13Þ

0

7 8

/

where H/ is a set of activations such that the kth row of H/
corresponds to the activations of the kth spectro-temporal
template Wk pitch-shifted by / frequency bins down, and#/ denotes the down-shift operator which moves each element in a matrix / lines down. Up-shift and left-shift operators can be introduced in the same way. Applying NMF,
NMFD or NMF2D to audio signals implies making asparseness assumption on the signal, as the additivity ofmagnitudes in the spectral domain is only true if the underlying components of the signal are sparse enough to minimize overlaps.
Lee and Seung (1999) introduced eﬃcient algorithms forcomputing the NMF of a matrix V based on both the

664

J. Le Roux et al. / Speech Communication 53 (2011) 658–676

least-squares error and the I -divergence, which have been
extended by Smaragdis for NMFD (Smaragdis, 2004)
and Schmidt and Mørup for NMF2D (Mørup andSchmidt, 2006; Schmidt and Mørup, 2006). These algorithms are based on multiplicative updates. If K is deﬁned
as in (15), we deﬁne the objective function asJ ðW ; H jV Þ ¼ 12 kV  KkF j for the least-squares error,
squares
where kkF denotes the Frobenius norm (sum of the PV
of all the elements), or J ðW ; H jV Þ ¼ m;n V m;n log Km;nm;nðV m;n  Km;n Þ for the I -divergence. For the least-squares
error, the updates can be written as
P
W

s

s

W 

"/ !s >//V

H

P "/ !s />
/K H

;

H

/

P #/ >s s
W V
H  s #/ > ;
s
P
s
K
s W
/

ð16Þ

while for the I -divergence they become
 !s />P "/
V s
P #/ >
V
s
H
/ K
s W
s
s
/
/
K
W
W 
; H
H 
;
P #/ >
P !s />
s
s W 1
/1 H

in the underlying acoustical scene, assuming the scene ischaracterized by the same kind of regularity.
A sparseness-promoting penalty function can be added
to the NMF2D objective function, in the form of the Lp(quasi-)norm, for 0 < p < 2, of the matrix H, or of thisquantity raised to the power p (Kameoka et al., 2009;
Mørup and Schmidt, 2006). Concretely, we deﬁne here a
new objective function as
p
X 

J s ðW ; H jV Þ ¼ J ðW ; H jV Þ þ k
ð18Þ
H /m;n  :/;m;nThe update equations for the minimization of this objectivefunction can be obtained similarly to the ones without thesparseness term through an auxiliary function approach(Lee and Seung, 2001), the sparseness term being dealt within the same way as in (Kameoka et al., 2009) for sparseness-based NMF and complex NMF. Although we shallskip here the derivation, the updates for H become
P #/ >s s
V
/
/
s W
H 
ðleast-squares errorÞ
H
P #/ >s s/ðp1Þ
W
K
þkpH
s
ð19Þ

ð17Þ
where  denotes the Hadamard product, i.e., element-wisematrix multiplication, matrix divisions are also performed
element-wise, > denotes the matrix transposition and 1 denotes a M  N matrix with all its elements set to 1.

H/

H/ 

P

P #/ >s V s
s W
K
#/ >
s
/ðp1ÞðI -divergenceÞ

s W 1 þ kpH

ð20Þ
/(p1)4.2. Sparseness as a key to global structure extraction
As pointed out by Mørup and Schmidt (2006), there isan intrinsic ambiguity in the decomposition (15): the structure of a factor in H can to some extent be put into the signature of the same factor in W and vice versa. One way to
alleviate this ambiguity is to impose sparseness on H, thus
forcing the structure to go into W. In the case of spectrograms with missing regions, this becomes even more critical
if one expects to retrieve a “meaningful” reconstruction of
the missing parts, and sparseness becomes compulsory, as
can be clearly seen in the particular case where some timeframes are completely missing: indeed, without a sparseness term on the activations, assuming that the spectralenvelopes were time- and pitch-invariant (which is only
approximately true), a perfect reconstruction of the spectrogram with gaps could be obtained with a single frame
representing the instantaneous spectral envelope templatein W modulated by the power envelope in the time direction (gaps included) in H. If sparseness of H is enforced,
then typical spectro-temporal templates would be learnt
in W, and seeing only a part of those templates in theincomplete spectrogram would give us information on theiractivations’ time, pitch and strength, which in turn would
enable us to reconstruct the unseen parts. The role ofsparseness is thus to ensure that global and recurring structures are extracted and used throughout the spectrogram,and it will be the key that will enable us to ﬁll in the gaps/

denotes the matrix H with all elements
where H
raised to the power (p  1). To take care of the fact thatthe sparseness term could be made arbitrarily small by scaling down H and correspondingly scaling up W, leaving
J ðW ; H jV Þ unchanged, a unit-norm constraint is introduced on W. This constraint could be introduced in the
update equations using Lagrange multipliers: in the I divergence case and using L1-norm normalization, thiswould simply amount to rescaling W and H appropriately
afterwards, while for L2-norm normalization as well as in
the least-squares case, as an analytical solution cannot be
obtained, one would need to resort to numerical computations. Here, we shall use instead update equations for W
which are derived in (Mørup and Schmidt, 2006) by replacing each Wk by Wk/kWkkF in the deﬁnition of the objectivefunction, for both the least-squares and I -divergence cases.
This idea was ﬁrst investigated by Eggert and Körner(2004) for sparse NMF based on least-squares minimization. Although the convergence of these updates is not proven, it is conjectured and has been observed on extensive
tests (Eggert and Körner, 2004; Mørup and Schmidt,
2006). We reproduce here the updates for the sake ofcompleteness:P"/ !s >/

P"/ !s >
/

/

Ws

Ws 

/

V H
K H

þW s diagPþW s diag

P

s0

1

s0 1

0
"/ !s >/!

0"/ !s >
/

V H

!!!
0

Ws

k H

!

!!! ;

ðleast-sq:Þ

0

Ws

ð21Þ

J. Le Roux et al. / Speech Communication 53 (2011) 658–676
P

/

W

s

!s0 >
"/
 !s />
P
0
s
/
V
H
þW
diag
ÞW s Þ0 1ðð1 H
s
K

s

W 

P

/

!s >

1 H / þW s diag

P

s0

1

!!

!
!!! ; ðI  div:Þ"/
 !s0/>
s0
V

W
H
K

ð22Þ

where diag denotes a diagonal matrix whose elements are
given by the argument. W is normalized at the beginningof each step, before performing the updates for H and W.
4.3. Use of prior distributions
The NMF framework can be considered in a Bayesian
way based on the correspondence between Bregman divergence-based optimization and ML estimation either for the
least-squares error or the I -divergence. Indeed, the NMF
objective function can be converted into a log-likelihood(Lee and Seung, 1999; Sajda et al., 2003), to which priorconstraints on the parameters can further be added(Cemgil, 2008; Févotte et al., 2009).
Sparseness terms evoked above involving Lp (quasi-)
norms of H can be considered as such, corresponding in
general to generalized Gaussian process priors, and theLaplace distribution for p = 1. But one can also introduce
Markovian constraints on the parameters to ensure smooth
solutions. Using Gamma chains on the coeﬃcients of W and
H in the time direction, one can show that analytical update
equations can still be obtained and the objective function
can be optimized based on the expectation–constrained
maximization (ECM) algorithm (Meng and Rubin, 1993).
As explained in Section 3, the auxiliary function method
we propose can be interpreted in an EM point of view inthe special case of Bregman divergences. The use of priordistributions on the parameters will thus be justiﬁed as wellfor the missing-data version of the NMF2D algorithm,
which we describe in more details in the next section, andlocal convergence will be guaranteed.
4.4. NMF2D on incomplete spectrograms
We consider the wavelet magnitude spectrogram of an
acoustical scene represented as a non-negative matrix Vm,n,deﬁned on a domain of deﬁnition D = [[1, M]]  [[1, N]] (corresponding for example to the time-frequency region
fx; t 2 RjX0 6 x 6 X1 ; T 0 6 t 6 T 0 þ T g, sampled in timeand frequency). We assume in general that the spectro-temporal patterns to be modeled are roughly pitch-invariant,and that the signals are sparse enough such that the additivity assumption on the magnitude spectrograms holds.
We assume that some regions of the magnitude spectrogram are degraded or missing and are interested in performing simultaneously an analysis of this acousticalscene with the NMF2D algorithm despite the presence ofgaps, and a reconstruction of the missing parts.
Even if the data matrix V is incomplete, i.e., if the values
Vm,n are missing or considered not reliable for some indices(m, n) 2 J  D, due to the fact that the NMF2D update665

equations (as well as the NMF and NMFD update equations) are in fact multiplicative versions of a gradient
update, it would actually be possible to still perform theminimization of the distance taken over the observed databy computing the gradient of this restricted objective function, in the same way as was done in (Virtanen et al., 2008)
for NMF. However, the formulation of the update equations would then become more intricate and less obvious
to interpret, and, although the updates could be originallycomputed simply and eﬃciently using FFT thanks to theirconvolutive nature, their missing-data version would
require an additional trick in order to compute them inthe same way (concretely, setting to zero the values of theterm against which H or W are convolved in the denominators of (16) and (17) where data is actually missing before
computing their FFT). In any case, using the method introduced in Section 2 is cleaner and easier to interpret, more
systematic and general. Finally, the simplicity and ease ofinterpretation of NMF2D make it a good example to illustrate the general principle we presented.
Applying the method introduced in Section 2.3 to
NMF2D leads to the following algorithm, which can beused to analyze incomplete spectrograms, with both objective functions: ðpÞKm;n if ðm; nÞ 2 J
ðpþ1Þ
Step 1 V m;n¼.
V m;n if ðm; nÞ R J
Step 2 Update W through (16) or (17) and H through (19)or (20)
We note that recent work by Smaragdis et al. (2009) can
be considered as another illustration of our framework. It
is based on a spectrogram model which is very close to
NMF. In the same way as we present here, this work relieson the EM algorithm to estimate the parameters on incomplete data, and it is thus very similar to our NMF2D example. As mentioned above for the NMF2D illustration, theoptimization could have been solved with multiplicative
updates as well, although the EM interpretation is moregeneral and elegant. Several applications to the reconstruction of missing audio data are considered there, investigating in particular the use of separate training data to
improve the estimation of the bases, which we did not consider here. Irregular repartitions of the missing data arealso considered, while we focus here on situations wheresome frames are entirely missing: such situations exploit
the convolutive nature of the NMF2D model, relying on
the time structure of the bases or on information fromother notes appearing at diﬀerent pitches, and are thus
appropriate to illustrate our framework.
5. Examples of application of missing-data NMF2D
5.1. Toy example: reconstructing a 2D image
We ﬁrst tested our algorithm on simulated data used in(Mørup and Schmidt, 2006). The data, shown in Fig. 2(a),666

J. Le Roux et al. / Speech Communication 53 (2011) 658–676

Fig. 2. NMF2D with missing-data on a toy problem. (a) Original simulated data. W consists of one cross in the ﬁrst factor and one circle in the second.
They are convolved with H given in the top of the ﬁgure to yield the full data matrix V. (b) Truncated data. The truncated areas are indicated in black. (c)Estimated factors and reconstructed image using the least-squares algorithm. (d) Estimated factors and reconstructed image using the I -divergencealgorithm.

were created with W consisting of one cross in the ﬁrst factor and one circle in the second, convolved with H given in
the top of the ﬁgure to yield the full data matrix V. The
NMF2D algorithm was used in the same conditions as in(Mørup and Schmidt, 2006), with s = {0, . . . , 16}, /= {0, . . . , 16} and an L1-norm sparseness penalty with coefﬁcient 1 for the least-squares algorithm and 0.001 for the I divergence algorithm. The circle and cross templates spanroughly 15 frames in both horizontal and vertical directions, while the whole data is 200 frames wide. To construct
the incomplete data, we erased 3 frames horizontally and 2

frames every 10 frames vertically, as shown in Fig. 2(b).
Note that none of the occurrences of the structures (circle
and cross) is fully available. However, in this ideal casewhere the original data is a strict convolution of the templates, the proposed algorithm is able to extract the original templates and their occurrences and to reconstruct
the original data, as can be seen in Fig. 2(c) (least-squares
update equations) and Fig. 2(d) (I -divergence update
equations). This shows that the reconstruction is basedon global features of the data learnt by gathering information from the whole domain.

J. Le Roux et al. / Speech Communication 53 (2011) 658–676

5.2. Audio example: reconstructing gaps in a sounderasing 80 ms of signal every 416 ms, leading to a signal
with about 20 % of data missing. Its spectrogram is shownin Fig. 3(b).
The mask indicating the region J to inpaint was built
according to the erased portions of the waveform. With aGabor wavelet transform, the inﬂuence of a local modiﬁcation of the signal theoretically spans the whole interval.
However, as the windows are Gaussian, one can consider
that the inﬂuence becomes almost null further than about
three times the standard deviation. This standard deviation
is inversely proportional with the frequency, and the inﬂuence should thus be considered to span a longer interval for
lower frequencies. Although it leaves some unreliable portions of the spectrogram out of the mask in the lower frequencies, for simplicity, we did not consider here thisdependence on frequency, and simply considered unreliable, for each 80 ms portion of waveform erased, 6 whole
spectrogram frames (corresponding to about 96 ms of signal in the highest frequencies). The incomplete spectrogramis shown in Fig. 3(c), with areas to inpaint in black.

8000

8000

3200

3200

1600

1600

Frequency (Hz)Frequency (Hz)5.2.1. Experimental setting
For auditory restoration experiments, contrary to what
is done in (Schmidt and Mørup, 2006), we did not use theshort-time Fourier transform afterwards converted into a
log-frequency magnitude spectrogram, but a wavelet transform, which directly gives a log-frequency spectrogram.
More precisely, the magnitude spectrogram was calculated
from the input signals digitized at a 16 kHz sampling rateusing a Gabor wavelet transform with a time resolutionof 16 ms for the lowest frequency subband. Higher subbands were downsampled to match the lowest subband resolution. The frequency range extended from 50 Hz to
8 kHz and was covered by 200 channels, for a frequency
resolution of 44 cent.
We used a 4.8 s piece of computer generated polyphonicmusic containing a trumpet and a piano, already used in(Schmidt and Mørup, 2006). Its spectrogram can beseen in Fig. 3(a). The incomplete waveform was built by

800
400
200

667

100

800
400
200
100

50
0

1.2

2.4

3.6

50

4.8

0

1.2

Time (s)

2.4

3.6

4.8

8000

τ

τ

φ φ

Time (s)Frequency (Hz)

3200
1600
800
400
200
100
50
0

1.2

2.4

3.6

4.8

8000

8000

3200

3200

1600

1600

Frequency (Hz)Frequency (Hz)Time (s)

800
400
200
100

800
400
200
100

50

50
0

1.2

2.4
Time (s)

3.6

4.8

0

1.2

2.4

3.6

4.8

Time (s)Fig. 3. NMF2D with missing-data on the spectrogram of a truncated waveform. (a) Spectrogram of the original waveform (a mixture of piano andtrumpet sounds). (b) Spectrogram of the truncated waveform. (c) Truncated spectrogram, with truncated regions indicated in black. (d) Estimated factors
and reconstructed spectrogram using the I-divergence algorithm. (e) Reconstructed and separated spectrogram of the piano part. (f) Reconstructed andseparated spectrogram of the trumpet part.

668

J. Le Roux et al. / Speech Communication 53 (2011) 658–676

The NMF2D parameters were as follows. As in(Schmidt and Mørup, 2006), we used two factors, d = 2,
since we are analyzing a scene with two instruments, andthe number of convolutive components in pitch was set
to / = {0, . . . , 11}, as the pitch of the notes in the dataspans three whole notes. For the convolutive components
in time, we used empirically s = {0, . . . , 31}, for a time
range of about 500 ms, thus roughly spanning the lengthof the eighth notes in the music sample. The I -divergence
was used as the distortion measure, and the L1 norm as
the sparseness term with the coeﬃcient k set to 0.001.
The algorithm was ran for 100 iterations.
5.2.2. Results and discussion
To evaluate the reconstruction accuracy of the spectrogram, we use two measures: Signal to Noise Ratio(SNR), deﬁned as 10log10 ðk b
S  Sk2 =kSk2 Þ where S denotes
the reference magnitude spectrogram and b
S the reconstructed magnitude spectrogram, and Segmental SNR(SSNR), computed as the median of the individual SNRsof all the frames. We note that computing the SNR directlyon the magnitude spectrogram amounts to assuming thatthe phase is perfectly reconstructed. The results are summarized in Table 1, where “in” refers to the measure computed inside the gaps (the inpainted part), “out” to themeasure computed outside the gaps (the part more classically reconstructed based on observed data), “M” refers
to the proposed Missing-data NMF2D, “O” to the original
NMF2D on the whole data with missing data (if any)
assumed to be zero, “C” to the magnitude spectrogramof the complete waveform, and “I” to the one of the incomplete waveform. Finally, “WX” refers to the spectrogram
reconstructed by applying algorithm W on spectrogram
X, and “Y/Z” to the comparison of spectrogram Y with
spectrogram Z as a reference. For example, the SNR of“MI/C” is the SNR of the spectrogram reconstructed using
our missing-data approach on the spectrogram of theincomplete data w.r.t. the spectrogram of the fullwaveform.

Table 1Results of the reconstruction experiment. “in” and “out” refer to themeasures computed inside and outside the gaps, respectively; “M” refers
to the proposed Missing-data NMF2D, “O” to the original NMF2D on
the whole data with missing-data (if any) assumed to be zero, “C” to themagnitude spectrogram of the complete waveform, and “I” to that of theincomplete waveform. Finally, “WX” refers to the spectrogram reconstructed by applying algorithm W on spectrogram X, and “Y/Z” to thecomparison of spectrogram Y with spectrogram Z as a reference.
SNR

OC/C
I/C
OI/C
OI/I
MI/C
MI/I

SSNR

In

Out

In

Out

13.2
2.5
4.3
5.8
10.6
3.8

13.1
21.3
10.6
10.8
13.1
13.1

12.6
2.7
4.1
7.3
10.5
3.3

12.4
28.2
10.1
10.5
12.1
12.1

The OC spectrogram is the spectrogram reconstructedby the original method applied to the complete spectrogram C. Comparing it to the complete spectrogram (OC/C results) thus gives us the modeling accuracy of NMF2D
in ideal conditions, and a reference on the performancethat we should aim for when trying to analyze the incomplete scene and to reconstruct its missing parts. On theother hand, comparing the spectrogram I of the incomplete
waveform to that of the complete waveform inside the gaps
(I/C “in” results) indicates the starting point before any
reconstruction is done. Let us now look at the performanceof the original NMF2D applied to the spectrogram of theincomplete waveform, under the crude assumption thatdata in the gaps are equal to zero. Comparing the reconstructed spectrogram OI to either that of the complete
waveform (OI/C) or the incomplete waveform (OI/I) showsﬁrst that a bias is introduced even in the reconstruction of
the observed data (“out”), and second that, as expected,
the missing data are not reconstructed (OI/C “in”). We
ﬁnally look at the performance of the proposed missingdata NMF2D also applied to the spectrogram of theincomplete waveform. Comparing the reconstructed spectrogram MI to that of the incomplete waveform (MI/I),
we see that the proposed algorithm correctly performs its
task of reconstructing the observed data (“out”). This
result is important in itself as it shows that the proposed
framework enables NMF2D, designed for complete data,
to be used on incomplete data without decrease of the performance measured on the observed part of the data, inparticular without letting the missing regions introduce abias in the analysis of the observed regions. One could
actually think of applications for which reconstruction of
the missing parts may be unnecessary, for example if only
the spectro-temporal templates or their activations themselves are desired. Comparing now MI to the spectrogramof the complete waveform (MI/C), we see that the formerly
erased regions (“in”) are correctly inpainted, with a great
improvement over the incomplete spectrogram, as seen earlier with the I/C results, and that our method performs closely to NMF2D applied on the complete spectrogram, as
we saw above with the OC/C results.
Graphical results are shown in Fig. 3(d)–(f), where one
can see in particular that the acoustical scene analysis(i.e., the learning of a spectro-temporal template for each
instrument and the estimation of the pitch and onset timeof each note) is performed correctly, and that blind source
separation is also performed in spite of the presence ofgaps.
Here again, it is interesting to note that the missing
information is not reconstructed from the neighboringparts as in classical interpolation techniques, but indirectlyfrom similar patterns in other regions of the spectrograms,
using the local information mainly to determine what similar patterns to use in the reconstruction. Although one
may think that in this particular simple example, where
pitches are ﬂat and spectral envelopes rather steady, interpolation techniques could be used as well, it is important toJ. Le Roux et al. / Speech Communication 53 (2011) 658–676

note that such techniques are helpless when higher-order
structure is necessary to reconstruct the missing regions:this is for example the case when these regions include
the beginning or the end of a note and typical note length
is thus a crucial key, when a complete harmonic is missing
and needs to be reconstructed by inferring it from the typical spectro-temporal envelope and the strength of the current note, or in general when the missing information is toocomplex to be inferred from the neighboring parts only.
One can relate this behavior of our model to the phonemic
illusion phenomenon mentioned earlier, where humans
arguably also use higher-order models (e.g., language models, production models) to infer the missing phonemes.
6. Missing-data harmonic-temporal clustering
The Harmonic-Temporal Clustering framework attempts
to perform the analysis of an acoustical scene by modeling its
power spectrogram as a constrained Gaussian mixture
model. It has been introduced by Kameoka et al. (2007) formusic signals and later extended to signals with continuously
varying pitch such as speech by Le Roux et al. (2007a). We
explain here how it can be further extended to deal withacoustical scenes with incomplete data, and how the inherent
continuity constraint on the fundamental frequency contour
imposed by the use of cubic spline functions and an extra
Markovian prior enables us to perform robust F0 estimation
on incomplete speech data.

669

where the weightPparameters
P wk, vkn and ukny are
P normalized such thatk wk ¼ 1;
n vkn ¼ 1; 8k andy ukny ¼ 1;
8k; n. The parameter sk gives the onset time of the source
model qk, wk its energy, vkn the ratio of energy inside its
nth partial, and the parameters ukny together with the duration parameter /k determine the shape of the temporalenvelope of this nth partial. A graphical representation of
an HTC source model qk(x, t; H) can be seen in Fig. 4.
The F0 contours lk(t) can be expressed using piece-wise ﬂatfunctions or cubic spline functions according to the signal
to be modeled.
The goal is toPminimize the diﬀerence between W(x, t)
K
and Qðx; t; HÞ ¼ k¼1 qk ðx; t; HÞ according to a certain criterion. We use the I -divergence (Csiszár, 1975) as a classical way to measure the distance between two non-negative
distributions:
Z Z
W ðx; tÞI ðW jQðHÞÞ ,W ðx; tÞ log
Qðx;
t; HÞ
D
 W ðx; tÞ  Qðx; t; HÞÞ dx dt;

ð25Þ

and we are thus looking for Hopt ¼ argmin I ðW jQðHÞÞ.
H

6.1. Overview of the HTC model
Consider the wavelet power spectrum W(x, t) of a signal
recorded from an acoustical scene, deﬁned on a domain of
deﬁnition
D ¼ fx; t 2 RjX0 6 x 6 X1 ; T 0 6 t 6 T 0 þ T g.
The problem considered is to approximate the power spectrum as well as possible as the sum of K parametric source
models qk(x, t; H) modeling the power spectrum of K
“objects” each with its own F0 contour lk(t) and its own
harmonic-temporal structure. As described in (Kameoka
et al., 2007; Le Roux et al., 2007a), the source models
qk(x, t; H) are expressed as a Gaussian mixture model(GMM) with constraints on the characteristics of the kernel distributions: supposing that there is harmonicity with
N partials modeled in the frequency direction, and thatthe power envelope is described using Y kernel functionsin the time direction, we can rewrite each source model inthe formqk ðx; t; HÞ ¼

N X
Y 1
X

S kny ðx; t; HÞ;ð23Þ

n¼1 y¼0

where H is the set of all parameters and with kernel densities Skny(x, t; H) which are assumed to have the following
shape:
wk vkn ukny eS kny ðx; t; HÞ ,2prk /kðxlk ðtÞlog nÞ2 ðtsk y/k Þ2

2r2
2/2
k
k

;

ð24Þ

Fig. 4. Graphical representation of an HTC source model. (a) shows thetime-frequency proﬁle of the model, while (b) shows a cross-section of themodel at constant time and (c) the evolution in time of the power envelope
function. The harmonic structure of the model can be seen in (b), and theapproximation of the power envelope in the time direction as a sum ofGaussian kernels can be seen in (c).

670

J. Le Roux et al. / Speech Communication 53 (2011) 658–676

6.2. Formulation of the model on incomplete data
The optimization process in HTC (Kameoka et al.,2007; Le Roux et al., 2007a) is nothing else than the ﬁttingof a model (in particular a constrained Gaussian mixture
model) to an observed distribution (the wavelet power
spectrum of an acoustical scene), using the I -divergence
as a measure of the goodness of ﬁt.
If some parts of the power spectrum are missing or corrupted, or if some parts of the HTC model are partially or
entirely lying outside the boundaries of the spectrogram (forexample if some harmonics of the model are above the maximum frequency and a prior is used to link the powers of theharmonics, ﬁtting the upper harmonics to zero will bias the
optimization), the estimation of the HTC model must be
performed under an incomplete-data framework, as in Section 2.3. In the same way as we showed there, optimization
can be performed in an iterative way by using the values ofthe model at the previous step as an estimation of the unobserved data. In the case of HTC, this results in a hierarchicalalgorithm with two levels. At the upper level is the iterativealgorithm described above. At the lower level, inside thestep 2 of the upper level, the EM algorithm is used as in
the classical formulation of the HTC optimization. Let W
be the observed part of the spectrogram and I  D its
domain of deﬁnition. The objective function to minimize
here is the same as (25) but restricted to the domain where
the spectrogram is observed:
Z Z
W ðx; tÞI ðW ; QðHÞÞ ,W ðx; tÞ log
Qðx;
t; HÞ
I
 ðW ðx; tÞ  Qðx; t; HÞÞ dx dt:

ð26Þ

We deﬁne the auxiliary function asI þ ðW ; V ; QðHÞÞ , I ðW ; QðHÞÞZ Z
V ðx; tÞV ðx; tÞ logþQðx; t; HÞ
DnI ðV ðx; tÞ  Qðx; t; HÞÞ dx dt:

ð27Þ

Then membership functions m can be further introduced asin the classical formulation of HTC to build the ﬁnal auxiliary function I þþ ðW ; V ; QðHÞ; mÞ. These membershipfunctions
P are non-negative and sum up to 1 for each(x, t) : k;n;y mkny ðx; tÞ ¼ 1. If we noteW ðx; tÞ if ðx; tÞ 2 I
Zðx; tÞ ¼:
V ðx; tÞ if ðx; tÞ 2 D n I
we deﬁneI þþ ðW ; V ; QðHÞ; mÞZ Z X
S kny ðx; t; HÞ,mkny ðx; tÞZðx; tÞ log
mkny ðx; tÞZðx; tÞDk;n;y ðZðx; tÞ  Qðx; t; HÞÞ dx dt:

ð28ÞUsing the concavity of the logarithm, one can see thatð29ÞI þ ðW ; V ; QðHÞÞ 6 I þþ ðW ; V ; QðHÞ; mÞwith equality for
S kny ðx; t; HÞ^ kny ðx; tÞ ¼ PK PN PY 1
:
ð30Þ
m
k¼1
n¼1
y¼0 S kny ðx; t; HÞAltogether, the optimization process can be formulated asfollows.
Step 1 Estimate V such that I ðW ; QðHÞÞ ¼ I þ ðW ; V ; HÞ:Vb ðx; tÞ ¼ Qðx; t; HÞ; 8ðx; tÞ 2 D n I:
ð31ÞStep 2 Update H with Vb ﬁxed:
b ¼ argmin I þ ðW ; Vb ; HÞ:
H

ð32Þ

H

To do so, perform one iteration of the classical formulationof HTC:
E-Step
S kny ðx; t; HÞ^ kny ðx; tÞ ¼ PK PN PY 1
m
;
k¼1
n¼1
y¼0 S kny ðx; t; HÞð33Þ

M-Step
b ¼ argmin I þþ ðW ; Vb ; H; mÞ^  log P ðHÞ
H

ð34Þ

H

where P(H) is a prior distribution on the parameters.

6.3. Optimization of the model
Analytical update equations for the M-step are derivedin (Kameoka et al., 2007; Le Roux et al., 2007a). However,when the F0 contour is modeled using cubic spline functions, which is relevant for speech or musical instruments
whose pitch can vary continuously, the spline parameters
were updated in (Le Roux et al., 2007a) not globally but
one after the other. The corresponding optimization procedure, called the expectation–constrained maximization
algorithm (ECM) (Meng and Rubin, 1993), does not
ensure the minimization in the M-step but nonethelessguarantees the decrease of the objective function. This
spline parameter update was thus not optimal but yet led
to very good results in F0 estimation accuracy. However,it suﬀered from some instability problems in long regionswith low harmonic energy (silence or unvoiced parts ofspeech for example). When dealing with missing-data problems, such issues become critical, and we thus need to use
better update equations for the spline parameters, brieﬂy
introduced in (Le Roux et al., 2007b) and which we present
in detail in Appendix B. Contrary to the update equations
previously described in (Le Roux et al., 2007a), the onespresented here are global analytical update equations
which lead to the minimum of the auxiliary function in
the M-step. They ensure a greater stability of the spline
model and a better F0 estimation accuracy, as shown in
the next section.

J. Le Roux et al. / Speech Communication 53 (2011) 658–676

7. F0 estimation on incomplete data with HTC
7.1. Importance of F0 estimation accuracy
Ensuring a very good accuracy for the F0 estimation is
not only important as a necessary step for computational
auditory induction by HTC, but it is also in itself a primaryissue. Indeed, being able to estimate the F0 accurately isimportant as well for some previous audio interpolation
methods such as Maher’s sinusoidal model based method(Maher, 1994), in which the harmonics before and after
the gap need to be linked, or Vaseghi and Rayner’s
extended AR model (Vaseghi and Rayner, 1990), whichtakes advantage of the long-term correlation structure of
the signals by introducing extra predictor parameters
around the pitch period.
7.2. Relevance of HTC’s F0 contour model
When applied to speech, HTC is based on a spline F0
contour. A Markovian prior, presented in Appendix B, isused on the parameters of the contour to ensure that it will
not move too abruptly. This Markovian prior penalizes thedeviation of a spline parameter from the linear interpolation of its neighbors. Altogether, HTC’s F0 contour model
is somewhere between a spline interpolation and a linear
interpolation, depending on the strength of the matching
between the HTC source model and the observed data.
Attempting to use this model for reconstruction ofincomplete data implies that the F0 contour inside thegap is close to an interpolation based on the values of
the contour outside the gap. To conﬁrm the relevance of
this interpolation, we thus need to ensure that, assumingthat the F0 estimation on complete parts of the data is
accurately performed, the F0 of the missing parts of thedata is accurately performed as well.
We thus evaluated as a preliminary experiment the accuracy of the F0 contour obtained by interpolating the reference F0 values outside the gaps on the whole interval using
both natural splines and linear interpolation. This can beconsidered as an evaluation of what can be expected by
HTC.
We then conducted experiments to conﬁrm the accuracyof the proposed method for F0 estimation. The goal here
was ﬁrst to conﬁrm that the new spline update equationsindeed outperform the former update equations on single-speaker F0 estimation in clean environment for complete data, then to evaluate the F0 accuracy with parts ofthe data missing.
7.3. Experimental setting
The general conditions of the experiments were exactly
the same as in (Le Roux et al., 2007a), and we shall brieﬂyreview them here, all the details being given there.
We used a database of speech recorded together with alaryngograph signal (Bagshaw et al., 1993), consisting of671

one male and one female speaker who each spoke 50 English sentences for a total of 5 min and 37 s of speech, for the
purpose of evaluation of F0-estimation algorithms. The
power spectrum W(x, t) was calculated from an input signal
digitized at a 16 kHz sampling rate (the original data of thedatabase was converted from 20 kHz to 16 kHz) using a
Gabor wavelet transform with a time resolution of 16 ms
for the lowest frequency subband. Higher subbands were
downsampled to match the lowest subband resolution.
The lower bound of the frequency range and the frequency
resolution were respectively 50 Hz and 14 cent. The spline
contour was initially ﬂat and set to 132 Hz for the malespeaker and 296 Hz for the female speaker. The length of
the interpolation intervals was ﬁxed to 4 frames. For
HTC, we used K = 10 source models, each of them with
N = 10 harmonics. We used as ground truth the F0 estimates and the reliability mask derived by de Cheveigné
and Kawahara (2002). As the spline function gives an analytical expression for the F0 contour, we compare our resultwith the reference values at a sampling rate of 20 kHz
although all the analysis was performed with a time resolution of 16 ms. Deviations over 20% from the reference weredeemed to be gross errors.
For the incomplete data, we prepared four sets of data
by replacing segments of the utterances of diﬀerent lengths
by silence. The sets are prepared such that approximately
20% of the data is lost, erasing segments of length L msevery 5L ms of data (the obtained utterances would thenbe successions of L ms of silence, 4L ms of speech, L msof silence, etc.). The four lengths we selected are 25 ms,
50 ms, 75 ms and 100 ms, and the corresponding data setsare denoted by Erase-25 ms, Erase-50 ms, Erase-75 ms
and Erase-100 ms, respectively. For example, the utterances in Erase-50 ms were produced by erasing 50 ms every
250 ms of data, leading to utterances which are successionsof 50 ms of silence, 200 ms of speech, 50 ms of silence, etc.
We shall note that gaps from 30 ms to 50 ms are alreadyconsidered very long gaps in the audio restoration literature (Esquef and Biscainho, 2006; Godsill and Rayner,
1998; Maher, 1994).
7.4. Preliminary experiment on F0 interpolation
We ﬁrst performed a preliminary experiment based on
the reference F0 and the reliability mask derived in (deCheveigné and Kawahara, 2002). The reliability mask
was used to determine the voiced regions of the speechutterance, and a global contour over the whole utterance
was derived by interpolating the values of the F0 reference
which were both inside the reliability mask and outside theerased segments of the data. We used both linear interpolation and cubic spline interpolation. We then computed
the gross error rates of the interpolated F0 values inside
the gaps (by construction the values outside the gaps areequal to the reference and thus no error can occur there).
Results for the incomplete data sets can be seen in Table 2.
Spline interpolation does not perform as well as linear

672

J. Le Roux et al. / Speech Communication 53 (2011) 658–676

Table 2Gross error rates (%) for F0 interpolation inside the gaps based onreference F0.
Data set
Erase-25 m s

Erase-50 m s

Erase-75 m s

Erase-100 m s

Cubic splines

Linear interpolation

Male
Female

0.9
0.2

0.0
0.0

Total

0.5

0.0

Male
Female

1.5
1.0

0.5
0.5

Total

1.2

0.5

Male
Female

4.0
3.0

1.3
0.5

Total

3.5

0.9

Male
Female

7.6
6.2

4.2
1.6

Total

6.9

2.9

interpolation due to the large variations that can occur
depending on the slope of the contour at the beginningor end of a voiced region. This is precisely what the Markovian prior in HTC’s F0 contour model aims to avoid.
7.5. Accuracy on complete data
We ﬁrst used the classical HTC formulation on complete
data, using the new spline update equations. Here, only
step 2 of the algorithm devised in Subsection 6.2 is thus
used (iteration of equations (33) and (34)).
The results can be seen in Table 3, with for comparison
the results obtained with HTC using the former spline
update equations as well as the ones obtained with thestate-of-the-art algorithm YIN (de Cheveigné and Kawahara, 2002). We note that we obtained 2.1% gross error
rate for YIN using the code made available by its authors,as opposed to 1.3% reported in the original paper. We cansee that HTC with the newly proposed spline update equations now performs comparably to YIN.

Table 3Gross error rates for F0 estimation on complete data (clean single-speakerspeech).
MethodYIN

HTC (former spline update)

HTC (proposed spline update)

Gross error (%)
Male
Female

3.2
1.0

Total

2.1

Male
Female

3.2
3.7

Total

3.5

Male
Female

1.1
1.3

Total

1.2

7.6. Accuracy on incomplete data
The wavelet transforms were performed on the truncated waveforms of the data sets introduced above. The
regions DnI which are to be considered missing in the spectrogram were deﬁned as the frames corresponding to theerased parts of the waveform. The inﬂuence of the erased
portion is larger for low frequencies, but we neglect thisand consider missing a whole frame regardless of the frequency bin.
In such situations where part of the data is irrelevant,one might think that algorithms which perform F0 estimation more locally should be used, using interpolation
between the preceding and following voiced portions to
obtain F0 values inside the gap. If the estimation can be
accurately performed outside the gaps, such a methodshould lead to very good results, as we saw in Section
7.4. However, one needs to note that if such algorithms
are used, a robust voice activity detection (VAD) must be
performed as well to determine which points should be
used in the interpolation. A poor VAD accuracy could leadto very bad results in the interpolation process, as unreliable values for the F0 could be used as base points forthe interpolation, leading to wrong results on the wholeinterpolation region. To illustrate this and as a comparison
with HTC, we used the algorithm YIN to perform F0 estimation outside the gaps, and used a linear interpolation to
obtain values inside the gaps, using the closest voiced
regions outside the gaps as boundaries. The positions of
the gaps were given, and the voiced regions were determined using the aperiodicity measure given by YIN, witha threshold of 0.2. The results given here were obtainedusing linear interpolation, but cubic spline interpolation
gave similar results.
Results for HTC and YIN are given in Table 4, withgross error rates for the whole ﬁle as well as for the erased
segments only. We can see that the performance of HTC
degrades as the gaps become longer. HTC performs better
than the algorithm based on YIN for the total accuracy aswell as for the accuracy inside the gaps with 25 ms and50 ms erased segments, while the algorithm based on
YIN performs better inside the gaps with 75 ms and
100 ms erased segments but is still outperformed on the
total error.
These results raise several remarks. We note ﬁrst thatthe accuracy of YIN on the whole waveform is stable asgaps become longer while decreasing inside the gaps, meaning that it tends to increase outside the gaps. Although thismay ﬁrst sound surprising, it is related to the fact thatinterpolation is performed by using as anchors the closest
regions outside the gaps with a suﬃciently low aperiodicity
measure; the presence of gaps will thus inﬂuence the interpolation process and the F0 estimates around them. But asgaps become longer, the results inside and outside becomeless dependent: each gap inducing a loss of information
which inﬂuences a small neighborhood through interpolation, many small gaps are likely to harm the estimates

J. Le Roux et al. / Speech Communication 53 (2011) 658–676
Table 4Gross error rates for F0 estimation on incomplete data with HTC and YIN(clean single-speaker speech). The results for YIN are indicated inparentheses.
Data set
Erase-25 ms

Erase-50 ms

Erase-75 ms

Erase-100 ms

673

simple example how the proposed framework was able tosimultaneously perform acoustical scene analysis and gap
interpolation in a musical piece with NMF2D, and how arobust F0 estimation could be performed on incomplete
data with HTC. While we assumed here that the gap locations were known, future work will investigate their joint
estimation together with the model parameters and themissing data, in a similar way to Barker et al. (2005) for
missing-feature speech recognition.

Error in the gaps (%)

Total error (%)

Male
Female6.0 (12.0)
4.6 (4.7)

3.9 (10.0)
3.1 (3.3)

Total

5.3 (8.3)

3.5 (6.5)

Male
Female8.8 (14.2)
6.9 (5.7)

4.2 (9.5)
2.7 (3.2)

Total

7.9 (9.9)

3.4 (6.3)

Acknowledgments

Male
Female14.1 (15.5)
13.4 (6.4)

4.5 (9.5)
4.1 (3.3)

Total

13.7 (10.9)

4.3 (6.3)

We would like to thank the anonymous reviewers fortheir insightful comments which greatly helped improvethe clarity and the general quality of this article.

Male
Female22.1 (19.2)
19.5 (7.4)

7.0 (10.3)
5.6 (3.6)

Appendix A. Derivation of the Q-function of Section 3.2

Total

20.9 (13.5)

6.3 (6.9)

outside the gaps more than a few long ones. On the other
hand, HTC is inﬂuenced in a diﬀerent way by the lengthof the gaps. Results outside the gaps are stable, while accuracy inside the gaps decreases faster than for YIN whengaps become too long. One reaches here the limits of theabove-mentioned “prediction power” of the model inside
the gaps, and its degree of freedom becomes too high tosystematically converge to a relevant solution. This problem could be coped with by investigating the introductionof more complex pitch contour models which try to encompass the long-term dynamics of the F0, such as the Fujisaki
model (Fujisaki and Nagashima, 1969) for example, ormore complex priors on the spline contours, although inboth cases one faces the risk of making the optimization
intractable.
Altogether, the results show that HTC’s F0 estimation
accuracy, while degrading in extreme cases, is very goodeven in the presence of long gaps, and that, although other
F0 estimation algorithms could be used as well, it is not
obvious, regardless of their performance on complete data,
whether they can be turned into eﬀective algorithms on
incomplete data, due in particular to the importance of arobust VAD for the interpolation to be eﬀective.

We compute here the explicit form of the Q-function
QðH; HÞ involved in the application of the EM algorithm
to the ML problem introduced in Section 3.2 in which
the data are supposed to have been generated independently at each x from a probability distribution of an exponential family F / with expectation parameter g(x, H). In
the following, we will denote by mx,/,H(z) the density of thisprobability distribution, which can be written as explainedin Section 3.1 directly using the corresponding Bregman
divergence:mx;/;H ðzÞ ¼ ed / ðz;gðx;HÞÞ b/ ðzÞ:As the data are supposed to have been generated independently at each x from the probability distribution withdensity mx,/,H(z), observed and unobserved data are in particular independent conditionally to H, and the Q-function
can be written as follows:QðH; HÞ ¼ Eðlog P ðhjHÞÞP ðhjf ;HÞ þ Eðlog P ðf jHÞÞP ðhjf ;HÞ
Z¼
Eðlog P ðhðxÞjHÞÞP ðhðxÞjHÞ dx
Rn nI

¼

We presented a computational framework to model
auditory induction, i.e., the human auditory system’s ability to estimate the missing parts of a continuous auditory
stream brieﬂy covered by noise, by extending acousticalscene analysis methods based on global statistical models
such as HTC and NMF2D to handle unobserved data.
We related the method to the EM algorithm, enabling
the use of priors on the parameters. We illustrated on a

Z

þ
Z

Rn nI

¼

Z
P ðhjf ; HÞlog P ðf ðxÞjHÞdx
IZ
Z
mx;/;H ðzÞ log mx;/;H ðzÞ dz dx þlog P ðf ðxÞjHÞdx
I

Z

Z
Rn nI



mx;/;H ðzÞ log b/ ðzÞ  d / ðz; gðx; HÞÞ dz dx

ðA:2Þ

Z



log b/ ðf ðxÞÞ  d / ðf ðxÞ; gðx; HÞÞ dxZZmx;/;H ðzÞd / ðz; gðx; HÞÞdz dx
¼
Rn nI
Z
 d / ðf ðxÞ; gðx; HÞÞdx þ C 1 ðf ; HÞ;þ

8. Conclusion

ðA:1Þ

I

I

where C 1 ðf ; HÞ does not depend on H. If we now rewrited/(z, g(x; H)) asd / ðz; gðx; HÞÞ ¼ d / ðgðx; HÞ; gðx; HÞÞ þ /ðzÞ  /ðgðx; HÞÞ hz  gðx; HÞ; r/ðgðx; HÞÞi;

ðA:3Þ

674

J. Le Roux et al. / Speech Communication 53 (2011) 658–676

we can simplify the ﬁrst term in Eq. (A.2):
Zmx;/;H ðzÞd / ðz; gðx; HÞÞdz
Z¼

1
ððtiþ1  tÞei þ ðt  ti Þeiþ1
tiþ1  ti
mx;/;H ðzÞdz d / ðgðx; HÞ; gðx; HÞÞZAðtÞ ¼ðt  ti Þðtiþ1  tÞ>
ðtiþ2  tÞM >i þ ðt  t i1 ÞM iþ1
6
ðB:3Þ

ðz  gðx; HÞÞmx;/;H ðzÞdz; r/ðgðx; HÞÞ

where Mj denotes the jth row of the matrix M and ej denotes the jth vector of the canonical basis. We note furthermore that A(t) = $zl(t; z).

þ C 2 ðf ; HÞ ¼ d / ðgðx; HÞ; gðx; HÞÞ þ C 2 ðf ; HÞ;
where C 2 ðf ; HÞ does not depend on H. To lead the calculation above, we used the fact that the mass of a probability
distribution of an exponential family with expectation
parameter gðx; HÞ is 1 and its mean is gðx; HÞ:Zmx;/;H ðzÞdz ¼ 1;
ðA:4ÞZzmx;/;H ðzÞdz ¼ gðx; HÞ:
ðA:5Þ
We then obtain for the Q-functionZd / ðgðx; HÞ; gðx; HÞÞdx
QðH; HÞ ¼ 

B.2. Optimization of the objective function
During the M-step of the EM algorithm, one wants to^  log P ðHÞ
minimizeJ ðHÞ ¼ I þþ ðW ; Vb ; H; mÞwithrespect to H. We can compute the gradient with respect
to z:
Z Z X
‘kny ðx; tÞðx  lðt; zÞ  log nÞAðtÞdx dt
rz J ¼ 
r2k
D k;n;y rz log P ðHÞZ Z X‘kny ðx; tÞ>¼ðx  AðtÞ z  log nÞAðtÞdx dt
r2k
D k;n;y

Rn nI



Z

d / ðf ðxÞ; gðx; HÞÞdx þ Cðf ; HÞIþ

¼ L ðH; gðx; HÞÞ þ Cðf ; HÞ;
where Cðf ; HÞ again does not depend on H.

Appendix B. Derivation of the spline parameter update
equations in HTC
B.1. Spline contour
The analysis interval is divided into subintervals [ti, ti+1)
of equal length . The parameters of the spline contour
model are then the values zi of the F0 at each bounding
point ti. Assuming that the second derivative vanishes atthe bounds of the analysis interval leads to the so-called
natural splines. Under this assumption, one can explicitly
compute oﬄine a matrix M linking the values z00i of the second derivative of the contour at ti with the values zi, suchthat z00 = Mz. An analytical expression for the contourl(t; z) as a concatenation of third order polynomials canthen be classically obtained. For t 2 [ti, ti+1):1ðzi ðtiþ1  tÞ þ ziþ1 ðt  ti Þtiþ1  ti
1 ðt  ti Þðtiþ1  tÞ ðtiþ2  tÞz00i þ ðt  ti1 Þz00iþ1 :6lðt; zÞ ,ðB:1Þ
One can notice that the expression of l(t; z) is actually linear in z:>lðt; zÞ ¼ AðtÞ z

 rz log P ðHÞ

ðA:6Þ

ðB:2Þ

where A(t) is a column vector such that, for t 2 [ti, ti+1),‘kny(x, t) = mkny(x, t)Z(x, t). Note that the term
RwhereRQðx;t; HÞdx dt in (28) does not contribute to the gradiD
ent w.r.t. z as the spline parameters do not inﬂuence thenormalization of the model. Let
Z X‘kny ðx; tÞ/ðtÞ ¼
ðx  log nÞdx;
r2k
k;n;y
Z X
‘kny ðx; tÞdx:cðtÞ ¼
r2k
k;n;y
Then
rz J ¼ 

ZZ/ðtÞAðtÞdt þcðtÞAðtÞAðtÞ> dt z

 rz log P ðHÞ:One can then obtain the Hessian matrix:Z
H z J ¼ cðtÞAðtÞAðtÞ> dt  H z log P ðHÞ:ðB:4Þ

If one uses a Markov assumption on the spline parameterswith Gaussian distributions for the state transitions, theprior distribution becomesP ðzÞ ¼ P ðz0 Þ

jzjYP ðzj jzj1 Þ;

j¼1

with z0 following a uniform distribution and
ðz z

Þ2

j j1

1
P ðzj jzj1 Þ ¼ pﬃﬃﬃﬃﬃﬃ e 2r2s :2prs

J. Le Roux et al. / Speech Communication 53 (2011) 658–676

Then
0

1
B 1
B
1B
rz log P ðHÞ ¼  2 B
rs B
B
@

1
2
..
.

0
¼ H z log P ðHÞz:

0

1

1
..
.

..

1

2

C
C
C
Cz
C
C
1 A

1

1

.

ðB:5Þ

Putting to 0 the gradient w.r.t. z, one can ﬁnd the update
equation for z:Z1/ðtÞAðtÞ dt:
ðB:6Þ
z ¼ ðH z J Þ
The convexity can be studied by looking at H z J in Eq.
(B.4). The ﬁrst term is indeed non-negative, as c(t) P 0, "t.
For the second term, coming from the prior distribution,
we recognize a tridiagonal matrix, for which the principal
minors can be easily calculated. If T = (tij) is a tridiagonal
matrix and an its nth principal minor, thenan ¼ tn;n an1 þ tn;n1 tn1;n an2 :ðB:7ÞIn our case, we see that the principal minors of Hzlog P(H)
are all non-positive. The matrix Hzlog P(H) is thus positive
semi-deﬁnite. Altogether, H z J is at least positive semi-deﬁnite, and the update (B.6) thus corresponds to a minimum.
References
Achan, K., Roweis, S., Hertzmann, A., Frey, B., 2005. A segment-based
probabilistic generative model of speech. In: Proc. Internat. Conf. on
Acoust. Speech Signal Process. (ICASSP), vol. 5. pp. 221–224.
Bagshaw, P.C., Hiller, S.M., Jack, M.A., 1993. Enhanced pitch tracking
and the processing of F0 contours for computer and intonationteaching. In: Proc. Eur. Conf. on Speech Comm. Technol. (Eurospeech). pp. 1003–1006.
Banerjee, A., Merugu, S., Dhillon, I.S., Ghosh, J., 2005. Clustering with
Bregman divergences. J. Machine Learn. Res. 6, 1705–1749.
Barker, J., Cooke, M.P., Ellis, D.P.W., 2005. Decoding speech in thepresence of other sources. Speech Comm. 45 (1), 5–25.
Barker, J.P., 2006. Robust automatic speech recognition. In: Wang, D.-L.,
Brown, G.J. (Eds.), Computational Auditory Scene Analysis: Principles Algorithms, and Applications. IEEE Press/Wiley, pp. 297–350.
Barlow, H.B., 2001. The exploitation of regularities in the environment by
the brain. Behav. Brain Sci. 24 (4), 602–607.
Bertalmio, M., Sapiro, G., Caselles, V., Ballester, C., 2000. Image
inpainting. In: SIGGRAPH 2000, Proc. Comput. Graphics. pp. 417–
424.
Bregman, A.S., 1990. Auditory Scene Analysis. The MIT Press.
Bregman, L.M., 1967. The relaxation method of ﬁnding the common
points of convex sets and its application to the solution of problems
in convex programming. Comput. Math. Math. Phys. 7 (3), 620–
631.
Cemgil, A.T., 2008. Bayesian Inference in Non-Negative Matrix Factorization Models. Technical Report CUED/F-INFENG/TR.609, University of Cambridge.
Cemgil, A.T., Godsill, S.J., 2005. Probabilistic phase vocoder and its
application to interpolation of missing values in audio signals. In:
Proc. Eur. Signal Process. Conf. (EUSIPCO).
Cherry, E.C., 1953. Some experiments on the recognition of speech, with
one and two ears. J. Acoust. Soc. Amer. 25 (5), 975–979.

675

de Cheveigné, A., Kawahara, H., 2002. YIN, a fundamental frequency
estimator for speech and music. J. Acoust. Soc. Amer. 111 (4), 1917–
1930.
Clark, P., Atlas, L., 2008. Modulation decompositions for the interpolation of long gaps in acoustic signals. In: Proc. Internat. Conf. on
Acoust. Speech Signal Process. (ICASSP). pp. 3741–3744.
Cooke, M.P., Green, P., Josifovski, L., Vizinho, A., 2001. Robust
automatic speech recognition with missing and unreliable acoustic
data. Speech Comm. 34 (3), 267–285.
Criminisi, A., Pérez, P., Toyama, K., 2004. Region ﬁlling and object
removal by exemplar-based image inpainting. IEEE Trans. Image
Process. 13 (9), 1200–1212.
Csiszár, I., 1975. I-divergence geometry of probability distributions andminimization problems. Ann. Probab. 3 (1), 146–158.
Eggert, J., Körner, E., 2004. Sparse coding and NMF. In: Proc. IEEE
Internat. Joint Conf. on Neural Networks, vol. 4. pp. 2529–2533.
Ellis, D.P.W., 1993. Hierarchic models of sound for separation andrestoration. In: Proc. IEEE Workshop Appl. Signal Process. Audio
Acoust. (WASPAA).
Ellis, D.P.W., 1996. Prediction-Driven Computational Auditory Scene
Analysis. Ph.D. Thesis, Massachusetts Institute of Technology.
Esquef, P.A.A., Biscainho, L.W.P., 2006. An eﬃcient model-based
multirate method for reconstruction of audio signals across long gaps.
IEEE Trans. Audio Speech Lang. Process. 14 (4), 1391–1400.
Févotte, C., Bertin, N., Durrieu, J.-L., 2009. Nonnegative matrix
factorization with the Itakura–Saito divergence with application tomusic analysis. Neural Comput. 21 (3), 793–830.
Fujisaki, H., Nagashima, S., 1969. A Model for Synthesis of Pitch
Contours of Connected Speech. Annual Report of the Engineering
Research Institute, Faculty of Engineering, University of Tokyo 28.
pp. 53–60.
Godsill, S.J., Rayner, P.J.W., 1998. Digital Audio Restoration: A
Statistical Model Based Approach. Springer-Verlag New York, Inc.,
Secaucus, NJ, USA.
Griﬃn, D.W., Lim, J.S., 1984. Signal estimation from modiﬁed short-time
Fourier transform. IEEE Trans. Acoust. Speech Signal Process. 32 (2),236–243.
Grünwald, P.D., 2007. The Minimum Description Length Principle. The
MIT Press.
von Helmholtz, H., 1954. On the Sensations of Tone as a Physiological
Basis for the Theory of Music (second English ed. A.J. Ellis, Trans.
1885, ﬁrst German ed. 1863). Reprinted by Dover Publications.
Janssen, A.J.E.M., Veldhuis, R., Vries, L.B., 1986. Adaptive interpolationof discrete-time signals that can be modeled as AR processes. IEEE
Trans. Acoust. Speech Signal Process. 34 (2), 317–330.
Kameoka, H., 2007. Statistical Approach to Multipitch Analysis. Ph.D.
Thesis, The University of Tokyo.
Kameoka, H., Nishimoto, T., Sagayama, S., 2007. A multipitch analyzer
based on harmonic temporal structured clustering. IEEE Trans. Audio
Speech Lang. Process. 15 (3), 982–994.
Kameoka, H., Ono, N., Kashino, K., Sagayama, S., 2009. Complex
NMF: a new sparse representation for acoustic signals. In: Proc.
Internat. Conf. on Acoust. Speech Signal Process. (ICASSP). pp.
3437–3440.
Kashino, M., 2006. Phonemic restoration: the brain creates missing speechsounds. Acoust. Sci. Technol. 27 (6), 318–321.
Le Roux, J., Kameoka, H., Ono, N., de Cheveigné, A., Sagayama, S.,
2007a. Single and multiple F 0 contour estimation through parametric
spectrogram modeling of speech in noisy environments. IEEE Trans.
Audio Speech Lang. Process. 15 (4), 1135–1145.
Le Roux, J., Kameoka, H., Ono, N., de Cheveigné, A., Sagayama, S.,
2007b. Monaural speech separation through harmonic-temporalclustering of the power spectrum. In: Proc. Acoust. Soc. Jpn. Autumn
Meeting. No. 3-4-3. pp. 351–352.
Le Roux, J., Kameoka, H., Ono, N., Sagayama, S., 2008a. On theInterpretation of I-Divergence-Based Distribution-Fitting as a Maximum-Likelihood Estimation Problem. Technical Report METR 200811, The University of Tokyo.

676

J. Le Roux et al. / Speech Communication 53 (2011) 658–676

Le Roux, J., Ono, N., Sagayama, S., 2008b. Explicit consistency
constraints for STFT spectrograms and their application to phasereconstruction. In: Proc. ISCA Workshop Statist. Perceptual Audition(SAPA). pp. 23–28.
Lee, D.D., Seung, H.S., 1999. Learning of the parts of objects by nonnegative matrix factorization. Nature 401, 788–791.
Lee, D.D., Seung, H.S., 2001. Algorithms for non-negative matrix
factorization. Advances in Neural Information Processing Systems
13 (Proc. NIPS*2000). The MIT Press, Cambridge, MA, pp. 556–562.
Lu, L., Mao, Y., Liu, W., Zhang, H.-J., 2003. Audio restoration by
constrained audio texture synthesis. In: Proc. Internat. Conf. on
Acoust. Speech Signal Process. (ICASSP), vol. 5. pp. 636–639.
Maher, R.C., 1993. A method for extrapolation of missing digital audio
data. In: 95th AES Convention, New York. pp. 1–19.
Maher, R.C., 1994. A method for extrapolation of missing digital audio
data. J. Audio Eng. Soc. 42 (5), 350–357.
McAulay, R.J., Quatieri, T.F., 1986. Speech analysis/synthesis based on asinusoidal representation. IEEE Trans. Acoust. Speech Signal Process.
34 (4), 744–754.
Meng, X.-L., Rubin, D.B., 1993. Maximum likelihood estimation via the
ECM algorithm: a general framework. Biometrika 80 (2), 267–278.
Mørup, M., Schmidt, M.N., 2006. Sparse Non-Negative Matrix Factor 2-D
Deconvolution. Technical Report, Technical University of Denmark.
Raj, B., Seltzer, M.L., Stern, R.M., 2004. Reconstruction of missing
features for robust speech recognition. Speech Comm. 43 (4), 275–296.
Raj, B., Stern, R.M., 2005. Missing-feature approaches in speech
recognition. IEEE Signal Process. Mag. 22 (5), 101–116.
Rajan, J.J., Rayner, P.J.W., Godsill, S.J., 1997. A Bayesian approach toparameter estimation and interpolation of time-varying autoregressive
processes using the Gibbs sampler. IEE Proc. Vision Image Signal
Process. 144 (4), 249–256.
Rayner, P.J.W., Godsill, S.J., 1991. The detection and correction ofartefacts in archived grammophone recordings. In: Proc. IEEE
Workshop Appl. Signal Process. Audio Acoust. (WASPAA).
Reyes-Gomez, M., Jojic, N., Ellis, D.P.W., 2004. Towards single-channel
unsupervised source separation of speech mixtures: the layered

harmonics/formants separation-tracking model. In: Proc. ISCA
Workshop Statist. Perceptual Audition (SAPA). pp. 25–30.
Sajda, P., Du, S., Parra, L.C., 2003. Recovery of constituent spectra using
non-negative matrix factorization. In: Proc. SPIE Wavelets X. pp.
321–331.
Schmidt, M.N., Mørup, M., 2006. Nonnegative matrix factor 2-D
deconvolution for blind single channel source separation. In: Proc.
Internat. Conf. on Independent Compon. Anal. Blind Source Sep.
(ICA 2006). pp. 700–707.
Smaragdis, P., 2004. Non-negative matrix factor deconvolution; extraction of multiple sound sources from monophonic inputs. In: Proc.
Internat. Conf. on Independent Compon. Anal. Blind Signal Sep.
(ICA 2004). pp. 494–499.
Smaragdis, P., Raj, B., Shashanka, M., 2009. Missing data imputation forspectral audio signals. In: Proc. IEEE Workshop Machine Learn.
Signal Process. (MLSP).
Vaseghi, S.V., Rayner, P.J.W., 1990. Detection and suppression of
impulsive noise in speech communication systems. Comm. Speech
Vision IEE Proc. I 137 (1), 38–46.
Veldhuis, R., 1990. Restoration of Lost Samples in Digital Audio Signals.
Prentice-Hall, Englewood Cliﬀs, NJ.
Virtanen, T., Mesaros, A., Ryynänen, M., 2008. Combining pitch-based
inference and non-negative spectrogram factorization in separating
vocals from polyphonic music. In: Proc. ISCA Workshop Statist.
Perceptual Audition (SAPA). pp. 17–22.
Wang, D.-L., Brown, G.J. (Eds.), 2006. Computational Auditory
Scene Analysis: Principles, Algorithms, and Applications. IEEE
Press/Wiley.
Warren, R.M., 1970. Perceptual restoration of missing speech sounds.
Science 167 (3917), 392–393.
Warren, R.M., 1982. Auditory Perception: A New Synthesis. Pergamon,
New York, NY.
Wolfe, P.J., Godsill, S.J., 2005. Interpolation of missing data values foraudio signal restoration using a Gabor regression model. In: Proc.
Internat. Conf. on Acoust. Speech Signal Process. (ICASSP), vol. V.
pp. 517–520.