982

IEEE TRANSACTIONS ON AUDIO, SPEECH, AND LANGUAGE PROCESSING, VOL. 15, NO. 3, MARCH 2007

A Multipitch Analyzer Based on Harmonic
Temporal Structured Clustering
Hirokazu Kameoka, Student Member, IEEE, Takuya Nishimoto, and Shigeki Sagayama, Member, IEEE

Abstract—This paper proposes a multipitch analyzer called theharmonic temporal structured clustering (HTC) method, thatjointly estimates pitch, intensity, onset, duration, etc., of eachunderlying source in a multipitch audio signal. HTC decomposes
the energy patterns diffused in time-frequency space, i.e., thepower spectrum time series, into distinct clusters such that eachhas originated from a single source. The problem is equivalent toapproximating the observed power spectrum time series by superimposed HTC source models, whose parameters are associatedwith the acoustic features that we wish to extract. The update
equations of the HTC are explicitly derived by formulating the
HTC source model with a Gaussian kernel representation. We
verified through experiments the potential of the HTC method.
Index Terms—Computational acoustic scene analysis, harmonic
temporal structured clustering (HTC), multipitch analyzer.

Fig. 1. Wavelet spectrogram of real performed music signal ranging from T
to T in time direction and from X to X in frequency direction.

I. INTRODUCTION

along time. In auditory scene analysis (ASA), these two processes in human audition are generally called “segregation”
and “integration,” respectively.
It is quite obvious that the more accurate the segregation
process, the more reliable the result of the integration process.
On the other hand, we hope to know, if possible, the pitches
and the spectral component powers at preceding and succeeding
frames to estimate a high-precision result of the segregation
process at the current frame assuming they change gradually
over time. Therefore, these two processes should be done essentially in a cooperative way and not independently with successive estimations for even more reliable results. This belief
has led us to formulate a unified estimation framework for thetwo-dimensional structure of time-frequency power spectra, incontrast to the conventional strategy. The method that is presented in this article formulates the problem as a localization andshape detection of distinct spectrogram portions in the acoustic
scene (time-frequency space), which essentially amounts to theestimation of pitch contour, onset, duration, and timbre featureof each sound source.

W

E HAVE been working on a new method for computational acoustic scene analysis having in mind, forexample, a music content description system, a new equalizer
system enabling volume and bass/treble controls for eachseparate sound source, a music information retrieval (MIR)
system, and a precise acoustic environment [background music(BGM), phone ringing, interfering speech, etc.] detector fora wide range of speech applications. This paper describes a
new multipitch analyzer that jointly estimate acoustic features
such as pitch, onset, duration, energy, spectral, and temporalenvelopes of each underlying source in a multipitch acoustic
signal.
While the standard level of the numerous conventional multipitch analyzers has been considered to be far from practicaluse, recent pioneering ideas, e.g., graphical model-based [1],filterbank-based [2], nonparametric Kalman filtering-based[3], [4], multiagent-based [5], cochlear filtering-based [6], andparametric signal and spectrum modelings-based approaches[7]–[11] have brought remarkable progress to the practicalstep. Most of these methods take two major steps to resolve
the problem: they start with a separation or a pitch feature
extraction of concurrent sources at each short-time segment(frame) and then find the most likely overall pitch trajectories

Manuscript received May 6, 2005; revised July 11, 2006. The associate editor
coordinating the review of this manuscript and approving it for publication was
Dr. Michael Davies.
The authors are with Graduate School of Information Science and Technology, University of Tokyo, Tokyo 113-8656, Japan (e-mail: kameoka@
hil.t.u-tokyo.ac.jp; nishi@hil.t.u-tokyo.ac.jp; sagayama@hil.t.u-tokyo.ac.jp).
Color versions of Figs. 5–9 are available online at http://ieeexplore.ieee.org.
Digital Object Identifier 10.1109/TASL.2006.885248

II. FORMULATION
Consider an observed power spectrum time series(see Fig. 1), where and are log-frequency and time, definedon a domain of definition(1)
The problem of interest is to decompose this observed pattern
into sequential spectral streams, i.e., clusters, each of whichis assumed to have originated from a single distinct source activation. This problem is an unsupervised categorization of theat each coordinate.
energy density

1558-7916/$25.00 © 2006 IEEE

KAMEOKA et al.: MULTIPITCH ANALYZER BASED ON HTC

