COMPARISON OF AUTO-REGRESSIVE,
NON-STATIONARY EXCITED SIGNAL
PARAMETER ESTIMATION METHODS
1

Akira Sasou1 , Masataka Goto 1 , Satoru Hayamizu2, Kazuyo Tanaka3
National Institute of Advanced Industrial Science and Technology (AIST)
2
3{a-sasou,m.goto}@aist.go.jp

Faculty of Engineering, Gifu University
hayamizu@info.gifu-u.ac.jp

Institute of Library and Information Science, University of Tsukuba
ktanaka@ulis.ac.jp

Abstract. Previously, we proposed an Auto-Regressive Hidden
Markov Model (AR-HMM) and an accompanying parameter estimation method. An AR-HMM was obtained by combining an
AR process with an HMM introduced as a non-stationary excitation model. We demonstrated that the AR-HMM can accurately
estimate the characteristics of both articulatory systems and excitation signals from high-pitched speech. As the parameter estimation method iteratively executes learning processes of HMM
parameters, the proposed method was calculation-intensive. Here,
we propose two novel kinds of auto-regressive, non-stationary excited signal parameter estimation methods to reduce the amountof calculation required.

INTRODUCTION
The linear prediction (LP) method is widely used for the analysis of speech
signals [1, 2]. However, several problems remain to be resolved. For example, (1) local peaks of LP spectral estimates are strongly biased towardharmonics, especially for high-pitched speech [3], and (2) addition of white
noise to the Auto-Regressive (AR) process markedly alters the spectral estimate [4]. These phenomena result in deterioration of the perceived qualityof re-synthesized speech and can also cause speech recognition errors.
LP methods assume that the excitation signal conforms to an Identically
Independent Distributed (IID) normal distribution. However, the actual excitation signal indicates non-stationary properties especially in the case ofhigh fundamental frequency. As a result, local peaks in the LP spectralenvelope estimated from high-pitched speech are strongly biased toward harmonics. To correct this, we proposed an Auto-Regressive Hidden Markov

Model (AR-HMM) and an accompanying parameter estimation method [5]in which the HMM was introduced as a non-stationary excitation model.
We also demonstrated that the proposed method can accurately estimate
the characteristics of both articulatory systems and excitation signals fromhigh-pitched speech.
As learning processes of HMM parameters are executed iteratively, theparameter estimation method proposed in our previous study is calculationintensive. Here, we propose two novel auto-regressive, non-stationary excited
signal parameter estimation methods to reduce the amounts of calculation
required.

AUTO-REGRESSIVE HIDDEN MARKOV MODEL
Previously, we proposed an AR-HMM that was obtained by combining an AR
process with an HMM introduced as a non-stationary excitation model. Figure 1 shows examples of this AR-HMM. The output probability distributionof each node in the excitation HMM is assumed to be a single normal distribution. The nodes of the ﬁrst AR-HMM are concatenated in a ring state,
so the state transition occurs in order. Therefore, this type of AR-HMM
can be used to represent periodically excited signals. An ergodic HMM, as
shown in the lower part of Fig. 1, can be used to represent an aperiodically
excited signal. The AR-HMM can represent various types of signals throughappropriate design of the network topology. The number of nodes and theprediction order are determined according to the signal. Usually, an Akaike
Information Criterion (AIC) is employed to determine the model [6].

AR(p)

AR(p)

Figure 1: Examples of AR-HMMs.

ITERATIVE PARAMETER ESTIMATION METHOD USINGVITERBI ALIGNMENTThe AR-HMM parameters are the AR coeﬃcients and the parameters of the
HMM. Previously, we presented an algorithm that iteratively estimates theseparameters from a signal x(t), t = 0, · · ·, T −1 [5]. In the following, P denotes
the prediction order of the AR process. Let a(i) = [a(i) (1), · · · , a (i) (P )]T
represent the ith estimate of the AR coeﬃcients. The ith estimate of theexcitation signal e(i) (t), t = P, · · ·, T − 1 is given by:(i)eP = xP − Ωa(i)
where

