1264

IEEE
TRANSACTIONSON ACOUSTICS, SPEECH,
AND
SIGNALPROCESSING,

VOL. ASSP-35, NO. 9, SEPTEMBER 1987

Accurate Autoregressive Spectrum Estimation at Low
Signal-to-Noise Ratio Using a Phase Matching
Technique
Abstract-This paper describes a new method of accurately estimat111, accurate AR parameters are obtainedby a least squareing the parametersof an autoregressive (AR) process contaminated by
fit with the phase error even in the case of low signal-tohigh-level white noise. Based on the phase matching technique, it minnoise ratio (SNR). The proposed method involves nonlinimizes the difference between the phase of the all-zero model and theear
optimization.However,itscomputationalrequirephase of the maximum phase signal reconstructed from the powerments are quite modest as shownby the simulation exspectrum of the observed signal. The parametersof the AR model are
obtained from the finite length sequenceof the estimated all-zero model. perimentsinSection
V . Finally,inSection
VI, other
The proposed method works only when the order of the A R model isadvantages of the proposed methodare also shownby exknown a priori at present. However, since the phase matching techperiments involving resonant vibration generated in a menique satisfies the conditions needed to apply the least mean-squarechanical system.
lowmethod, the AR parameters are estimated accurately even at asignal-to-noise ratio. With the iterativeor noniterative methods as dis11. DISADVANTAGESOF RECENTMETHODSFOR
cussed in the recent literature, it is not possible to reconstruct the allESTIMATINGAR PARAMETERS
zero model from the power spectrum when there
are dips and peaks
having no correlation with the polesof original AR signal in the power
If a signal x ( n ) is contaminated by white noise w( n ),spectrum. The method proposed in this paper allows one to accurately
theobserved signal y ( n ) is described as follows:
reconstruct the phase from the power spectrum in such cases. Finally,
it is confirmed with computer simulations and experiments that they ( n ) = x(.)
W(.)'
(1)
proposed method is useful for accurate estimation of the AR parameIf x ( n ) is the output signal of a pth-order AR model exters.

