IEEE TRANSACTIONS ON AUDIO, SPEECH, AND LANGUAGE PROCESSING, VOL. 15, NO. 4, MAY 2007

1135

Single and Multiple F0 Contour Estimation ThroughParametric Spectrogram Modeling of Speechin Noisy Environments
Jonathan Le Roux, Student Member, IEEE, Hirokazu Kameoka, Student Member, IEEE,
Nobutaka Ono, Member, IEEE, Alain de Cheveigné, Member, IEEE, and Shigeki Sagayama, Member, IEEE

Abstract—This paper proposes a novel 0 contour estimation algorithm based on a precise parametric description of the voiced
parts of speech derived from the power spectrum. The algorithm
is able to perform in a wide variety of noisy environments as wellas to estimate the 0 s of cochannel concurrent speech. The speech
spectrum is modeled as a sequence of spectral clusters governed by
a common 0 contour expressed as a spline curve. These clusters
are obtained by an unsupervised 2-D time-frequency clustering of
the power density using a new formulation of the EM algorithm,
and their common 0 contour is estimated at the same time. A
smooth 0 contour is extracted for the whole utterance, linking
together its voiced parts. A noise model is used to cope with nonharmonic background noise, which would otherwise interfere withthe clustering of the harmonic portions of speech. We evaluate ouralgorithm in comparison with existing methods on several tasks,
and show 1) that it is competitive on clean single-speaker speech,2) that it outperforms existing methods in the presence of noise,and 3) that it outperforms existing methods for the estimation ofmultiple 0 contours of cochannel concurrent speech.
Index Terms—Acoustic scene analysis, expectation-maximization (EM) algorithm, harmonic-temporal structured clustering(HTC), multipitch estimation, noisy speech, spline 0 contour.

I. INTRODUCTIONT

HE analysis of complex and varied acoustic scenes is afundamental and challenging problem for today’s acoustic
signal processing. At its core is the need for a method which cancontour of harmonic
determine precisely and robustly the
analysis has a very broad
signals such as speech. Robust
range of applications in computational auditory scene analysis(CASA), speech recognition, prosody analysis, speech enhancement, or speaker identification. So far, many pitch determination
algorithms (PDAs) have been proposed [1], some of them withvery good accuracy [2]. However, they are usually limited to theclean speech of a single speaker and fail in moderate amounts ofbackground noise or the presence of other speakers. Ideally, theManuscript received July 24, 2006; revised November 22, 2006. The associate
editor coordinating the review of this manuscript and approving it for publication was Dr. Simon King.
J. Le Roux, H. Kameoka, N. Ono, and S. Sagayama are with the Graduate School of Information Science and Technology, University of Tokyo,
Tokyo 113–8656, Japan (e-mail: leroux@hil.t.u-tokyo.ac.jp; kameoka@hil.t.utokyo.ac.jp; onono@hil.t.u-tokyo.ac.jp; sagayama@hil.t.u-tokyo.ac.jp).
A. de Cheveigné is with the Centre National de la Recherche Scientifique,
Université Paris V, and Ecole Normale Supérieure, 75230 Paris Cedex 05,
France (e-mail: Alain.de.Cheveigne@ens.fr).
Color versions of one or more of the figures in this paper are available online
at http://ieeexplore.ieee.org.
Digital Object Identifier 10.1109/TASL.2007.894510

performance of a PDA should stay high in as wide a range ofbackground noises as possible (white noise, pink noise, noise
bursts, music, other speech, etc.). Furthermore, the possibility
contours of several concurrent
to extract simultaneously thevoices is also a desirable feature. Several PDAs already exist
that deal with the tracking of multiple s [3]–[7]. Several ofthese algorithms rely on an initial frame-by-frame analysis followed by post-processing to reduce errors and obtain a smooth
contour, for example using hidden Markov models (HMMs)(see [3] for a review). Here, we propose to perform estimationand model-based interpolation simultaneously, based on a parametric model of the time and frequency shape of the spectralenvelope of the voiced parts of speech.
To this end, we rely on a parametric description of the wavelet
power spectrum that accounts for its structure simultaneously
in time and frequency directions. In brief, we model the power
as a sum ofparametric source models
spectrum, where is log-frequency, is time, andis the. In thisset of model parameters:
paper, this model is applied to the voiced parts of speech, but
in other contexts it could be applied to other acoustic events,
provided a mathematical description of those events. Here, weapply two constraints to the spectro-temporal model. Alongthe frequency axis, it consists of a series of harmonics withfrequencies multiple of a common . Along the temporal axis,
follows a contour modeled as a cubic spline, and thethisamplitude of each partial follows a smooth temporal envelopemodeled as a sum of Gaussian functions. The cubic spline
contour was preferred to other options such as the Fujisaki
model [8] as it can be applied to a wider range of stimuli
in addition to speech. Moreover, as we will explain later, thischoice enabled us to obtain analytic update equations during the
optimization process. To handle the slowly-varying harmonic
structure of the short-term wavelet spectrum, we used a multipitch analysis method initially developed for feature extractionof music signals, the harmonic-temporal structured clustering(HTC) method [9], which gives a parametric representation of
the harmonic parts of the power spectrum.
In the original formulation of HTC, each source model represents a sound stream with constant . However, by considering speech as a succession of models that each corresponds
to a phoneme (or more generally to a segment of the speechutterance with steady acoustic characteristics), we can use the
HTC method to model the spectrum as a sequence of speccontour. Contrary totral cluster models with a continuous

1558-7916/$25.00 © 2007 IEEE

1136

IEEE TRANSACTIONS ON AUDIO, SPEECH, AND LANGUAGE PROCESSING, VOL. 15, NO. 4, MAY 2007