(1)

(i)

eP = [e(i) (P ) e(i) (P + 1) · · · e(i) (T − 1)]T ∈ RT −P ,xt = [x(t) x(t + 1) · · · x(t + T − P − 1)] T ∈ RT −P ,
Ω = [xP −1 xP −2 · · · x 0 ] ∈ R (T −P )×P

We allocate a unique number from S = {1, · · ·, N } to each node of the excitation HMM to distinguish them from other nodes, where N is the number
(i)
2 (i)
of nodes. Let µn , σn , n ∈ S represent the ith estimates of the output
distribution population parameters in each node. Given a state-transition
sequence s(t) ∈ S, t = P, · · · , T − 1, the population parameters of an excita(i)2 (i)tion signal at time t are given by m (i) (t) = µs(t) , v (i) (t) = σs(t) . Hence, theexpectation vector of the excitation signal vector is represented by:(i)mP = [m(i) (P ) m(i) (P + 1) · · · m(i) (T − 1)]T

(2)

Based on the assumption that the samples of the excitation signal at diﬀerent
instants are mutually independent, the covariance matrix of the excitation
signal vector is deﬁned as a diagonal matrix given by:(i)ΣP = diag(v (i) (P ), v (i) (P + 1), · · · , v (i) (T − 1))(3)

The algorithm for parameter estimation consists of the following processes.
1. The initial population parameters of the excitation signal are prepared(0)(0)
as mP = 0, ΣP = I. Repeat the following processes from i = 0.
(i+1)

are estimated
2. The AR coeﬃcients a(i+1) and the excitation signal eP
(i+1)
(i)
(i)
by maximizing the likelihood given by L(e P ; mP , ΣP ).
(i+1)

(i+1)

, Σp
of the excitation signal vec3. The population parameters mp
(i+1)
(i+1)
tor are estimated by maximizing the likelihood given by L(e P ; mP ,
(i+1)
Σp
).
4. If the likelihood has converged, the algorithm stops. Otherwise, repeat
the above processes for i ← i + 1 from step 2.
By repeating the above processes, the likelihood increases almost monotonically in practical situations and converges to the optimum or to a local
optimum value.

The details of each step are as follows. In step 2, the AR coeﬃcient vector
can be obtained by solving

∂
(i)
(i) 
= 0.
log L(xP − Ωa; mP , ΣP )
∂a
(i+1)
a=a

The solution is given by:
(i)

(i)

(i)

a(i+1) = [ΩT (ΣP )−1 Ω]−1 ΩT (ΣP )−1 (xP − mP ).

(4)

(i+1)
The excitation signal vector eP
is derived from (1).

In step 3, the population parameters of the excitation signal vector are
estimated according to the following processes.
3.1 The Baum-Welch algorithm [7] estimates the population parameters(i+1)
2 (i+1)(i+1)µm , σm, m ∈ S of each output distribution using eP .
3.2 The Viterbi algorithm [8] estimates a state transition sequence s(t), t =P, P + 1, · · · , T − 1.
(i+1)

and the diagonal covariance matrix
3.3 The expectation vector mP(i+1)
ΣPof the excitation signal vector are estimated using (2) and (3).
PARAMETER ESTIMATION METHOD BASED ON
THE EXPECTATION-MAXIMIZATION ALGORITHM
The method described in the previous section estimates AR coeﬃcients and
HMM parameters separately and iteratively; HMM parameter estimation
requires the largest amount of calculation in the algorithm. When the method
is applied to real-time processing, it is necessary to reduce the number of
HMM parameter estimation iterations. To do this, we adopt an expectationmaximization (EM) algorithm. That is, using an EM algorithm, the AR
coeﬃcient estimation process can be embedded into the HMM parameter
estimation process, and all the parameters can be estimated during a single
HMM parameter estimation.
In the following, we describe the EM-based parameter estimation method
for generalized situations where several signals are available for learning andthe output distribution of each node is a Gaussian mixture. Let xn (t), n =1, . . . , N, t = 0, . . . , T n − 1 and a(k), k = 1, . . ., P represent observed signals
and AR coeﬃcients, respectively. The excitation signal en (t) emitted fromthe HMM is given by:en (t) = xn (t) +P


