Proc. ICASSP, pp. 5066-5069, 2010.

NON-NEGATIVE TEMPORAL DECOMPOSITION OF SPEECH PARAMETERS
Sadao Hiroya
NTT Communication Science Laboratories, NTT Corporation
3-1 Morinosato-Wakamiya, Atsugi-shi, Kanagawa, 243-0198 Japan

1ABSTRACTWe present a non-negative temporal decomposition method for line
spectrum pair and articulatory parameters. Based on the multiplicative update rules derived from a non-negative matrix factorization
algorithm, these parameters decompose into a set of temporally overlapped event functions that are restricted to the range [0, 1] and corresponding event vectors. With the proposed method, the RMS errorof the measured and estimated articulatory parameters is 0.16 mm
and the spectral distance of the measured and estimated line spectrum pair parameters is 1.97 dB. These results also show that theseestimation errors of proposed method are significantly smaller thanthose of the conventional method. This technique will be useful formany applications, such as speech coding and speech modification.

k-1

k

k+1

k0
tk-1

tkIndex Terms— Non-negative matrix factorization, Temporaldecomposition, articulatory parameters

tk+1Fig. 1. Example of event functions.

1. INTRODUCTION
In speech processing, it is important to efficiently represent multidimensional time-varying speech parameters, such as line spectrum
pair (LSP) and articulatory parameters. A hidden Markov model(HMM) [1] is one of the models that can represent dynamical behavior as well as the trajectory smoothness of speech parameters.
Temporal decomposition (TD) [2] can represent speech parameters
as a set of temporally overlapped event functions and corresponding event vectors. TD has been used for many applications: speechcoding [2, 3], segmentation of speech signals [4, 5], analysis of articulatory parameters [6, 7] and modification for the speech spectrum[8] as well as for the speaking rhythm [9]. For these purposes, theevent functions should be restricted to the range [0, 1]. This restriction contributes to reducing the quantization error in speech codingand modifying speech parameters easily. However, they have beenclipped at the range [0, 1] because it is difficult to apply the restriction to the conventional TD algorithm directly. Thus, a novel TDalgorithm with the restriction is required.
In this paper, we present a non-negative temporal decomposition (NTD) method to overcome the problem. For given speech parameters, this method can optimize the event functions, which arerestricted to the range [0, 1], and the event vectors using the multiplicative update rules derived from a non-negative matrix factorization (NMF) [10] algorithm.
2. CONVENTIONAL TEMPORAL DECOMPOSITION
The original temporal decomposition [2] approximates the ithspeech parameter yi (t) of time t toŷi (t) =m∑k=1

where ai,k is the kth event vector, φk (t) is the kth event function,p is the dimension of the speech parameter, T is the length of theparameter sequence and m is the number of event function. Shiraki[11] has assumed that φk (t) is zero for t < tk−1 and t > tk+1 (seeFig. 1). Thus, ŷi (t) can be represented as a linear combination ofai,k and ai,k−1 :ŷi (t) = ai,k φk (t) + ai,k−1 φk−1 (t),1 ≤ t ≤ T, 1 ≤ i ≤ p(1)(2)

Kim [3] has proposed that the sum of all event functions is one at
any time t:ŷi (t) = ai,k φk (t) + ai,k−1 (1 − φk (t)).

(3)

The event functions φk (t) are determined by minimizing the leastsquares error between yi (t) and ŷi (t):p∑φk (t) =(ai,k − ai,k−1 )(yi (t) − ai,k−1 )i=1p∑(ai,k − ai,k−1 )2.

(4)

i=1

However, the obtained event functions are not restricted to the nonnegative value, nor the range [0, 1]. Thus, Kim clipped event functions at the range [0, 1] [i.e. φk (t) = min(1, max(0, φk (t)))] andthen updated event vectors ai,k using the least-squares method asfollows.
m∑ai,k φk (t),tk−1 ≤ t ≤ tk .

k=1ai,kT∑
t=1

φk (t)φj (t) =T∑
t=1yi (t)φj (t),

1≤j≤m

(5)

p
∑

φk (t) ←ai,k yi (t) + α

