AUXILIARY FUNCTION APPROACH TO PARAMETER ESTIMATION OF CONSTRAINEDSINUSOIDAL MODEL FOR MONAURAL SPEECH SEPARATION
Hirokazu Kameoka ∗Nobutaka Ono, and Shigeki Sagayama

NTT Communication Science Laboratories,
3-1 Morinosato Wakamiya, Atsugi,
Kanagawa, 243-0198, Japan

Graduate School of Information Scienceand Technology, The University of Tokyo
7-3-1 Hongo, Bunkyo-ku, Tokyo, 113-8656, Japan

ABSTRACT

While McAulay et al. used as the model a mixture of K
pure tone signals, one can also consider the use of the superposition of K harmonic signals (source signal composed of
N harmonic components, where the frequency of n-th harmonic component is n times the F0 ). This model is often
used for single channel source separation especially when thetarget mixed signal is assumed to be composed of harmonic
signals [3, 4, 5]. Most of the source separation approaches
using this model are based on gradient search or sampling
methods[3, 4, 5]. However, this kind of numerical computation is often beset with local optimum problems and computational costs. A global optimum is not guaranteed to beobtained unless, in the case of the gradient search methods
the iterative computation is led to convergence for an inﬁnity of initial points, or in the case of the sampling methods
an inﬁnite number of trial is performed. For that reason, theproblem is to know if the search for the solution can be performed with a low computational cost or if it has an ability to
avoid local optima.
As explained above, albeit the sinusoidal signal model
represents extremely well acoustic signals such as speech andmusic, room was left for discussion on how to estimate its parameters. Against this background, we describe in this papera new optimization algorithm to obtain the maximum likelihood parameter of the constrained sinusoidal model.

We introduce in this paper an auxiliary function approach
to parameter estimation of the constrained sinusoidal model,
which enables us to derive a complex-spectrum-domain EMlike multiple F0 estimation algorithm. Through simulations,
we evaluated the performance of the presented method in theability to avoid locally optimal solutions. We implemented
a monaural speech separation system based on the presented
method and conﬁrmed its performance on compound signalsof real speech.
Index Terms— Acoustic signal analysis, Speech enhancement, Optimization methods, Parameter estimation
1. INTRODUCTION
Many conventional methods for multiple F0 estimation are
based on the power spectrum domain approach, in which theinﬂuence of the interferences between frequency components
of different sources are assumed negligible. However, it becomes usually difﬁcult to infer F0 s only from power spectrum
when two voices are mixed with close F0 s. In such a situation not only the harmonic structure (powers of harmonics)
but also the phase of each component is an important cue for
precise estimation of F0 .
Generally speaking, if the compound signal were separated into single voices, then it would be a simple matter toobtain an F0 estimate and phase estimates of the harmonics
from each stream. On the other hand, if the F0 s and the phases
were known in advance, these information could be very useful for any separation algorithms. This leads to a “chicken andegg” situation: F0 /phase estimation and source separation areeach a prerequisite of the other. This fact leads us naturally
to formulate F0 /phase estimation and source separation as a
joint optimization problem. The method described in this paper performs F0 /phase estimation step and source separation
step iteratively using a constrained sinusoidal model.
The range of application of sinusoidal model has widened
to Text-To-Speech synthesis, speech modiﬁcation, coding,etc. since McAulay et al. [1] showed that the sinusoidal signal model could be applied to Analysis-by-Synthesis systems
to obtain high-quality synthesized speech. In particular, as
the possibility to generate high-quality synthesized speechshows that the sinusoidal signal model represents extremely
well acoustic signals such as speech and music, we can havehigh expectations for its application to source separation.

2. PROBLEM SETTING
Consider as the time-varying acoustic signal the sum of pseudoperiodic signals given in an analytic signal representation:s(t) =

Ak,n (t)ejnθk (t) , t ∈ (−∞, ∞),(1)

k=1 n=1

where the instantaneous phase θk (t) of the fundamental component, and the instantaneous complex amplitude Ak,n (t)of the n-th partial component are the unknown parameters. μk (t) = θ̇k (t) amounts to the instantaneous F0 andak,n (t) = |Ak,n (t)| the instantaneous amplitude, which areboth assumed here to change gradually over time. These are
the free parameters that one wants to estimate, which wedenote for convenience by Θ(t). Now letting y(t) be theobserved signal of interest, we assume the following model:y(t) = s(t) + n(t), t ∈ (−∞, ∞),(2)

where n(t) is a Gaussian white noise. The maximum likelihood estimate of Θ(t) can thus be obtained by minimizing∗ The author performed the work while at the University of Tokyo.

1-4244-1484-9/08/$25.00 ©2008 IEEE

K 
N


29

ICASSP 2008

the L2 norm of the error signal in t ∈ R(−∞, ∞):
 