a(k)xn (t − k).

(5)

k=1

The output probability distribution os (e) of the node s is a Gaussian mixture
given by:os (e) =

M

m=1

2
λs,m N (e; µs,m , σs,m
),

M

m=1

λs,m = 1.

(6)

If we let πs and qs1,s2 represent initial probabilities and transition probabilities, respectively, AR-HMM parameter θ is given by θ = {a, π, q, λ, µ, σ 2}.
When observed signals xn (t), state-transition sequence sn (t) and Gaussiandistribution sequence mn (t) are given, the likelihood L(θ) of parameter θ isgiven by:L(θ | xn , sn , mn ) = πsn (0)

T
n −1

qsn(t),sn (t+1) λsn (t),mn (t)t=0

×N (en (t); µsn (t),mn (t) , σs2n (t),mn (t) ).

(7)

In the above equation, because the variables sn (t) and mn (t) are actuallyunobservable, it is necessary to estimate parameter θ from incomplete dataxn (t). This can be achieved using an EM algorithm, which consists of two
steps as follows.
• E-Step
We represent the probability of the unobservable variables sn , mn using the currently estimated parameter θ, then evaluate the expectation
Q(θ̂ | θ) of logarithmic likelihood log(L( θ̂)) with respect to the unobservable variables, where θ̂ represents the updated parameter. This
expectation is given according to the following equation:M

N T −1

Q(θ̂ | θ) =

n
 
1  
γn (t, i, j, m) log(L(θ̂))
N n=1 t=0m=1

(8)i∈S j∈S

where the parameter γn (t, i, j, m) represents the probability that themth Gaussian of the ith node emitted the excitation signal en (t) and atransition from the ith node to the jth node occurred. This parameteris evaluated by:γn (t, i, j, m) =2α(i, t − 1)qi,j λi,m N (en (t); µi,m , σi,m)β(j, t)L(θ | xn )(9)where α and β are obtained by applying a Forward-Backward algorithm
to the excitation signal.
• M-Step
In the M-step, the expectation Q(θ̂ | θ) is maximized with respect to θ̂.
The updated AR coeﬃcients â(k) are given as a solution of ∂Q/∂â(k) =0, k = 1, . . ., P . The solution is represented by the following equation:
ck,l =MN T
n −1   
γn (t, i, j, m)

n=1 t=0 i∈S j∈S m=1
2×xn (t − k)xn (t − l)/σi,m

dk =

M
N T
n −1   
n=1 t=0 i∈S j∈S m=1

γn (t, i, j, m)

(10)

2
×{xn (t) − µi,m }xn (t − l)/σi,mT[â(1), · · · , â(P )] = −C−1 dwhere C = (ck,l )k,l=1,···,P and d = [d1 , · · · , d P ]T .
The updated HMM parameters are given according to the following
equations:
M
N 
j∈S
n=1
m=1 γn (0, s, j, m)
π̂s = N 
(11)

M
i∈S
j∈S
n=1
m=1 γn (0, i, j, m)Nq̂s1 ,s2 = N

n=1

n=1

Tn −1 M
t=0

m=1 γn (t, s1 , s2 , m)Mj∈S
m=1 γn (t, s1 , j, m)

Tn −1 
t=0(12)

Tn −1 
n=1t=0j∈S γn (t, s, j, m)
λ̂s,m = N Tn −1 
M

j∈S
n=1
t=0
m =1 γn (t, s, j, m )(13)µ̂s,m =

Tn −1 j∈S γn (t, s, j, m)en (t)n=1t=0N Tn −1 j∈S γn (t, s, j, m)
n=1
t=0

(14)

2
j∈S γn (t, s, j, m) (en (t) − µs,m )t=0
N Tn −1 j∈S γn (t, s, j, m)
n=1
t=0

(15)

