Trajectory Modeling based on HMMs with the Explicit
Relationship between Static and Dynamic Features
Keiichi Tokuda, Heiga Zen, Tadashi Kitamura
Department of Computer Science
Nagoya Institute of Technology, Japan{tokuda,zen,kitamura}@ics.nitech.ac.jp

Abstract
This paper shows that the HMM whose state output vector includes static and dynamic feature parameters can be reformulated as a trajectory model by imposing the explicit relationship between the static and dynamic features. The derived
model, named trajectory HMM, can alleviate the limitations of
HMMs: i) constant statistics within an HMM state and ii) independence assumption of state output probabilities. We also derive a Viterbi-type training algorithm for the trajectory HMM.
A preliminary speech recognition experiment based on N -bestrescoring demonstrates that the training algorithm can improve
the recognition performance significantly even though the trajectory HMM has the same parameterization as the standard
HMM.

1. Introduction
The HMM has successfully been applied to modeling the sequence of speech spectra in speech recognition systems. Its
tractable and efficient implementations are achieved by assumptions: i) constant statistics within an HMM state and ii) independence assumption of state output probabilities. To overcome
these limitations in the standard HMM framework, alternativemodels have been proposed, e.g., [1]–[10]. Although thesemodels can improve the speech recognition performance, theygenerally require an increase in model parameters and computational complexity. The use of the dynamic features (delta
and delta-delta features) [11] also improves the performance of
HMM-based speech recognizers. However, it has been thoughtof as an ad hoc rather than an essential solution.
This paper derives a trajectory model by reformulating thestandard HMM whose state output vector includes static anddynamic feature parameters. The standard HMM allows inconsistent static and dynamic features. By imposing the explicit
relationship between static features and dynamic features, thestandard HMM is naturally translated into a trajectory model,
referred to as “trajectory HMM” in this paper. We also derivea Viterbi-type training algorithm for the trajectory model, andevaluate it in a continuous speech recognition task.
The formulation of the trajectory HMM is closely related
to a technique for parameter generation from HMM [12]–[14],in which the speech parameter sequence is determined so as tomaximize its output probability for the HMM under the constraints between static and dynamic features. While we derived the speech parameter generation algorithm in order toconstruct HMM-based speech synthesizers [15] which can synthesize speech with various voice characteristics [16]–[18], thegeneration algorithm was also applied to speech recognition in[19]. The proposed training algorithm can be viewed as a train-

ing algorithm for the method in [19]. The model derived in[20] based on the maximum entropy framework is also closely
related to the proposed trajectory model.
The rest of this paper is organized as follows. Section 2 defines the trajectory model. In Section 3 the training algorithm
for the model is derived. Results of a continuous speech recognition experiment are shown in Section 4. Concluding remarks
and future plans are presented in the final section.

2. Reformulating HMM as trajectory model
The output probability of a speech parameter vector sequence



 
o = o
for the standard HMM is given by
1 , o2 , . . . , oTP (o | λ) =P (o | q, λ)P (q | λ),
(1)
all q
where q = {q1 , q2 , . . . , qT } is a state sequence. In mostof speech recognition systems, the speech parameter vector
ot is assumed to consist of the static feature vector ct =[ct (1), ct (2), . . . , ct (M )] (e.g., cepstral coefficients), and dynamic feature vectors ∆ct , ∆2 ct (e.g., delta and delta-delta2  
cepstral coefficients), that is ot = [ct , ∆ct , ∆ ct ] . The
dynamic features calculated by
(1)

L+

∆ct



=

w(1) (τ )ct+τ

(2)

w(2) (τ )ct+τ .

(3)

(1)
τ =−L−(2)∆2 c t

L+



=

(2)
τ =−L−

correspond to the first and second time-derivative of the staticfeature ct , respectively. Conditions (2) and (3) can be arranged
in a matrix form:o = W c,
(4)
where
c

=

W

=

wt

=

(n)
wt

=





c
1 , c2 , . . . , cT



(5)

[w 1 , w 2 , . . . , w T ]
(6)


(0)
(1)
(2)
(7)
wt , wt , wt(n)0M ×M , . . . , 0M ×M , w(n) (−L− )I M ×M ,
1st