2


y(t) − s(t) dt.

Lemma 2. Denoting by Φ(Θ) the objective function, and by
Φ+ (Θ, m) the auxiliary function of Φ(Θ), then the objectivefunction Φ(Θ) can be decreased monotonically by minimizing Φ+ (Θ, m) iteratively with respect to m = (m1 , · · · , mJ )and with respect to Θ:
m = argmin Φ+ Θ, m
(10)
m


Θ = argmin Φ+ Θ, m
(11)
Θ

(3)

R

We now show that this time domain objective can beequivalently formalized in the time-frequency domain deﬁned by the Gabor transform (STFT).
Lemma 1. The time-frequency components of y(t) and s(t)by Gabor transform is by deﬁnition given by


,
(4)
Gy (ω, t)  y(u), ψω,t (u)

 u∈R
,
(5)
Gs (ω, t)  s(u), ψω,t (u)If Φ(Θ) is bounded below, then the parameter Θ converges
to a stationary point.
Proof. Suppose we set Θ to an arbitrary value Θ. We
will prove that Φ(Θ) is non-increasing after the update
Eq. (10) and Eq. (11). From Eq. (10), one obtains
Φ(Θ) = Φ+ (Θ, m), and it is obvious from Eq. (11) that
Φ+ (Θ, m)  Φ+ (Θ, m). By deﬁnition, one sees from Eq.
(9) that Φ+ (Θ, m)  Φ(Θ). Therefore, we can immediately
prove that Φ(Θ) = Φ+ (Θ, m)  Φ+ (Θ, m)  Φ(Θ).

u∈Rwhere ψω,t (u) is the Gabor function, which is a nonorthogonal basis used to measure the component of frequency ω attime t, and deﬁned as the product of the complex sinusoid withfrequency of ω and the Gaussian window centered at time t:2
ω02

ψω,t (u) = e− 2 (u−t) +jω(u−t) ,(6)

where ω0 is a time spread parameter of the Gaussian window,that can be chosen arbitrarily. Though trivial, we then have
 
 
2
2




−
s(t)
dt
=
η

y(t)
Gy (ω, t) − Gs (ω, t) dωdt,

It should be emphasized here that the EM algorithm canbe considered as a special case of this method.
One possible auxiliary function of Φ(Θ) can be derived
using the lemma suggested for example in [2].

where η is a constant.

Lemma 3. If a complex function mi (x) satisﬁes∗
i mi (x) = 1, then for x ∈ R(−∞, ∞)R2

R

If we now assume that θk (u) and Ak,n (u) are respectively
piecewise linear and piecewise constant, thenGs (ω, t) =

K 
N


−

Ak,n (t)e(ω−nμk (t))2
2
2ω0

2
 



y(x) −si (x)
 dx

R

.

i

(7)



k=1 n=1

2 1  mi (x)y(x) − si (x) dx, (12)


βi
i

In the case of discrete-time observations, the problem is tominimize2 (ω−nμk )2



−
2
2ω0
Y (ω) − dω, (8)AeΦ(Θ) =k,nRi mi (x) =R

where βi is a constant such thati βi = 1, βi ∈ (0, 1).
2

2

Putting Sk,n (ω)  Ak,n e−(ω−nμk ) /2ω0 for simplicity ofnotation, then by the Lemma 3 and from Eq. (8) we have thefollowing inequality:k,nat each discrete time point with respect to Ak,n and μk . For
clarity of notations the time index of Ak,n (t), μk (t) and Θ(t)is omitted and Y (ω) is a simpliﬁed notation of Gy (ω, t).

2
 





Φ(Θ) =
S
(ω)
Y
(ω)
−
 dω
k,n


R
k,n

2

 1



mk,n (ω)Y (ω) − Sk,n (ω) dω,βk,n R

3. PARAMETER OPTIMIZATION ALGORITHM
The parameter optimization algorithm we propose in thispaper is based on a principle called the auxiliary function
method. We ﬁrst deﬁne the auxiliary function and then show
how the iterative algorithm is performed.

(13)

k,n

where βk,n ∈ (0, 1),k,n βk,n = 1 and equality holds if1Sk,n (ω) .
mk,n (ω) =Sk,n (ω) + βk,n Y (ω) −Y (ω)

Deﬁnition 1. (Auxiliary function) Let Φ(Θ) be the objectivefunction that one wants to minimize with respect to the parameter Θ = (Θ1 , · · · , ΘI ). We then deﬁne Φ+ (Θ, m) as
the auxiliary function of Φ(Θ), and m = (m1 , · · · , mJ ) as
the auxiliary parameter if Φ+ (Θ, m) satisﬁes


 


 
Φ Θ = min Φ+ Θ, m ⇒ Φ Θ  Φ+ Θ, m . (9)

k,n

