IIEEE TRANSACTIONS ON SIGNAL PROCESSING, VOL. 39, NO. 2, FEBRUARY 1991

41 1

Discrete All-Pole ModelingAmro El-Jaroudi, Member, IEEE, and John Makhoul, Fellow, IEEE

Abstract-A new method is introduced for parametric modeling ofspectral envelopes when only a discrete set of spectral points is given.
This method, which we call discrete all-pole (DAP) modeling, uses a
discrete version of the Itakura-Saito distortion measure as its error
criterion. One result is a new autocorrelation matching condition that
overcomes the limitations of linear prediction and produces better fitting spectral envelopes for spectra that are representable by a relatively small discrete set of values, such as in voiced speech.
We present an iterative algorithm for DAP modeling that is shownto converge to a unique global minimum. We also present results ofapplying DAP modeling to real and synthetic speech. DAP modeling is
extended to allow frequency-dependent weighting of the error measure, so that spectral accuracy can be enhanced in certain frequency
regions relative to others.

1. INTRODUCTIONPOPULAR approach to speech modeling is to use a timevarying all-pole spectral shaping filter with a variable excitation source. This all-pole filter is typically estimated using
linear prediction methods [ 111. It has been known for some timethat linear prediction (LP) suffers from drawbacks that are especially manifested during voiced segments of speech. Specifically, the peaks of LP spectral estimates during these segments
are highly biased towards the pitch harmonics, especially forhigh-pitched sounds and voices [ 111. It has been pointed outthat the drawbacks of LP are inherent to its error criterion [ 1 11.
To overcome these drawbacks, we propose in this paper a new
all-pole modeling method based on a discrete form of the Itakura-Saito distance measure.
The new method, which we call discrete all-pole (DAP) modeling, overcomes the well-known limitations of LP and generally gives better all-pole spectral envelopes that are less biased
towards the pitch harmonics. In DAP, we approximate the spectrum of voiced speech, which has its energy located approximately at the harmonics of the fundamental pitch, by a line
(discrete) spectrum. The problem of finding the spectral envelope is then reduced to fitting an all-pole spectrum to a finite setof spectral points so as to minimize the discrete form of the
Itakura-Saito (I-S) distance measure. We present an algorithm
to compute the optimal envelopes and show that DAP modelingproduces generally better fitting spectral envelopes than LP.
Another advantage of DAP modeling is that the method canbe extended to allow error weighting as a function of frequency.
This weighting can be used to emphasize the spectral fit in sections of the spectrum at the expense of the fit at other sections.
This method has known applications in various aspects of speechprocessing where better spectral accuracy is often desired atlower frequencies relative to high frequencies.

A

In the next section, we present the known limitations of linear
prediction. In Section 111, we introduce DAP modeling and derive the minimization conditions for the discrete form of the I-S
distance measure. We also examine the possible forms of solutions one obtains using DAP modeling. In Section IV, wederive an iterative algorithm to compute the optimal all-polemodel and discuss its convergence properties. The results ofapplying DAP modeling to various signals are presented in Section V. In Section VI, we modify DAP modeling to allow
weighting of the error measure as a function of frequency, andthe results of weighted DAP are compared to those without
weighting.
OF LINEARPREDICTIONS
11. LIMITATIONSIn this section, we review linear prediction spectral analysis
and demonstrate, using examples, some of its inherent limitations.
The basic concept of linear prediction (LP) is to predict, subject to an error criterion, the present value of a signal based on
its previousp values, wherep is the prediction order. Normally,
the error criterion used is a least squares distance measure between the actual and predicted values. In the frequency domain,for discrete spectra defined at a set of frequencies U , E Q, theLP error criterion is equivalent to minimizing [ 121where P( 0 ) is the spectrum of th,e given signal, N is the numberof discrete frequencies U,, and P ( U ) is the spectrum of the allpole envelope defined as

Note that, in the definition of P ( w ) , we have incorporated thegain factor in the coefficients of the denominator, i.e., a. is not
restricted to 1. We emphasize here that the frequencies U , inthe set Q, which includes both positive and negative frequencies, can be arbitrary and need not be equally spaced.
The minimization of ELpwith respect to the predictor coefficients, ak, 1 5 k Ip , yields a well-known set of normal
equations:P

C R ( i - k)ak = 0,k=O

15 i 5p(3)

andManuscript received November 21, 1988; revised December 11, 1989.
This work was supported in part by the Department of Defense.
A. El-Jaroudi is with the Department of Electrical Engineering, University of Pittsburgh, Pittsburgh, PA 15261.
J. Makhoul is with Bolt Beranek and Newman Inc., Cambridge, MA
02138.
IEEE Log Number 9041 149.

P

1

R(k)Uk = k=O1053-587X/91/0200-0411$01.00
0 1991 IEEE

(4)

a0

(5)

412

IEEE TRANSACTIONS ON SIGNAL PROCESSING, VOL. 39, NO. 2, FEBRUARY 1991

is the autocorrelation of the discrete signal spectrum P ( U,).
Equation (4) was obtained by substituting (3) in (1) and setting
the energy of the continuous LP spectrum equal to the energy
R ( 0 ) of the given discrete spectrum.
It is well known [ 111 that by minimizing ELP,we are matching the autocorrelation of the continuous LP envelope to that of
the given discrete spectrum. In other words, we are effectively
setting

where

and R ( i ) is the autocorrelation of the discrete spectrum in (5).
In the following example, we show a typical behavior of LP
spectral analysis of discrete or harmonic spectra. We will use
LP to attempt to recover the original all-pole envelope of a periodic signal. The spectral envelope shown in solid line in'Fig.
l(a) represents a synthesis all-pole ( p = 12) filter which, when
excited with a periodic pulse sequence having a discrete spectrum with N = 30 frequency points, will give a periodic signal
whose line spectrum is shown in the figure. If we perform LP
analysis on the periodic signal, we obtain the dashed envelope
shown in Fig. l(a). Clearly, the LP envelope does not matchthe original envelope nor does it fit the discrete spectrum. While
there is a unique all-pole envelope (the original) that perfectly
fits the discrete spectrum, LP fails to recover that all-pole envelope. It has been shown [12] that for discrete spectra the LP
error measure (1) possessss an error cancellation property whichmakes it select an envelope other than the only one which passes
through all the spectral points.
We now show that it would be unreasonable to expect LP torecover the original envelope from the discrete spectral samples. We define Rorgto be the autocorrelation corresponding tothe original all-pole filter with spectrum P( U ) . Rorg( i ) andP ( w ) are related by the standard transform pair

i

1

and
m

C

P ( O ) = I = - m R , ~ ~ (e-]"'.
I)(9)

The autocorrelation R corresponding to the discrete samples ofthe synthesis envelope is defined in (5). By substituting (9) in(3,we find that R is related to R,, by

l

N

c

R(i) =R o r g ( I )e - J u t f 1 ( ' -,' )N m = !I=-m

for all i.

(IO)

This equation shows the aliasing that occurs in the autocorrelation domain whenever a spectral envelope is sampled at a discrete set of frequencies. For the periodic excitation case, thefrequencies w, will be equally spaced at U, = 27r(m - l ) / N ,and (10) reduces tomR(i) =[= -m

R o r g ( i- I N ) ,for all i.

(11)

For the example above, we show, in Figs. l(b) and (c), Rorg(i )and R ( i ) for lags 0 Ii 5 75 and excitation period N = 30.
Note that R ( i ) as obtained by aliasing Rorg( i ) as given by (1 1).

Fig. 1. (a) Example of the limitations of linear prediction spectral analysis. The solid line is the original 12-pole envelope which goes through all
the points. The dashed line is the 12-pole LP model for N = 30 spectral
lines. (b) Rorg( i ) for lags 0 5 i 5 75 and excitation period N = 30. (c)Autocorrelation R ( i ), 0 c i 5 75, corresponding to the discrete spectrum
0 5 i I75, corresponding toin Fig. l(a). (d) Autocorrelation RLP(i),the LP envelope in Fig. I(a). ( e ) Discrete frequency sampled impulse response h( i ), -37 Ii I37 corresponding to the discrete spectrum in Fig.
I@).

I

EL-JAROUDI AND MAKHOUL: DISCRETE ALL-POLE MODELING413

As shown in ( 6 ) , LP matches the autocorrelation of the continuous model spectrum to that of the given spectrum. By applying this matching to the example shown above, we obtainIn other words, since the autocorrelation corresponding to the
LP envelope will always equal an aliased version of Rorg(forthe discrete spectrum case), the LP envelope will not equal theoriginal envelope. It is also important to note that LP producesa unique all-pole model given a set of autocorrelations, whichmeans that the original all pole is not a possible solution to (3).
The autocorrelation R,, ( i ), 0 Ii s 75, corresponding to the
LP envelope in the example above is shown in Fig. l(d). We
note that RLp(i ) in Fig. l(d) is equal to R( i ) in Fig. l(c) for
the first 13 lags (since p = 12) and, very importantly, that
RLP( i ) is very different from the original autocorrelation Rorg
shown in Fig. l(b). The LP error criterion matches the autocorrelation of the continuous all-pole model to the autocorrelation of the given signal without taking into account the aliasingthat has occurred in the discrete spectrum. This interpretationof LP spectral estimation explains the sensitivity of LP estimates to high-pitched sounds. As the pitch increases, we havefewer and fewer harmonics (spectral samples) and the autocorrelation aliasing becomes more and more severe, which in turn,leads to worse LP models.
In Fig. 2(a), we show the same example as in Fig. 1 for N =30 and nonperiodic excitation. The discrete spectral lines shownin the figure are not equidistant and are not multiples of somefundamental frequency. The original envelope is shown in solidline and the LP envelope in dashed. Once again, LP fails torecover the original envelope. It is easy to understand why the
LP estimate was so poor by examining the autocorrelations corresponding to the original envelope and the discrete spectrum.
These autocorrelations are shown in Fig. 2(b) and (c), respectively. We note how very different the aliased autocorrelation
is from the original one. Consequently, the LP envelope whose
autocorrelation is shown in Fig. 2(d) and is equal to the aliased
autocorrelation for the first 13 lags, will also be very different
from the original one. These examples show that LP is thewrong approach to envelope estimation for discrete spectra since
it does not account for the aliasing caused by spectral sampling(periodic or nonperiodic).
This disregard to autocorrelation aliasing is demonstrated
again in the following one-pole example. Given a single polefilter, with a real pole at z = p , excited by a periodic sequencewith period N , one can show [2] that the resulting signal will
have a normalized first autocorrelation lag R ( 1 ), given by

When we perform LP analysis on the periodic signal, the resulting one-pole LP envelope will have a pole at z = pLp =R ( l ) = ( p N - ' + p)/(lp N ) . F o r p = 0.95 and N = 10, forexample, we have pLp = 0.99. Clearly, the LP envelope will
have a much narrower bandwidth than the original one. We notefrom (13) that, as N gets larger, R( l ) , and hence pLp,will approach p , which as N + 03 corresponds to the continuous spectrum case. However, for many applications, especially for highpitched speech, as N decreases, LP autocorrelation matching+

Fig. 2. (a) Example of the limitations of linear prediction spectral analysis
for nonperiodic but discrete spectra. The solid line is the original 12-poleenvelope; t& dashed line is the 12-pole LP model for N = 30 spectral
lines. (b) R o c g ( i) for lags 0 5 i 5 75 corresponding to the solid envelopein Fig. 2(a). (c) Autocorrelation R ( i ), 0 5 i I75, corresponding to thediscrete spectrum-in Fig. 2(a) with N = 30 nonharmonic spectral lines. (d)Autocorrelation RLP(i), 0 5 i I75, corresponding to the LP envelope in
Fig. 2(a).

often produces suboptimal and inadequate results. This matching is forced upon us by the error measure used for LP.
To improve upon the LP estimate, researchers have devised
methods with either a different error criterion or with added
constraints to regular LP (e.g., [4], 151, [8], [9], [15]-1171,

IEEE TRANSACTIONS ON SIGNAL PROCESSING, VOL. 39, NO. 2, FEBRUARY 1991414[20]). These methods can be divided into time domain techniques and frequency domain techniques. In his time domain
approach, Lee [9] weights the residual error so that time segments that are well predicted influence the all-pole estimate
more than the other segments which contain large residual values due to excitation pulses. A similar but less general method
was proposed earlier by Mizoguchi [15]. In another time domain approach, Rose and Clements [19] use a residual peakiness measure and select the all-pole model that maximizes it.
Also, there are some all-pole estimation methods [81, [ 161, 1201,[ 171, especially in coding applications, where, given the excitation, one selects the predictor that best matches the synthetic
to the original speech.
In his frequency domain approach, Hermansky [5] attempts
to generate more frequency samples of the original envelope byinterpolating between the measured harmonic peaks, and thenfits an all pole to the new set of frequency points using LP. Ina different method, Hermansky [4] improves the LP fit by firsttransforming the measured spectrum to modify its dynamic
range, then fits an LP envelope to the transformed spectrum. In
another approach, McAulay [ 141 makes a Gaussian assumption
about the observed time sequence then uses a maximum likelihood approach and some approximations to estimate the parameters of an all-pole filter given the samples of the discrete
spectrum.
The method we have selected employs the discrete form of
the Itakura-Saito distance measure. The reasons for this choice
and the details of this modeling technique are presented below.

111. DISCRETE
ALL-POLEMODELINGIn this section, we introduce a new method for all-pole envelope estimation of discrete spectra, which we call discrete allpole (DAP) modeling. We present the error criterion used in theestimation and derive its minimization conditions. We then derive an expression for the minimum error and show that there
is a unique model spectrum that achieves it.

A. Properties of the Error Measure
The Itakura-Saito (I-S) error measure was defined originallyfor continuous spectra [6], [7]. However, it can be adapted tothe discrete case as follows:

the spectral flatness of the error spectrum P ( w , ) / p ( w , ) , where
the spectral flatness is defined as the geometric mean of thespectral samples divided by their arithmetic mean. The proof ofthis property parallels the one for the continuous case given in[ 131. The major consequence of this property is that our optimalmodel is the one which makes the residual (error) spectrum as
flat as possible.
It has been shown that, for small values of EIS,the I-S error
approximates the mean-squared distance between the log spectra [6]. Based on this property, we shall define
EdB= 6.142(15)

for small EIS.
In this paper, we shall use EdB when comparing error values
since it provides an approximate estimate of the spectral error
in decibels.
It is important to note that, for the continuous case of the
Itakura-Saito measure, the optimal all-pole model is the sameas the one produced by LP. Therefore, by using this error measure, we do not sacrifice any of the advantages or performanceof LP in unvoiced segments of speech.
B. Error Minimization

For pedagogical reasons, we shall first minimize the error
measure in (14) with P ( o ) expressed ask=O

where { d k } can be shown to be equal to
P

do =

c a:

k=O

P-'

di = 2c akak+ir 1

Ii 5 p .

k=O

(19)

Note that d , is twice the autocorrelation of { a k } for 1 Ii 5
p , and do is equal to the zero-lag autocorrelation. We then set
a E , , / a d i = 0 for i = 0, *, p . The result can be shown toyield a set of correlation matching conditions, given by-

where P ( w , ) is the given discrete spectrum defined at N frequency points a, E Q, and P ( w , ) is the all-pole model spectrum, defined in (2), evaluated at the same frequencies. Thiserror measure is always nonnegative and is equal to zero only
when P(w,) = P(w,) for all w, E Q (i.e., P(w,) = P ( w , )gives a minimum for E,, but, as shown in the previous section,
not necessarily for E L p ) .
The continuous form of this error measure was originally presented as part of a maximum likelihood approach to linear prediction and was shown to produce the same result as LP for
continuous spectra. Later, the discrete version shown in (14)was derived by McAulay for the maximum likelihood spectral
modeling of periodic speech signals with Gaussian statistics~41.
A spectral flatness interpretation of this discrete error measure makes it a very reasonable choice for the problem of fitting
an envelope to a set of discrete spectral values. It can be shownthat minimizing the error in (14) is equivalent to maximizingk ( i ) = R(i),0 Ii 5 p(20)

where R ( i ) is the autocorrelation corresponding to the given
discrete spectrum defined in (5) and R ( i ) is the autocorrelation
corresponding to the all-pole model sampled at the same discrete frequencies as the given spectrum.

N

Equation (20) looks deceptively similar to the usual LP autocorrelation matching condition in (6). The major difference,
however, is that in LP, RLpJi ) is the autocorrelation of the continuous all-pole spectrum P( a),while here, R ( i ) in (21) is the
autocorrelation of a discrete sampling of the all-pole spectrum.
From (20), we see that DAP requires matching the given aliased
autocorrelation to the autocorrelation of the all-pole model
aliased in the same manner. It is this improved correlation

I

415

EL-JAROUDI AND MAKHOUL: DISCRETE ALL-POLE MODELING

matching condition, which incorporates the autocorrelation aliasing, that makes DAP better suited than LP for analyzing
voiced segments of speech, and discrete spectra in general.
While the matching condition (20) provides us with insight
into the solution of the modeling problem, it does not give us away of computing the parameters of the optimal all-pole model.
The all-pole model is obtained by using the definition of P ( CO)in (2) and setting a E , , / a a , = 0, i = 0 , * . . , p. This yields
the following set of equations relating the predictor coefficients{ ak} to the autocorrelations of the given discrete spectrum andthe sampled all-pole modelP2

c u k [ R ( i- k ) - k ( i - k ) ]0 5 i 5 p.

= 0,

k=O

(22)

The conditions in (22) can be expressed in matrix notation as2(R - h ) a = 0

(23)

or
Ra = Ra(24)

whey U is the column vector of predictor coeficients, and R
and R are symmetric Toeplitz matrices with elemfnts R ( i - j )and R ( i - j ), 0 5 i , j 5 p, respectively. Since R is a function
1 nonlinear equations in p f
of a, (23) constitutes a set of p
1 unknowns. Before we present a solution to these equations,
we derive the minimum error.

+

C. Minimum Error
The minimum error is obtained by substituting the condition
for minimization (24) in (14). We begin by simplifying the firstterm in the error measure (14):PP

c -P(w,)- - , = o C a,u,R(k - j )
N
-

m=l

j=o

= aTRa= a’Ra.

We obtain (27) by substituting (24) into (26). But
Pa’Ra =

P

1 c -c c u,ujR(k - j ) = -Nm=lP(w,)
P(mm) -

k=O j=O

Therefore, at the minimum, we have

Consequently, from (14) we haveof the ratio of the geometric means of the model spectrum andthe given spectrum. Both these properties have their equivalent
in continuous spectrum LP [ 111, Based on the similarities between the two methods (DAP and LP) and the‘fact that DAP
reduces to LP for the continuous spectrum case while LP does
not reduce to DAP for the discrete spectrum case, we conclude
that LP is just a special case of DAP modeling where the number of spectral points N goes to infinity.

D. The Solution and its UniquenessWe now focus on the solution of the minimization conditionsin (23). These equations allow one of two possible solutions:
1) A “matching” solution in which R = R and the modelsatisfies the conditions in (20).
2) A “singular” solution in which# R, and therefore,
the predictor vector a will be an eigenvector of the difference
matrix ( R - R ) corresponding to an eigenvalue equal to 0.
(Note that the trivial solution a, = 0 is not possible since it
produces unbounded values for R.)
Consequently, the optimal all-pole model will belong to oneof these two classes of solutions. It will either have an aliased
autocorrelation equal to that of the given discrete signal (matching solution), or it will not (singular solution). In Appendix A,
we show that the error function is convex and that the optimal
all-pole model, independent of which class of solutions it belongs to, will be unique under some mild conditions on thenumber of spectral points N . We also show that if the optimalmodel is the singular solution (i.e., it does not satisfy thematching conditions in (20)), it will have some roots on the unit
circle. This statement implies that there are cases where thegiven autocorrelation R(i ) cannot be generated (matched) usinga sampled all-pole spectrum, and consequently, the optimal allpole model will be unstable. It is important to point out thatthese cases rarely occur for real signals and, when they do occur, there are remedies. Below, we will present some examplesof the two types of solutions for all-pole models of orders 1 and
2.

E. Examples
In the examples below, the given autocorrelations are assumed to correspond to harmonic spectra (i.e., the frequency
points are equidistant). This simplifying assumption reduces themathematical complexity of the derivations without affecting
any of the conclusions. The results can be generalized to arbitrary discrete spectra. In these examples, we show that the rangeof autocorrelations generated by periodically exciting an all-polefilter is smaller than the one spanned by autocorrelations corresponding to the given general discrete spectra. (In contrast,in the continuous case, every given autocorrelation sequence of
1 can be generated with an all-pole of order p ) .
length p1) One-Pole Models: It can be shown [2] that the normalized first autocorrelation lag R ( 1) of an arbitrary harmonic
spectrum with N spectral lines has the range

+

Lm= I

L

J

_I

We conclude from (29) and (30) that, at, the minimum, the energy in the residual spectrum P ( w,)/P( CO,) is automaticallynormalized to 1 and the minimum error is equal to the logarithmcoss(N - 1 ) 5 R ( l ) 5 1~

forNodd(31)

forNeven.

(32)

N-1 5 R(l) 5 1

However, as shown earlier in (13), for a single-pole model witha pole at z = p , sampled at N harmonic frequency points, thep )/normalized first autocorrelation lag R ( 1 ) equals ( p - I+

IEEE TRANSACTIONS ON SIGNAL PROCESSING, VOL. 39, NO. 2, FEBRUARY 1991

416

(1 + p N ) . k ( 1 ) has the range
-1

L F
T
___---+ -N2 5 k ( l ) 5 1

forNodd(33)-1 5 k(1) 5 1

forNeven.

(34)

We note from (32) and (34) that, for N even, all possible valuesof R ( 1) can be generated (matched) using a one-pole filter.
Therefore, for N even, we will always have a matching solution. However, if N is odd, there may be cases where we may
not be able to find an all-pole that satisfies the correlation
matching condition (20). For example, for N = 7, we see from(31) that for an arbitrary spectrum we must have R ( 1 ) >-0.901, while ffom (33) we see that for a sampled single-polemodel we haveR( 1) > -0.714. Therefore, any value o f R ( 1 )that is in the range -0.901 < R ( l ) < -0.714 cannot bematched by a sampled one-pole model. In this range, the optimal all-pole will have a pole at z = - 1 and will be unstable.
2 ) Two-Pole Models: We now give an example of the possible solutions to (22) for second-order filters. It can be shown[2] that for the first two normalized autocorrelation lags R ( 1 ),R ( 2 ) corresponding to an N-point harmonic spectrum arebounded by a polygon with ( N / 2 ) + 1 sides for N even and( N + 1)/2 sides for N odd. The vertices of the polygon lie on
the parabola enclosing the autocorrelations corresponding to arbitrary continuous positive-definite spectra. Each vertex corresponds to a discrete spectrum consisting of a nonzero spectral
value at only one of the harmonic frequencies. An example ofthese regions is given in Fig. 3(a) for N = 11. The dashedpolygon bounds the region of autocorrelations corresponding toarbitrary 1 1-point harmonic spectra. The solid curve inside thepolygon bounds the autocorrelations generated by a 2-pole filter
sampled at the same 11 harmonic points. If the given autocorrelations lie within this nonconvex region, they can be matchedusing a 2-pole filter and we will have the matching solution.
However, if the given autocorrelations lie between this nonconvex region and the polygon, they cannot be matched using a2-pole filter and we will have the singular solution. In this case,as mentioned earlier, the optimal all-pole filter will have rootson the unit circle. Fig. 3(b) shows the same regions for N =10. In general, as N
CO, the number of vertices of the polygon
increases until the polygon coincides with the parabola. Also,
the region of all-pole generated autocorrelations will coincidewith the parabola and consequently, we will always have thematching solution.
In summary, the optimal all-pole model which satisfies (22)may or may not. satisfy the conditions in (20). If it does, we will
have the autocorrelation matching solution. If the optimal allpole model does not satisfy (20), we will have the singular solution and the optimal model will be unstable. It is necessary toemphasize that the optimal all-pole model in both types of solutions minimizes the spectral distance measure and is unique.
+Iv. ALGORITHM
FOR DISCRETE
ALL-POLEMODELINGIn this section, we present a simple iterative algorithm thatsolves the set of nonlinear equations shown in (22). We thenapply the algorithm to some simple envelope estimation problems to examine its accuracy and convergence behavior. We
also demonstrate that this algorithm will always converge to the
optimal solution and we present possible modifications for improving its speed of convergence.

c-

U P o s i t i v a Definite(b)
Fig. 3. (a) Autocorrelation boundaries for N = 1 1 . The parabola shows
the boundary for the autocorrelations R ( 1 ), R ( 2 ) of positive-definite continuous spectra. The dashed polygon is the boundary for the autocorrelations of general harmonic spectra with N = 11. The nonconvex curve inside
the polygon is the boundary for the autocorrelations of all-pole harmonic
spectra with N = 11 and p = 2. (b) Autocorrelation boundaries for N =IO. The parabola shows the boundary for the autocorrelations R ( 1 ), R ( 2 )of positive-definite continuous spectra. The dashed polygon is the boundary for the autocorrelations of general harmonic spectra with N = 10. The
nonconvex curve inside the polygon is the boundary for the autocorrelations of all-pole harmonic spectra with N = 10 and p = 2.

A . Sirnplijication
To simplify the solution of the nonlinear problem in (22), we
use the following property of sampled all-pole filters:
P

C akk(i- k ) = k ( - i ) ,k=O

for all i(35)where h^( -i ) is the (time-reversed) impulse response of the discrete frequency sampled all-pole model, given by

To prove this property, we start with the identity&(U,)

A(Um) = 1.

(37)

I

417

EL-JAROUDI AND MAKHOUL: DISCRETE ALL-POLE MODELING

Multiplying the H*(a,),and substituting (2), we get

p(w,)

*

A ( w , ) = fi*(w,).

(38)

By expanding A ( a,), we haveP

C u k P ( u me-jwfnk)= A*(o,).

k=O

(39)

Multiplying both sides by eJwmr,
averaging over U , E a, andusing the definition of 8 (i ) in (21), we obtain.

DN

for all i .
By substituting the all-pole property (35) into the minimization condition (22), we obtain the following set of equations
relating the all-pole predictor coefficients to the given autocorrelation sequence:
P

c

U$(i

k=O

- k ) = &( - j ) ,0 Ii 5 p .

(41)

The equations in (41) are written in vector form as

Ra = h(42 1

where h is a column vector with elements 6(- i ), 0 s i s p .
For the continuous spectrum case, we have &( i ) = 0 for i <0 (or h( - i ) = 0 for i > 0) and this set of equations reduces
to that of-regular linear prediction in (3) and (4). In fact, LP
assumes h ( i ) = 0 for i < 0 for both discrete and continuous
spectra. However, for the discrete spectrum case, h ( - i ) isnonzero in general. For the example in Fig. 1, h ( i ) is shownfor -37 s i r: 37 in Fig. l(e). It is clear that assuming h ( i )= 0 for i < 0 is a gross approximation which produces large
errors in the LP envelope estimates. In DAP modeling, we areafter the exact solution to the set ofequations in (41). However,the set in (41) is nonlinear since h ( - i ) depends on the valuesof the all-pole coefficients { uk}. This nonlinear set of equations
can be solved iteratively using the algorithm presented below.

B. m e Algorithm
The algofithm used for determining the predictors ak isstraightforward; it involves two steps repeated iteratively:
Given an estimate of the predictor, evaluate & ( - i using(36).
Given the new estimate of &(), solve the now
equations (41) for a new estimate of the predictors.
The algorithm in detail is as follows.
1) Perform peak picking on the spectrum of the speech signal. Obtain the locations w,, the magnitudes P(w,), and thenumber N of the peaks. (Note that the peak frequencies do nothave to be intege; multiples of some fundamental pitch. However, to minimize the number of spurious peaks, one could apply quasi-harmonic constraints when dealing with voiced speecht141.1
2) Given w, and P(U,) for 1 Irn IN , compute R( i )from ( 5 ) .
3) Using ordinary linear prediction, find an initial estimateof ( a k } , 0 Ik Ip .
4) Compute A ( u , ) for 1 Im IN .
5) Evaluate h( - i ) , 0 5 i Ip , using (36).

6) Solve (41) for the new estimate of { a k } .
7) Evaluate E,, using (2) and (14).
8) If the reduction in E I Sfrom previous iterations is greaterthan some threshold T, go to step 4; else continue.
9) Normalize the coefficients to satisfy (29).
10) stop.
Computationally, the DAP algorithm is more intensive thanLP. Each iteration of the above algorithm requires two real
DFT's of size N and the solution of a set of p
1 linear equations. It is important to note that the autocorrelation matrix R
in (41) is constant throughout the iterations, and therefore need
only be inverted once. Also, R is Toeplitz symmetric, whichallows the use of efficient algorithms in the solution of (41).