(n)

(t−L− )-th(n)

. . . , w(n) (0)I M ×M , . . . , w(n) (L+ )I M ×M ,
t-th

(n)

(t+L+ )-th

0M ×M , . . . , 0M ×M



T -th

(0)

,

n = 0, 1, 2(0)

L− = L+ = 0, and w(0) (0) = 1. The above model is improper in the sense of statistical modeling: it allows inconsistentstatic and dynamic feature vector sequences even though they
are constrained by (4). The statistical model should be defined
as a function of c because the original observation is c rather
than the augmented variable o.
When each state output probability distribution is assumed
to be single Gaussian, P (o | q, λ) is given by
P (o | q, λ) =

T

t=1

N (ot | µqt , U qt ) = N (o | µq , U q ), (9)where µqt and U qt are the 3M × 1 mean vector and the 3M ×
3M covariance matrix, respectively, associated with qt -th state,
and




µq =
µ
(10)
q1 , µq2 , ..., µqT
Uq=

diag [U q1 , U q2 , ..., U qT ] .

3. Training Algorithm

(8)

(11)

In this section, we derive a training algorithm for the trajectorymodel. It should be noted that although the model has the same
parameterization as the standard HMM, the output probability
is defined by (18) rather than by (1). Accordingly, the model
parameters should be trained based on (18).
An auxiliary function of current parameter set λ and newparameter set λ is defined byQ(λ, λ ) =P (q | c, λ) log P (c, q | λ ).
(21)
all q

It can be shown that by substituting λ which maximizes
Q(λ, λ ) for λ, the likelihood increases unless λ is a criticalpoint of the likelihood. However, it is not tractable since wehave to evaluate all possible state sequences. To avoid this difficulty, we apply the single Viterbi path approximation1 . As
a result, the problem is broken down into the following maximization problems:
q max

=

arg max P (c, q | λ)

(22)

λ

=

arg max P (c, q max | λ)(23)

By substituting (4) for (9), we can rewrite (9) as follows:
P (W c | q, λ)==N (W c | µq , U q )Kq · N (c | cq , P q ),(13)

First, we solve the maximization problem of (23). The problem
is equivalent to maximizing1
log P (c | q, λ) = − M T log(2π) − log |Rq |2 + c P −1(24)
q c + rq P qrq − 2 rq c

(14)with respect to
m

and
Rq

W  U −1
q W

=



−1

rq

=

W U q µq

(15)

Pq

=

(16)

Kq

=

R−1
q

(2π)M T |P q |(2π)3M T |U q |

1
−1

µ
· exp −q U q µq − r q P q r q2. (17)

By omitting the normalization constant Kq in the above expression (12), a new statistical model can be defined:P (c | λ) =P (c | q, λ)P (q | λ),
(18)
all q
where

(19)
P (c | q, λ) = N (c | cq , P q ).
It is interesting to note that the mean cq is exactly the same as
the speech parameter trajectory obtained by the speech parameter generation technique (Case 1 in [14], see Appendix also),that is,cq = arg max P (o | q, λ) = arg max P (W c | q, λ). (20)c

λ3.1. Model parameter update(12)

where cq is given by
R q cq = r qqcThis means that by assuming the parameter trajectory cq as
the mean for the spectral parameter vector sequence c corresponding to an utterance, the standard HMM can naturally betranslated into a trajectory model: the state output probabilityof observing the static part of the output vector changes during
a state, and is affected by statistics of neighboring states. Notethat the spectral parameter vector sequence c is modeled by a
mixture of Gaussians whose dimensionality is T M , and theircovariances P q are generally full.

σ

=
=






µ
1 , µ2 , . . . , µN


−1 

−1
U −1
1 ,U2 ,...,UN

(25)
,

(26)

where N is the total number of HMM states. In this paper, we
refer to P (c | q, λ) as “trajectory likelihood.”
By setting ∂ log P (c | q, λ)/∂m = 0, we obtain a set of
linear equations

−1

S
q W P qW SqΣ m = Sq W c(27)

for determination of m which maximizes log P (c | q, λ),
where
Σ−1
µq

=
=

diag (σ)
Sqm

(28)
(29)

U −1
q

=