983

clustering to minimize w.r.t.
the above measure

Fig. 2. Profile of the k th HTC source model q( ;22) ((19)).

The unknown variables being
are solving is summarized asx; t

The observed energy density
at each coordinate
has not necessarily originated from a single component
in a single source but is in the general case rather a mixtureof the energy patterns generated from other components ordifferent sources. Thus, it is natural to assume that the energy
has an unknown fuzzy membership to the thdensity.
source, introduced as a spectral masking function
Approximately assuming that observed power spectral densities
are the sum of the independent power densities generated fromunderlying sources, which is true in the expectation sense (if we
assume the phase of every component is uniformly random),must satisfy

and

and

, the sum overof(5)
, the optimization we(6)
In order to find the optimal and , we shall find it most convenient to recursively optimize and while keeping the othervariable fixed as there is no analytical solution to (6). As each
iteration necessarily decreases , which is bounded below,
and gradually converge to a stationary point. Note, however,that with this procedure and do not always converge to a
global optimum but possibly to a local optimum.
when is fixed can be obtained by finding
The optimal
the extreme value of

(2)

(7)

Therefore,
denotes the decomposed spectral
density of the th source, which we will refer to as the thsuch that allcluster. We may wish to decompose
clusters are temporally continuous and also have a harmonic
structure, if we assume all sources are periodic signals andtheir components evolve gradually over time. Let us define, governed by parameter vector , that charachereterizes the typical spectrogram of a single source (a graphicalrepresentation can be seen in Fig. 2). Note that the class of thesources mentioned here, in general, includes not only periodic
signals but also nonperiodic signals such as white/pink noises,
drum sounds or any others, as far as those spectrograms can, webe modeled with a suitable expression. Using
are now able to assess a “goodness” of the partitioned clusterwith the Kullback–Leibler (KL) divergenceofandwhere is a Lagrange undetermined multiplier. The partialis thus given as
derivative of the integrand w.r.t.

(3)

(8)
Setting it to zero, one obtains(9)
From (2) and (9), we have(10)
Substituting (10) in (9), one finally obtains the optimal spectral
masking function for a fixed as
(11)

with
(4)
The condition given as (4) makes the distortion measure (3) nonnegative. In the case of discrete-time observations, the distortion
measure is still given by (3) where the integral is replaced by thesum over all discrete points of and . Now one notices that as
andbecome closer, (3) approaches
zero. Hence, one can choose as the global cost function of the

Substituting (11) in (5), it becomes clear that during this clustering we are also iteratively decreasing the KL divergence beand the sum of
tween the whole observed spectrogram
over
(12)

Therefore, this clustering is understood as a geometric optimal approximation ofusing the model.

984

IEEE TRANSACTIONS ON AUDIO, SPEECH, AND LANGUAGE PROCESSING, VOL. 15, NO. 3, MARCH 2007

Note that this result proves without using Bayes rules the convergence of the expectation-maximization (EM) algorithm, as
as a complete
one regards as missing data anddata pdf. The correspondence to the EM algorithm
becomes much clearer by comparing the result of (11) and theterms in that depends on(13)
Fig. 3. Cutting plane of qwith the( ;22) at time .
x; ttfunction, given byThis discussion implies that practically the same procedureas the EM algorithm can be used even thoughandare not pdf’s.
updated by theWith fixed membership degree
E-step (11), the parameter , on the other hand, should be
updated by the M-step

Given the pitch contourin th HTC source model, the. Now let the frecontour of the th partial isquency spread of each harmonic component be approximated by
a Gaussian distribution function when the spectra are obtainedby1 the wavelet transform (constant transform) using Gaborwavelet basis function. The accuracy of this approximation isinvestigated through a numerical computation in Appendix I.
the power envelope of the th partial (preDenoting by
sumed to be a function that is normalizable sincehas
to satisfy (4)), such that
(16)
then the normalized energy density of the th partial in the th
and a
HTC source model is given as a multiplication ofGaussian distribution centered at(14)
depending on the specific form of. The update equais introtion for will be given in Section V after
duced. We will call this approach “harmonic temporal structured
is asclustering (HTC)” in the particular case where
sumed to have a harmonic structure.
The HTC procedure is summarized as follows. The input
to the system is an observed (known) signal, characterized by, where and are log-frequency andits spectrogramof th source/stream istime. The membership degree
unknown (spectral masking function). On the other hand, thespectrogram of the th source can be modeled by a function, whereis the set of model parameters. These
are the unknown variables that we want to estimate. The HTC
method works by iteratively updating the estimates of: 1)
with fixed by (11) and 2) with
fixed.
III. HTC MODEL
A. Model Representation
In this section, we will introduce the HTC source model. Let us assume through the rest of this paper that all
sources are periodic signals having smooth power envelopes.
Supposing the pitch contour during a single source activation isexpressed with a polynomial (imagine vibrato or glissando)(15)
a cutting plane ofat particular time represents aharmonic structure of pitch frequency(see Fig. 3).

