111 J. J. Thomas, âSelected digital signal processing algorithms over a. Cornput., vol. 25, Apr. 1971. York: Academic, 1975. finite field,â Ph.D. dissertation, Univ.
IEEE TRANSACTIONS ON ACOUSTICS, SPEECH, AND SIGNAL PROCESSING, VOL. ASSP-31, NO. 1, FEBRUARY 1983
[g] J. M. Pollard, “The fast Fourier transform in a fiiite field,” Math. Cornput., vol. 25, Apr. 1971. [lo] K. G . Beauchamp, Walsh Functionsand Their Applications. New York: Academic, 1975. [ 111 J. J. Thomas, “Selected digital signal processing algorithms overa finite field,” Ph.D. dissertation, Univ. Missouri, 1980.
On the Zeros of the Linear Prediction-Error Filter for DeterministicSignals RAMDASKUMARESAN Abstract-The zeros of a linear prediction-error filter polynomial for a class of deterministic signals, that are a sum of samples ofM exponentially damped/undampedsinusoids,isstudied. It isassumed that N samples are available for processing and that they are uncorrupted by noise. It is shown that the exponent parameters of the M signals can be determined from M zeros (called “signal zeros”) of an Lth degree prediction-error filter polynomial (L > M ) if L lies in between M and N - M (N - M / 2 in a special case). The rest of the L - M zeros of the filter polynomial switch are called extraneous zeros, are shown to be approximately uniformly distributed with in the unit circle, regardless of the type of exponentials in the signal, if the prediction filter coefficientsarechosen to haveminimumEuclideanlength. The results obtained provide insightinto theestimation problem with noisy data. I. INTRODUCTION The N samples of a squence y ( n ) , which are available for processing are assumed t o be composed of a sum of samples of M exponentially damped/undamped sinusoidal signals. That is
217
11. THE LOCATIONSOF SIGNALZEROS OF c ( Z )
The noiseless datasequence y ( n ) isgiven by (1). The g‘ be a exponentparameters sk areassumeddistinct.Let ( L + 1 ) X (1) vector of coefficients ofa polynomial C(z).
k=O
“T” denotes matrix transpose. Let us arrange the data samples y ( n ) in the form of an ( N - L ) X ( L + 1 ) Toeplitz matrix A; (forward data matrix).
Proposition 1 : If the coefficient vector g‘ satisfies A i g ‘ = 0, and if L satisfies the inequality M < L d ( N - M ) , then G ( z ) has M of its L zeros at esk, K = 1, 2 , . . * ,M . This fact is easy t o see from the following observation. Any row of A ; can be written asa linear combination of the following M linearly independent vectors.
That is, the rank of A; is M as long as A; has at leastM rows, i.e., if ( N - I,) > M or L d ( N - M). The null space of A ’ h;fs dim,ension L + 1 - M . Since g ’ liesin the null space o r Af, f k g = 0, k = 1, 2, . . * ,M , if L d ( N - M ) . Therefore
In other words, the polynomial G(z) = zk=o L g k z - k has zeros a t esk: k = 1, 2;. . . ,M . However, L has to be, equal to or where ak and s k , k = 1 , 2 , . . . ,M are unknown complex num- greater than M , in order that the null space of A f has at least bers. The sequence y ( n ) could, for example, be the autocordimension one. If L = M , the null space of A ’ has a unique relation sequence of an AR/ARMA process. The original idea vector which is the same (with in a scale factor) as that given that the Sk’S can be obtained from the zeros of an Mth degree by Prony metbod. If L > ( N - M ) , A ; has less than M rows polynomial was due to Prony [ 11. In practice, however, the (the rank of A f < M ) and the property in (6) will not be true observed signal is often noise corrupted anda polynomial C ( z ) in general.Hence, L shouldsatisfytheinequality M dL < of degree L > M (assuming N > 2 M ) is needed to accurately ( N - M ) in order that G ( z ) may have zeros at esk, k = 1 , 2 , [ 21). estimatetheunknownparameters(seeforexample ..-,M. TuftsandKumaresan [3] -[ 61 havepresentedtechniques Now let us define a ( N - L ) X ( L + 1 ) Hankel matrix, called based on singular value decomposition (SVD) of linear predica backward data matrix A ; , in which theN data samples, y ( n ) tionequations to estimatethe sk’s accurately.Inthiscorare complex conjugated and arranged in reversed order comrespondence we study the properties of the zeros of an Lth pared to the forward data matrix defined in (4). “*” denotes G ( z ) . We assume complex conjugate. degreeprediction-errorfilterpolynomial that the data sequence y ( n ) is not corrupted by noise. The Y * ( L + 1>1 behavior of the zeros of the prediction-error filter in the noisy data casecan be better understood using the results of this investigation, It is shown that the parameters s k , k = 1 , 2 , . ‘ . , M can be obtained from the M zeros (signal zeros) of G ( z ) ,if I, satisfies the inequality M d L d (N- M ) ( ( N - M / 2 ) in a special case). Proposition 2: If ,a yoefficient vector g ’ satisfies the homogThe 1, - M extraneous zeros of C ( z ) are shown to be approximately uniformly distributed with in the unit circle, if the sum eneousequation A b g = 0, andif I. satisfiestheinequality M d L d ( N - M ) , the polynomial G ( z ) haszerosat e-Si, of G ( z ) are chosen of magnitude squares of the coefficients k=l,2;-.,M. to be minimum (with the first coefficient held at unity). The above is easy to show along the samelines as Proposition 1. Note that these signal zeros are the complex conjugate of ManuscriptreceivedAugust 5, 1981; revised January 7, 1982and July 29, 1982. This research was supported in part by a grant from the the reciprocal of the corresponding M signal zeros in Proposition 1. Office of Naval Research. If the sequence y ( n ) is composed of undamped sinusoidal The author was with the Department of Electrical Engineering, UniversityofRhodeIsland,Kingston, RI 02881. Heisnow withGTE sk s arepurelyimaginarynumbers)thenthe signal(i.e.,all Business Communication Systems, Inc., McLean, VA 22102. polynomials in Proposition 1 and 2 will have the same signal 0096-35 18/83/0200-0217$0 1.OO O 1983 IEEE
218
VOL. ASSP-31, NO. 1, FEBRUARY
IEEE TRANSACTIONS ON ACOUSTICS, SPEECH, AND SIGNAL PROCESSING,
zeros, since esk = e-’I. In this case, the rows of A b and A ; are spanned by the same vectors. Thus the two data matrices can be combined to f y m a (2(N-L) X (L + 1))forward-backward data matrix A f b .
Q = 1 + Ig1 I2 + 182 I2 + . . . + IgL I 2
1983
(1 1)
Note that the first coefficientis chosen t o be unity. Since the polynomialG(z) is the productof G I (z) and G2 (z), its coefficient are the convolution of the coefficients of G1 (z) and Gz ( z ) . That is, m
Then the followingis true. Proposition 3: Ifa vector g‘ satisfies A;bg’ = 0 and if L satisfies the inequality M < L < (N - M/2),then the polynomial G ( z ) has M of its zeros on the unit circle at esk = KSi, k = 1, 2 , . . . ,M. The data sequence y ( n ) is assumed t o be composed of undamped sinusoids. M / 2 stands for an integer if M / 2 is notan rounded to thenext largerintegervalue integer. The arguments in this case are i!entical t o those :i Proposition 1. The number of rows in A f b is 2(N-L). A f b should have at least M rows or 2 (N- L ) 2 M ,if G (z) should have the property given in (6) or (9). Hence, L should satisfy M < L < (N- M / 2 ) . This extreme value ofL = N - M/2 was used in [71 for estimating frequencies of closely space sinusoidal signals in noise. 111. T H E EXTRANEOUS Z E R O S OF G ( Z )
The prediction-error filter polynomial G ( z ) has I , - M other of the extraneous zeros called extraneous zeros. The lqcation zeros depends OF the clfoice of t h e g vector, which lies in the null space of A f (or A b or A i b as the caTe may ,be). Recall that the dimension of the null space of A f (or A b o r A i b ) is L + 1 - M. Hence,thy:are$tleast L + 1 - M different choices for 8‘. Since A f g = 0, g will also lie in the null space well. stlands for matrix complex conjugate of A ’ ~ A ’as transpose.Therefore, g canbechosenasaneigenvector I, + 1 - M zeroeigenvalues of cqrrFsponding to one of the A A . Thisapproach is relatedtoPisarenko’smethod[SI. The problem with these choices of g is that the extraneous zeros of the corresponding polynomial could be located any where in the Z plane and may not be identifiable from the M signal zeros without prior information or additional computation. I? thefollowingproposition y e shallsqecif;yaunique vector g (satisfying A f g = 0 or A b g = 0 or A f b g = 0 as the case may be) which results in a polynomial G(z) with desirable properties. Proposition 4: The L - M extraneous zeros of the Lth degree polynomial G(Z) = Z i = og k z - k fall inside the unit circle Iz I = 1, regarqless of the IM signal zero locations if its cqefficient vecforg = (go,, glf, g, , . . . , gL)T which satisfies A,ig 0 (or A b g = 0 or A f b g = 0) is chosentohavethefollowing constraints: g o = 1 and lgl 1’ + ( g , ,1 + lg, j 2 + . . . + / g L I* is minimum. It is assumed that L satisfies the inequalities specified in Proposit,ion 1 (or 2 o r 3 as the case may be). As long as g satisfies the homogeneous equation it follows from Proposition 1 (or 2 o r that G ( z ) has M signal zeros at e s k , k = 1, 2, . . . , M (or e--sk, k = 1, 2, . . . , M ) . We can factor the polynomial G ( z ) into
7
2)
where
The factor G 1 ( 2 ) has the M signal zeros which fall at known locationsasperPropositions1 o r 2 o r 3. We are,however, interested in the locations of the extraneous zeros, which are C,(z), when the coefficients of the zeros of the polynomial G ( z ) are chosen such that Q , defined below, is minimum.
grz =
k = .m
(12)
cn-kbk
where go = 1, bo = 1, co = 1, bi=O for i > M and i < O , and c j = 0 for i > L - M and i < 0. One can now imagine the above convolution, as performing linear prediction. This method of linear prediction iscalled“autocorrelation”method.Seefor example,[9] on the“datasequence”1, b l , b,, . . . , b~ [which are the coefficients of G1 (z)] using a “prediction-error filter” defined by the vector of coefficients of the polynomial C,(z), C, C = ( I , c1, c 2 , . . . , c~ - M ) ~ .Then the coefficients of the polynomial G(z), 1,g l ,g, , . . . , gL is the “error sequence” o u t of the filter c . Thus the problem of minimizing Q in (1 1) is the same as one of finding a “prediction-error filter” c that minimizes the error energy at its output, operating on a “data sequence” 1, b l , b 2 , . . . , b M . It is well known that this problem results in a set of Toeplitz system of normal equations (Yule-Walker type; see [ 91 for example). These equations are R(l)...R(L-MR(O)...R(Z,-M-
1) 2)
The elements of the Toeplitz matrix and the vector right side are given by the formula
on t h e
where 6 , are the coefficients of GI ( z ) ( b o = 1). The R (i) will be zero beyond lag M since the “data sequence”itself has only M + 1 samples ( b k k ) . The solution to the above linear set of equations minimizes Q, and determines the polynomial G , (2). It is well known [9] that such a “prediction-error filter” c has all its zeros inside the unit circle. Thus the I, - M extraneous zeros of G(z) will lie with in the unitcircle. The above discussion is valid for any two polynomials G 1 (z) (with known zeros) and G 2 ( z )(the zeros of which are chosen). Hence, the following statement: If an Lth degree polynomial with leading Coefficient unity has IM known zeros and the rest L - M zerosarechosen t o minimizethesum of magnitude squares of its L + 1 coefficients, then the L - M chosen zeros will have magnitude less than unity. This, of course, is a restatement of the autocorre!ation method of linear prediction. The coefficient vector g of the polynomial C(z), satisfying the constraints mentioned in Proposition 4, can be found by simply solving a linear system of equations. Siqce g ’ satisfies the homogeneous equations A;g‘ = 0 (or A b g = 0 or A;,g’ = S) and go = 1, we can rewrite the homogeneous equations asfollows. A f g = - hf
(15)
where g T = ( g l , g,, . . . , g~),g‘ = (1, g T f , and A; is partitioned as A>= (hfIAf). It is easy t o see that the columns of
IEEE TRANSACTIONS ON
PROCESSING, VOL. ASSP-31, NO. 1, FEBRUARY 1983 219
ACOUSTICS, SPEECH, AND SIGNAL
Af and the column vector hf lie in the span of the M linearly e-sk, e-2sk, * * . , e - ( L - l ) s k ) T , k = independent vectors, (1, 1 , 2, . . . ,M . The rank of Af is also M . Although there are many vectors which satisfy (1 6 ) , the one which has minimum Euclidean length (or minimum value for lgl l2 + lg2 l2 + * * + lgL 12) can be found as follows:
IO0 0
,
-
(16)
g = -A#h
where A # is the pseudoinverse of A [ 101. Nowcynsiderthesituationwhenwefindtwocpefficient vectors gf and gh using the data matrices A; and A b , respectively, on, the same signal sequence ,defined in formula 1. Let gi and gb satisfy A;gj = o and Abgb = 0, respectively, and 4. Let the corhave the constraints specified in Proposition responding polynomials be called Gf(z) and G b (z). Proposition 5: The Lth degree ( ( N - M ) > L > M ) polyG b (z), have the same L - M extraneous nomials Gf(z) and zeros. This result can be used in identifying the M signal zeros from the L - M extraneous zeros in the presence of noise in the data. From Propositions 1 and 2 it follows that theM signal zeros of Gf(z)arethereflections of the M correspondingsignal zeros of G b ( z ) about the unit circle. ‘f, e s k , k = 1, 2, * . * ,M are the signal zeros of Gf(z), then e - s k , k = 1, 2, . . . , M are the signalzeros of G b (z).Letthepolynomialfactorscorres onding to the signal zeros of Gf(z) and C b ( z ) be named G&z) and G f ( z ) , respectively. Due to this special property of the signalzeros thecoefficients of Cf(z) arethereversed, complex conjugated, scaled versionof the coefficients ofG f (z). That is, if the coefficients of Cf(z) are 1, b l , b z , . * . , b M , then the coefficients of G f ( z )are b & - l , . . . ,b : , I)/&$. Now, recall that the locationof the L - M extraneous zeros are determined by the solution of the equations in ( 13). The solution vector c depends on the values of R(i)’s. However, looking at the way in which the R (i)’s are computed in (14), it is clear thattheyareunaffectedbythe reversal andcomplex conjugation of bk’s. Andthescalefactorl/b&occuringin the coefficients of G f (z) affects both sides of (13) equally, leaving the vector c unchanged. Therefore we conclude that the L - M extraneous zeros of Gf(z) and G b(z) remain the same while the signal zeros are reflected in radius about the unitcircleasthe predictio: direction fs :hanged (i.e.,the homogeneous equations A;gf = O and A b g b = o are used) on a given data sequence. esk and e-’E When sk are purely imaginary (sinusoid case), are the same. Hence, the polynomials Gf(z) and Gf ( z ) (the factors of Gf(z) and C b ( z ) , respectively) are the same. Thus the polynomials Gf(z), G b ( z ) [andalso G fb (z)]havethe same signal and extraneous zeros and are, in fact, identical.
(bs,
I v . SPREADOF THE EXTRANEOUSZEROS In the last section we showed that the L - M extraneous zeros of G ( z ) fall,insidetheunitcircle if thepolynomial coefficient vector g obeyed certain constraints. We shall now give aqualitativeargument to showthattheseextraneous zerosareapproximatelyuniformlydistributedaroundthe inside of the unit circle. We also give an example to illustrate this point. In ( 1 3) the R (i)’s are zero beyond lag M . Therefore one can imagine the R ( i ) sequence to be a legitimate autocorrelation sequence of a moving average (MA) process. This hypothetical MA process would be the output of all zero filter excited by M signal whitenoise,thezerosofwhichcoincidewiththe zeros. Further, wecanimagine thatthe“prediction-error process,invoking the filter” c is attemptingtowhitenthis well known whitening property of the autocorrelation method of linear prediction [ 91. In fact, an infinite order filter c will c atberequiredtoperfectlywhitenan MA process.Since tempts to whiten the MA spectrum (which is the Z transform
(C>
Fig. 1. (a) The spectrum of a hypothetical MA process generated by a filter Gf(z) with three (signal)zerosat eS1, esz and es3. (b) The reciprocal of the magnitude of the filter transfer function G%(z) which attemptsto whiten the spectrum shown in (a) L = 5 0 . (c) The zeros of Gf(z) ( = G f ( z ) G%(z)). The locations of the signal zeros are shown by “x” and the extraneous zeros are shown by dots. L = 5 0 , L - M =41. of the sequence R ( i ) evaluated on the unit circle) it has to havezeros[theextraneouszeros of G(z)l aroundtheunit circle wherever there is concentration of power, i.e., wherever the Z transform of the sequence R ( i ) does not have a (signal) zero close to the unit circle. We now give anexampleinwhichthesignalsequence is given by the formula y ( n ) = e S l n + eSzn + e S 3 * , n = 1, 2, . . . where s1 = (0.1 +j2n(0.35)), s2 =j2.rr(0.5 and s3 = (-0.3 + j2n(0.65)). The signal zeros of G f ( z ) (which can be found by solving A ; g j = 0) are located at e s l , esz and e s 3 . The factor Gf(z) of G f ( z ) is (1 - e s k z - ’ , ) . Theplot of corresponding “MA spectrum” ( I Gf(z), z = e l w 1’) is shown in Fig. l(a). L was chosen as 50. The factor G{(z) of Gf(z) corresponding to the L - 3 ( ~ 4 7 extraneous ) zeros was computed. This was done by finding the vector c (which gives the coefficients of G.f(z))bysolving(13).Themagnitude of thereciprocal of G [ ( z ) evaluatedontheunit circleis plottedinFig. l(b). Fig. l(c) shows the 47 extraneous zero and 3 signal zeros of G f ( z ) . Notethattheextraneouszerosareapproximately
220
IEEE TRANSACTIONS ON ACOUSTICS, SPEECH, AND SIGNAL PROCESSING,
uniformlydistributedwithinthecircle.However, if there is noise in the data the extraneous zeros tend to fallcloser to or outside the unit circle and are usually observed as spu[ 111 ). rious spikes in spectral estimate plots (see for example However, in the caseofSVDbasedmethods(3)-(6),the to emulatethebehaviourpredicted extraneouszerostend here. This is due to the effectivesignal to noise ratio enhancement achieved by SVD methods. If the data sequence is composed of expon:ntially damped sign+, using thebackwarddatamatrix Ab (i.e., solving A6gb. = 0 to find the polynomial G b ( z ) ) gives a mechanism to discriminate the signal zerosfromtheextraneousones, since the signal zeros fall outside the circle (see Proposition 2) andtheextraneouszeros fallinside t h e circle. Thisfact is utilized in [61.
VOL. ASSP-31, NO.
i~
1, FEBRUARY 1983
I
TAPE ~~~
~
PRINTER
~
~ -
~
~
~
~
M
I
1
Fig. 1. A block diagram for the recognition system used in this work.
I. INTRODUCTION Thisworkreportsastudyofmachinerecognition of the Cantonese digits using bandpass filters. The study of human speech and its information pattern is a complicated and someREFERENCES what difficult process. One of the main difficulties seems to lie in the fact that human speech is not stationary and that [ 11 F. B. Hildebrand, Introduction to Numerical Analysis. New variations exist even for a given speaker when the utterances York: McGraw-Hill, 1956, ch. 9. aredeliveredatdifferenttimes.Averyaccurateprocessing [2] S. M. Kay and S. L.Marple, “Spectrum analysis-A modern perspective,”Proc. ZEEE, vol. 69, pp. 1380-1419, Nov. 1981. methodneednotalwaysbeanadvantage,andforpractical [3] R. Kumaresanand D. W. Tufts,“Singularvaluedecomposition purposes, speed could be an equally important factor. For the and spectral analysis,” in Proc. 1st IEEE Acoust., Speech, Signal of bandpass filters has recognition of isolated words the use Processing Workshop Spectral Estimation, Aug. 1981, Hamilton, been fairly successful for a small vocabulary [ 11, [21. This Ont., Canada, pp. 6.4.1-6.4.12. work concerns the recognition of the Cantonese digits, a Chi[4] D. W. Tuftsand R. Kumaresan,“Frequencyestimationofmultiple sinusoids: Making linear prediction perform like maximum nese dialect used by most of the populace in Hong Kong. This Chinese dialect, unlike English, is monophonemic and is a true likelihood,” Proc. IEEE, SpecialIssue on SpectralEstimation, tonal language [ 3 ] . Thisstudywas divided intotwoparts. to bepublished. The first part examined the spectral composition of the Can[5] -, “Singularvaluedecompositionandfrequencyestimation by linear prediction,” ZEEETrans. Acoust., Speech, Signal Pro- tonese digits and assessed the suitability for machine recognicessing, to bepublished. tion. Acomparisonwasmadeontherecognitionscoresof [6] R. Kumaresanand D. W. Tufts,“Estimating the parameters of these digits with their English counterparts.. The second part exponentially damped sinusoidal signals and pole-zero modeling,” extended the recognition to consider two weighting schemes. ZEEE Trans. Acoust., Speech, Signal Processing,to be published. On-line/off-line and speaker-dependent/speaker-independent “Improvedspectralresolution 111: Efficientrealization,” [7] -, recognition were performed. Proc. IEEE, vol. 68, pp. 1354-1355, Oct. 1980. [8] V. F. Pisarenko,“The retrieval of harmonics from a covariance function,” Geophys. J. Roy. Astron, SOC.,vol. 33, pp. 247-266. 11. EXPERIMENTAL PROCEDURES [9] J. Makhoul,“Linearprediction: A tutorial review,” Proc. IEEE, Fig. 1 shows a block diagram for the hardware used in the vol. 63, pp. 561-580, Apr. 1975. recognitionprocess [ 4 ] . Thespeechsignalobtainedfroma [ 101 C. L. Lawson and R. J. Hanson, Solving Least Squares Problems. high-sensitivity microphone (JVC M5 10) triggers an amplifier Englewood Cliffs, NJ: Prentice-Hall, 1974. [ 11] T. J. Ulrych and R. W. Clayton, “Time series modelling and maxi- with 60 dBgainoccurringwhenthesignalexceedsapreset mum entropy,” Phys. Earth Planetary Interiors, vol. 12, pp. threshold slightly above the peak noise level. The signal then 188-200, Aug. 1976. passes through a set of twelve &pole bandpass filters spanning the voice band from 100-1 06 kHz. The filter outputs arehalfwave rectifiedandintegratedforapresetperiod of 0.5 s. These outputs are multiplexed, A/D converted to 256 quantization levels andfedto aM6800microcomputer.Forrealtimerecognition,anautomatic signal gain-controlcircuit is Machine Recognition of the Cantonese Digits used. Recognition was done in the microcomputer. Software includes provisions for the creation of the reference templates, Using Bandpass Filters storage in the microcomputer and the subsequent retrieval for H. L. KWOK, L. C. TAI, A N D Y. M. FUNG comparison t o test inputs. No pretest training was given to the sixmalespeakers.Theywereselectedrandomlyandwere aged between 21-35. All werenativespeakers of Cantonese. Abstract-This work presentsthe results of an experiment on machine The spoken digits could be either prerecorded or real-time proor recognition of the Cantonese digits (0-10) using spectral matching. A cessed,Thereferencetemplateswerecreatedusingthree two-step procedure was used in the normalization and weightingof the more repetitions of the digits. The filter outputs are not suitable for direct comparison as fiter outputs. It wasobserved thatthe scoresvariedfrom 55-100 For thispurpose,normalization is they are not normalized. percent. done by dividing each of the filter outputs by their sums. This removestheoverallamplitudevariationsforthedifferent Manuscript received May 22, 1981; revised November 19, 1981 and speakers. To account for the relative importance of the freMay 11,1982. The authors are with the Department of Electronics, Chinese Univer- quencyspectrum,twoweightingschemeswereconsidered. A ) considered the amplitude variaThe first scheme (scheme sity of Hong Kong,Shatin, New Territories, Hong Kong.
0096-35 18/83/0200-0220$01.OO 0 1983 IEEE
J