(14)
Let us denote by Φ+ (Θ, m) the right-hand side of Eq. (13).
By Deﬁnition 1, Φ+ (Θ, m) is an auxiliary function of theobjective Φ(Θ), and mk,n (ω) is an auxiliary parameter. By
Lemma 2, we consider next to minimize Φ+ (Θ, m) withm30

respectto Θ. Using√ the result 2of the Gaussian integral:
2
S
(ω)
dω
=
πω0 Ak,n  , one obtainsk,nR+Φ (Θ, m) =
+

√

πω0

 1
βk,n
k,n

−

− 2e



 Ak,n 2(ω−nμk )2
2
2ω0

k,n



R

βk,n


mk,n (ω)Y (ω)2
Re Ak,n m∗k,n (ω)Y ∗ (ω) dω.

(15)

However, one notices from Eq. (15) that one still cannot
obtain analytically the update equation for μk because of the22nonlinear part e−(ω−nμk ) /2ω0 in Eq. (15). One may want toderive another auxiliary function such that the update equation for μk can be obtained analytically. In order to derive
such an auxiliary function, we focused on two points: one isthat −e−x is a continuously differentiable concave function,and second is that we have the following theorem about continuously differentiable concave function.

Fig. 1. An illustration of the iterative algorithm
Setting to 0 the partial derivative of Φ+ (Θ, m, α) with respect
to μk , one obtains the F0 parameter update rule:


 n 
e−αk,n (ω) Re Ak,n m∗k,n (ω)Y ∗ (ω) ωdωβk,n Rn.
μk = 
n2−αk,n (ω)∗∗
e
Re Ak,n mk,n (ω)Y (ω) dω
βk,n Rn(20)
Setting to 0 the partial derivative of Φ+ (Θ, m, α) with respect
to A∗k,n , the update equation for Ak,n can also be derived as

Lemma 4. Let f (x) be a real function of x that is continuously differentiable and concave. Then, for any point α ∈ R,
f (x)  f (α) + (x − α)f  (α).

(16)

Since −e−x is a differentiable concave function, using
Lemma 4 we have the inequality −e−x  (x − α − 1)e−α ,2for any point α ∈ R. Replacing x with (ω − nμk ) /2ω02 andα with a real function αk,n (ω), then from Eq. (15),Φ+ (Θ, m) ≤

√

πω0



 Ak,n 2
k,n

βk,n

Ak,n = √

+

1
πω0


R

mk,n (ω)Y (ω)e−αk,n (ω)
2(ω − nμk )− αk,n (ω) − 1 dω.
2ω02

(21)

The presented algorithm is summarized as follows (Fig. 1).

 1  
mk,n (ω)Y (ω)2 +2Re Ak,n m∗k,n (ω)Y ∗ (ω)βk,n Rk,n
2−αk,n (ω) (ω − nμk )− αk,n (ω) − 1 dω. (17)e2ω02Step 0 Initial setting of {μk , {Ak,n }1≤n≤N }1≤k≤K .
Step 1 Update mk,n (ω) by Eq. (14).
Step 2 Update e−αk,n (ω) by Eq. (19).
Step 3 Update Ak,n by Eq. (21).
Step 4 Update μk by Eq. (20) and return to Step 1.

Denoting by Φ+ (Θ, m, α) the right-hand side of this inequation, Φ+ (Θ, m, α) can also be considered as an auxiliary
function of Φ(Θ) because

The algorithm described above generates the suboptimal estimates not only of μk and Ak,n but also of the separate signal
mk,n (ω)Y (ω). This fact shows that the algorithm performs
parameter estimation and source separation simultaneously.

Φ(Θ)  Φ+ (Θ, m)  Φ+ (Θ, m, α).

4. EXPERIMENTAL EVALUATION(18)

We ﬁrst compared the dependency on the initial parameterand the convergence speed of the gradient search method
and the proposed method. In this comparative experiment,
we used a synthetic signal as test data which was created byadding together two periodic signals (with F0 s of 200Hz and
270Hz) composed of 10 harmonic components, with eachcomponent’s amplitude and phase determined by random
generation. In the sinusoidal signal model, we set K = 2
and N = 10. A Gabor transform with diffusion parameterω0 = 0.366 was performed on this synthetic signal (16kHz
sampling frequency) to obtain Y (ω). The courses of the
update of the F0 parameters μ1 , μ2 at each step through the

In this case both mk,n (ω) and αk,n (ω) are the corresponding
auxiliary parameters. The equality Φ(Θ) = Φ+ (Θ, m, α)
holds when mk,n (ω) is given by Eq. (14) and αk,n (ω) by
αk,n (ω) =2(ω − nμk )
.
2ω02

(19)

The advantage worth mentioning of deriving this auxiliary
function is that it enables the analytical expression of the update equation for the F0 parameter μk , allowing us a complexspectrum-domain EM-like multiple F0 estimation algorithm.