(17)
where denotes the frequency spread of every partial, andis the relative energy of the th partial, satisfying(18)
Therefore, the power density of the th HTC source model
as a whole (see Fig. 2) can be written as

(19)
where
indicates the total energy of the th source. The sum
isof a number of the HTC source models
supposed to be a model for the observed overall spectrogram.
Let us now discuss how we should model the power enve. In the general case where one does notlope function
know in advance what the sources are, it is perhaps wise not tobased only upon a
use a physically oriented model forparticular sound production mechanism. It is thus important toas generic a model as possible.
introduce foris a function that is
What we hope to choose fortemporally continuous, nonnegative, having a time spread from
1Having music applications in mind, since usually the musical notes areequally spaced on a log-frequency scale, there is a need to use a transformthat has a variable window length characteristic. In particular, it is difficult
to distinguish the sources of lower notes with short-time Fourier transform
spectra, for example, that has a fixed window length across all frequencies. Thisis one of the reasons we chose to use wavelet transform as a front-end for HTC.

KAMEOKA et al.: MULTIPITCH ANALYZER BASED ON HTC

985

where

(23)

Fig. 4. Power envelope function U

the problem could be equivalently simplified by further breaking
each cluster down into-labeled subclusters, associated.
with the subsource modelthat
Introducing now another masking function
into thedecomposes the th partitioned cluster
th subcluster and satisfies for all andt( ) [(20)].

minus to plus infinity (assuming the Gabor-wavelet basis as themother wavelet) and adaptable to various curves. Furthermore,
it should satisfy (16). We came up with a function satisfying allthese requirements, given as(20)
is the center of the forefront Gaussian, that could be considered as an onset time estimate,
the weight parameter foreach kernel, that allows the function to have variable shapes. To
must only be normalized to unity
satisfy (16),(21)
The particularity of this function is that the centers of theGaussian function kernels are spaced by a distance proportionalwith a proportionto the common diffusion parameter
ality coefficient , which we henceforth set to 1 (see Fig. 4).
This tying ensures the smoothness of the curve by preventing
also
adjacent kernels to be separated from each other.
inworks as a parameter to make a linear stretch of
the time direction allowing to express various durations ofsources. Moreover, by forbidding switches in the position of
the kernels, it reduces the singularity of the system, improving
the optimization perspectives.
All the parameters of the HTC model are listed in Table I.
B. Subclustering
While the representation of, the HTC source
model, has been introduced, we are not yet able to obtain theanalytical solution for (14) for each parameter in . However,as the HTC source model is specified as the sum of the subsource model

we have the Jensen’s inequality for all , as shown in (24) at thebottom of the page, and equality holds when(25)
The proof is omitted since it can be easily obtained by following
the same way as in Section II. This means that if and only if
is given by (25),is minimized to. One notices that obtaining
the global cost function
parameter update equations throughseems much easier than(14) since the summation over , , and no longer appears
with theinside the logarithm function. After makingby (25), one can indirectly decrease theupdate ofobjective function by decreasing
through an update of .
This is because is always guaranteed by the inequation (24) to. Decreasing
w.r.t. can be done
be even smaller thanby increasing (26), as shown at the bottom of the next page.
The HTC procedure in the particular case where the HTC
source model is given by (22) is again summarized as follows.
, among the total enThe masking function
for each subcluster is unknown. On theergy density
other hand, the HTC subsource model for the subcluster energy, where is the set of model paramedensity isters (pitch, spectral envelope, temporal envelope, intensity, harmonicity, cf. Table I). The HTC works by iteratively updating
the estimates ofby1) (E-step)
(27)

(22)
2) (M-step)