HMMs, which assume discrete phonetic states, we aim here tomodel smooth transitions in the temporal succession of the spectral structures.
Finally, we introduce also a noise model to deal with thenonharmonic power coming from background noise. While thespectrogram of voiced speech is characterized by harmonic parts
with strong relative power, noise tends to have a more flat spectrogram. Extracting voiced speech from such noise corresponds
to searching for local, harmonically structured “islands” withina “sea” of unstructured noise.
The parametric model we obtain in this way is optimized
using a new formulation of the EM algorithm. The spectral clusters are obtained by an unsupervised 2-D clustering of the power
density, performed simultaneously with the estimation of thecontour of the whole utterance. We describe our method more
precisely in the following section.
II. FORMULATION OF THE MODEL
A. General HTC Method
Consider the wavelet power spectrumof a signal
recorded from an acoustic scene, defined on a domain of def.
initionThe problem considered is to approximate the power spectrum
parametric source models
as well as possible as the sum ofmodeling the power spectrum of“objects” eachcontourand its own harmonic-temporalwith its own
structure. We note that the formulation described hereafter is not
limited to any particular time-frequency representation, and it
could work with a standard linear-frequency short-time Fouriertransform (STFT). However, we found better performance withwavelets, possibly because it offers relatively better spectral resolution for the low-frequency harmonics of the voice. We willthus use the wavelet power spectrum in the following.
are exAs described in [9], the source models
pressed as a Gaussian mixture model (GMM) with constraints
on the characteristics of the kernel distributions: supposing thatpartials modeled in the frequency
there is harmonicity withdirection, and that the power envelope is described using
kernel functions in the time direction, we can rewrite eachsource model in the form
(1)
where

Fig. 1. Graphical representation of a HTC source model. (a) shows the timefrequency profile of the model, while (b) shows a cross section of the model atconstant time, and (c) the evolution in time of the power envelope function. The
harmonic structure of the model can be seen in (b), and the approximation of
the power envelope in the time direction as a sum of Gaussian kernels can beseen in (c).

We use the -divergence [10] as a classical way to measure the“distance” between two distributions(3)and we are thus looking for.
Keeping only the terms depending on and reversing the signof this expression, one defines the following function to maximize with respect to :is the set of all parameters and with kernel densities
which are assumed to have the following shape:
Using this function
rameter :(4), one can derive the likelihood of the pa-(2)
where the parameters ,andare normalized to unity.
A graphical representation of a HTC source model
can be seen in Fig. 1.
and
Our goal is to minimize the difference between
according to a certain criterion.

(5)is the Gamma function, and the second part of thewhereexponent ensures that we obtain a probability measure. One canindeed see this probability as the joint probability of all the variables
independently following Poisson distributions of

LE ROUX et al.: SINGLE AND MULTIPLE

CONTOUR ESTIMATION

1137

parameter. This way of presenting the problem enables
us to interpret it as a maximum a posteriori (MAP) estimation
problem and to introduce prior functions on the parameters as
follows, using Bayes theorem:(6)
Our goal is now equivalent to the maximization with respect
to of. In the following, we will writefor. The problem is that in the term
simply, there is a sum inside the logarithm, and we thus cannot obtain an analytical solution. However, if we introduce nonnegative membership desumming to 1 for each, one can write,
grees
using the concavity of the logarithm:

(7)
where
denotes the convex combination with coefficients. Moreover, the inequality (7) becomes an equality for

We can thus use an EM-like algorithm to maximize the likelihood by alternately updating and the membership degrees ,which act as auxiliary parameters, while keeping the other fixed:
E-stepM-step
with

(8)
where
have from (7) that

. For all, we indeed(9)Fig. 2. Optimization through the EM algorithm. During the E-step, the auxil; m . Then, during theiary parameter m is updated to m so that
M-step,; m is optimized w.r.t.
, ensuring that;m >. The maximization of
can thus be performed through;m
the maximization of the auxiliary function; m alternately w.r.t. m and.

J (22 ^ )J (22 ^ ) = J (22)^2

J (22) = J (22 ^ )J (2^ )  J (2^ ^ )J (22)
2
J (22 )

can be used as an auxiliary function to maximize,andenabling us to obtain analytical update equations. The optimization process is illustrated in Fig. 2.
A slightly different approach was presented in [9], but leads
to the same optimization process as the one we present here.
The E-step is straightforward and is dealt with in exactly thesame way as in [9]. However, the possibility to obtain analytical
update equations during the M-step depends on the actual ex, and it could be performed there in the specialpression ofcontour. More gencase of a piece-wise flat
erally, if the other HTC parameters ( , ,,, , ), then the update equado not enter in the expression oftions obtained in [9] for these parameters can be used as is, and
we only need to obtain the M-step update equations for thecontour parameters. We will explain how we proceed in the following subsection.
B. Speech ModelingIn the following, in order to model the spectrum of a speechutterance, we will make several assumptions. First, we assume
contour is smooth and defined on the whole interval:that thevalues are
we will not make voiced/unvoiced decisions, andassumed continuous. Second, we fix the harmonic structure ofeach HTC source model so that it corresponds to segments ofspeech with steady acoustic characteristics, and assume that aspeech segment is a succession of such steady segments sharing
contour.
a common1) Spline
Contour: In our previous work [9], the HTC
contours,
method has only been applied to piece-wise flat
which is relevant for certain instruments like the piano for example, but of course not in speech. Looking for a smooth
contour, we chose to use cubic spline functions as a general classof smooth functions, so as to be able to deal in the future with awide variety of acoustic phenomena (background music, phoneringing, etc). Moreover, their simple algebraic formulation enables us to optimize the parameters through the EM algorithm,
as update equations can be obtained analytically. It may happen
assumption is invalid, such as at the end ofthat the smoothutterances where -halving can occur, but this problem is facedby all algorithms that exploit continuity of , and the assumption is justified empirically in that including it tends to reduceoverall error rates. Moreover, failure to track -halving at theend of utterances is perhaps not too serious, as the halving is1138