N

N

N
2
=
σ̂s,m

n=1

Tn −1 

The above two steps are iterated while replacing the current parameter θwith the updated parameter θ̂ until the likelihood converges.
POPULATION PARAMETER ESTIMATION BASED ON SAMPLE AVERAGE AND SAMPLE VARIANCE
The method described above can reduce the HMM parameter estimation iterations by adopting the EM algorithm. However, even if the HMM parameter
estimation is executed only once, the estimation process still requires a large
number of calculations. In this section, we propose a simpliﬁed method that
does not use the HMM to estimate the population parameter of the excitation
signal. In this method, the AR coeﬃcients and the population parametersof the excitation signal are estimated separately and iteratively as in themethod described in section 3, but the population parameters are evaluatedwith only the sample average and sample variance of the excitation signal.
The processes of this method are the same as processes 1 to 4 describedin section 3. The diﬀerence is that the excitation signal vector population
(i+1)
(i+1)
parameters mp, Σp
are estimated as follows. First, the sample aver(i+1)
age m(t) and sample variance v(i+1) (t) are estimated with the following
equations:Ts1m(i+1) (t) =e(i+1) (t + k)
(16)
2Ts + 1k=−Ts

v(i+1) (t) =

Ts 
2
1 
e(i+1) (t + k) − m(i+1) (t)
2Ts

(17)

k=−Ts

where e(i+1) (t) is the (i + 1)th estimate of the excitation signal obtained from
(1) and Ts is a parameter that determines the number of samples used for(i+1)
sample average and variance evaluations. The population parameters mp
,
(i+1)
Σp
are then obtained by substituting the evaluated sample average andvariance into (2) and (3).