31

330
320
310
300
290
280
270
260
250
240
230
220
210
300
290
280
270
260
250
240
230
220
210
200
190
180

(a)

(b)

(c)

0

10

20

30

40

50

iteration #

60

70

80

90

Fig. 4. Utterance by a female speaker (a), a male speaker (b)and their mixed signal (c).

100Fig. 2. Course of the update of μ1 , μ2 (the proposed method)
330
320
310
300
290
280
270
260
250
240
230
220
210
300
290
280
270
260
250
240
230
220
210
200
190
180

Fig. 5. Separated signals corresponding to the female (top)and the male speaker (bottom).

0

10

20

30

40

50

iteration #

60

70

80

90

sults of the separation of the mixed signal shown in Fig. 4
is shown in Fig. 5. After separation performed on the mixed
signal of the male speaker A and the female speaker B (witha SNR of -0.3dB seen from the male speaker A), the SNRs
for the speakers were 7.2dB and 6.4dB. On the mixed signal of the female speaker A and the female speaker B (with a
SNR of 1.5dB seen from the female speaker A), we obtained
the SNRs of 6.0dB and 4.8dB, and on the mixed signal ofthe male speaker A and the male speaker B (with a SNR of-0.3dB seen from the male speaker A), we obtained the SNRsof 4.8dB and 4.3dB after the separation was performed.
Although the method presented in this paper estimates theparameter independently for each frame, one may expect a
substantial reduction of the musical noise and an improvement of the SNR if a coordinated parameter estimation acrossseveral adjacent frames could be performed. This shall be oneof our future works.

100

Fig. 3. Course of the update of μ1 , μ2 (the gradient method)
proposed method and the gradient method, starting from various initial parameter conditions, are shown in Fig. 2 andFig. 3. The transitions of the update values of μ1 and μ2
corresponding to the same iterative computation are shown ineach ﬁgure respectively in the upper and lower part with thesame color and same line type. The initial value for the amplitude Ak,n was set to 0. One sees from Fig. 2 and Fig. 3 thatthe gradient method often gets trapped into stationary points
different from the true values for initial values of μ1 , μ2
which are not sufﬁciently close to the true values (200Hz,
270Hz), while the proposed method converges quickly fromany initial points in a large interval to the true values.
Next, we conﬁrm here the performance for single channel
source separation. We use the ATR B-set speech database tobuild the mixed signals by adding together the waveforms ofutterances from two male speakers, two female speakers, ora male speaker and a female speaker. For all the speech data
the sampling rate was 16kHz, and the frequency analysis wasdone using a Gabor transform with a frame interval of 10ms.
ω0 was set to 0.366 and N was set to 30. The initial values for μk were obtained by ﬁnding all the frequencies givinga minimum or a maximum of the real part or the imaginary
part of Y (ω), and selecting the frequencies with the 10 largest
powers. K was initially set to 10 and after the parameter converged, the source models with the 2 largest total-powers werechosen as the ﬁnal source signal estimates. In this experiment,in order to conﬁrm the basic source separation performance inthe situation where the permutation problem would be dealt
with, we determine to which source the separated signals correspond by looking at their proximity to each signal prior tothe mixing. Under the above conditions, an example of re-

5. SUMMARY
In this paper, we introduced the auxiliary function method
for parameter optimization of constrained sinusoidal model,
which enabled us to derive a complex-spectrum-domain EMlike multiple F0 estimation algorithm. We conﬁrmed the ability to avoid local solutions and the convergence speed of thepresented method and the performance on speech separation.
6. REFERENCES[1] R. J. McAulay and T. F. Quatieri, “Speech Analysis/Synthesis Basedon a Sinusoidal Representation,” IEEE Trans. Acoust., Speech, Signal
Process., Vol. ASSP-34, No. 4, pp. 744–754, 1986.
[2] M. Feder and E. Weinstein, “Parameter Estimation of Superimposed
Signals Using the EM Algorithm,” IEEE Trans. Acoust., Speech, Signal Process., Vol. ASSP-36, No. 4, pp. 477–489, 1988.
[3] D. Chazan, Y. Stettiner and D. Malah, “Optimal Multi-Pitch Estimation Using the EM Algorithm for Co-Channel Speech Separation,” In
Proc. ICASSP’93, Vol. 2, pp. 728–731, 1993.
[4] P. Jinachitra, “Constrained EM Estimates for Harmonic Source Separation,” In Proc. ICASSP2003, Vol. 6, pp. 609–612, 2003.
[5] S. Godsill and M. Davy, “Bayesian Harmonic Models for Musical
Pitch Estimation and Analysis,” In Proc. IEEE ICASSP2002, Vol. 2,
pp. 1769–1772, 2002.

32