i=1p∑(ai,k−1 ai,k φk−1 (t) + a2i,k φk (t)) + α(φk−1 (t) + φk (t))φk (t)(7)

φk−1 (t)

(8)

i=1p∑

φk−1 (t) ←ai,k−1 yi (t) + α

i=1p∑(ai,k−1 ai,k φk (t) + a2i,k−1 φk−1 (t)) + α(φk−1 (t) + φk (t))
i=1∑

tk+1ai,k ←yi (t)φk (t)t=tk−1tk −1∑∑ai,ktk+1(ai,k−1 φk−1 (t)φk (t) + ai,k φ2k (t)) +t=tk−1

t=tk

3. NON-NEGATIVE TEMPORAL DECOMPOSITION
It is difficult to determine the event functions that are restricted
to the range [0, 1] by the least-squares method without a clipping,but the non-negative matrix factorization (NMF) [10, 12] algorithm
can iteratively determine the non-negative value of event functions. Moreover, if we assume the model in Eq. (2) for yi (t) andφk (t)+φk−1 (t) = 1 at any time t, we can determine the event functions that are restricted to the range [0, 1]. The problem is defined
astkpm∑∑∑min(yi (t) − ai,k φk (t) − ai,k−1 φk−1 (t))2k=2 t=tk−1 i=1that minimizes total distortion for the whole interval 1 ≤ t ≤ Tis derived efficiently by utilizing the dynamic programming (DP)method [11]. That is, we have(D(tk ) =tkpm∑∑∑(yi (t) − ai,k φk (t) − ai,k−1 φk−1 (t))2k=2 t=tk−1 i=1(φk (t) + φk−1 (t) − 1)2 ,(6)k=2 t=tk−1

where 1 = t1 < t2 < · · · < tm = T and α is the weight. In linewith a previous idea [12], we can obtain the multiplicative update
rule for the event function and vector [Eqs. (7)-(9)].
4. OPTIMIZATION OF EVENT TIMINGThe distortion d(y(t), ŷ(t)) of the cost function for each interval
tk−1 ≤ t ≤ tk only depends on time tk−1 and tk . Therefore, theevent timing
tk = arg

mint2 ,··· ,tm−1T∑
t=1

d(y(t), ŷ(t)),

2≤k ≤m−1

tk−1 ∈Rk−1

D(tk−1 ) +tk∑)d(y(t), ŷ(t))(11)

t=tk−1

Rk−1 = {t|tk−1 − δ ≤ t ≤ tk−1 + δ},(12)

D(tk ) is an accumulated minimal distortion at tk and δ is a searchrange.
5. ALGORITHMS

Thus, we consider minimizing the following cost function by NMF
algorithm:tkm∑∑

minwheres.t. φk (t) ≥ 0, φk−1 (t) ≥ 0, ai,k ≥ 0, ai,k−1 ≥ 0φk (t) + φk−1 (t) = 1 for ∀t, i, k+α(9)(ai,k+1 φk (t)φk+1 (t) + ai,k φ2k (t))(10)

5.1. Restricted TD (RTD)
In RTD [3], we first set the initial event vectors as ai,k = yi (tk )for the initial event timing tk . Then, the event functions and theevent timings are determined using Eq. (4) and DP. Finally, the eventfunctions are clipped at the range [0, 1] and only the event vectors areupdated by Eq. (5) for given event functions and event timings. Notethat Kim [3] has proposed a method for holding an ordering propertyof LSP parameters. But this is irrelevant to the main subject, thus wedon’t use this in the paper.
5.2. NTDIn NTD, the initial event vectors and event timings are the same as
RTD. Then, the event functions and the event timings are determinedby DP and the NMF update rules of Eqs. (7)-(8). Note that the eventvectors are not updated at this stage. Finally, for given event timings,
the event functions and vectors are updated using Eqs. (7)-(9).
6. EXPERIMENTAL CONDITIONS
We compare the proposed method (NTD) with the conventional
method (RTD). We evaluate the proposed method using LSP andarticulatory parameters. Articulatory parameters and speech signal

2.3
Spectral distance (dB)

2.25
Mic.
16 kHz
UL UI
T1

LL

32-ms Hamming window

LI

T2 T3

T

250 Hz

2.2
2.15
2.1

RTD