IEEE TRANSACTIONS ON AUDIO, SPEECH, AND LANGUAGE PROCESSING, VOL. 15, NO. 4, MAY 2007

usually not salient perceptively, other than as a roughness of indeterminate pitch. Furthermore, it is one of a wider class of iris
regular voicing phenomena (diplophony, creak) for which
hard to define [11].
The analysis interval is divided into subintervals
which are assumed of equal length. Following [12], the parameters of the spline contour model are then the values of theat each bounding point . The values of the second derivative
for a cerat those points are given by the expression
which can be explicitly computed offline, undertain matrix
the hypothesis that the first-order derivative is 0 at the boundsof the analysis interval. We can assume so if we set the boundsof the interval outside the region where there is speech. One canthen classically obtain a simple algebraic formula for the conon the whole interval. For
,
tour

(10)
2) Optimization of the Model: To design our model, we further make the following assumptions. For simplicity, we firstdescribe here the case of a single speaker, and the multiple
speakers case will be presented in Section II-B.4.
We make all the source models inside the HTC method sharecontour:, , by plugging the anathe same
contour into (2), suchlytical expression (10) of the spline
expresthat all the source models are driven by the same
sion. Our intention is to have a succession in time of slightly
overlapping source models which correspond if possible to successive phonemes, or at least to segments of the speech utterance with steady acoustic characteristics. As the structure is assumed harmonic, the model takes advantage of the voiced partsof the speech utterance, which it uses as anchors. When usedon natural speech, if the unvoiced/silent parts are too long, it
may happen that the spline contour becomes unstable, which
contour extraction immecan deteriorate the accuracy of thediately before or after a section of unvoiced speech, especially
if the neighboring voiced parts are not strongly enough voiced.
If they are not, we believe there is no particular incidence on thecontour near unvoiced parts of the speech.
accuracy of the
The results of the experimental evaluations will show that thisassumption is justified.
We also assume that inside a source model the same powerenvelope is used for all harmonics, as we want to isolate structures in the speech flow with stable acoustic characteristics. The
model allows source models to overlap, so a given spectral shape
can merge progressively into another, which allows it to fit arbitrary spectro-temporal shapes. The subscript can thus be ex. We note however that the following discussion
cluded inon the optimization of the model is independent of this particularassumption and that the algorithm is general enough to allow
separate power envelope functions.
The optimization process goes as follows: we start by updating the HTC parameters which do not enter in the spline,, , ) through analytical upmodel (namely , ,date equations which are the same as in [9]. Once these parame-ters have been updated, we compute the derivatives ofwith respect to the spline parameters . A global maximum isdifficult to obtain, but we can update the one after the other,starting for example from and using the already updated parameters for the update of the next one. This way of performing
the updates is referred to as the coordinate descent method [13],and can be summarized as

..
.
(11)
where
denotes the set of whole parameters except the .
The corresponding optimization procedure, called the expectation-constrained maximization algorithm (ECM) [14], does not
ensure the maximization in the M-step but guarantees the in. Putting the derivatives with recrease of the function
spect to to 0, one then finds update equations analytically atstep :

(12)
where
and
does not depend on andonly depends on andthe fixed matrix .
The partial derivatives with respect to the other parameters,, , ) are the same as in [9], as mentioned( , ,above.
3) Prior Distribution: As seen in Section II-A, the optimization of our model can be naturally extended to a MAP estimaon the parameters,
tion by introducing prior distributions
which work as penalty functions that try to keep the parameters
within a specified range. The parameters which are the best compromise with empirical constraints are then obtained throughequation (6).
, it becomes
By introducing such a prior distribution on
possible to prevent half-pitch errors, as the resulting source
model would usually have a harmonic structure with zero power
for all the odd order harmonics, which is abnormal for speech.
We apply the Dirichlet distribution, which is explicitly given by
(13)
where

is the most preferred “expected” value ofsuch that, the contribution degree of the prior andtheis taken whenGamma function. The maximum value forfor all . Whenis zero,become uniformLE ROUX et al.: SINGLE AND MULTIPLE

CONTOUR ESTIMATION

1139

Fig. 3. Comparison of observed and modeled spectra (“Tsuuyaku denwa kokusai kaigi jimukyoku desu,” female speaker). The estimated F contour is reproduced on both the observed and modeled spectrograms to show the precision of the algorithm. (a) Observed spectrogram and estimated F contour. (b) Modeled
spectrogram and estimated F contour.

distributions. The choice of this particular distribution allows
us to give an analytical form of the update equations of.
Although the spline model can be used as is, one can also
introduce in the same way a prior distribution on the parameof the spline
contour, in order to avoid an overfittersting problem with the spline function. Indeed, spline functions
have a tendency to take large variations, which is not naturalcontour of a speech utterance. Moreover, theconfor thetour might also be hard to obtain on voiced parts with relatively
lower power or poor harmonicity. The neighboring voiced portions with higher power help the estimation over these intervals
by providing a good prior distribution.
To build this prior distribution, we assume that the form aMarkov chain, such thatand assume furthermore that
follows a uniform distribution,
follows a Gaussian distriand that, conditionally tobution of centerand variancecorresponding to theweighting parameter of the prior distribution:

In the derivative with respect to
up two new terms

used above to obtain (12) add

and the update equation (12) then becomes

(14)

and
are, respectively, the numerator and dewhere
nominator of the right term of equation (12). The update equation for the boundary points is derived similarly.
An example is presented in Fig. 3, based on the Japanese sentence “Tsuuyaku denwa kokusai kaigi jimukyoku desu” utteredby a female speaker. Shown are 2-D representations of the observed and modeled spectra (after 30 iterations of the estimation
contour estimated through our method is realgorithm). The
produced on both the observed and modeled spectrograms toshow the precision of our algorithm. One can see that the model
contour is acapproximates well the spectrum and that the
curately estimated.
Estimation: The multiple-speakers case is a4) Multiple
simple extension of the single-speaker one. While for a single
estimation all the source models share the same
contour,
for multiple speakers, according to the number ofcontoursthat we want to estimate, we group together source models into
subsets such that source models inside a subset share a common
contour. For example, if we use
models in total fortwo speakers, the models with index
will be attributed to the first speaker, while the others will be attributed
to the second one. We thus have pools of source models driven
contour for each of the pools and corresponding
by a single
to one of the speakers, and we only need to introduce a set ofcontours. These sets can be
spline parameters for each of the
optimized independently and simultaneously in the exact same
way as in the single-speaker case: the E-step is unchanged, andthe M-step is performed by first updating the HTC parameterswhich do not enter in the spline model, and then updating thecontour through the ECM algospline parameters of eachrithm. This last update can be done independently as the derivawith regards to the parameters of one of thetives ofcontours do not include any term depending on the parameters of the other contours.
The method handles overlapping harmonics by, at each iteration of the algorithm, re-estimating simultaneously the contributions of each voice in the spectrum. It relies on the assumption that the spectra of the contributing sources combine additively, and neglects phase-dependent vector summation. This is1140

IEEE TRANSACTIONS ON AUDIO, SPEECH, AND LANGUAGE PROCESSING, VOL. 15, NO. 4, MAY 2007of course a rough approximation in the view of perfect source
separation. However, the speech spectrum is usually relatively
estimation can rely mostly on the composparse, and thenents that are isolated from each other. This is the reason why
we think this rough approximation is not a serious problem, as
estimation is our main objective in this paper.

and inIf we noteandas in Section II-A, thetroduce
EM algorithm can be applied in the same way as describedabove, just differing in the range of the summations on , , and. We only need to specify the update equations of the M-stepfor the noise model parameters

C. Noise ModelingWe introduce a noise model to cope with the backgroundnoise that can be a disturbance in the process of clustering
the harmonic portions of speech. Indeed, it would be more
rewarding, in the purpose of decreasing the -divergence,for the harmonic source models to take very large variancesin the log-frequency direction and have the centers of theGaussian distribution go on portions of the spectrogram withstrong broad power, even though there is no harmonic structure
corresponding to these portions, especially if the noise power
is comparable to or even larger than the speech signal power.
The spectrogram of a quasi-periodic signal such as voiced
speech consists of a large number of line spectrum components
and has spikes that are strongly marked, while the spectrogramof a white or pink noise has a tendency to be more flat, without
significant spiky regions. The idea to design the noise model
was thus that detecting the harmonic parts of the spectrogramin a noisy background corresponds to searching for thin andharmonically distributed “islands” which rise out of a “sea” ofnoise. We thus chose to model the noise using a mixture ofGaussian distributions with large fixed variance and with centers fixed on a grid, the only parameters of the model being therespective weights of the Gaussians in the model and the ratio
of noise power inside the whole spectrogram. This noise-cancelling approach can be considered quite close to spectral subtraction in the sense that the power spectrum is in both cases assumed additive. However, while spectral subtraction generally
needs voiced/unvoiced decision to obtain the power spectrum of
the background noise, which furthermore has to be assumed stationary, our approach estimates adaptively the noise part of thespectrum even when there is speech, taking advantage of the valleys of the spectral structure. The only assumption we make isthat the noise spectrum is smooth in both time and frequency direction, whereas speech spectrum is more spiky in the frequency
direction. Therefore, we can expect our model’s performance to
be close to the best performance that spectral subtraction could
reach.
be the number of columns of the grid (in the timeLet
the number of rows (in the log-frequency didirection) andrection). The noise model is defined as(15)
where is the proportion of the noise model in the total model,
are the weights of the Gaussian functions in the mixture
thethe variances of the Gaussian
and add up to 1,
the log-frequency index of the centers of the thfunctions,row, andthe time index of the centers of the th column.

(16)
(17)

where the superscript
refers to the iteration cycle. The update
equations for the other parameters remain the same, the onlychanges coming from the E-step where the noise model now
enters in the summation.
The noise cancelling performed through the introduction of
the noise model is illustrated in Fig. 4: Fig. 4(a) gives a 3-D viewof the original clean spectrogram of a part of the sentence “It was
important to be perfect since there were no prompts” uttered bya female speaker, and Fig. 4(b) shows the spectrogram of thesame part of the utterance to which white noise at a signal-tonoise ratio (SNR) of 2 dB has been added. The estimation of
the spectrogram where the noise has been cancelled, shown in
Fig. 4(c), is obtained through the optimized masking functions
as in the following formula:(18)
In the course of the optimization of the model, the contour estimation is performed on this “cleaned” part of the spectrogram,which enables our estimation algorithm to perform well evenin very noisy environments, as we will show in the next section.
D. Parametric Representation and Potential Applications
We would like to stress the fact that our algorithm not only
contour, but also gives a parametric repreestimates thesentation of the voiced parts of the spectrogram. This can be
useful especially in the analysis of cochannel speech by multiple speakers, as one can get a parametric representation of theharmonic parts of the separated spectrograms of each utterance,
as shown in Fig. 5: Fig. 5(b) represents the modeled spectrogram of the Japanese utterance “oi wo ou” by a male speaker,
and Fig. 5(c) the one of the utterance “aoi” by a female speaker,
conextracted from the mixed signal shown in Fig. 5(a) (thetour near the boundary is not relevant as there is no sound by thesecond speaker there). These parametric representations could
be used for example to cluster the spectrogram of the mixed
sound and separate the speakers, as well as for noise cancelling,
as we showed in Fig. 4. The reader might observe that the harmonic tracks at the end of the utterance in Fig. 5(b) possess
greater energy than their original counterparts. This might berelated to the fact that the power envelope functions of all theharmonics are linked together. We shall note however that thisLE ROUX et al.: SINGLE AND MULTIPLE

