1-Q-41
Vocal Tract Spectrogram Estimation withFormant Frequency Contour Factorization∗
Yunhan Zou1,2 , Li Li1,3 , Hirokazu Kameoka1
1
NTT Communication Science Laboratories
2
Georgia Institute of Technology 3 University of Tsukuba

1

Introduction

The spectral envelope of speech contains rich information about the voice characteristics of the
speaker. According to the source-filter model for
speech production, spectral envelopes correspond to
the resonance characteristics of the vocal tract. If
we can assume that the vocal tract spectra observed
at each time frame can be approximated as a linear sum of spectral templates scaled by time-varying
amplitudes, the vocal tract spectrogram, interpreted
as a non-negative matrix, can be approximated as
the product of two non-negative matrices, one containing the time-independent spectral templates arranged as column vectors and the other containing
the time-varying amplitudes arranged as row vectors. This way of representing a data matrix iscalled non-negative matrix factorization (NMF) andapplying it to a vocal tract spectrogram has several
attractive features.
First, it allows us to decouple vocal tract spectra
into time-independent and time-dependent factors,
namely the spectral templates and the temporal activations. This decomposition is noteworthy since
the former factor roughly describes the voice characteristics of the speaker whereas the latter contains
information about the transcription of the utteredsentence. Thus, if we convert the spectral templates
while keeping the temporal activations unchanged,
we can modify a speaker’s voice to sound like it werespoken by another speaker [1, 2]. This technique iscalled voice conversion.
It also provides a novel solution for tackling thesource-filter decomposition problem. The sourcefilter model assumes that a speech spectrum is given
as the product of the vocal tract and the glottal excitation spectra. Source-filter decomposition refers
to the problem of estimating the vocal tract andthe excitation source spectra solely from a speech
spectrum. Within a voiced segment where we assume the excitation source signal to be a periodic
pulse train, the spectrum of the produced speech is
given by the product of the vocal tract spectrum
and the equally spaced pulses with an interval equal
to the fundamental frequency (F0 ). This process
is shown in Fig. 1. Thus, a speech spectrum canbe seen as a sampled version of the vocal tract spectrum with missing information between the harmonics. Widely used vocoders such as STRAIGHT [3]and WORLD [4] provide ways of extracting a spectral envelope from a voiced spectrum by smoothlyinterpolating between the harmonics. However, it∗

Fig. 1

Source-filter model

is generally impossible to impute the missing data
and restore the true vocal tract spectrum only fromthe speech spectrum observed at a particular frame.
This illustrates the limitation of the frame-by-frame
analyzers that perform source-filter decompositionat each time frame independently. An attempt has
recently been made to overcome this limitation byleveraging multi-frame observations. The idea is toemploy a missing data imputation technique using
NMF to estimate the entire vocal tract spectrogram[5]. The underlying assumption behind this method
is that the vocal tract spectra are expected to be
represented fairly well as a linear sum of a limited
number of spectral templates. This expectation may
be supported by the fact that the number of vowel
types is usually limited in normal speech.
As these examples show, factorization of vocal
tract spectrogram has significant potential. In theprevious work, NMF provides a template-based representation of spectral magnitude. However, spectral magnitude is not the only important factor incharacterizing speech. As shown in Fig. 2, vocal
tract spectrograms are typically characterized byseveral dominant peaks varying continuously overtime. These peaks are often called formants. The
frequencies of these peaks correspond to the resonance frequencies of the vocal tract and the frequencies of the first two formants are known to be important in determining the quality of vowels. Fig. 3
shows an example of the spectrogram of Fig. 2 as aresult when NMF is applied to the decomposition ofspectral magnitude only. By looking at the flat formant trajectories in this example, we can confirm
that the NMF model is not well suited to express
the continuous trajectories of the formant frequencies. This is because the NMF model is only able toapply a linear sum of spectral templates in expressing a spectrum and does not have the ability to interpolate the peak frequencies of the templates. To
obtain a template-based representation better suited
to express vocal tract spectrogram, we propose yet
another model consisting of formant frequency set
templates.

フォルマント周波数行列分解による声道スペクトログラム推定, 鄒 雲漢 (NTT, ジョージア工科大), 李 莉 (NTT,筑波大), 亀岡 弘和 (NTT)