using

.

(24)

986

IEEE TRANSACTIONS ON AUDIO, SPEECH, AND LANGUAGE PROCESSING, VOL. 15, NO. 3, MARCH 2007

TABLE IPARAMETERS OF THE HTC SOURCE MODEL

With this procedure, converges to a stationary point as well aswith the procedure described in Section III.
The philosophy we adopt in the HTC is to approximate asource spectrogram in a real environment as closely as possible
by a mathematically simple and compact model characterized
by meaningful parameters rather than to try to make a perfect fit
by an extremely complex model whose parameters are no longer
meaningful.

nomial-like distribution of parameterof presenting the problem such that

. This way(33)
enables us to extend it as a maximum a posteriori (MAP) estimation problem and to introduce prior distributions on the parameters as follows, using Bayes theorem

IV. EXTENSION TO MAP AND THE USE OF PRIORS
Keeping only the terms depending on in (12) and taking the
opposite, one defines the following function to maximize w.r.t.
:
(28)
(34)
Using this function and letting

be
(29)

one can derive the likelihood of the parameter(30)
where Dirac delta

ensures that the density is zero if(31)

The parameter

This optimization can be performed by iteratively updatingby (27) and using.
Prior distribution works as a penalty function that tries to keep
the parameters within a specified range. By introducing such a, it becomes possible to prevent subharprior distribution on
monic (half-pitch) errors, as the resulting source model would
usually have a harmonic structure with zero power for all the oddorder harmonics, which is abnormal for speech and for many, on the other hand,
instruments. A prior distribution on
helps to avoid overfitting many source models to the observed
power envelope of a single source, as the resulting individual
source models in this case would often have abnormal powerenvelopes.
For this purpose, we apply the Dirichlet distribution as theprior distribution, which is explicitly given by

given by
(32)

(35)

is theensures that we obtain a probability measure whereGamma function. One can indeed see this probability as thefollowing a multijoint probability of all the variables

(36)

(26)

KAMEOKA et al.: MULTIPITCH ANALYZER BASED ON HTC

987

Fig. 5. Three-dimensional and top views of the observed power spectrum time series (top and bottom left) and the optimized HTC source models (top and bottom
right).

whereand
are the most preferred “expected” values ofandsuch thatand, , andregulate the strength of the priors. Whenandare zero,andbecome uniform distributions. The maxand
are taken respectively whenimum value forandfor all and if
and.
The choice of this particular distribution allows us to give an
and
. The
analytical form of the update equations ofjoint prior distribution is now given byfinding the extreme value of (38), as shown at the bottom of thepage.
V. PARAMETER UPDATE EQUATIONSHaving music applications in mind, let us assume here thatand
all pitch contours are zero-order polynomials:
each partial stream in the HTC source model has the same powerenvelope (only a single power envelope function is assumed in, for the purpose of rethe HTC source model):ducing the dimensionality of the features to extract. From (23),is given bylogarithmic subsource model(37)
Denoting by , andthe Lagrange multipliers for ,and, respectively, the update equation of (M-step)that is guaranteed to increase
can be derived by

(39)

(38)

988

IEEE TRANSACTIONS ON AUDIO, SPEECH, AND LANGUAGE PROCESSING, VOL. 15, NO. 3, MARCH 2007

Fig. 6. Estimates of  ,  , Y  (top), and the reference MIDI data displayedin piano-roll form (bottom).

Setting to zero the partial derivative of (38), the update equationof each parameter at M-step of the th iteration is derived asfollows (see appendix for its derivation):

2

Fig. 7. Accuracy =!
100 (percent) of the result for each time spread parameter d ranging from 1 to 40 when ! is set to 2 100 (top) and 2 800(bottom) rad/s.

22

the coordinate descent method [12] and the corresponding optimization procedure is called the expectation-constrained maximization (ECM) algorithm [13].

(40)
(41)
(42)

(43)

(44)
(45)

(46)