EXPERIMENTS
We conducted an experiment with synthetic speech to compare the estimation accuracy of the three proposed methods with the conventional autocorrelation linear prediction (ALP) method. In the following, we refer to the
proposed methods described in Sections 3 to 5 as Method-1, Method-2, andMethod-3, respectively. Speech was synthesized using the excitation signalsof an impulse train and AR coeﬃcients of order 16 extracted from a male’svowel /a/. The excitation signals were generated by adding white Gaussian
noise N (0, 0.1) to impulse trains of amplitude 50 in fundamental frequencies
ranging from 100 Hz to 900 Hz. The sampling frequency was set to 16 kHz,the analysis frame length was set to 30 ms, and the prediction order was set
to 16. The excitation HMM used for Method-1 and Method-2 consisted of
two nodes: one representing an impulse, and another representing a noise segment. Both methods iterate the respective algorithm until the following conθ̂))−log(L(θ))
dition is satisﬁed during 5 successive iterations: | log(L(log(L(θ))| < 1.0e−7.
The parameter Ts used for Method-3 was set to 5. The estimated AR coeﬃcients were transformed to LPC Mel-Cepstrum coeﬃcients, which play avery important role in applications such as speech recognition. The estimation accuracy of each method was then evaluated by the Euclidean distance
between the original LPC Mel-Cepstrum coeﬃcients and those obtained bythis method.
Figure 2 shows the vocal tract spectra estimated by Method-1, Method-2,
Method-3, and the ALP method. The original vocal tract spectrum that wasused for synthesizing speech signals is shown at the top of each ﬁgure. Figure
3 shows the evaluated Mel-Cepstrum distances of all of the methods. As
shown in these ﬁgures, the estimation accuracy of the ALP tended to degrade
as the fundamental frequency increased. In contrast, Method-1 can extractspeech features in the widest range of fundamental frequency; its estimation
accuracy was also the highest of all of the methods examined. The estimation
accuracy of Method-2 up to a fundamental frequency of 550 Hz was almost
equivalent to that of Method-1. In the fundamental frequency range above
600 Hz, however, the estimation accuracy of Method-2 deteriorated markedly
as compared to the other methods. In Method-3, the estimation accuracy
tended to asymptotically approach that of the ALP when the fundamental
frequency exceeded 600 Hz. However, the estimation accuracy was very closeto those of Method-1 and Method-2 up to a fundamental frequency of 550
Hz. Figure 4 shows the processing times required by each method given as
values relative to the processing time of Method-1.
These results indicated that all of the proposed methods are capable of extracting features when analyzing speech signals with fundamental frequencies
below 550 Hz. In addition, Method-3 can be applied to real-time processing
because it does not use the HMM to estimate the population parameter of theexcitation signal. However, Method-1 is still the best method if maximum
estimation accuracy is necessary and there is suﬃcient time.

CONCLUSIONS
Here, we proposed two novel auto-regressive, non-stationary excited signal
parameter estimation methods, in addition to the method proposed previously. We are currently planning to construct a singing voice recognition
system and apply the proposed method to feature extraction at the frontend of the system.

REFERENCES[1] F.Itakura and S.Saito, “A statistical method for estimation of speech spectral
density and formant frequencies,” Electronics and Communications in Japan,
Vol.53-A, No.1, pp.36-43, January 1970.
[2] B.S.Atal and S.L.Hanauer, “Speech analysis and synthesis by linear predictionof the speech wave,” J. Acoust. Soc. Am., Vol.50, pp.637-644, 1971.
[3] J.Makhoul, “Linear Prediction: A Tutorial Review,” in Proc.of IEEE, Vol.63,
No.4, pp.561-580, April 1975.
[4] S.M.Kay, “The Eﬀects of Noise on the Autoregressive Spectral Estimator,”
IEEE ASSP-27, No.5, pp.478-485, Oct. 1979.
[5] A.Sasou, K.Tanaka, “Glottal Source Modeling Using HMM and Robust Analysis of High Fundamental Frequency Speech,” in Trans. IEICE, Vol. J84-D-II,
No.9, pp.1960-1969, Sep. 2001(in Japanese).
[6] H.Akaike, “A new-look at the statistical model identiﬁcation,”
IEEE Trans. Autom. Control, Vol.AC-19, No.6, pp.716-723, 1974.
[7] L.E.Baum, “An inequality and associated maximization technique in statistical estimation of probabilistic function of a Markov process,” Inequalities,
Vol.3, pp.1-8, 1972.
[8] A.J.Viterbi, “Error bounds for convolutional codes and an asymptotically optimum decoding algorithm,” IEEE Trans. on Information Theory, Vol.IT-13,
No.2, pp.260-269, April 1967.

100

100

150

150

200

200

250

250

300

300

Fundamental Frequency[Hz]Original

Fundamental Frequency[Hz]Original

350
400
450
500
550
600

350
400
450
500
550
600

650

650

700

700

750

750

800

800

850

850

900

900

0

1000

2000

3000

4000
5000
Frequency[Hz]

6000

7000

8000

0

1000

2000

(a)Method-1

6000

7000

8000

6000

7000

8000

Original

100

100

150

150

200

200

250

250

300

300

350

Fundamental Frequency[Hz]Fundamental Frequency[Hz]

4000
5000
Frequency[Hz](b)Method-2

Original

400
450
500
550
600
650

350
400
450
500
550
600
650

700

700

750

750

800

800

850

850

900

900

0

3000

1000

2000

3000

4000
5000
Frequency[Hz](c)Method-3

6000

7000

8000

0

1000

2000

3000

4000
5000
Frequency[Hz](d)Auto-correlation LP

Figure 2: Vocal tract spectra estimated by (a)Method-1, (b)Method-2, (c)Method-3
and (d)Auto-Correlation LP

LPC Mel-Cepstrum Distance

4.5

4

Method-1
Method-2
Method-3
Auto-Correlation LP

3.5

3

2.5

2

1.5

1

0.5

0
100

200

300

400

500

600

700

800

Fundamental Frequency[Hz]Figure 3: Evaluated Mel-Cepstrum Distances

10

1

1.0

0.1

2.195e-2
0.01

0.001

3.605e-48.373e-50.0001

1e-05

Method-1

Method-2

Method-3

Figure 4: Relative processing times

ALP

900