CONTOUR ESTIMATION

Fig. 4. Estimation of the clean part of a noisy spectrogram. (a) shows a 3-Dview of the original clean spectrogram of a part of the sentence “It was important
to be perfect since there were no prompts” uttered by a female speaker. (b) shows
the spectrogram of the same part of the utterance to which white noise at an SNRof 2 dB has been added. (c) shows the estimated noise-cancelled part of thespectrogram of (b).

0

is not an important issue for the applications that we foresee,such as speech enhancement or clustering, because we plan to
use the parametric models only to obtain the proportions of eachaudio object inside the original spectrogram, and then multiply
them by the original spectrogram to extract the desired object,

1141

Fig. 5. Parametric representation of separated spectrograms. (a) shows the
spectrogram of a signal obtained by mixing the two Japanese utterances“oi wo ou” by a male speaker and “aoi” by a female speaker, together with
the F contours estimated by our method. (b) and (c) show the parametricrepresentations of the spectrograms of the utterances by the male and female
speaker, respectively, extracted from the mixed signal shown in (a).

as in (18). If the original spectrogram had low energy at a certain time-frequency point, the obtained spectrogram will have
low energy as well.

1142

IEEE TRANSACTIONS ON AUDIO, SPEECH, AND LANGUAGE PROCESSING, VOL. 15, NO. 4, MAY 2007

III. EXPERIMENTAL EVALUATION
A. Single-Speaker

Estimation in Clean Environment

We evaluated the accuracy of thecontour estimation of
our model on a database of speech recorded together with alaryngograph signal [15], consisting of one male and one female speaker who each spoke 50 English sentences for a totalof 0.12 h of speech, for the purpose of evaluation of -estimation algorithms.
was calculated from an input
The power spectrum
signal digitized at a 16-kHz sampling rate (the original data ofthe database was converted from 20 kHz to 16 kHz) using a
Gabor wavelet transform with a time resolution of 16 ms forthe lowest frequency subband. Higher subbands were downsampled to match the lowest subband resolution. The lower boundof the frequency range and the frequency resolution were, respectively, 50 Hz and 14 cent. The spline contour was initiallyflat and set to 132 Hz for the male speaker and 296 Hz for thefemale speaker. These values were tuned in frequency bins aftera few preliminary experiments. We shall note that, as we wereable to obtain very good performance with the same initial conditions on various data, we believe that the performance is notcontour parameters. Thevery sensitive to the priming of thelength of the interpolation intervals was fixed to four frames.
source models, each of them withFor HTC, we usedharmonics. This is enough forestimation. For abetter modeling of the spectrogram, one can use 40 or 60 harmonics for example. Temporal envelope functions were modGaussian kernels. The
were set to,eled using
to
, to
,
to 32 ms, and to 422
cents. For the prior functions,
was fixed to 0.4,
to 0.04,. The algorithm
andwas run on the utterances from which the initial and final silence parts were manually removed. We note that the minimum
time window on which the algorithm can work is one frame, inwhich case the algorithm works frame-by-frame. This gives acceptable results as described in [16], but perfomance improves
with longer windows, as documented in [9]. It is thus used hereon the whole time interval. The computational cost is also discussed in [9].
estimates and the reliability
We used as ground truth themask derived by de Cheveigné et al. [2] under the following
estimate was obvicriteria: 1) any estimate for which theously incorrect was excluded; and (2) any remaining estimate
for which there was evidence of vocal fold vibration was included. Frames outside the reliability mask were not taken into
account during our computation of the accuracy, although ouralgorithm gives values for every point of the analysis interval byconstruction. As the spline function gives an analytical expression for the contour, we compare our result with 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 were deemed to
be gross errors. The results can be seen in Table I, with forcomparison the results obtained by de Cheveigné et al. [2] forseveral other algorithms. Notations stand for the method used,
as follows. ac: Boersma’s autocorrelation method [17], [18], cc:
cross-correlation [18], shs: spectral subharmonic summation

TABLE IGROSS ERROR RATES FOR SEVERAL F ESTIMATION ALGORITHMS
ON CLEAN SINGLE-SPEAKER SPEECH[18], [19], pda: eSRPD algorithm [15], [20], fxac: autocorrelation function (ACF) of the cubed waveform [21], fxcep:
cepstrum [21], additive: probabilistic spectrum-based method[22], acf: ACF [2], nacf: normalized ACF [2], TEMPO: the
TEMPO algorithm [23], YIN: the YIN algorithm [2]. Moredetails concerning these algorithms can be found in [2]. We cansee that our model’s accuracy for clean speech is comparable
extraction algorithms
to the best existing single-speakerdesigned for that purpose.
B. Single
Pink Noise