+cited by white noise e ( n ) ,P
I. INTRODUCTIONx(.)= - C ai x ( n - i )e(.)(2)
HIS paper concerns a method for describing the resoi= 1nant vibration system using the parameters of an AR
where e ( n ) is the white noise uncorrelated withw ( n )andmodel. When noise is added to a signal under analysis,var ( e ( n ) } = 2
the signal is not described by the AR process, even if thesignal is the response of an
AR system because of thevar { w ( n ) } = ow.
2
introduction of the spectral zeros by the additive noise. A
number of methods for estimatingAR parameters of such Each z transform of (1) and (2) is represented as follows:
AR-plus-noise models have been proposed recently in the
Y ( z ) = X ( z ) + w(z)(1')
literature [1]-[5]. However,, accurate AR parameters ofandthe model are not obtained by such methods as describedin Section 11.
Therefore, we have developed a new method of estimating AR parameters accurately, even when the signal
where A ( z ) = 1
a,z-'
U , Z - ~ , and Y ( z ) ,is contaminated by white noise. The proposed method is
X ( z ) , W ( z > ,and E ( z ) denote the z transforms of y ( n ) ,based on a significant extension of signal reconstructionx ( n ), w ( n ) , and e ( n), respectively. The power spectrum
techniques using the phase estimated from the observedP),( z ) of the observed signaly ( n )is expressed as follows:
signal. Since the errorin the estimated phase is distributed
according to the normal distribution with constant mean
and variance for all frequencies as described in Section+T+--- ++

Manuscript received February 3, 1986; revised March 9, 1987.
H. Kanai was with the Research Center for Applied Information Sciences, Tohoku University, Sendai, Japan. He is now with the Education
Center for Information Processing, Tohoku University, Sendai, Japan.
M. Abe and K. Kid0 are with the Research Center for Applied Information Sciences, Tohoku University, Sendai, Japan.
IEEE Log Number 8715564.

- 0:
-

+

A ( z ) * A*( 1 / z * )A ( z ) * A*(l/z*)
02, *

.,"

(3)

where z* denotes the complex conjugate of z . Let the nuB ( z ) * B*( l / z * )merator of the last equation bewhere B ( z ) = 1 ET=, bj z P i ;then y ( n ) may be mod-+0096-3518/87/0900-1264$01.OO @ 1987 IEEE-

KANA1 et a l . : AR SPECTRUM ESTIMATIONAT

1265

LOW SNR

eled as the ARMA(p , p ) model B ( z ) / A ( z ) driven by thewhite noise u ( n ) of the power ai. The z transform of theoutput signal y ( n ) is expressed as follows:

where U ( z ) is the z transform of u ( n ) . Thus, the observation noise introduces spectral zeros, which are locatedbetween the true AR poles and the origin in the z plane.
Thezerosmovetowardthetrue
AR poles as the SNR
decreases [3].
A number of methods for estimating AR parameters in[ 11such cases have been proposed in the recent literature[5] and can be divided into the following five classes [l],121.
1) The larger orderAR spectral estimation in which themodel is approximated to a larger order AR model [3].
2) The method using modified Yule-Walker equations.
3) The noise compensation method for the correlation
function in which the estimated additive noise power
issubtracted from the main diagonal of the autocorrelation
Fig. 1 . Illustration showing the procedure of the proposed method for esmatrix [4].
timation of the parameters of the AR process in white noise.
4) The directpowerspectrummatchingtechniquewhich minimizes the average squared difference between
the nonparametric power spectrum ( = periodogram ) es- sumption that the real signal y ( n ) has a rational z transtimated from the observed signal and the parametric powerform of the AR-plus-noise signal. The five processes ofthe method are explained below using the simple example
spectrum of the model [5].
za and theof a first-order AR model which has a pole
5 ) The methodusingtheoverdeterminedmodified
zb,
(
I
z
b
1
<
1
za
I
<
1
) in the z
introduced
spectralzero
Yule-Walker equations [2].
plane.
By
using
the
same
example,theright-handside ofAll of these five methods estimate theAR parameters1
shows
schematicallytheprocessesofthemethod,Fig.
using the autocorrelation function estimate
R ( 7 ) or its
1 illustrates the processes
equivalent, the power spectrum estimatePy( w ), which is and the left-hand side of Fig.
byapplying
them
to
the
second-order
AR model used inequal to the periodogram in the fourth method
4) whereV.
Sectionboldface is used throughout the paper to distinguish esti1) The power spectrum estimate Py ( w ), which is themates from the true values. However, accurate
AR paaveraged
periodogram, is obtained using FFT of the obrameterscannot be obtained by thosemethods asdey ( n ) [see Fig. 1 , part ( l ) ] . Since P,,(w ) has
servedsignalscribed below. Each of the real and imaginary parts of the
the
pole
z,,
the
zero zb, and their complex conjugate respectrum Y ( w ) , estimated from the AR-plus-noise signal
1/
z $ and 1 / z $ , the power spectrump,, ( z )ciprocalpairsy ( n ) of (4),is distributed independently according to theisexpressedasfollows:
normaldistribution N [ O , u:(w)] wherethevarianceis.:(a)= J B ( w ) / A ( w ) *( ~ = P y ( w ) and P y ( w ) de2
B ( z ) B*(l/z*)
P y ( z ) = uu *notes the true power spectrumof y ( n). Therefore, the peA ( z ) A*( 1 / z * )riodogram estimate Py( w ) of y ( n ) is distributed according to the chi-squared distribution with 2 degrees of freedom where the mean andthe variance are ( w ) = Py( w )and a ; ( w ) = P ; ( w ) , respectively [6, ch. 11.31. Thus,
2 ) Basedon the noniterativeminimumphasereconthe periodogram Py ( w ) has a large variance, especially in
struction
algorithm [7], [8], theminimumphasesignal
the frequency band where the signal-to-noise ratio is high.
Xmin(z)is reconstructedfrom the periodogram Py( w )Therefore, accurate AR parameters are not obtained by
using
FFT
two times. Since this operation is equivalent
the methods described above, even
if the (generalized)
tomapping
maximum phase poles and zeros to their conleast square fitting technique is applied to the minimizajugate
symmetric
counterparts with the unit circle of thetion of the average squared errors involved in the estiz
plane,Xmin(z)isobtained as follows:
mated power spectrum or the estimated autocorrelation.

.$

-

.',

111. PHASE MATCHING BASEDMETHODFOR ESTIMATINGAR PARAMETERS[seeFig. 1 , part (2)]. Thisequation
Here, we propose a new method for estimating the AR where C, =parameters. The following discussion is based on the as- showsthattheminimumphasereconstructedsignal

a

1266

IEEE
TRANSACTIONSON ACOUSTICS.
SPEECH,AND SIGNAL
PROCESSING,

VOL. ASSP-35, NO. 9, SEPTEMBER 1987

Xmin(z) isthe
output of the ARMA( p , p ) model
B (z) / A ( z ) excited by the impulse Co * 6 (12).
3 ) The phase Omin(w ) of theminimum phase signal
Xmin(z)is calculated. By changing the sign of e m i n ( @ ) ,the new phasefunction 8 ( w ) is obtained. Since a) bymapping the poles and the zeros to their conjugate symmetric counterparts with the unit circle of the z plane, thesign of thephasefunction is changed, and b) thepole
1/ z $ or the zero 1 / z $ makes the same contribution to thephase as the zero 2, or the pole zb, respectively, the folFig. 2. Illustration showing the phase emor A@ ( w ) in the minimum phaselowing four models have the same phases as 8 (a):model X,,,,,( w ) which is reconstructed from the power spectrum of theobserved signal y(n).