日本音響学会講演論文集- 323 -

2017年3月

Frequency (kHz)

8
7
6
5
4
3
2
1
0

0

0.5
Time (s)Fig. 4

1

GMM representation for spectral envelope.

Frequency (kHz)Fig. 2 Vocal tract spectrogram.
8
7
6
5
4
3
2
1
0

0

0.5
Time (s)1Fig. 3 NMF applied to the spectrogram of Fig. 2.

2

Model

2.1 Motivation
Each type of vowels is characterized by a formant frequency set. During the transition froma phoneme to another, one formant frequency set
varies continuously towards another formant frequency set. Since the phoneme number is limited, itwould be reasonable to assume that the formant frequency set at each time frame can be represented asa linear combination of formant frequency set templates. This idea implies that a matrix that contains
the formant frequency values at diﬀerent time frame,
which we hereafter call “formant frequency matrix”,
can be modeled as the product of two matrices. In
the following, we propose a novel vocal tract spectrogram model by incorporating this formant frequency matrix model. We also derive a convergenceguaranteed algorithm for estimating the parametersfrom an observed spectrogram.
2.2 Formant frequency matrix factorization
We start by describing a spectral envelope using a
Gaussian mixture model (GMM) [6, 7], interpreted
as a function of frequency (see Fig. 4 for its graphical illustration):F (ω, t) = β(t)I∑αi (t)Gi (ω, t),

(1)

i=1
2

(ω−µi (t))−
1
2σ 2 (t)ie,
Gi (ω, t) = √2πσi (t)(2)

where ω and t denote frequency and time, respectively, and I is the number of Gaussian mixture components. This representation consists of parameters(roughly) corresponding to the frequency and powerof spectral peaks and is particularly convenient forincorporating the above idea. µi (t), σi2 (t) and αi (t)日本音響学会講演論文集

are the mean, variance and weight of each Gaussiancomponent, corresponding to the frequency, peakiness and power of a dominant peak in the spectralenvelope. To eliminate an indeterminacy in the scalings of β(t) and αi (t), we assume αi (t) to satisfy asum-to-one constraint so that β(t) represents the energy of the spectrum in each frame.
If we assume that each Gaussian corresponds to aformant, µi (t) can be seen as the frequency of the ith formant at time t. As discussed in 2.1, we assumeµi (t) to be represented as a weighted sum of K formant frequency set templates wi,1 , . . . , wi,K . K canthus be thought of as the number of vowel types.
In the same way, we assume the magnitude ai,k to
be represented as a weighted sum of K magnitude
set templates ai,1 , . . . , ai,K . Here, we constrain theweight coeﬃcients hk (t) for wi,k and ai,k to be equal:∑µi (t) =wi,k hk (t),(3)
kαi (t) =∑ai,k hk (t).

(4)

k

The way these parameters are tied to each other isextremely important since we may want to associatea {wi,k , ai,k } pair with a particular shape of the ith prototype spectral envelope. Similarly, for thevariance we assume∑σi2 (t) =ci,k dk (t).
(5)
kSince formant frequencies usually vary continuously over time, we may want hk (t) to be a smoothfunction. One convenient way of ensuring smoothness is to express hi (t) as a weighted sum of halfoverlapping Hanning functions g1 (t), . . . , gJ (t) thatcover the entire time range:∑hk (t) =gj (t)lj,k ,(6)jand treat lj,k as the parameter to estimate insteadof hk (t). In this way, the smoothness of hk (t) canbe easily controlled by window size settings.
We place the first Gaussian function to center at
ω = 0 to account for the large energy present at
lower frequencies in the vocal tract spectrogram.
Namely, we set w1,k = 0.

3

Parameter estimation
We use Θ to denote the set of the parameters:- 324 -

W = {wi,k }i,k ,(7)

2017年3月

A = {ai,k }i,k ,L = {lj,k }j,k ,C = {ci,k }i,k ,D = {dk (t)}k,t ,β = {β(t)}t .

(8)
(9)
(10)
(11)
(12)