Estimation on Speech Mixed With White andWe performed
estimation experiments on speech to which
a white noise, band-passed between 50 and 3300 Hz, was added,
with SNRs of 0, 2, and 10 dB. These SNRs were selected
because 0 dB corresponds to equal power, 2 dB is used in astudy [24] that used the same database as we use in III-C, and10 dB is a relatively noisy condition to illustrate the effectiveness of our algorithm in that case. The database mentioned
above [15] was again used, and the white noise added was generated independently for each utterance. We also performed experiments with pink noise, band-passed between 50 and 3300 Hz,
with an SNR of 2 dB. The spectrum of pink noise is closer
to that of speech than white noise. The noise model was initialand theall equal to, while theized withcent and,variances were fixed toof the Gaussian distributions of theand the centersnoise model were fixed on a grid such that the distances between
them in the time and log-frequency directions were all equal toand, respectively, to ensure a good overlap between theGaussian distributions. The determination of the variances wasmade after a few experiments while keeping in mind thatshould be significantly larger than the typical variance in the
log-frequency direction of the Gaussian of the harmonic model
but small enough to still be able to model fluctuations in thenoise power.
As a comparison, we present results obtained on the same
database using YIN [2] and the algorithm of Wu, Wang, andtracking in a noisy enBrown [4], specifically designed forvironment, and that can also handle the estimation of two simultaneous s. Their code is made available on the Internetby their respective authors. The algorithm of Wu, Wang, andBrown will be referred to as the WWB algorithm. According toLE ROUX et al.: SINGLE AND MULTIPLE

CONTOUR ESTIMATION

1143

TABLE II
ACCURACY (%) OF THE F ESTIMATION OF SINGLE-SPEAKER SPEECH MIXED WITH WHITE AND PINK NOISESTABLE III
CATEGORIZATION OF INTERFERENCE SIGNALSTABLE IV
ACCURACY (%) OF THE F ESTIMATION OF VOICED SPEECH
WITH SEVERAL KINDS OF INTERFERENCESTABLE V

its authors, the parameters of this algorithm could be tuned on a
new database to obtain the best performances, but they mention[4] that it is supposed to work fine in the version made available(trained on a corpus [24] that we will use later).
We obtained good results, presented in Table II, showing therobustness of our method on noisy speech, when noise is not harmonic. Note that the level of the noise added is greater than thatof the original signal, and that at 10 dB it is even difficult forhuman ears to follow the pitch of the original signal. The harmonic structure of our model is effective for detecting speechin the presence of background noise. YIN and the WWB algorithm were both outperformed, although we should note againthat their code was used as is, whereas ours was developed withthe task in mind. Thus, this comparison may not do them fulljustice.

F ESTIMATION OF CONCURRENT SPEECH BY MULTIPLE SPEAKERS,GROSS ERROR FOR A DIFFERENCE WITH THE REFERENCEHIGHER THAN 20% AND 10%

The experiments were performed in the same conditions as
described in Section III-B for HTC, and the results are presented
in Table IV. One can see that our algorithm again outperforms
YIN and the WWB algorithm in all the interference categories.
D. Multipitch Estimation

C. Validation on a Corpus of Speech Mixed With a Wide
Range of Interferences
In order to show the wide applicability of our method, we
also performed experiments using a corpus of 100 mixtures ofvoiced speech and interference [24], commonly used in CASA
research. In [4], half of the corpus is used for model parameter
estimation and the other half for system evaluation. As it is not
specified in that paper which part of the corpus was used forwhich purpose, we decided to use the full corpus as the evaluation set for comparison of the algorithms, which can only be
an advantage for the WWB algorithm. The results we present
for the WWB algorithm differ from the ones given in [4] as thecriterion we use is different. To be able to compare it with our
method, which does not perform a voiced/unvoiced decision, we
do not take into account errors on the estimation of the numberof s, but only look at the accuracy of the output of the pitch
determination algorithm. Moreover, we focus on theestimation of the main voiced speech, as we want here to show thatin a wide range of noisy
our algorithm robustly estimates theenvironments. The ten interferences are grouped into three categories: 1) those with no pitch; 2) those with some pitch qualities; and 3) other speech, as shown in Table III. The reference
contours for the ten voiced utterances were built using YIN
on the clean speech and manually corrected.

We present here results on the estimation of thecontourof the cochannel speech of two speakers speaking simultaneously with equal average power. We used again the databasementioned above [15] and produced a total of 150 mixed utterances, 50 for each of the “male–male,” “female–female,” and“male–female” patterns, using each utterance only once andmixing it with another such that two utterances of the same sentence were never mixed together. We used our algorithm in thesame experimental conditions as described in Section III-A for
contours.
clean single-speaker speech, but using two spline
The spline contours were initially flat and set to 155 and 296 Hz
in the male–female case, 112 and 168 Hz in the male–male case,and 252 and 378 Hz in the female–female case.
The evaluation was done in the following way: only times
inside the reliability mask of either of the two references werecounted; for each reference point, if either one of the two spline
contours lies within a criterion distance of the reference, we
considered the estimation correct. We present scores for two
criterion thresholds: 10% and 20%. For comparison, tests using
the WWB algorithm [4] introduced earlier were also performed,
using the code made available by its authors. YIN could not
be used as it does not perform multipitch estimation. Resultssummarized in Table V show that our algorithm outperforms
the WWB algorithm on this experiment.

1144IEEE TRANSACTIONS ON AUDIO, SPEECH, AND LANGUAGE PROCESSING, VOL. 15, NO. 4, MAY 2007