VI. RELATION WITH PREFEST [10] AND
OUR PREVIOUS WORK [11]A cutting plane of the HTC source model at a particular time
gives a similar representation to the harmonic structure model
introduced by Goto in [10] and followed independently by ourprevious work described in [11]. Goto used a tone model, a harmonic structure model represented by a Gaussian mixture, andtried to model an observed short-time power spectrum by the
mixture of a large enough number of the tone models densely
spaced by a fixed interval. He then used the EM algorithm frameby frame for estimating the MAP mixture weights of the tonemodels to measure pitch likelihoods from an observed power
spectrum. In our previous work described in [11], on the other
hand, we started with a formulation based on the clustering
principle using a harmonic structured Gaussian mixture clustermodel in order to try to directly estimate each mean parameter,
the pitch estimate itself. The source spectrogram model introduced in the HTC method designs the overall smooth temporal
evolution of the spectral structure model introduced by us and
Goto.
VII. EXPERIMENTAL EVALUATION

where denotes the subcluster densities

A. Conditions

The superscript
refers to the iteration cycle. Note that the
update equations of the higher order coefficientsof the pitch trajectory model can all be derived as well as,
if necessary.
Note, however, that the update (40)–(46) do not ensure themaximization in the M-step but guarantee the increase of. This way of performing the updates is referred to as

To verify the performance of the HTC method, we evaluatedpitch estimation accuracies with a set of real performed musicsignals excerpted from the RWC music database [14]. The first23 s of each music signal data were used for the evaluation.
We implemented a GUI editor to create a ground truth data
set of pitch sequences (a screenshot of the GUI editor can beseen in Fig. 9). Since the RWC database also includes a supplement MIDI file associated with each real-performed musicsignal data, we created each ground truth data set with the GUI

KAMEOKA et al.: MULTIPITCH ANALYZER BASED ON HTC

989