2.05

NTD

2
1.95

16 LSP coefficients

EMA system

Vocal-tract spectrum

Articulatory positions

1.9verticalpositions

0

Articulatory-acoustic pair data

1234 5 6 7
Number of delta

8910

horizontal positions

0.3

0.25

0.2
RTDNTD

0.15

0.1
012345678910

Number of delta

Fig. 4. Spectral distance of LSP parameters for the number of δ’s.
Vertical bars indicate the standard error of the mean.
Cross-correlation coef. of event functions

RMS error of articulatory parameters (mm)Fig. 2. Simultaneous observations to obtain articulatory-acousticdata.

1
0.98
0.96
0.94
articulatory
0.92

LSP

0.90.880.86
0-1 1-2 2-3 3-4 4-5 5-6 6-7 7-8 8-9 9-10

successive delta

Fig. 3. RMS error of articulatory parameters for the number of δ’s.
Vertical bars indicate the standard error of the mean.

Fig. 5. Cross-correlation coefficients of event functions. Verticalbars indicate the standard error of the mean.

data were obtained from simultaneous observations using an electromagnetic articulographic (EMA) system [1] and audio recording ofcontinuous speech utterances (Fig. 2). The articulatory data werecollected using the EMA at a sampling rate of 250 Hz. The articulatory parameters were represented by the vertical and horizontal
positions of six receiver coils, which were placed on the lower incisor (LI), the upper and lower lips (UL, LL), and the tongue (T1,T2, T3; three positions). The speech signal was recorded at a sampling rate of 16 kHz. Sixteen LSP parameter coefficients without
the 0-th coefficient were obtained using a 32-ms Hamming window
with a 4-ms frame shift.
In articulatory-acoustic recordings, one Japanese male subject
read 16 Japanese sentences. The average numbers of phonemes andframes per sentence were 37 and 748, respectively. We labeled thecentral point of each phoneme on the time axis, which we call initial
event timing tk . Thus, m is the number of phonemes. The labelingwas done manually by an expert. The α was set to 106 and 103 forarticulatory and LSP parameters, respectively.

was smaller than that of RTD. Figure 4 shows the spectral distanceof LSP parameters for the number of δ’s. For every δ, the spectral
distance of NTD was smaller than that of RTD. For both articulatoryand LSP parameters, the estimation error using NTD monotonicallydecreased for the number of δ’s. For NTD, φk (t) + φk−1 (t) at any
time t was one. This indicates that the proposed method efficiently
reduces the estimation error under the condition that the event functions are restricted to the range [0, 1].
The estimation error of articulatory parameters for NTD (0.16mm) was much smaller than 1.22 mm for an articulatory HMM [13]and 1.65 mm for a kinematic triphone model [14] because NTD parameters have not been statistically modeled in this study. But weexpect that the error of NTD with statistical modeling will be equivalent to or smaller than that of other articulatory modeling by using
redundancy in NTD parameters.
For the application of NTD to speech coding and speech modification, better stability of the event function would contribute to reducing the quantization error and improving the quality of modifiedspeech. To investigate this, we assessed the cross-correlation coefficients of event functions for the successive δ’s for articulatory and
LSP parameters (Fig. 5). This was done by concatenating the event
functions as [φ1 , φ3 , φ5 , · · · ]: we calculated the cross-correlation ofsolid lines in Fig. 6(d) and (e). The coefficient of the articulatory

7. RESULTS
Figure 3 shows the estimation error of articulatory parameters using
RTD and NTD for the number of δ’s. For every δ, the error of NTD

(a) Speech
articulatory parameters(b) Vertical positions ofarticulatory parameters.
9. ACKNOWLEDGEMENTS

LIUL

The author thanks Drs. H. Gomi and K. Kameoka for many usefuland helpful discussions.

LL
T1
T2

10. REFERENCES
10 mm

T3[1] Hiroya, S. and Honda, M., “Estimation of Articulatory Movements from Speech Acoustics Using an HMM-Based Speech
Production Model,” IEEE Trans. Speech Audio Process., pp.
175–185, 2004.

(c) The first six LSP parameters

1.4
1.2