diag (S q σ) ,(30)

and S q is a 3T ×3M N matrix whose elements are 0 or 1 determined according to the state sequence q. The dimensionality of(27) is 3M N : although it could be tens of thousands, it is still
small enough to solve the set of linear equations using currently
available computational resources.
For maximizing log P (c | q, λ) with respect to σ, we apply
a steepest descent algorithm using the first derivative∂ log P (c | q, λ)∂σ

1 S q diag−1 W P q W 2− W cc W  + 2µq c W =+ W cq c
q W − 2µq cq W

because (24) is not a quadratic function of σ.
1 We can also use the N -best approximation.

(31)

c(1)
2.0
1.0
training data
HMM mean sequence
trajectory model (w/o parameter update)trajectory model (mean & variance update)

0.0
-1.0
0.1

0.2

0.4

sil

j

0.6

i b

u

N n

0.8

o

1.0

j

i

1.2

ts u ry o

1.4

k u w

a

1.6 (sec)pau

Figure 1: Generated trajectories cq for the trajectory HMM (thin line) and the standard HMM (broken line), the state mean sequenceµq (dotted line), and one of training data c (thick line).

Table 1: Average log “trajectory likelihood” P (c | q, λ).

training data
test data

w/o update
6.86
6.93

trajectory HMMm
σ
7.91 12.0
7.83 12.1

m& σ
14.8
14.8

Table 2: Recognition error rates (%) for the trajectory HMM
and the standard HMM.

error rate
rel. imp.

w/o update
19.4
0.5%

trajectory HMM
m
σ
18.5
19.0
5%
2.6%

m& σ
18.1
7%

standard
HMM
19.5
ref

3.2. Optimal path search
Next, we discuss the optimal path search problem of (22).
Based on the approximation
q max

=

arg max P (q, c | λ)=arg maxq

1
P (q, o | λ)Kqarg max P (q, o | λ),
q
q

(32)
(33)

we can use the Viterbi algorithm for the standard HMM. To reduce the inaccuracy caused by the approximation, the obtained
state boundaries can be adjusted by a recursive algorithm (case2 in [14]) based on the “trajectory likelihood.” By incorporating time-recursive calculation of the normalization factor Kq in
the Viterbi algorithm, we can obtain an optimal path directly bya modified Viterbi algorithm. This would be possible by using
the speech parameter generation algorithm of [21] which canoperates in a time-recursive manner with some look-ahead.
Noted that extensions of the training algorithm to the multiple observations and multi-mixture distributions are straightforward.
3.3. Decoding
For decoding, we can use i) rescoring scheme, ii) recursive calculation of normalization factor Kq in the Viterbi algorithm.
In approach i), N -best candidates are rescored by the “trajectory likelihood.” The state boundaries of each candidate canalso be adjusted by the recursive algorithm (case 2 in [14]). On
the other hand, in approach ii), Kq is calculated in the Viterbi
decoding using the time-recursive algorithm of [21] with some
look-ahead.

4. Experiment
A preliminary experiment was carried out. We used phonetically balanced 450 sentences from ATR Japanese speech
database for training speaker-dependent monophone HMMs
and their trajectory versions. Diagonal covariances were used
for both types of models. Speech signals were sampled at 16