+

C . Examples
To demonstrate that the DAP theory and the algorithm given
above are basically sound, we applied the algorithm to the harmonic discrete spectrum in Fig. l(a) and the nonharmonic discrete spectrum in Fig. 2(a). As expected, in both cases DAP
was able to recover the original all-pole filter, with the resulting
spectrum in each case being almost identical to the original(solid) envelope shown in the figures.
In Section V, we present and discuss the results of applying
DAP to real and synthetic speech. But first, we present the convergence properties of the DAP algorithm.

D. Algorithm Convergence
To examine the convergence characteristics of the DAP algorithm, we show in Fig. 4 plots of the spectral error EdBas afunction of the iteration number for discrete spectra obtained
from the 12-pole example of Fig. 1. The lower line shows the
behavior of EdB for a spectral sampling of N = 60 while theupper curve corresponds to a sampling of N = 34. By comparing the two curves, we notice that the algorithm converges fasterfor a higher value of N , partly because the initial LP estimate
is closer to the final DAP result. For smaller values of N close
to 2p, as shown by the upper curve, we notice that the error
decreases every two iterations instead of every iteration. The
reason for this behavior will become evident when we analyze
the convergence properties of the algorithm.
While the algorithm is intuitive and straightforward and wasderived without rigorous mathematics, we will show that it isfundamentally sound and its behavior is not at all surprising.
We will use the vector notation introduced earlier to prove thatthe algorithm is equivalent to a well-known fast gradient technique and suggest an improvement to the algorithm which increases convergence speed.
The algorithm in Section IV-B can be given in vector notation
as follows. Given an estimate a, of the vector a at iteration m,
we compute the vector h,. Then, the new estimate a, + is given
from (42) as
U,+, = ~ - I i i , = R - I R m am.
(43),