The Power spectrum time series was analyzed by the wavelet
analysis) using Gabor-wavelet basis
transform (constantfunctions with a time resolution of 16 ms for the lowest
frequency subband on an input signal digitalized at a 16-kHz
sampling rate. To speed up the computation time, we set the time
resolution across all the subbands equally to 16 ms. The lower
bound of the frequency range and the frequency resolution were
60 Hz and 12 cents (where 12 cents amount to one octave),respectively. The initial parameters offor the HTC source models were automatically determined bypicking the 60 largest peaks in the observed spectrogram of 400
consecutive frames (6.4 s). After the parameters converged, thesource model, whose energy per unit time given by
was smaller than a threshold, was considered to be silent. The
experimental conditions are shown in detail in Table III.
We chose2 “PreFEst” [10] for comparison, as it is one of
the most frequently cited works dedicated to multipitch analysis. Since PreFEst extracts only the most dominant pitch trajectory and does not include a specific procedure of estimating
the number of sources, we included intensity thresholding aswell for the pitch candidate truncation.
As the HTC method generates pitch, onset time and offset
time with continuous values, we quantize them to the closest
note and the closest frame number in order to match with theformat of the reference. Using the hand-labeled ground truth
data as references, pitch accuracies were computed by

Fig. 8. Graphical representations of (62) with d = 10 (top) and d = 40(bottom) when ! is set to 2 100 rad/s together with the fitted Gaussians.

2

where
number of the total frames of the voiced part;number of deletion errors;
number of insertion errors;
number of substitution errors.
B. Results

Fig. 9. Screenshot of the GUI editor we implemented to create the ground truth
data set of note pitches, onsets, and durations. The note events of the supplement
MIDI data included in the RWC database, which are not temporally aligned withthe corresponding real performed signal data, are displayed as rectangular objects over the spectrogram of the real performed signal. We are then able to edit
the rectangular objects to align carefully the onset and offset times according to
the background spectrogram.

by adjusting the onset and offset times of each MIDI file witha spectrogram image for a background. The list of the experimental data sets and the time average of the number of concurrent sources, the total number of the frames for all pitches
divided by the number of the frames, are shown in Table II.

A typical example of the pitch, onset, and offset estimateson a particular test data is shown in Fig. 6 together with thehand-labeled ground truth data. The optimized model and theobserved power spectrum time series are shown with 3-D andgrayscale displays in Fig. 5.
To validate the performance of the proposed method, we compared the highest accuracy of the HTC method with that of thePreFEst among all the thresholds that were tested, which also
shows the limit of the potential capability. The highest accuracies of PreFEst and HTC among all the thresholds we tested areshown in Table IV together with the number of insertion, deletion, and substitution errors, respectively. Comparing these accuracies between PreFEst and HTC, HTC outperforms PreFEst
for most of the data, which verifies its potential.
The workstation used to perform the experiments had a Pentium IV processor with 3.2-GHz clock speed and 2-GB memory.
With our implementation with the conditions listed in Table III,
2Note that we implemented for the evaluation only the module called“PreFEst-core,” a frame-wise pitch likelihood estimation, and not included theone called “PreFEst-back-end,” a multiagent-based pitch tracking algorithm.
Refer to [10] for their details.

990IEEE TRANSACTIONS ON AUDIO, SPEECH, AND LANGUAGE PROCESSING, VOL. 15, NO. 3, MARCH 2007

TABLE II
LIST OF THE EXPERIMENTAL DATA EXCERPTED FROM RWC MUSIC DATABASE.[14]B. Audio Coding

TABLE III
EXPERIMENTAL CONDITIONS

must be a good approximation to theSince
is the optimal
observed power spectrum time series wheremodel parameter vector, the residual power spectrum(48)
may be effectively compressed by the well-known Huffman encoding. As the original power spectrum can be restored with theresidual and the optimized parameter vector, our method is thus
also expected to be used as an effective audio coding application
if acoustic signals can be intelligibly reconstructed from power
spectra. However, its effectiveness over traditional lossless andlossy coding techniques will need to be evaluated in future investigations.
IX. CONCLUSION

the computational time for analyzing an acoustic signal of 25.6-slength was about 2 min. In most cases, the parameters of the
HTC source models converged within less than 100 iteration
cycles.
We also compared the HTC performances with different conditions: the time range of an analyzing spectrogram segment of100, 200, and 400 frames, and the number of the HTC source
models of 15, 30, and 60, respectively. Comparative results areshown in Table V. From the results, one can see that the larger
the time range of a spectrogram segment, the higher the accuracies. This shows that the domain of definition of should be aslarge as possible for a higher performance of the HTC.

We established a new framework for multipitch analysis based upon two-dimensional geometric modeling andestimation of the distinct spectral streams localized in thetime-frequency space, the acoustic scene, and investigatedthrough experiments its effectiveness over the conventional
method, PreFEst.
The method described in this paper still has many interesting
issues to consider, e.g., estimation of the number of note events
without relying on heuristic thresholding, further precise modeling of pitch contour, and the inharmonicity factor parameterand a further investigation on its application.
APPENDIX IJUSTIFICATION FOR THE ASSUMPTION OF THE FREQUENCYSPREAD BEING A GAUSSIAN DISTRIBUTION

VIII. DISCUSSION OF THE POTENTIAL APPLICATIONS
A. Sound Source Segregation
It should be emphasized that the HTC method can also be
used to extract a portion of the power spectrum time series associated with a single source by

anda given signal and its Fourier
Denoting bytransform, we define the Fourier transform pair as follows:
(49)
(50)

(47)
Denoting by
where is the optimized model parameter vector. Using a decoding technique (such as the inverse wavelet transform, phase
vocoder, etc.) to reconstruct acoustic signals from the power
spectrum time series, the HTC also enables the separation ofsources.

the wavelet basis function
(51)

where
is the analyzing wavelet. and refer to the scaleand shift parameters, which have the dimensions of period andKAMEOKA et al.: MULTIPITCH ANALYZER BASED ON HTC

991

TABLE IV
ACCURACIES OF THE PREFEST [10] AND THE HTC

TABLE VCOMPARISON OF THE HTC PERFORMANCES WITH DIFFERENT RANGES OF A SPECTROGRAM SEGMENT AND THE NUMBER OF SOURCE MODELS

time, respectively. Based on the above definitions and the generalized Parseval’s theorem, the continuous wavelet transformis defined in the frequency domain as an inner product between
and, the Fourier transform ofdoes not depend
from which we obviously see thaton and is thus uniform over time.
a Gabor function of frequency 1
Now, if we choose forrad/s with time spread

(52)

(59)

Denoting by

the Fourier transform ofsuch thatis then a Gaussian function centered at 1 rad/s(53)

then the wavelet transform amounts to an inverse Fourier transform of the subband-filtered spectrum
(54)

(60)
is thus a Gaussian subband filter with center frerad s. Substituting (60) into (58), the power denquency ofsity in the period domain is given explicitly as(61)

is a sinuLet us now consider the very simple case where
soidal wave and see how the frequency spread at a certain time
be anslice of the wavelet scalogram is given by. Letting
analytic signal with an angular frequency of and with an amplitude ofin the log-frequency domain
The density
at a particular time slice ofassuch thatis thus given

(62)
(55)
its Fourier transform is then given by(56)
denotes the Dirac delta function. Substituting (56)
whereinto (54), one immediately obtains(57)
Its power is thus given by(58)

which is obviously not a Gaussian distribution function.
In order to justify the assumption made in the HTC source
model (Section III) that the frequency spread of the wavelet
power spectra is close to a Gaussian distribution, we will now
investigate through a numerical computation how much themean parameter of the Gaussian distribution that best approximates the above function with regard to the KL-divergencecriterion deviates from the true frequency . Denoting by
the Gaussian distribution model of the mean parameter
(63)

992

IEEE TRANSACTIONS ON AUDIO, SPEECH, AND LANGUAGE PROCESSING, VOL. 15, NO. 3, MARCH 2007

the KL divergence between

andis defined as

As we find that the denominator of the above is equal to (68),
one finally obtains
(71)

(64)
with

Onset Time:w.r.t.

: Setting to zero the partial derivative of(65)
The optimal , that is supposed to be here the estimate of , canthus be computed by finding the extreme value of . Setting tozero the partial derivative of with respect to , which is givenas(72)
one obtains
(73)

(66)
one obtains

Relative Energy of the Harmonics:
w.r.t.
the partial derivative of: Setting to zero

(67)
Since it is impossible to obtain analytically both the integrals
in the numerator and denominator in (67), we give numericalis set to
100 and
800 rad/s,
results for when
100 (%) of therespectively. Fig. 7 shows the accuracy
result for each time spread parameter ranging from 1 to 40.
, the value we chose in the experiment describedWhen
100 (%) within Section VII, the accuracy for both
100 and
800 were around 100.002 (%). This numerical result shows that the assumption of the frequency spreadbeing a Gaussian distribution may not necessarily affect critically the result of the pitch estimations. Graphical representaand
can be seen in Fig. 8,
tions of (62) withfrom which one is able to see how each of them is close to a
Gaussian distribution.

(74)
one obtains(75)

From (75) and (18), the Lagrange multiplierplicitly as

is given ex-(76)
which finally gives us the following:

APPENDIX IIDERIVATION OF THE PARAMETER UPDATE EQUATIONS
We will show in this section how the parameter update(40)–(46) were derived.
Energy: : From (4), one immediately obtains

Coefficients of the Power Envelope Function:
w.r.t.
to zero the partial derivative of(77): Setting(68)
Constant Term of the Pitch Contour:
the partial derivative of
w.r.t.

: Setting to zero(78)
one obtains
(79)

(69)
one obtains

From (79) and (21), the Lagrange multiplierplicitly as(70)

is given ex-(80)

KAMEOKA et al.: MULTIPITCH ANALYZER BASED ON HTC

993

which finally gives us the following:

Duration:w.r.t.

(81): Setting to zero the partial derivative of(82)
can be rewritten in the form (a quadratic equation)
(83)
where

with

, from which one finally obtains[4] M. Abe and S. Ando, “Auditory scene analysis based on time-frequency
integration of shared FM and AM (II): Optimum time-domain integration and stream sound reconstruction,” (in Japanese) Trans. IEICE, vol.
J83-D-II, no. 2, pp. 468–477, 2000.
[5] T. Nakatani, “Computational auditory scene analysis based on residuedriven architecture and its application to mixed speech recognition,”
Ph.D. dissertation, Kyoto Univ., Kyoto, Japan, 2002.
[6] M. Wu, D. Wang, and G. J. Brown, “A multipitch tracking algorithm
for noisy speech,” IEEE Trans. Speech Audio Process., vol. 11, no. 3,
pp. 229–241, May 2003.
[7] M. Feder and E. Weinstein, “Parameter estimation of superimposed
signals using the EM algorithm,” IEEE Trans. Acoust., Speech, Signal
Process., vol. 36, no. 4, pp. 477–489, Apr. 1988.
[8] D. Chazan, Y. Stettiner, and D. Malah, “Optimal multi-pitch estimation
using the EM algorithm for co-channel speech separation,” in Proc.
IEEE ICASSP, 1993, vol. 2, pp. 728–731.
[9] S. Godsill and M. Davy, “Baysian harmonic models for musical pitch
estimation and analysis,” in Proc. IEEE ICASSP, 2002, vol. 2, pp.
1769–1772.
[10] M. Goto, “A real-time music-scene-description system: Predominant-F0 estimation for detecting melody and bass lines in real-world
audio signals,” ISCA J., vol. 43, no. 4, pp. 311–329, 2004.
[11] 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. IEEE ICASSP, 2004, vol. 4, pp.
297–300.
[12] W. I. Zangwill, Nonlinear Programming: A Unified Approach. Englewood Cliffs, NJ: Prentice-Hall, 1969.
[13] X. L. Meng and D. B. Rubin, “Maximum likelihood estimation via the
ECM algorithm,” Biometrika, vol. 80, no. 2, pp. 267–278, 1993.
[14] M. Goto, H. Hashiguchi, T. Nishimura, and R. Oka, “RWC music database: Popular, classical, and jazz music database,” in Proc. ISMIR,
2002, pp. 287–288.
[15] H. Kameoka, T. Nishimoto, and S. Sagayama, “Audio stream segregation of multi-pitch music signals based on time-space clustering using
Gaussian kernel 2-dimensional model,” in Proc. IEEE ICASSP, 2005,
vol. 3, pp. 5–8.

(84)
Frequency Spread of Partial Distributions:zero the partial derivative of
w.r.t.

: Setting to(85)
one obtains(86)
ACKNOWLEDGMENT
The authors would like to thank the reviewers for their valuable comments, Dr. M. Goto, J. Le Roux, and Dr. N. Ono for
their fruitful discussions on this work, and Dr. M. Schuster forhis help with the English language.
REFERENCES[1] K. Kashino, K. Nakadai, and H. Tanaka, “Organization of hierarchicalperceptual sounds: Music scene analysis with autonomous processing
modules and a quantitive information integration mechanism,” in Proc.
IJCAI , 1995, vol. 1, pp. 158–164.
[2] A. Klapuri, T. Virtanen, and J. Holm, “Robust multipitch estimation for
the analysis and manipulation of polyphonic musical signals,” in Proc.
COST-G6 Conf. Digital Audio Effects, 2000, pp. 233–236.
[3] K. Nishi and S. Ando, “Optimum harmonics tracking filter for auditory
scene analysis,” in Proc. IEEE ICASSP ’96, 1996, vol. 1, pp. 573–576.

Hirokazu Kameoka (S’05) received the B.E. and
M.E. degrees from University of Tokyo, Tokyo,
Japan, in 2002 and 2004, respectively. He is currently pursuing the Ph.D. degree at the Graduate
School of Information Science and Technology,
University of Tokyo.
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), IPSJ, Acoustical Society ofJapan (ASJ), and International Speech Communication Association (ISCA). He
was awarded the Yamashita Memorial Research Award from the Information
Processing Society of Japan (IPSJ), Best Student Paper Award Finalist at theIEEE International Conference on Acoustics, Speech, and Signal Processing(ICASSP), and the 20th Telecom System Technology Student Award from the
Telecomunications Advancement Foundation (TAF), all in 2005.

Takuya Nishimoto received the B.E. and M.E. degrees from Waseda University, Tokyo, Japan, in 1993
and 1995, respectively.
He is a Research Associate at the Graduate Schoolof Information Science and Technology, Universityof Tokyo. His research interests include spoken dialogue systems and human–machine interfaces.
Mr. Nishimoto is a member of the Instituteof Electronics, Information, and Communication
Engineers (IEICE), Japan, Information Processing
Society of Japan (IPSJ), Acoustical Society of Japan(ASJ), Japanese Society for Artificial Intelligence (JSAI), and Human Interface
Society (HIS).

994

IEEE TRANSACTIONS ON AUDIO, SPEECH, AND LANGUAGE PROCESSING, VOL. 15, NO. 3, MARCH 2007

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 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
was appointed Professor of the Graduate School of Information Science andTechnology (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), IEICEJ, and IPSJ. He received the National Invention Award from the Institute ofInvention of Japan in 1991, the Chief Official’s Award for Research Achievement from the Science and Technology Agency of Japan in 1996, and otheracademic awards including Paper Awards from the Institute of Electronics, Information, and Communications Engineers (IEICEJ) Japan, in 1996 and fromthe Information Processing Society of Japan (IPSJ) in 1995.