Note that all these values must be non-negative.
Here we derive a parameter estimation algorithm
that fits the present model F (ω, t) to an observed
spectrogram Y (ω, t). For convenience of the derivation of the optimization algorithm, we use theKullback-Leibler (KL) divergence (also known as
the I divergence) to measure the diﬀerence between
the present model and an observed spectrogram:I(Θ) =∫∫ (Y (ω, t)Y (ω, t) log
F (ω, t))− Y (ω, t) + F (ω, t) dωdt. (13)Note that the integral taken over the interval [0, ∞)of the first Gaussian centered at ω = 0 is 1/2 andthat of each of the other Gaussians is approximately
equal to 1. Thus, to make them consistent, we rede2
2
1 (t))}. Since theexp{− (ω−µfine G1 (ω, t) as √2πσ
2σ12 (t)1 (t)∑Hanning functions are set to satisfy j gj (t) = 1,
we can show that the integral of F (ω, t) taken overthe interval [0, ∞) is β(t) whenI∑ai,k = 1,

(14)

i=1

∑

lj,k = 1,(15)

Although it is diﬃcult to analytically find theglobal minimum point of the current optimization
problem, we can search for a stationary point using a majorization-minimization (MM) algorithm.
An MM algorithm refers to an iterative algorithm
that consists of iteratively minimizing an auxiliary
function called a “majorizer”. Suppose J (Θ) is anobjective function that we wish to minimize with respect to Θ. A majorizer J + (Θ, Ξ) is then defined asa function satisfying J (Θ) = minΞ J + (Θ, Ξ), whereΞ is called an auxiliary variable. An algorithm thatconsists of iteratively minimizing J + (Θ, Ξ) with respect to Θ and Ξ is guaranteed to converge to astationary point of the objective function. For our
objective function, we can show that
J + (Θ, Ξ) =[∫∫Y (ω, t) log Y (ω, t) − log β(t){∑∑ai,k gj (t)lj,k−λi (ω, t)
θi,k,j (t) log
θi,k,j (t)
i
k,j
(∑
)
1
k ci,k dk (t) − ζi (t)−+ log ζi (t)2ζi (t)()
1 ∑ ρ2i,k (t)−×2
ci,k dk (t)
k
(
2 )
2
∑ wi,k∑gj2 (t)lj,k2
wi,k gj (t)lj,k +ω − 2ω
γi,k,j (t)
k,j
k,j
}]
+ log λi (ω, t) dωdt + η∫+ β(t)dt + R(A),(19)

kis an auxiliary function where Ξ is the set of theauxiliary variables

are satisfied. (13) can therefore be rewritten as
I(Θ) =)∫∫ (Y (ω, t)Y (ω, t) log− Y (ω, t) dωdt
F (ω, t)∫+ β(t)dt, (16)under the redefined G1 (ω, t) and the conditions (14)and (15).
In addition, since we want each Gaussian to be associated with a dominant peak in a spectrum, we do
not want the magnitude of any of the Gaussians to
be extremely small. Thus, to prevent the magnitude
matrix ai,k from becoming sparse in the process ofmodel fitting, we incorporate a penalty term givenby the negative logarithm of a symmetric Dirichlet
distribution. A symmetric Dirichlet distribution is
maximized when all the arguments become exactly
equal. The penalty term is thus given as(∏)
R(A) = − log
ai,k φ .
(17)i,kHence, our objective function to be minimized isJ(Θ) = I(Θ) + R(A).

日本音響学会講演論文集(18)

λ = {λi (ω, t)}i,ω,t ,θ = {θi,k,j (t)}i,k,j,t ,ζ = {ζi (t)}i,t ,ρ = {ρi,k (t)}i,k,t ,γ = {γi,k,j (t)}i,k,j,t ,

(20)
(21)
(22)
(23)
(24)

∑
that must satisfy 0 ≤ λi∑(ω, t) ≤ 1,i λi (ω, t) =1, 0 ≤ θi,k,j (t) ≤ 1,θ(t)=1, 0 ≤
i,k,j
k,j
∑
ρ
(t)
≤
1,
ρ
(t)
=
1,
0
≤
γ
(t)
≤
1, and
i,k
i,k,j
k i,k
∑
γ
(t)
=
1.
The
updateequation
for
each
of
k,j i,k,j
the parameters and the auxiliary variables can beobtained in closed form through (19). The details
are omitted owing to space limitations.

4

Experiments

4.1 Experimental Conditions
To show the ability of the proposed method toexpress a vocal tract spectrogram with natural formant trajectories using a compact template-based
representation, we tested the proposed method on
the spectrogram obtained with the STRAIGHT
analysis.

- 325 -

2017年3月

Frequency (kHz)

8
7
6
5
4
3
2
1
0
0

0.1

Fig. 5Frequency (kHz)The STRAIGHT analysis was executed with 5[ms] time shift and the sampling frequency of thetest signal was 16 [kHz]. In the experiment, we used
K = 10 templates in our model. wi,k describing thei-th formant location in the k-th template was initialized by observing the formant locations from the
STRAIGHT spectrogram. In addition, it was initialized to include I = 5 Gaussians. The Hanning
function gj (t) was designed to half-overlap with eachcovering approximately 100 [ms] time span. We used
speech samples from the ATR speech database. The
parameter update process is divided into a total of 5
parts: during the initial fitting, basis matrix wi,k isunchanged while other parameters are updated for
70 iterations. This is done under the expectationthat with other parameters initialized to proper values, the subsequent update on wi,k , the key to thisexperiment, will therefore be more eﬃcient and lesserror-prone. Then follows the part where formant
location set wi,k joins the update along with all others for 120 iterations. In the last 3 parts, the number of Gaussians is further doubled (9, 17 and 33) tomake the fitting more accurate and each part runsfor 120 iterations.

8
7
6
5
4
3
2
1
0
0

0.1

0.2

0.3
Time (s)0.4

0.5

STRAIGHT spectrogram Y (ω, t).

0.2

0.3
Time (s)0.4

0.5

Fig. 6 F (ω, t) fitted to Y (ω, t) along with µi (t).
To confirm the model fitting accuracy, we plotted
the spectrograms of specific segments obtained withthe proposed method along with the estimated formant trajectories µi (t). These graphs can be foundin Figs. 5 and 6. We can see from these graphsthat the estimated formant trajectory is not oversmoothed and follows the formant transition well.
To verify the reliability of extracted formant locations, columns of matrix wi,k and ai,k were swappedin a certain way to obtain a new speech signal. By
listening to this speech signal, the speech was converted in an utterance with a completely diﬀerent
sentence while it still sounded natural enough as if
it was spoken by the same speaker. Modified spectrogram is plotted to show the same segments as in
the previous step, as in Fig. 7.
As the final step of evaluation, the distance between first and second formant locations in wi,k is
expanded. By listening to the generated speech, we
concluded that speech element sounds more distinguished from each other. The modified spectrogram
can be seen in Fig. 8.

5Frequency (kHz)Experimental Results

0.1

0.2

0.3
Time (s)0.4

0.5Reconstructed F (ω, t) obtained by swappingcolumns of W and A.

Conclusion

8
7
6
5
4
3
2
1
0
0

0.1

0.2

0.3
Time (s)0.4

0.5

Fig. 8 Reconstructed F (ω, t) with modified W .

To obtain a template-based representation well
suited to express vocal tract spectrograms, this paper proposed a vocal spectrogram model that allows
us to co-factorize the formant frequency/magnitude
matrices in the process of model fitting. We developed an update algorithm for the model parameters
based on an auxiliary function approach. We confirmed through experiments that the formant contour factorization is able to represent speech formant
trajectory, which resulted in natural-sounding synthetic speech.
Acknowledgment This work was supported by JSPS
KAKENHI Grant Numbers 26730100 and 26280060.

日本音響学会講演論文集

8
7
6
5
4
3
2
1
0
0

Fig. 7Frequency (kHz)

4.2

参考文献
[1] Takashima et al., IEICE Trans. IS, E96-A(10),1946–1953, 2013.
[2] Wu et al., IEEE/ACM Trans. ASLP, 22(10), 1506–
1521, 2014.
[3] Kawahara et al., Proc. ICASSP, 3933–3936, 2008.
[4] Morise et al., IEICE Trans. IS, E99-D(7), 1877–
1884, 2016.
[5] 中村, 亀岡, 音講論 (春), 3-P-33, 393–396, 2016.
[6] Kameoka et al., IEEE Trans. ASLP, 18(6), 1507–
1516, 2010.
[7] Hojo et al., Proc. SSW8, 129–134, 2013.

- 326 -

2017年3月