With simple manipulation of the above equation, we obtain= U, - R - ~ ( R- R,)u,.

U,,,

(44)

But we see from the derivation of (23) thatg, = 2(R - Rm)am(45 1

is the gradient g at iteration m . Therefore, the update equationin (44) reduces to
%+I

= a,

-

(2R)-'g",.

(46)

418

IEEE TRANSACTIONS ON SIGNAL PROCESSING, VOL. 39, NO. 2, FEBRUARY 1991

a

4

040
I?lR

'0

:?U

Fig. 4. Spectral error E,,, as a function of the iteration number for the DAP
modeling example in Fig. l(a) and for excitation periods N = 34 (top) andN = 60 (bottom).

The update equation (46) is recognized as a special form of
the classical quasi-Newton method [ 101 which in general is expressed asa,,,= a, - as-lg,

where S is the Hessian (matrix of second derivatives) of the costfunction and a! is a scalar determining the length of the update
vector S-lg,. We note that, in the case of the DAP algorithm,
the 2R matrix is used as an estimate for the Hessian. This choice
is appropriate since 2R is a good estimate of the Hessian forlarge N (it is exact for N --t 03 ) and it is positive definite, which
is a requirement for the convergence of such algorithms. We
also note that, for the DAP algorithm, the update scalar is always set to l . While this value is optimal for the quadratic case,
it is suboptimal for nonquadratic problems such as this one. The
use of a suboptimal value of a! is the reason why the error does
not necessarily decrease every iteration. To improve the convergence speed, we incorporate an update scalar a! in our algorithm. The algorithm is then modified to
U,+l

= a,

= a,(l- Q!(2R)-lgm- a!)+ a!R-'h^,,(47)
0 5 Q! I 1 .

(48)

A proper choice of a! can guarantee a decrease in the error every
iteration, which would provide for faster convergence. Fig. 5
shows the convergence curves for the DAP algorithm for thesame example in Fig. 4 with N = 34 and for a! = 0.1, 0.5, 1.O.
We note that, for Q! = 1 (the standard algorithm in (43)), we
have slow convergence. Moreover, we do not have an error decrease every iteration. On the other hand, for a! = 0.5, we havea fast convergence rate and an error decrease every iteration.
The algorithm with a! = 0.: has a much slower convergence
rate and the curve flattens before attaining the minimum error.
This behavior is typical of gradient techniques where if the update parameter is larger than its optimal value, oscillations occur and if it is lower than its optimal value, slow convergenceoccurs. In most experiments, we have found that cx in the range
0.4 Ia! I0.8 was adequate and achieved good results in termsof convergence speed.

V. EXPERIMENTAL
RESULTS
We now present the results of applying DAP modeling toshort-term spectra obtained from synthetic speech and from realspeech.

Fig. 5 . E,, versus iteration for the DAP algorithm for the example in Fig.
l(a) using various values of a with N = 34 andp = 12.

A . Synthetic Speech

We applied DAP modeling to the estimation of formant frequencies of synthetic vowels. We synthesized one-second longsegments of each of three vowels, led, lit, and lul. Each vowel
was synthesized by exciting a fixed 12-pole filter with known
formant frequencies by a sequence of pulses with varying pitch.
The pitch values within each vowel were changed continuously
from 100 to 350 Hz, but never higher than the value of the firstformant. We then estimated the spectral envelope of the synthetic speech for every frame of 20 ms (for a total of 50 frames
per vowel) using LP and DAP modeling with p = 12 poles.
From the envelope in each frame, we located the first three local
maxima as estimates of the formant frequencies and computed
the percent error ofieach estimate. The percent error was defined as I % I = ( I F - { I / F ) x 100 where F is the known
value of the formant and F is the estimate. The first row in Table
1 shows the percent error for the LP formant estimates averaged
over all 50 frames for each of the three formants of each vowel.
The second row in the table shows the average percent error for
the DAP estimates. (The third row will be described in Section
VI). We note that DAP produced better estimates of all the formants for all three vowels. The improvement was the greatest
for the first formant of each vowel and constituted a reductionof a factor of two to three in the percent error. Also, for eachmethod we show the average value of E,, for each vowel. The
decrease in EdBwas on the order of 1 (dB) for all three vowels.
It would be misleading to claim that DAP modeling will always necessarily provide a better estimate of the vocal tract resonances for real speech than other modeling techniques. This
claim would be valid only if the vocal tract can be well modeledusing an all-pole envelope which, as is well known is not always the case. What DAP modeling always provides is an allpole envelope that is on the average closer to the harmonic peaks
than the LP envelope.

B. Real Speech
As an example of the performance of DAP modeling on real
segments of speech, we show in Fig. 6(a) the spectrum of a
20-ms frame of voiced speech and its line spectrum approximation. We also show the LP and DAP 12-pole estimates of thespectral envelope. The LP estimate is shown as a dashed line
while the DAP envelope is the solid line. We note that the DAP
envelope generally fits closer to the harmonic peaks than the LP
envelope. Also, the DAP formants are less biased toward thegiven spectral peaks than the LP formants. The difference be-

I

419

EL-JAROUDI AND MAKHOUL: DISCRETE ALL-POLE MODELINGTABLE IDAP VERSUSLP I N THE ESTIMATION
OF THE FORMANTS
OF SYNTHETIC
VOWELS

/er/

I

F2
1100

F3
1475

2.3
1.0
1.1

7.0
1.8
1.8

LPC
DAP
WDAP '

-

F1

430

til

2.6
0.9
0.9

F2

lul

E,,

F1
265

2296

F3
3007

E,,

F1
296

F2
867

F3
2240

Em

1 .o
0.2
0.1

6.6
3.0
2.3

0.8
0.1
0.2

0.7
0.2
0.7

1.1
0.3
0.1

6.3
2.4
2.2

2.3
0.8
0.8

1 .o
0.2
0.3

0.9
0.2
0.1

tween how LP and DAP spectra match the spectral peaks canbe seen by examining the spectral deviation plots in Fig. 6(b).
This figure shows the spectral error P(o,)/~(w,) in decibels
at the spectral peaks for the LP model (dashed lines) and the
DAP model (solid lines). As expected, DAP modeling provides
a better fit at most frequency points. The convergence behaviorof the DAP algorithm for this spectrum is shown in Fig. 6(c).
We note that the error EdBdecreased from 2.2 for the LP estimate to 1.4 for the DAP estimate in four iterations and that theerror decrease for subsequent iterations was very small. This
convergence behavior, which is typical for real speech segments, allows us to limit the algorithm to very few iterations.
This makes DAP modeling an attractive procedure with someextra computational load but not a prohibitive one. Early stopping of the iterations also has an advantage in the case where
the optimal solution is singular because the algorithm will notreach the unstable optimal all-pole model, but will still achieve
an improved spectral fit and a decrease in error using the firstfew iterations.
In order to compare the general behavior of DAP modelingwith LP modeling for real speech, we present in Fig. 7 a scatterplot of E d , for the two methods when applied to the voiced
segments of a sentence from a female speaker. We note that, as
expected, EdBfor the DAP estimate is always less than E d B forfor these segments ranged
the LP estimate. The decrease infrom 0.15 to 3.5, with an average of 0.65.
Fig. 8 shows the spectrum of a 20-111s frame of unvoiced
speech and the LP and DAP 12-pole estimates of the spectralenvelope. The LP envelope is shown in dashed line while the
DAP envelope is in solid line. We note that the DAP and LP
envelopes provide similar fits to the spectrum. This conclusion
becomes evident when we note that E , j B decreased from 2.5 to
2.4 in 10 iterations.
This result is expected since unvoiced speech has a large
number of spectral peaks and as shown earlier, DAP modeling
reduces to LP for large N .

A

VI. WEIGHTEDDAP

-0

1

2

3

5

4

6

7

(1

9

10

ITER

(C)
Fig. 6. (a) LP and DAP 12-pole estimates of the spectral envelope ofvoiced speech. The DAP envelope is the solid line and the LP envelope isthe dashed line. Number of spectral points N = 40. (b) Spectral e m r plot
for the LP (dashed line) and DAP (solid line) estimates shown in Fig. 6(a).
(c) Convergence behavior of the DAP algorithm for the spectrum shown in
Fig. 6(a).

As mentioned earlier, the use of the Itakura-Saito error measure allows us to weight the error as a function of frequency.
Frequency-dependent weighting has been used in many applications to obtain a better envelope fit in certain frequency regions relative to others. It has been shown that, in speechcoding, a better fitting envelope at the lower harmonics can improve the quality of the synthetic speech [22]. It has also beenshown that proper error weighting can improve the recognition
performance of speech in noise [ZI]. In this section, we modify
our DAP procedure to include a frequency-dependent weightingof the I-S error. We call this method weighed DAP, or WDAP.

IEEE TRANSACTIONS ON SIGNAL PROCESSING, VOL. 39, NO. 2, FEBRUARY 1991

420

Similar to DAP, the WDAP minimizaiton conditions admit two
types of solutions: the autocorrelation matching solution inwhich R , ( i ) = R w ( i ) , 0 5 i Ip , or the singular solution
where the autocorrelation matching is not possible and the optimal all-pole model will have roots on the unit circle.
Also similar to DAP, the minimization conditions (51) canbe simplified using an identity for sampled all-pole filters,
namely,P

C akRw(i - k ) = & ( - i ) ,for all i

(54)

k=O

where h,( - i ) is defined as
9 .D

1.0

2.0

9.0

3.0

5.0

(55)

Ed8 for LFC

Fig. 7 . Scatterplot of E,, for LP and DAP for a female speaker.

By substituting (54) into (51), we obtain a set of nonlinear equations

The set of equations in (56) is solved iteratively using a twostep procedure similar to the one presented earlier for the noweighting case. These two steps, which are repeated until convergence, areGiven an estimate of the predictor, evaluate k,( - i ) using(55).
Given the new estimate of 6,( - i ), solve the now "linear" equations (56) for a new estimate of the predictors.

=. .U

1000.02000.0:ooo.o:ooo.oEJUO.9w o (HZ)
Fig. 8. LP (dashed envelope) and DAP (solid envelope) 12-pole estimatesof the spectral envelope of unvoiced speech.

A . Minimization of the Weighted Error Measure

The weighted error measure for the discrete case is expressedaswhere W ( U , ) is a positive weighting function defined only at
the frequencies of interest U, E Q and normalized such thatIcN

W ( U , ) = 1.
Nm=1

(50)

For an all-pole model, the minimization of E w l s with respect tothe predictor coefficients yields the following conditions:PP

C akR,(i - k ) - C akk,(i - k ) = 0,

k=O

k=O

o5 i 5p
(51)

where

c

l N
R,(i) = W(o,)P(o,) cos ami,

(52)

I N
R,(i) = N mC
= \ ~ ( w , , ~ ) P ( wcos, ) w,i.

(53)

N m = lThe WDAP algorithm exhibits all the convergence properties
mentioned earlier for the DAP algorithm.

B. Results
We applied WDAP to the estimation of formant frequenciesof the same synthetic vowels described in Section V-A. The
only difference in this experiment is that the spectral error isweighted using the function W ( U ) = [ 1/ 1
( U / a c )] where
the cutoff frequency U , corresponds to 800 Hz. This weighting
function corresponds to a mel-scale based weighting. The third
row in Table I shows the average percent error for WDAP estimates for each of the formants for the three synthetic vowels.
We note that WDAP modeling improves the estimates for formants whose values are lower than the cutoff frequency at theexpense of the estimates for the formants that are higher thanthe cutoff frequencies. Estimates of the formants that are relatively near U , do not vary greatly when we use WDAP insteadof DAP modeling.
It is important to note that while a correlation matching solution may exist for a certain weighting function, it may notnecessarily exist for another, even though the spectrum P ( U , )is the same in both cases. This point is illustrated in Fig. 9.
Fig. 9(a) shows the admissible region of the normalized autocorrelation lags R,( 1 ) and &(, 2 ) that can be generated froma 2-pole model, excited with a pulse train with period N = 10
and weighted with W ( U ) in Fig. 9(b). The admissible region in
Fig. 9(a) is to be compared with that in Fig. 3(b) for the uniform
weighting case. It is clear from the two figures that the admissible autocorrelations depend on the weighting functions: therefore, the existence of correlation matching solutions will also

+

I

42 I

EL-JAROUDI AND MAKHOUL: DISCRETE ALL-POLE MODELING

R(2)

(b)

Fig. 9. (a) Normalized autocorrelation lags for two-pole filters weighted
using weighting function W ( w ) in Fig. 9(b) and period N = 10. (b)Weighting function W ( w ) for Fig. 9(a).

be weight dependent. Note that this change in the region ofadmissible autocorrelations with change in spectral weighting
only occurs for the discrete spectral case. For continuous spectra, the region of admissible autocorrelations is always the whole
positive definition region, with or without spectral weighting.

VII. CONCLUSIONSIn this paper, we introduced a new method for estimating theenvelope of discrete spectra. By matching the appropriate autocorrelations, discrete all-pole (DAP) modeling is able to produce better fitting envelopes to a set of spectral points than linear
prediction. We also presented a simple iterative algorithm which
we showed to be a fast gradient technique guaranteed to converge. The basic DAP method was extended to allow frequency-dependent weighting of the spectral error to improve thespectral fit in certain frequency regions of frequency at the expense of the fit at other regions. The basic method can also be
extended to allow pole-zero modeling of discrete spectra [2],[31.

APPENDIXA
In this Appendix, we shall prove the following properties ofthe error measure and the optimal all-pole model: 1 ) The allpole model satisfying the minimization condition (22) is uniqueas long as the number of spectral points N exceeds a certain

threshhold. 2) If the optimal all-pole model does not satisfy the
autocorrelation matching condition in (20), it will have roots on
the unit circle.
Many of the arguments in this Appendix are based on the
one-to-one correspondence between predictor sequences { ak1with poles inside the unit circle and the positive-definite correlation sequence { d k } defined in (18) and (19). From one sequence, we can always compute the other. We start by showingthat the minimum is unique with respect to the sequence { dk } .
Note that we do not intend to compute the optimal set of { dk} .
We are only exploiting the properties of the error function ofthe { dk } parameters to find its properties for the { uk1 parameters. We then show that, depending on whether the optimal setof correlation coefficients corresponds to a positive-definite
spectrum or not, the optimal all-pole model will or will not satisfy the correlation matching condition (20). If the optimal allpole does not satisfy (20), we show it to have roots on the unit
circle.
1. Uniqueness with Respect to { dk}A sufficient condition for the minimum with respect to { dk 1 tobe unique (i.e., for the local minimum to be global) and for EIS
to be a convex function of dk is that the Hessian matrix S ofsecond derivatives be positive definite or, equivalently, havethe product x T S x be positive for all vectors x # 0 . We adapt
a proof presented by Preuss [ 181 for the continuous case to thediscrete case. We first compute the matrix S with elements s,,422IEEE TRANSACTIONS ON SIGNAL PROCESSING, VOL. 39, NO. 2, FEBRUARY 1991

= d 2 E I s / a d j a d j0
, Ii, j Ip. Using the definition of p(o)in (17), one can show that

a2EIS
adiadj1

s.. = -= - C B 2 ( w , ) cos iw, cos jam

‘’

(A.I)

N m= I

and
P

.

I
=-

P

2. Autocorrelation Matching Solution

N

D

D

j=o

r=O

c ~*(w,) ,2,kxjxi cos jw, cos iw,.

Nm=1

Since the summations over i and j are separable, we can write(‘4.3)For the Hessian to be positive definite, we need to have at least
one term inside the summation be different from zero for all x# 0. In other words, we needI)

C xi cos w,i

# 0,

i=O

for x # 0 for at least one frequency point w, E 0.

(A.4)

This condition can be written in vector notation as
Cx # 0,forx # 0(A.5)+

where C is an N X ( p1 ) matrix with cos w, i being the entryof the mth row and the ith column. Since x # 0, the condition
1. The rankin (A.5) is satisfied if the matrix is of rank r >- por number of linearly independent rows of this matrix equals
the number of distinct frequency points. It is important to notethat w and - w do not qualify as distinct frequencies since they
produce identical rows in the matrix. Therefore, r >- p
1 can
be written as++

p + 1 5 No,u.

where Nosuis the number of frequency points in the range [0,
XI.This inequality translates into the following conditions depending on what the set Q of discrete frequencies includes:N
p I2

ifOEQandnEQ

N-1
P I -

Note that (A.6) gives the minimum number of spectral points
required for uniqueness for a given model order.
The form of the optimal all-pole model will depend onwhether the unique optimal value of { dk} occurs within theregion corresponding to positive definite spectra or outside it.
Below we consider both possibilities and derive the propertiesof the optimal all-pole predictors { a k } based on those derived
above for the set { d , } .

2N
p s - 2- 1(‘4.7)

if either 0 E Q or ?r E Q(A.8)

otherwise.

(A.9)

It is these conditions that guarantee at least one of the terms inthe summation in (A.3) to be different from zero, and therefore,
the uniqueness condition to be satisfied.
In conclusion, the error is convex in { d k } and therefore theoptimal set { d k } is unique if the condition in (A.6) is satisfied.

If the optimal set of { dk 3 corresponds to a positive-definite
spectrum, then the corresponding { a k } sequence will exist andwill have poles inside the unit circle. In this case, the optimal
all-pole model will satisfy the matching condition (20) and theminimization condition (22) (hence the name matching solution). Moreover, since the optimal sequence { dk } is unique (asshown by (A.3) above), the corresponding minimum-phase allpole madel { a k } will also be unique. In terms of tke matrix
equation (23), this solution makes the matrix ( R - R ) identically zero and, therefore, (23) will be satisfied by the corresponding all-pole model.

3. Singular Solution
For the case where the unique local minimum in { dk } is nonpositive-definite, there will be no local minimum among positive definite spectral models. However, given the convexity ofthe error in { dk } , the optimal positive-definite { dk } sequence(the one which produces the smallest error among positive-definite spectral models) will be on the boundary of the positivedefinite space. Therefore, the corresponding { a k } will also beon the boundary of positive-definite space. In other words, theoptimal all-pole will have roots on the unit circle. In termsof the matrix equation (23), the optimal all-pole model is thevector a that corresponds to the zero eigenvalue of the matrix(R- 8 ) .
REFERENCES
B. S. Atal and M . R. Schroeder, “Recent advances in predictive
coding-applications to speech synthesis,” in Proc. Speech Commun. symp. 1974, pp. 27-31.
A. El-Jaroudi, “Discrete spectral modeling with application to
speech analysis,” Ph.D. dissertation, Northeastern University,
Aug. 1988.
A. El-Jaroudi and J . Makhoul, “Discrete pole-zero modeling andapplications,” in Proc. IEEE Int. Conf Acoust., Speech, Signal
Processing (Glasgow, Scotland), May 1989, pp. 2162-2165.
H. Hermansky, “Analysis and synthesis of speech based on spectral transform linear predictive method,” in Proc. IEEE In?. Conf
Acoust., Speech, Signal Processing, 1983, pp. 777-780.
H. Hermansky, “Spectral envelope sampling and interpolation in
linear predictive analysis of speech,” in Proc. IEEE Int. ConJ
Acoust. , Speech, Signal Processing, 1984, pp. 2.2.1-2.2.4.
F. Itakura and S. Saito, “Analysis synthesis telephony based on
the maximum likelihood method,” in Proc. 6th Int. Congr.
Acoust. (Tokyo, Japan), 1968, pp. C17-C20.
F. Itakura and S . Saito, “A statistical method for estimation of
speech spectral density and formant frequencies,” Electron.
Comrnun. Japan, vol. 53-A, 1970, pp. 36-43.
V. Jain and R. Hangartner, “Efficient algorithm for multipulse
LPC analysis of speech,” in Proc. IEEE Int. Con5 Acoust.,
Speech, Signal Processing, 1984, pp. 1.4.1- I .4.4.
C. Lee, “Robust linear prediction for speech analysis,” in Proc.
IEEE Int. Conf Acoust., Speech, Signal Processing, Apr. 1987,
pp. 289-292.

I

423

EL-JAROUDI AND MAKHOUL: DISCRETE ALL-POLE MODELING

Luenberger, Linear and Nonlinear Programming. Reading,
MA: Addison-Wesley, 1965.
J. Makhoul, “Linear prediction: A tutorial review,” Proc. IEEE,
vol. 63, no. 4, pp. 561-580, Apr. 1975.
J. Makhoul, “Spectral linear prediction: Properties and applications,” ZEEE Trans. Acoust., Speech, Signal Processing, vol.
ASSP-23, no. 3, pp. 283-246, June 1975.
J . D Markel and A. H. Gray, Jr., Linear Prediction ofSpeech.
New York: Springer, 1976.
R. J. McAulay, “Maximum likelihood spectral estimation andits application to narrow-band speech coding,” ZEEE Trans.
Acoust., Speech, Signal Processing, vol. ASSP-32, no. 2, pp.
243-251, Apr. 1984.
R. Mizoguchi, “Speech analysis by selective linear prediction in
the time domain,” in Proc. ZEEE Znt. Conf. Acoust., Speech,
Signal Processing, vol. 3, pp. 1573-1576, 1982.
[16] Ay Parker, “Low-bit rate spkech enhancement using a new methodof multiple impulse excitation,” in Proc. IEEE Int. Conf. Acoust.,
Speech, Signal Processing, 1984, pp. 1.5.1-1.5.4.
[17] J. Picone, “Joint estimation of the LPC parameters and the multipulse excitation,” Speech Commun., 1986, pp. 253-260.
[18] R. D. Preuss, “Autoregressive spectral estimation in noise withapplication to speech analysis,” Ph.D. dissertation, Oklahoma
State University, Dec. 1983.
[19] R. C. Rose and M. A. Clements, “All-pole speech modeling witha maximally pulse-like residual,” in Proc. ZEEE Znt. Con$
Acousr., Speech, Signal Processing, 1985, pp. 481-484.
[20] S. Singhal and B. Atal, “0ptimizing.LPC filter parameters for
multipulse excitation,” in Proc. IEEE Znt. Conf. Acousr., Speech,
Signal Processing, 1983, pp. 781-784.
1211 F. K. Soong and M. M. Sondhi, “A frequency-weighted Itakura
spectral distortion measure and its application to speech recognition in noise,” in IEEE Trans. Acoust., Speech, Signal Processing, vol. 36, pp. 41-48, Jan. 1988.
[22] V. Viswanathan, M. Berouti, A. Higgins, and W. Russell, “A
harmonic deviations vocoder for improved narrow-band speechtransmission,” in IEEE Int. Conf. Acoust., Speech, Signal Processing (Paris, France), May 1982, pp. 610-613.

Amro El-Jaroudi (S’85-M’88) was born inCairo, Egypt, in 1963. He received the B.S.
degree in 1984, the M.S. degree in 1984, and
the Ph.D. degree in 1988, all in electrical engineering from Northeastern University, Boston, MA.
Currently he is an Assistant Professor at theElectrical Engineering Department at the University of Pittsburgh, Pittsburgh, PA. His research interests include speech processing,
spectral estimation, and digital signal processing algorithms.

John Makhoul (S’64-M’70-SM’78-F’80) wasborn in Deirmimas, Lebanon, on September 19,
1942. He received the B.E. degree from theAmerican University of Beirut, Beirut, Lebanon, in 1964, the M.Sc. degree from Ohio State
University, Columbus, in 1965, and the Ph.D.
degree from the Massachusetts Institute ofTechnology, Cambridge, in 1970, all in electrical engineering.
Since 1970 he has been with Bolt Beranek
and Newman Inc., Cambridge, MA, where heis now a Chief Scientist, working on various aspects of speech communication and digital signal processing. He is also an Adjunct Professor at Northeastern University and a Research Affiliate at the M.I.T.
Research Laboratory of Electronics.
Dr. Makhoul was a recipient of the 1978 Senior Award, the 1982
Technical Achievement Award, and the 1988 Society Award of theIEEE Acoustics, Speech, and Signal Processing Society. His 1975
IEEE PROCEEDINGS
paper on linear prediction has been named a “Citation Classic” by the Institute for Scientific Information. He is a Fellow of the Acoustical Society of America.