kHz and windowed by a 25.6-ms Blackman window with a 10ms shift, and then mel-cepstral coefficients were obtained by amel-cepstral analysis technique [22]. The feature vector consists of 19 mel-cepstral coefficients including the zeroth coefficient as well as their delta and delta-delta coefficients. We
used 3-state left-to-right HMM structure with no skips. For agiven Viterbi path obtained by the standard HMM, parametersof the trajectory HMM were optimized by the algorithm described in Section 3.1 (i.e., the proposed training algorithm was
not iterated in this experiment). The trained trajectory HMM
was evaluated in a continuous phoneme recognition task basedon the N -best rescoring approach: first, a 10-best list was generated for each test utterance by using the conventional Viterbi
algorithm with the standard HMM, and then each candidate wasrescored by the trajectory HMM.
Table 1 shows average log likelihoods of the trajectory
HMMs (“w/o update”: model parameters were not updated, i.e.,
those of the standard HMM were used, “m”: only means were
updated, “σ”: only variances were updated, “m & σ”: bothmeans and variances were updated). It is shown that the training algorithm improves the “trajectory likelihoods” considerably for both training and test data sets.
Figure 1 shows the mean trajectories (the first coefficientin cq ) for the trajectory HMM (“w/o parameter update”: model
parameter of the standard HMM were used, “mean and varianceupdate”: means and variances were updated). The state mean
sequence of the standard HMM (the first coefficient in static partof µq ) and the trajectory extracted from one of the training data(the first coefficient in c) are also plotted in the figure. It canbe seen that the trajectory generated from the trained trajectory
HMM is closer to the trajectory of the training data. It should
be noted that the trajectory corresponding to a phoneme varies
according to its neighboring phonemes (see phoneme /j/ in the
figure). This means that the trajectory HMM has the capability
to capture the coarticulation effects naturally.
Table 2 shows the phoneme recognition error rates for the
trajectory HMM and the standard HMM. The trajectory HMM
achieves a relative error reduction of 7% over the standard
HMM. Given the fact that without parameter update the perfor-mance of the trajectory HMM did not improve, model parameter training based on the “trajectory likelihood” is essential in
the trajectory HMM approach.

[14] K. Tokuda, T. Yoshimura, T. Masuko, T. Kobayashi, T. Kitamura,“Speech parameter generation algorithms for HMM-based speechsynthesis,” Proc. ICASSP, vol.3, pp.1315–1318, June 2000.

5. Conclusion[15] T. Yoshimura, K. Tokuda, T. Masuko, T. Kobayashi and T. Kitamura, “Simultaneous modeling of spectrum, pitch and duration in
HMM-based speech synthesis,” Proc. EUROSPEECH, pp.2347–
2350, 1999.

This paper defined a new kind of trajectory model (trajectoryHMM) for acoustic modeling of speech by reformulating thestandard HMM whose observation vector consists of static anddynamic features. We also derived a Viterbi-type training algorithm for the trajectory model. It was shown that the model canwork as a trajectory model while it maintains the basic structureof the standard HMM. In a preliminary speech recognition experiment, a relative error rate reduction of 7% over the standard
HMM was achieved. One of the advantages of the trajectory
HMM is that huge amounts of software resources for the standard HMM could be easily reused.
In the near future, results of a complete speech recognition
experiment and application to Text-to-Speech synthesis will bepresented. Future work also includes derivation of a fast algorithm for solving (27) and a Baum-Welch-type training algorithm.

6. Acknowledgment
Authors would like to thank Prof. Takao Kobayashi and Dr.
Takashi Masuko for discussions.

7. References[1] M. J. F. Gales and S. J. Young, “Segmental Hidden Markov Models,” Proc. EUROSPEECH, pp.1579–1582, 1993.
[2] H. Gish and K. Ng, “Parametric trajectory models for speech
recognition,” Proc. ICSLP, pp.I-466–I-469, 1996.
[3] W. J. Holmes and M. J. Russell, “Speech recognition using a linear
dynamic segmental HMM,” Proc. ICASSP, pp.1611–1614, 1995.
[4] M. Ostendorf, V. Digalakis, and O. A. Kimball, “From HMMs tosegment models: a unified view of stochastic modeling for speech
recognition,” IEEE Trans. on Speech and Audio Processing, vol.4,
no.5, pp.360–378, 1996.
[5] L. Deng, M. Aksmanovic, X. Sun, J. Wu, “Speech recgnition using hiden Markov models with polynomial regression functions
as nonstationary states,” IEEE Trans. on Speech and Audio Processing, vol.2, no.4, pp.507–520, Oct. 1994.
[6] Y. Gong and J. P. Haton, “Stochastic trajectory modeling forspeech recognition,” Proc. ICASSP, vol.I, pp57–60, 1994.
[7] S. Takahashi, “Phoneme HMM’s constrained by frame correlations,” Proc. ICASSP, pp. 219-222, 1993.
[8] K. K. Paliwal, “Use of temporal correlation between successiveframes in hidden Markov model based Speech recognizer,” Proc.
ICASSP, pp.215-218, 1993.
[9] M. Ostendorf and S.Roukos, “A stochastic segment model forphoneme-based continuous speech recognition,” IEEE Trans. On
Acoustics, Speech and Signal Processing, vol.37, no.12, pp.1857–
1869, 1989.
[10] C. J. Wellekens, “Explicit correlation in hidden Markov model for
Speech Recognition,” Proc. ICASSP, pp.383-386, 1987.
[11] S. Furui, “Speaker independent isolated word rcognition using dynamic features of speech spectrum,” IEEE Trans. Acoust., Speech,
Signal Processing, vol.ASSP-34, pp.52–59, 1986.
[12] K. Tokuda, T. Kobayashi and S. Imai, “Speech parameter generation from HMM using dynamic features,” Proc. ICASSP, pp.660–
663, 1995.
[13] K. Tokuda, T. Masuko, T. Yamada, T. Kobayashi and S. Imai,“An algorithm for speech parameter generation from continuous
mixture HMMs with dynamic features,” Proc. EUROSPEECH,
pp.757–760, 1995.