[2] Atal, B.S., “Efficient Coding of LPC Parameters by Temporal
Decomposition,” in ICASSP, 1983, pp. 81–84.

10.8[3] Kim, S-J. and Oh, Y-H., “Efficient Quantisation Method for
LSF Parameters Based on Restricted Temporal Decomposition,” Electronics Letters, pp. 962–964, 1999.

0.6
0.4

[4] Bimbot, F., Chollet, G., Deleglise, P., and Montacie, C.,“Temporal Decomposition and Acoustic-phonetic Decoding ofSpeech,” in ICASSP, 1988, pp. 445–448.

(d) Event function ofarticulatory parameters(e) Event function ofLSP parameters

0.201[5] Dix, P.J. and Bloothooft, G., “A Breakpoint Analysis Procedure Based on Temporal Decomposition,” IEEE Trans. Speech
Audio Process., pp. 9–17, 1994.

0.8
0.6
0.4

[6] Jung, T-P., Krishnamurthy, A.K., Ahalt, S.C., Beckman, M.E.,
and Lee, S-H., “Deriving Gestural Scores from ArticulatoryMovement Records Using Weighted Temporal Decomposition,” IEEE Trans. Speech Audio Process., pp. 2–18, 1996.

0.2
010.8
0.6
0.4
0.2050100

150200250300350400450500550frameFig. 6. (a) Speech waveform. (b) Measured (thin lines) and estimated (thick lines) articulatory parameters of vertical positions. (c)Measured (thin lines) and estimated (thick lines) first six LSP parameters. (d) Event functions of articulatory parameters. (e) Eventfunctions of LSP parameters. δ = 5.

parameters is almost 0.99, but that of LSP is smaller. Moreover, wefound the event functions of LSP parameters were largely changedfor δ of more than 7. This indicated that the stability of event functions of articulatory parameters using NTD is better than that of LSPparameters. One possible reason is that the temporal patterns of thearticulatory parameters are simpler and smoother than those of theLSP parameters, but a more detailed analysis is required.
Finally, the cross-correlation coefficient of event functions between LSP (Fig. 6(e)) and articulatory parameters (Fig. 6(d)) using
NTD was calculated. The coefficients were 0.65, 0.53, and 0.36 forδ = 0, 5, and 10, respectively. The lag was around 3 msec.
8. CONCLUSIONSThe paper presented a non-negative temporal decomposition methodfor speech parameters. The error of NTD is smaller than that ofRTD under the condition that the event functions are restricted to therange [0, 1]. Using NTD, we plan to analyze speaker variability of[7] Collins, M.J., Krishnamurthy, A.K., and Ahalt. S.C., “Generating Gestural Scores from Articulatory Data Using Temporal
Decomposition,” IEEE Trans. Speech Audio Process., pp. 230–233, 1999.
[8] Nguyen, B.H. and Akagi, M., “A Flexible Spectral Modification Method Based on Temporal Decomposition and Gaussian
Mixture Model,” Acoust. Sci. & Tech., pp. 170–179, 2009.
[9] Shiraki, Y., “Optimal Temporal Decomposition for Voice Morphing Preserving ∆ Cepstrum,” IEICE Trans. Fundamentals,
pp. 577–583, 2004.
[10] Lee, D.D. and Seung, H.S., “Learning the Parts of Objects
by Non-negative Matrix Factorization,” Nature, pp. 788–791,
1999.
[11] Shiraki, Y. and Honda, M., “Extraction of Temporal Pattern ofSpectral Sequence Based on Minimum Distortion Criterion,”in ASJ Spring Meeting, 1991, pp. 233–234.
[12] Virtanen, T., “Monaural Sound Source Separation by Nonnegative Matrix Factorization With Temporal Continuity andSparseness Criteria,” IEEE Trans. Audio Speech and Lang.
Process., pp. 1066–1074, 2007.
[13] Hiroya, S. and Mochida, T., “Multi-speaker Articulatory Trajectory Formation Based on Speaker-independent Articulatory
HMMs,” Speech Commun., pp. 1677–1690, 2006.
[14] Okadome, T. and Honda, M., “Generation of Articulatory
Movements by Using a Kinematic Triphone Model,” J. Acoust.
Soc. Am., pp. 453–463, 2001.