IV. CONCLUSION AND FUTURE WORK
We introduced a new model describing the spectrum as a sequence of spectral cluster models governed by a common
contour function, with smooth transitions in the temporal succession of the spectral structures. The model enables an accucontour on the whole utterance byrate estimation of thetaking advantage of its voiced parts in clean as well as noisy
environments. We explained how to optimize its parameters efficiently and performed several experiments to evaluate the accuracy of our model. On single-speaker clean speech, we obtained
good results which we compared with existing methods specifically designed for that task. On cochannel concurrent speech,
single-speaker speech mixed with white noise, pink noise, andon a corpus of single-speaker speech mixed with a variety ofinterfering sounds, we showed that our algorithm outperforms
existing methods.
We are currently working on using the precise parametricexpression and clustering of the spectrogram we obtained fornoise cancelling, speech enhancement, and speech separation.
We also intend to improve the modeling of the spectral structure using a new method that we describe in [25] in order toembed the formant structure into our model.
ACKNOWLEDGMENT
The authors would like to thank Dr. E. McDermott and P. Masurel for fruitful discussions and comments and the anonymous
reviewers for their valuable comments which greatly helped improve the quality of this article.
REFERENCES[1] W. J. Hess, Pitch Determination of Speech Signals. New York:
Springer, 1983.
[2] A. de Cheveigné and H. Kawahara, “Yin, a fundamental frequency estimator for speech and music,” J. Acoust. Soc. Amer., vol. 111, no. 4,
pp. 1917–1930, 2002.
[3] A. de Cheveigné, “Multiple F estimation,” in Computational Auditory
Scene Analysis: Principles, Algorithms and Applications, D.-L. Wang
and G. J. Brown, Eds. New York: IEEE Press/Wiley, 2006.
[4] M. Wu, D. Wang, and G. Brown, “A multipitch tracking algorithm fornoisy speech,” IEEE Trans. Speech Audio Process., vol. 11, no. 3, pp.
229–241, May 2003.
[5] T. Tolonen and M. Karjalainen, “A computationally efficient multipitch
analysis model,” IEEE Trans. Speech Audio Process., vol. 8, no. 6, pp.
708–716, Nov. 2000.
[6] Y. H. Gu and W. M. G. van Bokhoven, “Co-channel speech separation
using frequency bin nonlinear adaptive filter,” in Proc. ICASSP, 1991,
pp. 949–952.
[7] A. Klapuri and M. Davy, Signal Processing Methods for Music Transcription. New York: Springer, 2005.
[8] H. Fujisaki and S. Nagashima, “A model for synthesis of pitch contoursof connected speech,” Faculty of Eng., Univ. Tokyo, Tokyo, Japan,
1969, vol. 28, pp. 53–60, Annu. Rep. Eng. Res. Inst..
[9] H. Kameoka, T. Nishimoto, and S. Sagayama, “A multipitch analyzer
based on harmonic temporal structured clustering,” IEEE Trans. Audio,
Speech, Lang. Process., vol. 15, no. 3, pp. 982–994, Mar. 2007.
[10] I. Csiszar, “I -divergence geometry of probability distributions andminimization problems,” Ann. Probability, vol. 3, no. 1, pp. 146–158,
1975.
[11] P. Hedelin and D. Huber, “Pitch period determination of aperiodic
speech signals,” in Proc. ICASSP, 1990, pp. 361–364.
[12] W. H. Press, W. T. Vetterling, S. A. Teukolsky, and B. P. Flannery,
Numerical Recipes in C++: The Art of Scientific Computing, 2nd ed.
Cambridge, U.K.: Cambridge Univ. Press, 2002.
[13] W. I. Zangwill, Nonlinear Programming: A Unified Approach. Englewood Cliffs, NJ: Prentice-Hall, 1969.

[14] X. L. Meng and D. B. Rubin, “Maximum likelihood estimation via the
ECM algorithm,” Biometrika, vol. 80, no. 2, pp. 267–278, 1993.
[15] P. C. Bagshaw, S. M. Hiller, and M. A. Jack, “Enhanced pitch tracking
and the processing of F contours for computer and intonationteaching,” in Proc. Eurospeech, 1993, pp. 1003–1006.
[16] H. Kameoka, T. Nishimoto, and S. Sagayama, “Separation of harmonic
structures based on tied Gaussian mixture model and information
criterion for concurrent sounds,” in Proc. ICASSP, 2004, vol. 4, pp.
297–300.
[17] P. Boersma, “Accurate short-term analysis of the fundamental frequency and the harmonics-to-noise ratio of a sampled sound,” in Proc.
Inst. Phonetic Sci., 1993, vol. 17, pp. 97–110.
[18] P. Boersma and D. Weenin, “Praat System,” [Online]. Available: http://
www.fon.hum.uva.nl/praat/[19] D. J. Hermes, “Measurement of pitch by subharmonic summation,” J.
Acoust. Soc. Amer., vol. 83, pp. 257–264, 1988.
[20] Edinburgh Speech Tools Library. [Online]. Available: http://www.
cstr.ed.ac.uk/[21] Speech Filing System. [Online]. Available: http://www.phon.ucl.ac.uk/
resource/sfs/[22] B. Doval, “Estimation de la fréquence fondamentale des signaux
sonores,” Ph.D. dissertation, Univ. Pierre et Marie Curie, Paris,
France, 1994.
[23] H. Kawahara, H. Katayose, A. de Cheveigné, and R. D. Patterson,“Fixed point analysis of frequency to instantaneous frequency mapping
for accurate estimation of F and periodicity,” in Proc. Eurospeech,
1999, vol. 6, pp. 2781–2784.
[24] M. P. Cooke, “Modeling auditory processing and organisation,” Ph.D.
dissertation, Univ. Sheffield, Sheffield, U.K., 1993.
[25] H. Kameoka, J. Le Roux, N. Ono, and S. Sagayama, “Speech analyzer
using a joint estimation model of spectral envelope and fine structure,”in Proc. ICSLP, 2006, pp. 2502–2505.