of the original AR model in the estimated power spectrum
Py ( w ) . Thus, in a practical case, the minimum phase signal Xmin(z ) of (6) calculated from Py( w ) is expressed as
the sum of two components:

andxzero(z)= c, . ( 1 - zuz-')(l - z $ z ) .

(7)4) By restricting the signal reconstructed from the phase
8 (w ) to beingfinitelength,onlytheall-zeromodel
X,,,,( z),having the same phase as8 ( w ), is obtained [seeFig. 1, part (4)J. A method of obtaining the signal fromits phase is described in Section IV.
5 ) As the method ordinarily used in the linear prediction of speech todetermine the poles of the AR model [9,p. 291, the poles of the AR model are determined fromthe roots of the polynomial:

where A U ( z ) denotes approximately the z transform ofGauss random noise A u ( a ) . Then, the spectrum Xmin(w)
of the estimated minimum phase model is expressed asfollows:
X m i n ( w > = Xrnin(w) + A x ( w )(11)where Xmin(w ) and AX( w ) denote the spectrum of the true
minimumphasemodeland
the spectrum error, respec2. Using (lo), Xmin(w) andtively,asshowninFig.
A X ( o ) are expressed as follows:P
Xzero(z)=n= -pxzero(n)

.z

-n

= 0.

Since (,x,,- p ) is not zero, by dividing (8) by x,,,,,(
ZP, the following polynomial is obtained:

-

(8)

and

-p)

2p f Ic x:,m(n) . z

-n

n =.o

= 0

= Xmin(0) * A U ( w ) .

(9)where x&( n ) = (,x,,n - p ) /x,,(- p ) . The roots ofthe polynomial contain the pole z, ( 1 z, 1 < 1) and thezero 1 / z t ( I 1/ z g I > 1 ) in ( 5 ) . If the roots whose magnitudes are less than one are selected, only the poles ofthe original AR model are obtained [see Fig. 1, part (5j].
As described in Section 11, almost all the methods inthe literature have estimated AR parameters using the estimated power spectrum ( = periodogram j or its equivalent, the estimated autocorrelation function. However, thevariance of the errors in the periodogram estimate is equal
to P : ( w ) and is large, especially in the frequency band
where the poles makepeaks.However,thevarianceofthe error A 8 ( w ) of the phase 0 ( w ) estimated in the third
process of the proposed method is constant for all the frequencies as described in the following.
Since the observed signal is the finite length, there aremany peaks and dips having no correlation with the poles(12b)

Thus, the estimated phase Omin(w ) of the minimum phasesignal Xmin(z)is expressed as the sum of the true phaseemin(w ) and the phase error A 9 (a):Ornin(w)= O,~,(U) + A ~ ( u ) .

(13)

The left-hand side of Fig. 1, part ( 3 ) shows that truephase
B ( w ) = - 8 m i n ( w ) and its estimate O ( o j = --emin(w).
When 1 A 9 ( w ) 1 is significantly less than 1.O, A 8 ( w ) is
approximated by tan A 0 ( a ) .Then, A 8 ( w ) is expressed
as follows:
A O ( w ) = tan AO(w){ amin(0) - emin(w>}tan ernin(w ) - tan Omin(w )

= tan

-

I + tan Omin(w ) . tan Om,,(w ) .

(14)

Each spectrum of (1 I ) is decomposed into real and imag-

KANA1 et al.: AR SPECTRUMESTIMATIONAT

efore,
ts
follows:
inary

1267

LOW SNRas
Xmin(w) = Rmin(w) + jlmin(0)A X ( w ) = AR(w)+ jAZ(w),= xzero(-U)

exp { +je(w)}.

P+ j[Zmin(w)+ AI(w)].(1%)

x,,,,(n)n = -p*sin { n u + 0 ( w ) } = 0(for all w ) .

error (14) is expressed as

I f I A R ( w ) l << l R m i n ( U ) l and I A I ( w ) ) << l I m i n ( W ) I ,the phase error is approximatedby the following equation:
A l ( w ) Rmin(w) - AR(w) Imin(W)

Rmin(w)2 + I m i n ( w >2- Im [ A X ( w ) * x$in(w)]I

2

IXmin(@)= ~m [AX(w)/ x m i n ((19)

Using the definition of the Fourier transform, the imaginary part of the equation is expressed as follows:

xmin(w) = [Rmin(w) + A R ( ~ ) ]A0(w)xzero( +a) ~ X P - j ~ ( w ) }

(15b)

and

Using these spectra, the phasefollows:

{

(154

3

w)]

= Im
(17)
[AU(w)].

(20 1+

Since 0 ( w ) is selected from ( 2 p
1 ) distinct values ofw, (2p1 ) linear homogeneous equations are obtained
from (20).
When 0 (w ) expresses the true phaseof the all-zeromodel, x,,,,( n ) can be determined from the linear equations. However, thereis a phase errorof A 0 ( w ) in phase
0 ( a )calculated from the estimated power spectrum
asshown in Fig. 1, part (3), and this method does not compensate for the phase error at all. Therefore, an accurate
all-zeromodeland
AR parameterscannotbeobtainedusing this method.

+That is, the phase error is equal to the imaginary part ofthe spectrum A U ( w ) of the Gauss random noise in (10). B. Iterative Method (71, [lO]-[l I ]Thus, the phase errors { A 0 ( w ) 1 shown in Fig. 1, partLet x,( n : i ) be the estimate of the all-zero model at(4) aredistributedaccordingtothenormaldistribution
the ith iteration. In the frequency domain, the phaseofN [ O , 0 2 ] where the variance is constant for all the fre(,x,,n : i ) is replaced by the estimated phase 0 ( w ), andquencies.Therefore,moreaccurate
AR parametersare
in the time domain, the finite duration constraint
is imobtained by the least mean-square fit of the phase of the posed. Then, the ( i1)th estimate is obtained, and theall-zero model Xzero(z)to the phase 0 ( 0 ) = -Omin(0 ) unknownspectrummagnitudeisretrievedgraduallyinestimatedfrom theobservedsignal.
A phasematching
successive iterations.
method is proposed in Section IV-C.
of the allWhen 0 ( w ) can accurately express the phasezero model, it is known that the iteration procedure alIV. RECONSTRUCTION
OF AN ALL-ZERO
MODELFROMn ) where A
ways converges to a unique limit ofA . xZem(ITS PHASE
is apositiveconstant [lo]-[ll]. However, as shown in
In the literature, there are two methods for reconstruct- Fig. 1, part (3), there is a difference of { A 0 ( w ) } being a finite duration mixed phase signal from its phase. In tween the obtained phase@ ( w ) and the true phase0 ( w )thissection,the difficulties inreconstructing the signal of the all-zero model as follows:from the noisy phase 0 ( w ) using these methods are exA 0 ( w ) = @(LO) - 0 ( ~ )(21)plained,
and
then
an
alternative
new reconstruction
where A 0 ( w ) of this equation is equal to -A 0 ( w ) of
method is proposed.
(13). Thus,thesatisfyingall-zeromodel
is notreconx,,,(n)structed. This may be explained as follows. Let
A . Noniterative Method (71, [lO]-[I 11and x,,,(n) bethesignalsofall-zeromodels
reconThis method reconstructs afinite length ( = 2p
1)
structed from the true phase0 ( w ) and its estimate@ ( w ) ,signal x,(n ) , ( -p In 5 + p ) of the all-zero model
respectively. Substituting 0 ( w ) of (21) into
(19), theXzero(z)from the phase estimate @(a).Since the signal
spectrum X,,( w) of x,,,( n ) is expressed as follows:
xzer0(n)of the all-zero model is real, the following relation holds well for the discrete Fouriertransform X,,( w ):~ X [P- j { ~ ( 0 -) ~ e ( w ) ) x]z e r o ( w >++= ~ X [j{@(W,P- ~ e ( w ) } ] Xzero(-w).

X,,,( w ) of x,,,,(
Comparing (22) to (19), the spectrum

(22)

n)

1268

IEEE
TRANSACTIONSACOUSTICS,
ON

SPEECH, AND SIGNAL
PROCESSING,

VOL. ASSP-35.
NO.

9, SEPTEMBER 1987

The finite durationsignal xZero(n)is obtainedfromthe
given phase 0 ( k ) by minimizing the following sum CY:

is expressed by [lo, sect. VI]

M- 1

where X is theconstant ensuring that x,,( 0) = 1; assuming that I A 8 ( w ) I << 1, X,,,( w ) of (23) is approximated byThus, x,,,(n) is described by

CY

1A8(k)12

=
k=O

+

MIN.

(29)

This nonlinear minimization is performed using
a standard optimization technique such as the Marquardt method[12]. Since the spectral zeros introduced by the observed
noise are closeto the poles of the AR model in the z plane,0 ) of thethe z transform of the initial estimate xzero(n:iteration is selected as follows:P

where e ( n ) is the real odd signal of the inverse Fouriertransform of the imaginary odd sequence j A 0 ( 0 ) .Since
the phase errors { A 8 ( w ) ] are distributed randomly according to the normal distribution as described in Section
111, e ( n ) has an infinite duration, and consequently, theduration of the reconstructed signal xzero(n)is also infinite. However, the finite durationconstraint is imposedin each iteration. Therefore, an accurate all-zero modelcannot be obtained by this method.

X,,(Z:O) = II ( 1 - z ~ z - ' ) (1 - z ~ z ) , ( 3 0 )i= 1

where { zi } are the poles obtained by the method using
the overdetermined modified Yule-Walker equations.
If the denominator of theright-hand side in (28) is equal
to zero, Ae ( w ) is evaluated as being unstable. Thus, the(i1)th estimate x,,( n : i1) is determined from thei ) under the condition that the deith estimate xZero(n:nominator of the right-hand side in (28) can only be positive.
C. A Proposed Reconstruction Method Using the Least
When the SNR is low, 1 A 8 ( w ) 1 is less than a / 2 for
Square Method Which Fits the Phase Error
w . Thus,the initial estimate of theall-zero modelall
Substitution of 8 ( w ) of (21) into (20) results in
Xzero(z)satisfies the conditionbecause every zero inP
volved in X,,,( z : 0) of (30) has its complex reciprocalxZem(n)* sin { n u + O ( w ) - A 8 ( w ) } = 0
pair,
and then the phase of X,,,,( z : 0 ) is zero for all w .
n = -pThe optimal estimate of the all-zero model also satisfies(for all a).
( 2 6 ) the condition because the phase of the denominator of (28)0( w ) =Using the formula sin ( Q l - Q,) = sin Q , cos f12 - cos is described by { phase of X,,,,(- w ) }
A
e
(
w
)
,
and
1
Ae(o)
1
is
less
than
a
/
2
for
all
w.
Q l sin Q,, and assuming that 1 A 8 ( w ) 1 << 1, this equaHowever,
when
the
SNR
is
high,
1
8
(
w
)
I
is
not
alwaystion is approximated by
less
than
n/2.
In
such
cases,
thepositive
constant
P isP
added to the power spectrum estimate P,,( z ) in (3) or ( 5 )so that the absolute value of 0 ( 0 ) calculated from themodified
power spectrum P;,(0 ) = P , (z) + p is less than* cos { n uO(w)}] = 0(for all a). (27)z / 2 for all w . By adding 0,the initial estimate satisfies
Then, the phase error A 8 ( w ) is expressed as follows:
the above condition, and each Marquardt iteration deterP
mines the successive estimate of the all-zero model underC
xzero(n)* sin (no ~ ( w ) }theabove condition; the optimum estimate
also satisfies
n = -pA8(w) =
the above condition. Sincethe variance of the phase errors
A e ( w ) does not increase by the addition of /3, and thexzem(n). cos {no + O ( w ) }n = -ppositions of all the poles involved in Pi (z) are equal tothose of the poles in P,(z); the AR parameters(for all w )(28aare obtained after convergence in the Marquardt method.
where xzero(0 ) = 1. When 0 ( w ) is selected from equally
Since the phase reconstructionis nonlinear, the positive
distributed distinct M values with intervals of 2 a / M , the constant /3 is obtained in the following way. The calcudiscrete phase error is expressed as follows:
lation is in two stages, the first of which determines theminimumvalue PO of the real part of Xmln(w)reconPC
xZero(n). sin { 2 a ( n k / ~ )+ ~ ( k ) ) structed from P,(w) for all w . If Po is positive, j @ ( @ ) I
n = -p= I - Omin(w ) I is less than n / 2 for all w and both the allA8(k)=zeroand AR parameters are estimated by minimizC x z e m ( n ). cos { 2 x ( n k / M ) + O ( k ) ) ing CYmodeln = -pof (29) using 8 ( 0 )without adding PO to P y (0 ) . If
Po is not positive, in the second stage, Xmin(w ) is recon~m [Xzero(-k) * ~ X P{ j ~ ( k ) } ]structed from Pl.(w) = P , ( o ) + P where = &, 2P0,
4/30, 8/30, . . . ,-&,where $6 is the first of these values at
Re [Xzeco( - k ) * exp { j @ ( k ) ) ]which the real part of Xmin(w ) reconstructed from Pi.( w )(fork = 0 , 1, 2, . * * , M - 1 ) . (28b) = PJ,(w )Ph is positive for all w . It then follows that

+

+

+

-

+

+

c

+

KANA1 et al.: AR SPECTRUMESTIMATIONATLOW SNR

1269

.\

I

-

.\',

VI

5-30

OVER-DETERMINED MODIFIED

YULE-WALKER METHOD'\POWER SPECTRUM

.'

Q:

-7r

-7r

0 ANGULAR FREQUENCY w

0 ANGULAR FREQUENCY

w

,

I

0

-5

QRe

'

+5

I

SIGNAL TO NOISE RAT10 (SNR) ( d B )Fig. 4. The comparison of the AR parameters estimated by the two methods for various SNR (for the first example).

Z-PLANE(b-2)
Fig. 3. The characteristics of the AR model (for the first example). (a) The
fourth-order AR model of (31). (b-1) AR model buriedin the white noise( S N R = - 5 dB). (b-2) The polesof the AR model and the spectral zeros
introduced by the additive noise.

[2, eq. (2.23)] are evaluated for 20 distinct valuesof 7satisfying 7 > 4. Estimated AR parameters { ai} obtainedby the methods are evaluatedby.the following normalized
mean-square error (NMSE) [14] as
NMSE =

the all-zero model and AR parameters can be estimatedusing 0 ( w ) obtained from Xmin(w ).
V. SIMULATIONRESULTS
In order to illustrate the characteristicsof the proposed
method and to compare the proposed method to the overdetermined modified Yule-Walker (OMYM) approach [2,eq. (2.27)], we choose themostpopulartwoexamples
used in the literature [4], [ 131.
First Example: Consider the following four-order AR
process:4~ ( z =) C ai - z -i
i=O

(314

where

a. = 1,al = -2.7607,

a3 = -2.6535,

anda2 = 3.8106,
a4 = 0.9238(3W

Fig.3(a)and(b)showsthecharacteristicsoftheARmodel and the AR-plus-noise model ( SNR = - 5 dB),
respectively. This AR process has four poles at both 0.98
exp ( *jO.69) and 0.98 exp ( k j O . 8 8 ) , which are close toeach other as shown in Fig. 3(b-2). The two methods are
implemented on a MELCOM-COSMO 700 (1 word = 32bit) computer using single precision arithmetic. The total
length of thesynthesizedAR-plus-noisesignalis8192
points. The nonparametricpowerspectrumusedinthe
proposed method is estimated by summation of 16 periodograms,each of which is computedusingaHanningwindowwithalengthof 512points.IntheOMYW
method, by usingthesameAR-plus-noiseprocessof8191-point length, the extended Yule-Walker equations

2(ai - k

ai)

C a:
I

where k is the scaling constant chosen
to minimize the
NMSE. Typically, by using the proposed method, about
10-50 Marquardt iterations were needed to achieve convergence in the case considered here.
Fig. 4 shows the
NMSE as a function of the SNR. When the SNR is about
5 dB, both methods estimate the AR parameters accurately. In this case,the evaluation of thephaseerror
A 0 ( w ) of (28) is unstable since 1 Omin(w ) 1 of the minimum phasesignalobtainedfrom
the powerspectrum
Py( w ) has a value greater than n / 2 for the frequencies
around the central frequencies of the poles as shown Fig.
5(a), and the denominator of the right-hand side of (28)
has small values for the frequencies. Then, as described
previously in Sectioh IV, by adding the positive constant
P to the power spectrum PY.(w ) so that 1 €Imin( w ) 1 calculated from P$( w ) = Py(w ) fl is less than 7r/2 for all'w
as shown in Fig. 5(bj, a satisfactory AR model is estimated by the proposed method as shown in Fig. 5(c).
When the additive noise level was increasedby 5 dB,
the error in estimated parameters obtainedby the OMYW
4.
methodincreased by about15dB asshowninFig.
Nevertheless, the proposed method estimates the AR parameters accurately, even in the lower SNR cases when
-5 dB as shown in Fig. 6.
the SNR
Secorid Example: Considerthefourth-order AR process with parameters which are given as follows:

+

+

-

a. = 1,al = -1.352,

a3 = -0.662,

anda2 = +1.338,
a4 = +0.240.

(33 1

Fig. 7 shows the characteristicsof the AR process and the1270IEEE TRANSACTIONSONACOUSTICS,SPEECH,AND SIGNAL PROCESSING, VOL. ASSP-35, NO. 9, SEPTEMBER 1987

ANGLiLAR FREQUENCY

W

(a)
(b)
Fig. 7. The characteristics of the AR model (for the second example). (a)The fourth-orderAR model of (33) and the AR model buried in the white
noise ( S N R = f 4 dB). (b) The poles of the AR model and the spectrgl
zeros introduced due to the additive noise.

(C)
Fig. 5. In case of high SNR, by adding the positive constant to the power
spectrum, the accurate AR parameters are estimated from the pro'posed
method (for the first example). (a) In case of SNR = + 5 dB, the evaluation of the phase error A 0 ( a )of (28) is unstable since I Om,,( w ) 1 of
the minimum phase signal obtained from the power spectrum estimate
Pj( w ) is larger than7r/2 for the frequencies near the central frequenciesof the poles. (b) By adding the positive constant6 to the power spectrum
P, ( w ), the absolute value of Omin( w ) calculated from P;( w ) = P, ( w )p is less than 7r/2 for all w . (c) A satisfactory AR model is estimated(a)(b)
by the proposed method.
andphase obtained by the two
Fig. 8. Theestimatedpowerspectrum
methods (for the second example ( S N R = +4 dB). (a) OMYW method.
(b) The proposed method.
rn
1

+

(a)
(b)
Fig. 6. Theestimatedpowerspectrum
and phase obtained by the two
methods (for the first example) ( S N R = -5 dB). (a) OMYW method.
(b) The proposed method.

AR-plus-noiseprocess (SNR = +4 dB).Fig. 8 showsthe estimated spectrum of the AR process obtained by theabove two methods. The SNR is about +4 dB. The AR
model estimated by the OMYW method has a large degreeof en,or as shown in Fig. 8(a). However, by using the
proposed method, a satisfactory spectrum is obtained.
'

VI. EXPERIMENTAND THE RESULTS
We applied the proposed method to the estimation ofcharacteristics of the resonant vibration caused by many
flaws on the rough race of a small-sized ball bearing.

The resonant vibration signal is detected by the following procedure [ 151, [ 151. The outer ring i s fixed by imposing an axial pressure, and the inner one revolves at aconstant speed o f 1800 RPM. Under such conditions, aflaw on the race causes radial movement of the outer ring,and the resulting signal from the movement is sensed bya vibration pickup attached to the outer ring in an Anderson meter [ 151. The signal is amplified and filtered througha high-pass filter to eliminate the primary frequency component ( = 30 Hz) corresponding to the rotation of theinner ring. The filtered signal is A / D converted with a
12 bit A / D converter at a sampling period of 30 ps.
Since theresonantvibrationsignaldetectedaboveis
buried in the.high level noise, the detected signal z ( n ) is
approximately expressed as the AR-plus-noise process in(1) and ( 2 ) . To estimate the characteristics of the resonant
vibration x ( i t ) using z ( n), the following five processes
are cairied out (see Fig. 9).
1 ) After multiplying each input sequence of 83 normal
bearingshavingsmoothracesby the Hanningwindow,
2048-point FFT's are carried out. The power spectra areaveraged for the 32 successive nonoverlapping sections.
Then, by summingup the powerspectra of 83 normal
bearings,theaveragedpowerspectrum
P,( w ) isobtained. Fig. 9(a) shows the inverse of P, ( w ). If the ad-

1271

KANA1 et 0 1 . : AR SPECTRUM ESTIMATION AT LOW SNR

posed method is set to be the unit impulse as follows:

h wv
PKE-WHITENED

Using the initial value, the denominator of the right-handside in (28) is always positive for all frequencies. Then,
as described in Section IV-C, by confining the denominator to a positive value in each iteration, 12 iterations
wereneeded toachieveconvergenceforthisexample.
Since the inputof the proposed method is the power spectrum of the.observed signal, the prewhitening process is
simple as described above.

VII. CONCLUDINGREMARKSA new spectrumestimatemethodwhich
is based on
ESTIMATED AR MODEL
phase matching is proposed for the accurate estimation of
AR processburiedinhighlevel
theparametersofthe
1 POWER SPECTRUM I/(e)I
noise.
From
experiments
withnumericalexamples,the
2.44k
6.6k
FREQUENCY(Hz)AR parameters of the AR-plus-noise process at low SNR
are estimated accurately using our new method.
Fig. 9. The estimated AR spectrum of the resonant vibration of a ball bearing obtained by the proposed method. (a) The inverse spectrum P,( m)-’Two issues remain for future research as follows.
of the average power spectrum of the vibration signals of the 83 normal
AR model should be derived by
1) The poles of theball bearings. (b) The power spectrum P , ( w ) of the vibration signal ofcomputing the roots of the polynomial obtained from thea defective ball bearing with rough races. (c) The prewhitened power
spectrum P , ( w ) obtained by multiplying the power spectrum P , ( w ) in
finite length all-zero model estimate as described in Sec(b) by P,(w)-’ of (a). (d) The power spectrum P; ( w ) obtained by setion11.
lecting and splitting the important frequency band of the power spectrum
2) In the experiments, the order of theAR model is
P , ( w ) of (c). (e) Estimated second-order AR model obtained using the
proposed method from Pi(w ) of (d).
alreadyknown.However,foranunknownsignal,
it isdifficult to choose the order of the AR model. In such acase, the ordinary spectrum estimation methods use comditive noise is white, P, ( w ) denotes the power spectrum putationally efficient order recursive techniques tofind ARof the transfer function of the signal measurement proparameters for the various ordersof the model. However,cess.
since -the method proposed in this paper uses nonlinear
2) After multiplying the sequence of a defective bearing
optimization for phase matching, it is difficult to use suchhaving a rough race by the Hanning window, 2048-pointrecursive techniques.
FFT’s are carried out. Thepower spectra are summed up
These important issues are currently under investigaforthe32nonoverlappingsuccessivesections.
The re- tion.
sulting averaged power spectrumP, ( o ) denotes the prodREFERENCES
uct of P, ( w ) and the power spectrum Py( 0 ) of the ARplus-noise signal.
S. M. Kay and S. L. Marple, Jr., “Spectral analysis-A modern perspective,” Proc. IEEE, vol. 69, pp. 1380-1419, Nov. 1981.
3) By dividing the power spectrum P, ( 0 ) by P, (a),J . A. Cadzow,“Spectralestimation:
An over-determined rational
prewhitened power spectrum P,( w ) is obtained as shown
model equation approach,” Proc. IEEE, vol. 70, pp. 907-939, Sept.
in Fig. 9(c).
1982.
S. M. Kay, “The effects of noise on the autoregressive spectral esti4) We .selected the important frequency band ofPy( w )mator,” IEEE Trans. Acoust., Speech, Signal Processing, vol. ASSPand clipped it from Py ( w ) , and then the power spectrum
27, pp. 478-485, Oct. 1979.
Pi ( w ) is obtained as shown in Fig. 9(d). The number of-,“Noise compensation for autoregressive spectral estimators,”points of the power spectrum Py ( a )is 512.
IEEE Trans. Acoust., Speech, Signal Processing, vol. ASSP-28, pp.
292-303, June 1980.
5) It is difficult for theproposed method to estimate theB. Friedlander and B. Porat,“Aspectral
matching technique for
AR parameters when the order ofthe AR model is not
ARMA parameter estimation,” IEEE Trans. Acoust., Speech, Signal
known a priori. However, it is known that the central freProcessing, vol. ASSP-32, pp. 338-343, Apr. 1984.
A . V. Oppenheim et a l . , Digital Signal Processing. Englewood
quency of the resonant vibration of the outer ring in ball
Cliffs, NJ: Prentice-Hall, 1975.
bearing is about 4 kHz [ 151, [ 161. Thus, we applied ourT. F. Quatieri, Jr. and A. V. Oppenheim, “Iterative techniques formethod to estimate the poles of the real vibration underminimum phase signal reconstruction from phase ormagnitude,”IEEE Trans. Acoust., Speech, Signal Processing, vol. ASS€’-29, pp.
the assumption that the vibration is described by the sec1187-1193, Oct. 1981.
ond-order AR model in the frequency band around 4 kHz.
B. Yegnanarayana and D. K.Saika, ‘,‘Significance of group delay
Using the proposed method, the AR parameters of the vifunctions in signal reconstruction from spectral magnitude or phase,”
IEEE Trans. Acoust., Speech, Signal Processing, vol. ASSP-32, pp.
bration signal are obtained from P i ( w ) . Fig. 9(e) shows
610-622,June 1984.
the characteristics of the estimated resonant vibration.The
J. D. Markel and A. H.Gray, Jr., Linear Prediction ofSpeech. New
initial value of the Marquardt interation used in the proYork: Springer-Verlag, 1976.

1272

IEEE
TRANSACTIONSACOUSTICS,
ON

SPEECH, AND SIGNAL
PROCESSING,[lo] M. H. Hayes, J. S. Lim, and A. V. Oppenheim, “Signal reconstruction from phase or magnitude,” IEEE Trans. Acoust., Speech,Signal
Processing, vol. ASSP-28, pp. 672-680, Dec. 1980. ’[ l l ] G. A. Merchant and T. W. Parks, “Reconstruction of signals fromphase: Efficient algorithms, segmentation, and generalizations,” IEEE
Trans. Acoust., Speech, Signal Processing, vol. ASSP-31, pp. 11351147, Oct. 1983.
[12] S. L. S. Jacoby, J. S. Kowalik, and J. T. Pizzo, Iterative Methods
for Nonlinear Optimization Problems. Englewood Cliffs, NJ: Prentice-Hall,1972.
[13] T. J. Ulrych and T. N. Bishop, “Maximum entropy spectral analysis
and autoregressive decomposition,”Rev. Geophys. Space Phys., vol.
13, pp. 183-200, Feb. 1975.
[14] C. Y. Espy and J . S. Lim, “Effects of additive noise on signal reconstruction from Fourier transform phase,” IEEE Trans. Acoust.,
Speech, Signal Processing, vol. ASSP-31, pp. 894-898, Aug. 1983.
[15] H. Kanai, M. Abe, and K. Kido, “Detection and discrimination offlaws in ball bearings by vibration analysis,” J . Acoust. SOC.Japan,
vol. 7, pp. 121-131, Mar. 1986.
“Detection of slight defects in ball bearings by nonperiodic[16] -,
analysis,” J . Acoust. SOC. Japan, vol. 7 , pp. 219-228, July 1986.

VOL. ASSP-35, NO. 9, SEPTEMBER 1987

Masato Abe (”85)

was born in Sendai, Japan,
on May 3, 1952. He received the B.S., M.S., andPh.D. degrees from Tohoku University, Sendai,
Japan, in 1976, 1978, and 1981, respectively.
He currently works as a Research Associatewith the Research Center forApplied Information
Sciences, Tohoku University. His research interest is in digital signal processing for acoustics andin computer architecture.