[16] T. Yoshimura, K. Tokuda, T. Masuko, T. Kobayashi and T. Kitamura, “Speaker Interpolation in HMM-Based Speech Synthesis
System,” Proc. EUROSPEECH, pp.2523–2526, 1997.
[17] M. Tamura, T. Masuko, K. Tokuda and T. Kobayashi, “Adaptation of pitch and spectrum for HMM-based speech synthesis using
MLLR,” Proc. ICASSP, vol.2, pp.805–808, May 2001.
[18] K. Shichiri, A. Sawabe, K. Tokuda, T. Masuko, T. Kobayashi, T.
Kitamura, “Eigenvoices for HMM-based speech synthesis,” Proc.
ICSLP, pp.1269–1272, 2002.
[19] Y. Minami, E. McDermott, A. Nakamura, S. Katagiri, “A recognition method with parametric trajectory synthesized using direct
relations between static and dynamic feature vector time series,”
Proc. ICASSP, vol.I, pp.957–960, 2002.
[20] K. S. Van Horn, “A Maximum-entropy solution to the framedependency problem in speech recognition,” Tech. Rep., Dept. ofComputer Science, North Dakota State University, Nov. 2001.
[21] K. Koishida, K. Tokuda, T. Masuko and T. Kobayashi, “Vector
quantization of speech spectral parameters using statistics of dynamic features,” Proc. ICSP, pp.247–252, 1997.
[22] T. Fukada, K. Tokuda, T. Kobayashi and S. Imai, “An adaptivealgorithm for mel-cepstral analysis of speech,” in Proc. ICASSP,
1992, pp.137–140.
[23] A. Acero, “Formant analysis and synthesis using hidden Markov
models,” Proc. EUROSPEECH, pp.1047–1050, 1999.

A. Appendix
A.1. Maximizing P (o | q, λ) with respect to o under o =Wc
The logarithm of P (o | q, λ) can be written as1
 −1
log P (o | q, λ) = − o U −1
q o + o U q µq + Const, (34)2where Const is independent of o.
It is obvious that P (o | q, λ) is maximized when o = µq
without the condition (4), that is, the speech parameter vector sequence becomes a sequence of the mean vectors. With
the condition (4), maximizing P (o | q, λ) with respect to o isequivalent to that with respect to c. By setting∂ log P (W c | q, λ)
= 0,
∂c

(35)

we obtain a set of equations
Rq c = r q .

(36)
3

3

For direct solution of (36), we need O(T M ) operations because Rq is a T M × T M matrix. By utilizing the special
structure of Rq , (36) can be solved by the Cholesky decomposition or the QR decomposition with O(T M 3 L2 ) operations3 ,(n)where L = maxn∈{1,2},s∈{−,+} Ls . Equation (36) can also
be solved by an algorithm derived in [12]–[14], which can operate in a time-recursive manner [21].
22 When U
3q,i is diagonal, it is reduced to O(T M ) since each of the
M -dimensions can be calculated independently.
3 When U
2q,i is diagonal, it is reduced to O(T M L ). Furthermore,
(1)

(1)

when L− = −1, L+ = 0, and w(2) (i) ≡ 0, it is reduced toO(T M ) as described in [23].