Jonathan Le Roux (S’06) received the degree inmathematics from the Ecole Normale Supérieure,
Paris, France, the M.Sc. degree in partial differential
equations from the University of Paris XI, Paris, in
2001, and the M.Sc. degree in stochastic processes
from the University of Paris VI, Paris, in 2003. He iscurrently pursuing the Ph.D. degree at the Graduate
School of Computer Science, Telecommunications
and Electronics of Paris, University of Paris VI, andat the Graduate School of Information Science andTechnology, Department of Information Physics andComputing, University of Tokyo, Tokyo, Japan.
His research interests include acoustic signal processing, speech processing,computational acoustic scene analysis and language acquisition modeling.
Mr. Le Roux is a student member of the Acoustical Society of Japan (ASJ)
and the International Speech Communication Association (ISCA).

Hirokazu Kameoka (S’05) received the B.E. and
M.E. degrees from the University of Tokyo, Tokyo,
Japan, in 2002 and 2004, respectively. He is currently
pursuing the Ph.D. degree at the Graduate School ofInformation Science and Technology, University ofTokyo.
His research interests include acoustic signal processing, speech processing, and music processing.
Mr. Kameoka is a student member of the Instituteof Electronics, Information and Communication Engineers (IEICE), Information Processing Society ofJapan (IPSJ), Acoustical Society of Japan (ASJ), and International Speech Communication Association (ISCA). He was awarded the Yamashita Memorial Research Award from IPSJ, Best Student Paper Award Finalist at IEEE International Conference on Acoustics, Speech, and Signal Processing (ICASSP), andthe 20th Telecom System Technology Student Award from the Telecomunications Advancement Foundation (TAF), all in 2005.

LE ROUX et al.: SINGLE AND MULTIPLE

CONTOUR ESTIMATION

Nobutaka Ono (M’02) received the B.E., M.S., andPh.D degrees in mathematical engineering and information physics from the University of Tokyo, Tokyo,
Japan, in 1996, 1998, 2001, respectively.
He joined the Graduate School of Information Science and Technology, University of Tokyo, in 2001
as a Research Associate and became a Lecturer in2005. His research interests include acoustic signal
processing, speech processing, music processing,sensing and measurement, and auditory modeling.
He has been a Secretary of the technical committeeof psychological and physiological acoustics in Japan since 2005.
He is a member of the Institute of Electrical Engineers of Japan (IEEJ), theAcoustical Society of America (ASA), the Acoustical Society of Japan (ASJ),
the Society of Instrument and Control Engineers (SICE), and the Japan Societyof Applied Physics (JSAP). He received the Sato Paper Award from the ASJ in
2000 and the Igarashi Award at the Sensor Symposium on Sensors, Micromachines, and Applied Systems from the IEEJ in 2004.

Alain de Cheveigné (M’06) received the M.Sc. degree in electronics from the University of Paris VI,
Paris, France, in 1977, the Ph.D. degree in linguistics
from the University of Paris VII in 1982, and the Habilitation in Neurosciences from Paris VI in 2000.
He is currently a CNRS Senior Researcher with
the Audition Group, Laboratoire de Psychologie de laPerception, affiliated with CNRS, University of Paris
V and Ecole Normale Supérieure. His research interests include psychophysics, modeling of hearing
processes such as pitch and auditory scene analysis,
speech and music signal processing such as F0 and spectral estimation, audio indexing and retrieval, and signal processing of magnetoencephalography signals.
Dr. de Cheveigné is member of the Acoustical Society of America (ASA), theAcoustical Society of Japan (ASJ), the Society for Neuroscience (SFN), and the
Association for Research in Otolaryngology (ARO).

1145

Shigeki Sagayama (M’82) was born in Hyogo,
Japan, in 1948. He received the B.E., M.E., andPh.D. degrees from the University of Tokyo, Tokyo,
Japan, in 1972, 1974, and 1998, respectively, all inmathematical engineering and information physics.
He joined Nippon Telegraph and Telephone Public
Corporation (currently, NTT) in 1974 and started his
career in speech analysis, synthesis, and recognition
at NTT Laboratories, Musashino, Japan. From 1990
to 1993, he was Head of the Speech Processing Department, ATR Interpreting Telephony Laboratories,
Kyoto, Japan, pursuing an automatic speech translation project. From 1993 to1998, he was responsible for speech recognition, synthesis, and dialog systems
at NTT Human Interface Laboratories, Yokosuka, Japan. In 1998, he became a
Professor of the Graduate School of Information Science, Japan Advanced Institute of Science and Technology (JAIST), Ishikawa, Japan. In 2000, he wasappointed Professor of the Graduate School of Information Science and Technology (formerly Graduate School of Engineering), University of Tokyo. His
major research interests include processing and recognition of speech, music,
acoustic signals, hand writing, and images. He was the leader of anthropomorphic spoken dialog agent project (Galatea Project) from 2000 to 2003.
Prof. Sagayama is a member of the Acoustical Society of Japan (ASJ), Institute of Electronics, Information and Communications Engineers, Japan (IEICEJ), and Information Processing Society of Japan (IPSJ). He received the National Invention Award from the Institute of Invention of Japan in 1991, the Chief
Official’s Award for Research Achievement from the Science and Technology
Agency of Japan in 1996, and other academic awards including Paper Awards
from IEICEJ in 1996 and from IPSJ in 1995.