MULTI-CHANNEL HOMOMORPHIC WAVELET ESTIMATION by MARK CHRISTOPHER B.Sc, The U n i v e r s i t y A THESIS SUBMITTED LANE Of V i c t o r i a , 1980 IN PARTIAL FULFILMENT OF THE REQUIREMENTS FOR THE DEGREE OF MASTER OF SCIENCE in THE FACULTY OF GRADUATE STUDIES D e p a r t m e n t Of G e o p h y s i c s And Astronomy We a c c e p t to this thesis the required as conforming standard THE UNIVERSITY OF BRITISH COLUMBIA November 1983 © Mark C h r i s t o p h e r Lane, 1983 In p r e s e n t i n g this thesis r e q u i r e m e n t s f o r an of British it freely available agree that for that Library s h a l l make for reference and study. I f o r extensive copying of h i s or be her g r a n t e d by shall not be Geophysics The U n i v e r s i t y o f B r i t i s h 1956 Main M a l l V a n c o u v e r , Canada V6T 1Y3 Date DE-6 (3/81) J/^. /?, /403 and of further this Astronomy Columbia thesis head o f this my It is thesis a l l o w e d w i t h o u t my permission. Department O f the representatives. copying or p u b l i c a t i o n f i n a n c i a l gain University the s c h o l a r l y p u r p o s e s may understood the the I agree that permission by f u l f i l m e n t of advanced degree a t Columbia, department or for in partial written ii Abstract Wavelet e s t i m a t i o n c a n be p o s e d a s a m u l t i - c h a n n e l information problem. convolution of homomorphic additive of be wavelet channel with of data an t r a n s f o r m maps t h e d a t a space. wavelet a Each The mapping may a l s o and impulses. estimated using sequence. A a c o n v o l u t i o n a l t o an effect In t h e a d d i t i v e averaging. i s modeled as t h e impulse from common partial space This is separation the wavelet can termed cepstral averaging. This thesis provides a available for synthesis components method the and homomorphic comparison i t s realization. for alternative reviews wavelet to cepstral The inferior solution to estimation is i s proposed the o r i g i n a l homomorphic produce averaging. c o n v o l u t i o n a l space. separation may t h e n this of e s t i m a t e s . for be u s e d principal data a as an o f n o i s e on t h i s cases which is a r e used i n estimated channel. best are an a l t e r n a t e components A wavelet to define of shows t h a t n o i s e may For these each techniques estimates i n which p r i n c i p a l components suite to and the proposed a v e r a g i n g . The e f f e c t components cepstral of method i s i n v e s t i g a t e d . The i n v e s t i g a t i o n cause p r i n c i p a l transform by Principal estimate from T a b l e of C o n t e n t s Abstract i i List of T a b l e s List of F i g u r e s v vi Acknowledgements viii I. INTRODUCTION 1 II. HOMOMORPHIC TRANSFORMS 3 2.1 Introduction 2.2 G e n e r a l i z e d 3 Superposition 4 2.3 The C h a r a c t e r i s t i c System 7 2.4 The I n v e r s e C h a r a c t e r i s t i c System 12 2.5 The L i n e a r 12 System 2.6 P r o p e r t i e s Of The Complex C e p s t r u m 13 2.7 C o m p u t a t i o n a l R e a l i z a t i o n Of The C h a r a c t e r i s t i c System 17 2.8 Summary III. 3.1 18 PHASE UNWRAPPING 19 Introduction 3.2 P r i n c i p a l 19 Value 3.3 I n t e g r a t i o n 23 Of The D e r i v a t i v e 24 3.4 F a c t o r i z a t i o n 28 3.5 Number T h e o r y ....28 3.5.1 C o m p u t a t i o n a l D e t a i l s 3.6 D i f f i c u l t i e s Of Phase U n w r a p p i n g 3.7 C o m p a r i s o n Of Phase U n w r a p p i n g 3.8 Summary 39 Techniques 40 42 43 iv IV. PRINCIPAL COMPONENT ANALYSIS 44 4.1 I n t r o d u c t i o n 44 4.2 P r i n c i p a l 45 Components 4.3 Some P r o p e r t i e s Of P r i n c i p a l Components 4.4 Summary V. WAVELET 55 57 ESTIMATION 58 5.1 The P r o b l e m 58 5.2 A S o l u t i o n 58 5.2.1 The S m o o t h i n g Function 62 5.2.2 The W a v e l e t 64 5.2.3 66 Principal Components v s A v e r a g i n g 5.2.4 M u l t i p l e Sequence 5.2.5 Single 5.2.6 Low P a s s I n p u t s 5.3 An A l t e r n a t e 5.4 P r a c t i c a l Sequence Data 69 Data 77 86 Solution 87 Considerations 90 5.5 Summary 90 VI. INVERTIBILITY OF THE HOMOMORPHIC TRANSFORM VII. DISCUSSION AND CONCLUSIONS 98 APPENDIX A - THE Z-TRANSFORM APPENDIX B - FOURIER TRANSFORMS 99 AS CHEBYSHEV POLYNOMIALS APPENDIX C - GENERATION OF A STURM SEQUENCE BIBLIOGRAPHY 92 .. 100 104 109 V List 1 . Wavelet Misfits of Tables 85 vi List 1. C a n o n i c Representation of Figures of a Homomorphic System 2. Representation of t h e C h a r a c t e r i s t i c 3. Representation of t h e 4. Complex F u n c t i o n 5. Principal 6. B r a n c h e s of a r c t a n g e n t 7. Smoothing F u n c t i o n 8. Complex C e p s t r a Due 9. An 10. Value System D?? i n t h e Complex and System D 11. W a v e l e t Plane Unwrapped Phase Function 20 22 30 64 to Impulses 65 66 Principal Components Estimates 12. M u l t i p l e 13. W a v e l e t by 8 12 Assumed W a v e l e t Discrimination 6 67 68 Sequence D a t a 70 Estimates ..71 14. N o i s y Data 72 15. N o i s e Inputs 73 16. Cepstra of N o i s e 17. Eigenvalues 18. Covariance 19. C o v a r i a n c e of 74 Input Matrix and Output of N o i s e Matrices Sequence D a t a 75 ..75 76 20. Single 21. S e g m e n t a t i o n and 22. Cepstra 23. Noisy 24. Wavelet Estimates 83 25. Segment's C e p s t r a 84 Estimation of Segments Single Sequence D a t a 79 80 81 82 vii 26. Low Pass S i g n a l s 87 27. Out o f Band E n e r g y 28. Time Domain 29. Principal 88 Components 89 Inversion of Cepstra 93 30. I n v e r s i o n of C e p s t r a 95 31. Inversion of C e p s t r a 96 32. I n v e r s i o n of C e p s t r a 97 vi i i Acknowledgement At for this point in a thesis, panegyrization. Therefore, of an u n t o l d m u l t i t u d e , My s u p e r v i s o r , managed to academic This he has expertise I will me by given who make r e s e a r c h the has labyrinth of of t h i s freely I am thesis. h i s wisdom, grateful f o r the me, by T a d , t o meet many what Levy, me w i t h topics myriad Ulrych, is people i s . particularly for sharing knowledge o f c o m p u t e r s and h i s p r o g r a m s . Shlomo provided the sharing I thank Brad P r a g e r , his vast Tadeusz through and p h i l o s o p h i e s . opportunities following the course t o the completion done i s allocated do s o . Professor guide research space Iain Jones and Kerry Stinson many s t i m u l a t i n g d i s c u s s i o n s dealt with in this about thesis, as w e l l as accompanies students others. The comraderie progressing which through the important. Accordingly, (Bandito) Cabrera, J e f f r e y . Without them system together and g l a d l y , I Don White, life thank is Julian Lynda F i s k and Bob would have been much my many duller. Christian Sutherland l u c u b r a t i o n s and p r o v i d e d has much survived needed balance and ix frame. I the I am deeply thank grateful. the U n i v e r s i t y Natural Engineering their financial and of B r i t i s h C o l u m b i a Research s u p p o r t i n t h e form o f Council and for fellowships. 1 I. INTRODUCTION In any i n v e s t i g a t i o n , gained which "Let The wavelet and i n s i g h t i s paramount. a l l t h i n g s be done u n t o original purpose estimation optimal i t i s the understanding linear by of of t h i s thesis combining system. investigation edifying." However, the 1 was to investigate a n o n - l i n e a r mapping and an i t has mapping, become the primarily system and an their interaction. Wavelet information e s t i m a t i o n c a n be p o s e d problem. convolution wavelet from can remains channel be u s e d additive In linear constant map the additive i s modeled impulse impulse this as the sequence. The sequence A n o n - l i n e a r homomorphic data from a varies transform convolutional partial to an or t o t a l and i m p u l s e s . space t h e problem We p r o p o s e Resolution to a c l e a r e r an the separate the wavelet the homomorphic St. Paul, Bible. while of the wavelet practice results. with of data T h i s mapping may a l s o e f f e c t techniques. further wavelet to channel. to the In 1 a space. separation to of Each channel a s a m u l t i - c h a n n e l common t h e use o f p r i n c i p a l yielded components some unexpected n e c e s s i t a t e d r e - e x a m i n a t i o n of t r a n s f o r m and p r i n c i p a l understanding with and i m p u l s e s . approach of t h e s e c a n be a n a l y z e d of t h e i r components. This led interaction. 1 C o r i n t h i a n s 14,26, K i n g James V e r s i o n o f t h e H o l y 2 Chapter II of t h i s t r a n s f o r m : i t s p r a c t i c a l r e a l i z a t i o n . of the phase c o n t i n o u s phase unwrapping unwrapping a l t h o u g h new u n d e r l y i n g i s i t Having a v e r a g i n g i s of a i s from be p r o v i d e d i n of used by a i s c o n v o l u t i o n a l wavelet IV in terms a n a l y s i s problem an r e l a t i o n s . i n a d d i t i v e An o p t i m a l A d e t a i l . space, a l t e r n a t e a n a l y s i s . of termed Phase e s t i m a t i o n . component i s III. examined to and c a l c u l a t i o n T h i s mathematical i n t e r e s t , for the Chapter n u m e r i c a l d i v e r s e p r i n c i p a l Chapter i n c l u d e s i n homomorphic p r o p e r t i e s t r a n s f o r m . d i s c u s s e d a the development, F o u r i e r t h e o r e t i c a l mapped examines R e a l i z a t i o n i n c o r p o r a t e can developed theory, e s s e n t i a l l y may a p p r o a c h , method and t h e s i s This i s i n f o r m a t i o n e x t r a c t i o n . The a p p l i c a t i o n e s t i m a t i o n c o r r u p t e d poorer space. then by d e a l t by be Recent t r a n s f o r m ' s to This components to i n V. case Chapter For p r i n c i p a l a v e r a g i n g . To a v o i d be used i n the c h a n n e l of data a wavelet d e f i n e s e p a r a t i o n . a best procedure l i t e r a t u r e i s has i n v e r t i b i l i t y . c a s t Chapter from d i s c u s s e d doubts VI wavelet b r i e f l y y i e l d problem, may be components may these on data c o n v o l u t i o n a l e s t i m a t e i n of may t h i s o r i g i n a l P r i n c i p a l e s t i m a t e a l s o the components may homomorphic used e s t i m a t e s . c o n c e r n . each p r i n c i p a l n o i s e , than components For w i t h a d d i t i v e e s t i m a t e s p r i n c i p a l found i s of i n d i v i d u a l Chapter the V. homomorphic a d d r e s s e s t h i s 3 II. 2.1 HOMOMORPHIC TRANSFORMS Introduction A problem which a r i s e s i n s i g n a l p r o c e s s i n g separating signals The of theory solution. signals be w h i c h have been combined mappings Consider a into different separated and space, y i e l d i n g signals are proven to transformation regions separate input separated, useful which mapped maps techniques original The p r o b l e m t h e n becomes t h a t o f such a transform individually. transform The exists One theory of (1965). may transformations thus c a l l e d The a combining must to homomorphic of between input was finding an vector (Lipschutz, space. case formalized algebraically and output that systems. superposition. by an appropriate c l a s s e s of nonlinear by systems These linear spaces and a r e 1974, p. 123). t h e r e l a t i o n s h i p between s i g n a l s i n the input them i n t h e o u t p u t i n the each o f homomorphic superposition defines f o r combining finding generalized represented i f the may be more space than consider systems He c o n s i d e r e d homomorphic s y s t e m s generalized rule be and may do n o t know, a p r i o r i , by t h e t h e o r y w h i c h obey a p r i n c i p l e systems We approach i s provided Oppenheim i n the transform transformation. combined Even separation appropriate i t s back t o t h e o r i g i n a l amenable t o t h i s space. of in These r e g i o n s signals. various that i n a known way. be of a space. individually the not has is space and t h a t f o r Linear systems are a 4 special 2.2 c a s e of Generalized For are homomorphic the Superposition purposes considered multidimensional vectors. as indicates that x Consider a x (n) and a l l input t h u s may be of the discrete system d e f i n e d by the linear input the signals system and T, by a signals considered The function Then superposition discrete is a be 2 scalars. this discussion, be or and of to sequences x,(n) systems. notation variable x(n) n. t r a n s f o r m T. and b be Let constant definition, satisfies 2 and the system (1) It and relation the for will be, input signals. input signals with and the of a from a rule combined Then a by In signal rule for a for to rule this combine scalar for a case transforms the input with an combining combining satisfies linear combining combining a This addition. scalars. those used system H a for the rule of rules f o r combining a addition is irrelevant. l e t + represent j_ r e p r e s e n t scalars. with combining * represent signal suitability signals different Similarly, signals signals signals for : represent signal. output rules in general, and of t o d e f i n e more g e n e r a l combining Let multiplication, shows separation corresponding signals scalar transformation is possible signals the to 2 o r d e r of superposition transform the relation T[a-x,<n) + b - x ( n ) ] = a ' T [ x , ( n ) ] + b - T [ x ( n ) ] Note t h a t as output generalized 5 p r i n c i p l e of s u p e r p o s i t i o n i f H[a:x, (n)*b:x (n) ] = a ^ H t x ^ n ) ] + bj_H[x (n)] 2 Comparison (2) 2 of equations (1) and (2) shows that l i n e a r systems s a t i s f y a g e n e r a l i z e d s u p e r p o s i t i o n with the o p e r a t i o n s as m u l t i p l i c a t i o n and There are the o p e r a t i o n s * and (1975) and of such The systems the d e t a i l s w i l l The class vector + as a d d i t i o n . of mathematical r e s t r i c t i o n s were developed not be given systems by space to an output here. (2) can transformations v e c t o r space. from The to v e c t o r a d d i t i o n and those combining signals to multiplication. as a cascade representation. homomorphic and for with A l l systems of t h i s c l a s s can be of In three this combining the itself m thesis known we will The operation as for scalar represented the canonic consider that : represents a s c a l a r s with input s i g n a l s and output multiplication. i n t e g e r m, systems an system i n which * r e p r e s e n t s d i s c r e t e c o n v o l u t i o n + represents a d d i t i o n . through correspond be rules for combining s i g n a l s correspond scalars and Oppenheim s p e c i f i e d by equation i n t e r p r e t e d as a l g e b r a i c a l l y l i n e a r input and a v a r i e t y of systems which obey a g e n e r a l i z e d p r i n c i p l e of s u p e r p o s i t i o n . the formalism : In operation the : corresponds _^ special which rule i s best d e f i n e d represents scalar case that the s c a l a r i s an to the c o n v o l u t i o n of the s i g n a l with times. In the i n t e r e s t of r e a d a b i l i t y , l e t us denote the output operations + and j_ by + and • respectively. Then e q u a t i o n (2) becomes H [ a : x , ( n ) * b : x ( n ) ] = a«H[x,(n)] 2 The c a n o n i c r e p r e s e n t a t i o n (3) + b-H[x (n)] 2 o f t h i s homomorphic system i s shown 1. in Figure x(n) Figure 1 - Canonic R e p r e s e n t a t i o n of a Homomorphic 0 maps from a c o n v o l u t i o n a l P is a convolutional the homomorphic space to an a d d i t i v e space. system s p a c e t o an a d d i t i v e which space. transforms It is from defined a by relationship D [ a : x , ( n ) * b : x ( n ) ] = a«D[x,(n)] 2 L (1). System is a linear The s y s t e m D~ relationship system, 1 defined i s the i n v e r s e + (4) b«D[x (n)] 2 by t h e r e l a t i o n in equation of D and i s d e f i n e d by t h e 7 D" {a-D[x,(n)] + b«D[x (n)]} = a.x, ( n ) * b : x ( n ) 1 2 In general representation canonic solution, The is system is called separating convolution 2.3 and representation p r o b l e m of system the shows signals reduced therefore, lies to the characteristic that, once which that i n the have of D canonic system. is been The fixed, the combined by linear filtering. s p e c i f i c a t i o n of the characteristic decomposed space into to system an i t s constituent the product the z-transforms That is of (Appendix A ) . convolution Z[x,(n)*x (n)] 2 D serves additive z - t r a n s f o r m of linear of to transform space. transforms. two signals of the It from a can multiplicative that the i s equal to the individual signals = z[x,(n) ] - Z [ x ( n ) ] representation space. (6) 2 shown i n F i g u r e The 2 and system may be D be Recall i s a homomorphic s y s t e m w h i c h maps from a c o n v o l u t i o n a l steps. The C h a r a c t e r i s t i c System convolutional a specifies L. The It D the (5) 2 has a r e a l i z e d in to canonic three 8 x(n) Figure First, the z-transform convolutional space multiplicative function from function, into a continuous space. Next, maps this Tukey along [they] space t h e power i s defined spectrum Finally , the "although reduce by word strange confusion on a s t h e power s p e c t r u m of a signal. this continuous the inverse space. This cepstrum. proposed of a from t h e a d d i t i v e B o g e r t , H e a l y and spectrum. o t h e r words, was p r o p o s e d because considerably cepstrum 1 was X(z), in i n t o another space. t h e complex (1963) a s a p a r a p h r a s e additive space x ( n ) , from a l o g a r i t h m maps x ( n ) , i n another a d d i t i v e cepstrum with fourteen function, continuous function output, x(n), i s c a l l e d word sequence, t h e complex the m u l t i p l i c a t i v e i n t o a sequence, The maps an i n p u t X ( z ) , i n an a d d i t i v e z-transform space 2 - R e p r e s e n t a t i o n of the C h a r a c t e r i s t i c System D This, f o r use i n t h i s at first balance." sight, 1 The of the l o g a r i t h m of Oppenheim, Schafer and B o g e r t e t a l , 1963, "The Q u e f r e n c y A l a n y s i s o f Time S e r i e s f o r Echoes: Cepstrum, Pseudo-Autocovariance, C r o s s - C e p s t r u m and Saphe C r a c k i n g , " , i n P r o c . Symp. on Time S e r i e s A n a l y s i s (N.Y.: John W i l e y and S o n s ) , p . 209. 9 Stockham the (1968) added complex t h e word complex t o e m p h a s i z e t r a n s f o r m s and t h e complex logarithm t h e use i n the of system D. We have implicitly assumed maps from a m u l t i p l i c a t i v e s p a c e means t h a t i s must be d e f i n e d log[X,(z)-X (z)] have also assumed t h a t Thus X ( z ) must be x(n). unique The specification restrict is, X(z) inverse of of and X ( z ) In t h e c a s e t h a t circle, and (Appendix us d e n o t e z=exp(jw), the e v a l u a t i o n by X ( j u ) . X(JCJ) where j i s t h e Xi(jw) be sequence, of X(z). the We may the r e g i o n s the u n i t c i r c l e . That i n c l u d i n g the unit i s n o t a n a l y t i c on i n t e g r a t i o n may of the be used A). Let x(n) X ( z ) or X(z) z-transform. requires such t h a t include an a l t e r n a t e c o n t o u r o f (7) some x(n) convergence real This 2 z-transform definition t o be space. + log[X (z)] X ( z ) has an of logarithm that X ( z ) and X ( z ) a r e a n a l y t i c i n a r e g i o n circle. unit of x(n) such valid of a r e g i o n x ( n ) and convergence the t h e complex t o an a d d i t i v e = log[X,(z)] 2 We that real In terms implies i s an odd of r e a l = Xr(jw) imaginary u n i t , that function o f X ( z ) on and imaginary circle, parts + j-Xi(jw) j 2 = -1. The Xr(jo>) i s an even of the u n i t It further (8) requirement function implies that of w that and X(JCJ) 10 is i n co w i t h p e r i o d 2ir. periodic unit a). circle We on the f u n c t i o n of have X(jw) Comparison provided X(z) analyticity X(jco), the d i s c r e t e The X i (jw) (10b) - a r g [ x ( jco) ] continuous C o n t i n u i t y of Xr(jco) = l o g | X ( j w ) | i s assured, and no z-transform time t h e argument Fourier the on unit (phase) may be unit circle. of t h e complex e v a l u a t e d on circle, Note t h a t z e r o s Xi(jw) unit the of = arg[X(jco)] logarithm. the by Note t h a t circle, i s just transform. result of the must C o n t i n u i t y of problem with logarithm 2n arg[X(jcj)] zeros X(jco) on essential complex implies that (10a) the d e f i n i t i o n the (8) X r ( j w ) = l o g | X ( jo>)| become p o l e s of X ( z ) . d e p e n d s on (9b) be has of (9a) (9a,b) w i t h e q u a t i o n log|X(ja)| f u n c t i o n s of a. = log[X(jo>)] = l o g | X ( j < j ) | + j - a r g [ X ( j<j) ] of e q u a t i o n s Therefore, the of X ( z ) i m p l i e s t h a t X(jco) must be a c o n t i n o u s X(jw) X(z) Analyticity a multiple of affecting the phase's b e i n g uniqueness from complex added and analyticity the a m b i g u i t y function. t o the p r i n c i p a l of in defining Any integer phase without r e p r e s e n t a t i v e of t h a t function. 11 However, only one c o n t i n u o u s the function. will satisfy (7). The Also, in general, the property procedure that, principal wrapped in the value back is obtained. Phase u n w r a p p i n g in of v i e w , The d e f i n i t i o n dealt with elsewhere The requirements periodic x(n). implicitly former follows directly Xi(jw) the oddness value be n e g a t i v e . because X(z) of X(jw) which for of the Xi(jw) the be from t h e above continuity at therefore the origin having theoretical defined. phase are The input X(jo) odd a n d sequences has a phase mean v a l u e . requirements, latter requirements. i s zero These imply at the The noting follows of a z e r o magnitude on t h e u n i t from that origin. circle. is precluded Now, t h e i s e q u a l t o t h e mean v a l u e must be p o s i t i v e . a sequence x(n) factor i s e q u a l t o i t s magnitude and cannot The p o s s i b i l i t y h a s no z e r o s have (-n,n). continuous, must a l s o have a p o s i t i v e of X(joj) a unwrapped allowable and hence t h e phase o f X(ja>), the to range of From the a phase, the p h a s e c a n be u n i q u e l y i s the phase of X ( j w ) . and from t o b e , a major x ( n ) must be s u c h t h a t It that Xi(jw), that a t u=it. i s commonly t h e phase i s s a i d value's phase equation results of of thesis. constrain In p a r t i c u l a r , zero Thus in this in unique phase cepstra. and computation continuous specified and c o n t i n u e s the continuous this computation principal has been, representative terminology Thus t h e c a l c u l a t i o n o f complex point of This the only this numerical into be of a d d i t i o n of f i n d i n g termed phase unwrapping. fact phase can It a negative is worth mean v a l u e , noting value of x(n) that, t h e phase o f 12 X(jo>) a t CJ=0 i s an 2.4 The The Inverse (5). convolutional in Figure transform, complex unit multiple characteristic It transforms s p a c e and 3. The has the transform complex e x p o n e n t i a l exponential circle, so ir and of is system from The Linear linear unspecified the and defined in a d d i t i v e space to representation inverse i f X(z) by a shown cascading a z- z-transform. The i s a n a l y t i c on the exp[X(z)]. 3 - Representation of the System D" 1 maps from an additive System system L once D outcome of i s the is fixed. characteristics intended is 1 is realized The i n v e r s e c h a r a c t e r s t 1 c system D"' space to a c o n v o l u t i o n a l space. The discontinuous. D" an canonic i s u n i q u e and, Figure 2.5 is C h a r a c t e r i s t i c System inverse equation odd of the the only system in Its specification inputs, procedure. the H will which is depend on mapping D and the 13 2.6 P r o p e r t i e s Of The Complex We w i l l x(n) c o n s i d e r here with p o s i t i v e with index no s i n g u l a r i t i e s (Schafer, Cepstrum only n. finite length input Such s e q u e n c e s have i n the z plane sequences z-transforms and c a n be r e p r e s e n t e d a s 1969) X(z) where the (11) |a(k)|<1 unit and |b(k)|<1. circle The a ( k ) a r e t h e m. z e r o s and t h e b ( k ) a r e t h e m circle. The f i r s t factor represents factor, a A, i s a shift of zeros outside the unit 0 constant the inside input and the sequence second by m 0 positions. Let For us c o n s i d e r t h e f i r s t real sequences contributes is be the sign input is real, and i t s contribution shown (Tribolet, 1979) t h a t t h e s i g n of x(n). multiplication by +1 consistent if A with (11). i s positive i t 1975). If t o x ( n ) i s more c o m p l i c a t e d . o f t h e mean v a l u e by i n equation o n l y t o x ( 0 ) (Oppenheim and S c h a f e r , negative, can A two f a c t o r s satisfying the Thus, A It of A i s equal t o normalizing the o r -1 t o make A p o s i t i v e i s continuity requirement of Xi(jcj) at the o r i g i n . In to a term the second factor, t h e exponent x , ( n ) i n t h e complex cepstrum of o f z, m , 0 gives rise 14 Xi(n) (Ulrych, 1971). circle, m cepstrum. Thus remove t h i s be p h a s e component examined most 1. This zeros outside the unit a n d ic, (n) may d o m i n a t e t h e complex to shift i s equivalent the input t o removing with sequence t o the linear s a t i s f y i n g the X i ( j w ) be c o n t i n u o u s a n d p e r i o d i c . properties (Oppenheim relevant large many o f X(jco) a n d i s c o n s i s t e n t requirements that Several has i t i s important term. (12) 0 I f X(z) can 0 = - [ m «cos(irn) ]/n of the and S c h a f e r , are outlined The complex complex c e p s t r u m have 1975; T r i b o l e t , 1979). been Those below. cepstrum decays a t l e a s t as f a s t a s 1/n. That i s |x(n)| < C-|d /n| n where C i s a c o n s t a n t and d is - oo the <n< ao maximum of |a(k)| and |b(k)|. 2. circle) If x(n) i s minimum p h a s e x ( n ) c a n be c a l c u l a t e d 3. circle) the unit then x(n) and (no z e r o s o u t s i d e = 0 recursively, I f x ( n ) i s maximum p h a s e then n<0 directly (no z e r o s from x ( n ) . inside the unit 15 x(n) = 0 , and x ( n ) c a n be c a l c u l a t e d 4. be o f I f x(n) infinite 5. property results a cepstrum is A. elsewhere. Since Note Recall unit of that unit duration, sequence can nevertheless we a sequence fact that This a smooth s m a l l and p r o p e r t y 1 rapidly scales origin only at to at maps the to factor In t h e g e n e r a l be c h a n g e d case Consider a origin origin i s always and that x(n) case the c o n t o u r zero and this property of having zeros has finite includes the integration w e i g h t i n g the becomes a positive. sequences = x(n)-expfa-n] has transform yields of c o n v e r g e n c e w h i c h i s , x(n) is complex and by e x p [ A ] , t h e by e x p o n e n t i a l l y That unit cepstrum convolution have e x c l u d e d i n p u t the a the the inverse the sequence In the w i t h n. exp[A] the s c a l e from spectrum are |b(k)| addition w(n) n=0. and is (Appendix A ) . whose far the time sequence. circle. near x(n) are X ( z ) has a r e g i o n circle. X(z) that the applying which c o n v o l u t i o n merely follows. x(n) w i l l concentrate non-zero Then sequence of impulse at the scaling which magnitude from x ( n ) . whose z e r o s x(n) decays to to of |a(k)| A d d i n g an equivalent on t h e tends result Thus that directly duration, cepstrum from a s e q u e n c e 6. time smooth is circle. implies i s of f i n i t e complex is recursively, duration. The spectrum n>0 input 16 where a i s a r e a l constant. The sequence w(n) has the z- transform W(z) = X(z-exp[-a]) or W(z) = X ( z ' b ) where b=exp[-a]. radially Thus, t h e z e r o s by t h e f a c t o r x(n)=x (n)*x (n), t 2 b. is exponential equivalent sequences. phase all outside the t o the input, 2 of the c o n v o l u t i o n a l of sequences convolution Exponential of weighting exponentially weighted c a n be u s e d t o make a m i x e d minimum p h a s e o r maximum p h a s e by moving and p o l e s unit convolutional x,(n)•exp[a-n]*x (n)«exp[a«n] weighting sequence e i t h e r i t s zeros I f we have a o f X ( z ) a r e moved then x(n)*exp[a*n] = Thus, and p o l e s either circle. minimum and maximum phase Then inside the the s p e c i a l unit circle or p r o p e r t i e s of the s e q u e n c e s c a n be e x p l o i t e d . 17 2.7 Computational We of have d i s c u s s e d continuous evaluated a on leads That the an is, if In x(n) 0 an is will the x(n) i n c r e a s e N and reduce the has the (1969) cepstrum, causing effect may The be phase of infinite x(n). X(ja>). and must computation Schafer, the to 1975). calculated sequence duration, the aliasing shown original However, as input x (n) x(n). will 0 x(n) be decays sequence will error. that one the exponentially may input reduce c e p s t r a l sequence. weighting i t t o d e c a y more r a p i d l y . the This complex E x a m p l e s of this i n S t o f f a et a l . ( l 9 7 4 ) . most d i f f i c u l t complex c e p s t r u m This cepstrum, i n the zeros has of found transform. discrete Fourier (Oppenheim and exponentially weighting effect the when = "2.x(n+k«N) of appending by Fourier transforms use i n terms z-transform, X(jw). true i s of exponentially, aliasing We s y s t e m , D, be version Schafer the number o f p o i n t s in general, aliased i s just of cepstrum 0 Since, the samples x (n) i s the fact, representation. aliased cepstrum, x ( n ) , where N C h a r a c t e r i s t i c System c a n n o t compute c o n t i n o u s to evaluate to The characteristic unit c i r c l e , discrete transform the transforms. In p r a c t i c e we use R e a l i z a t i o n Of lies We part, i n the must in practice, in calculating computation of the compute samples of p h a s e of X(jo>) t o d e f i n e X i ( j c j ) as in equation argument the (10b). the or continous If we 18 use the Xi(jw) relation and X(jcu), we general, thus ARG[X(jco)] = t a n " [ X i ( j w ) / X r (jco) ] , Xr(jco) obtain represent the the p r i n c i p a l will not continuous where 1 the imaginary principal value suffice. value and r e a l of o f t h e phase arg[X(ja>)]. In i s d i s c o n t i n u o u s and Methods of computing phase a r e d e a l t w i t h parts of i n t h e next samples o f t h e chapter. i For the s p e c i a l minimum or case maximum that the input phase, calculated, without and (1975) p r o v i d e Schafer the aliasing, complex directly details sequence is cepstrum from x(n). o f t h e method either can be Oppenheim used. 2.8 Summary We have discussed w h i c h maps c o n v o l u t i o n s effect the properties complex this which w i l l of and problems be d e a l t w i t h homomorphic were input inverse in This convolved transform on t h e a l l o w a b l e logarithms computational non-linear into additions. separation of restrictions a transform signals. presented. may Several There are s i g n a l s due t o t h e u s e o f z-transforms. computing i n t h e next transform the chapter. There complex are logarithm 19 III. 3.1 UNWRAPPING PHASE Introduction The many problem of p h a s e u n w r a p p i n g c a n directions. elucidate Say X(w). the we by Xi(co) represent them d e f i n e the as a path the real the with The angle the a X(co) crosses the We axis. -it. wish changes to This of phase can defined 2TT. principal straight real from approach line axis. to will a as real real That is be L e t Xr(cp) and in and the , axis. ARG[X(CJ)], origin to is to X(a>) this angle The angles I f the path ARG[X(co)] of undergoes a jump i s jump i s from -it t o +ir. phase, arg[(X(w)], including generally As complex i s m o v i n g down, t h e the let 4. the downward. axis, unique adding may Traditionally be continuously across by and in Figure path from transform X(co) p a r t s of negative everywhere co of p h a s e o f X(co) I f X(co) a Fourier plane. plane negative 27r. define or decrease of imaginary I t i t i s m o v i n g up, increase multiples complex t r a c e s out the continuously be and negative jump d i s c o n t i n u i t y to geometrical function i n the upward and -ir c o i n c i d e on from +IT valued positive i s measured p o s i t i v e +ir and a a x e s o f a complex shown. as use sequence, x ( n ) , with f r o m z e r o , X(co) increases makes real i s a complex represented defined us approached concepts. have a X(co) plane Let be t o ARG[X(co)] the negative unbounded this an axis. and which real will arg[X(co)] i n t e g e r number of 20 Xi(uj) Figure As plane. u 4 - Complex F u n c t i o n increases, F(u) traces i n the Complex Plane out a path i n the complex arg[X(o>)] = ARG[X(a>)] + L(G>)-2TT where L(a>) i s an i n t e g e r can change only when Xr(w)=0 and phase, argtX(cj)], t h e a c t of if only finding at is w h i c h makes a r g [ « ] discontinuities i n ARG[•]. Xi(cj) sign. definition problem of is the complex is related to occur continuous the are.multiplied, Note that c a n be no u n i q u e p h a s e . intimately i s the imaginary functions The L(o>) i s c a l l e d p h a s e u n w r a p p i n g . arg[X(w)] also These c a l l e d t h e unwrapped p h a s e o f X ( C J ) and X(o>) = 0 f o r some u> t h e r e This changes L(u) continuous. part related logarithm. of the This unique follows the l o g a r i t h m of X ( C J ) . requirement their to phases that, will if be as It two complex added. In 21 general ARGtX.tw)] * + ARG[X (G>)] 2 A R G [ X . ( w ) + X (o>)] 2 but arg[X,(w)] + arg[X (o>)] = arg[X, 2 Note t h a t then ARG[ • ] occur at points ratio its plane. Then axis tan ['3 follows L(o>) i s projected _ 1 can approaches phase. searched for removed A technique integral (1977). This discontinuities d o e s n o t have t h e s i g n (Xr(w)=0). unwrapped the (-7r/2,7r/2) and t h e change An when X(w) because i n f o r m a t i o n of onto the crosses example o f t h e p r i n c i p a l unwrapped right the value p h a s e a r e shown in 5. Various and on a phase and t h e c o r r e s p o n d i n g Figure 2 [ X i (w)/Xr (OJ) ] 1 on t h e i m a g i n a r y a x i s . X i (w)/Xr(CJ) imaginary of i s defined c o n s t i t u e n t s and t h u s half X (w)] i f one u s e s t h e d e f i n i t i o n ARG[X(w)] = t a n " the + Schafer by a d d i t i o n to e x p l o i t Bhanu (1969) discontinuities them of have been put f o r t h its used algorithm i n the p r i n c i p a l o f n-27r u s i n g the a n a l y t i c derivative (1977) an t o determine the which phase, A R G ( C J ) , the a p p r o p r i a t e n. d e f i n i t i o n of the phase as was developed used a d i f f e r e n t by technique Tribolet to exploit 22 ARG[X(CJU)] —n _ -2n _ Figure 5 - Principal Value and Unwrapped Phase The p r i n c i p a l v a l u e of a phase (ARG) i s wrapped ( - T T ) , T T ) w h i l e t h e unwrapped phase ( a r g ) Is c o n t i n u o u s . this same d e f i n i t i o n . polynomial novel factorization approach, developed S t e i g l i t z and D i c k i n s o n of b y . McGowan into (1982) p r o p o s e d a s a means o f c o m p u t i n g t h e p h a s e . some and theoretical Kuc (1982) interest, using A has been number theory 23 applied to polynomials. These methods McGowan and Kuc the concepts will will Principal The involved. between The method of to elucidate Some d i s c u s s i o n o f t h e p r o p e r t i e s o f phase unwrapping w i l l be presented. Value method of phase value is due principal briefly. be d i s c u s s e d i n more d e t a i l f u n c t i o n s which a f f e c t 3.2 be p r e s e n t e d the p r i n c i p a l unwrapping to by Schafer processing (1969). The the relation v a l u e , ARG(w), and t h e c o n t i n u o u s phase, arg(OJ) , i s ARG(CJ) That i s , the continuous phase. Schafer's ARG(a>) c a u s e d appropriate starts be p h a s e , modulo algorithm addition discontinuities to {arg(w) }M0D2TT 2TT, searches or which at subtraction i s 0 for real increasing w i t h i n which a d j a c e n t considered continuous. CJ. is removes 2ir. of The reduced. Any m i s - i d e n t i f i e d the unwrapping a t l a r g e r More another recently, must the algorithm searches user for specify a l i e f o r t h e phase general, becomes, more a c c u r a t e a s t h e s p a c i n g between is principal them by The x ( n ) , and samples must In the for discontinuities in by t h e modulo o p e r a t i o n and at arg(O), threshold = the algorithm adjacent discontinuities will samples affect w. Poggiagliolmi et a l . (1982) a l g o r i t h m f o r p r o c e s s i n g the p r i n c i p a l phase. proposed It also 24 searches for discontinuities, phase v a l u e s between transform. They its first moment This is an suggest may numerical these to The methods solution methods shifting remove recalculates values the input of sequence such calculating part number of of Fourier the the that transform. linear phase discontinuities. Derivative to compute of an are adaptively the before reduce the I n t e g r a t i o n Of Various original i s zero attempt component and 3.3 the but the analytic attempts unwrapped phase use expression. t o do numerical a Essentially, integration with constraints. Let Then, as x(n) be i n the a s e q u e n c e and previous X(co) be chapter, X(co) we i t s Fourier transform. define = log[X(co)] (2a) = log|X(cu)| + j »arg[X(co) ] (2b) or X(co) where t h e complex chapter transform. given and logarithm equation The derivative is (2b) defined as represents of X(co) is well in a the proper defined previous Fourier and is by X* (co) = d{log[X(co) ]}/dco (3a) 25 where the X'(CJ) = log'|X(cj)| X' = X ' (CJ)/X(CJ) (w) prime + j«arg'[X(«)] (3c) denotes d i f f e r e n t i a t i o n Comparison of equations (3b) derivative of arg[X(cj)] is X'(CJ)/X(CJ) . I t c a n be shown (3b) and equal with (3c) to CJ. respect shows that the t o the imaginary part of (Oppenheim and Schafer, 1975) that arg'[X(cj)] The = {Xr (CJ)-Xi phase, a r g [ X ( c j ) ] , ' (CJ) - X i ( C J ) - X r ' (CJ) }/| X ( C J ) | c a n be u n a m b i g u o u s l y (4) 2 d e f i n e d as (5) o with initial condition previous chapter, Equation the (5) arg[X(0)]=0. arg[X(to)] must be establishes continuity. mean p h a s e d e r i v a t i v e i s zero. As mentioned continuous i nthe and Oddness w i l l odd. result i f That i s (6) o If this will i s not t r u e , make i t s o . independent computed removal of Implicity, the from phase i t h a s been assumed knowledge o f a r g ' [ X ( c j ) ] . directly linear x(n) using In fact, component that this the r e l a t i o n we have can be (Tribolet, 26 1977) Xr'(co) + where F{«} j-Xi'(co) d e n o t e s the Fourier arg'[X(co)] by using in (5) y i e l d s the equation The (5). trapezoidal as the that rule most for spacing decreases. in The This an error and and integrating as numerical integration of solution Again, adjacent approach a l s o suffers i s to t h i s method samples of use a improves arg'[X(co)] from e r r o r a t u>=u, w i l l integration computing phase. i s the integration. Thus, (4) straight-forward between in (7) unwrapped problem (7) transform. equations computational equation = -j.F{n-x(n)} propagation a f f e c t the result a t co>co, . A more integration adaptive the was developed integration. phase derivative rule sophisticated derivative, is then between, say, that these multiple of The are co . carrying (1977) and the must be as in (5) using computed p r i n c i p a l value. close i s added. to That some i n t e g e r L ( c o ) , where E a well as phase trapezoidal The each o t h e r the called The value of the constraint to within an is | (ARG[X(co) ] + 2 T T - L ( C J ) } - a r g [ (X(co) ] | < for out is calculated. This 2 to 2ir of p r i n c i p a l p h a s e ARG(co) as a r g ' (co) co, and i s t h e n compared integer Tribolet integrated phase values by method i s some s m a l l E positive number 27 supplied satisfy the trapezoidal procedure procedure obtain This compared (1977) it investigated phase to computed from F { n « x ( n ) } as w e l l of contained the (-7T,<T) There f o r increased Note t h a t , is value. in terms 2 n«x(n) starts propagation. An f o r co>co, . algorithm and p r o p o s e d and interpolation r u l e of i n t e g r a t i o n . as the f i r s t i s a trade-off accuracy to as He a l s o polynomial The points of the second d e r i v a t i v e 2 computation through Tribolet's as a b e t t e r the evaluation principal from e r r o r carry t h e use o f p i e c e w i s e derivative recalculating techniques. requires the co, investigated alternate co . is satisfied. Then t h e i n t e g r a t i o n suffers L(co,) w i l l of 2 a DFT o f x ( n ) and higher method a l s o in determining Bhanu to point between co, and (co -co,)/2 a n d the i n e q u a l i t y by t a k i n g another midway between co, and ARG[X(co)] and a r g ' [ X ( c o ) ] . needed. the i s applied iterated until is initiated c h e c k s v a l u e s o f L(co) t o I f no L(co) s u f f i c e s , i s calculated rule a t co=0 and p r o c e e d s error The a l g o r i t h m the i n e q u a l i t y . phase d e r i v a t i v e The the by t h e u s e r . of This of the phase, derivative o f more and initial i n the i n t e g r a t i o n . although the p r i n c i p a l value of the arctangent ( - T T / 2 , TT/2 in ) , by r e t a i n i n g t a n " [ X i (co)/Xr (co) ] c a n be a s s i g n e d . 1 a the sign value information unique within 28 3.4 Factorization Steiglitz polynomial acurate" are on p h a s e of thus the the proposed z-transform r o o t s and unit circle be i s equal algorithm. saddle t o the of "reliable and r o o t s of a z- the product of may Fourier sum phase of of finding Dickinson They ignored then transform. the phases each d i p o l e , unambiguously. and the a use These d i p o l e s continuous p o i n t s on the Once t h e t o form t h e calculated Steiglitz as factored into (dipoles). The can be p r o b l e m becomes t h a t o f root-finding With r o o t s of a of and this complex u s e d a Newton-Raphson t h e p r o b l e m of grounds that these repeated rarely occur practice. 3 . 5 Number Theory McGowan continuous discussed In 1 the i t can the p r o d u c t product, polynomial. is the dipoles. method, t h e in known, of d e g r e e one evaluated e a c h of of (1982) method of p h a s e c o m p u t a t i o n . 1 polynomials The Dickinson factorization transform be and and Kuc (1982) p r o p o s e d a method of d e f i n i n g t h e p h a s e b a s e d on in detail linear important. and There are of the i s unstable. Steiglitz and Factorization," p. 984. theory. relevant proofs system c o n t r o l stability system number theory methods a system without theory will be presented. t h e p r o b l e m of which actually A particular The will determine finding method stability of which p a r t the of determining Dickinson, 1982, "Phase Unwrapping I . E . E . E . T r a n s . A.S.S.P., v . ASSP-30, No. by 6, 29 stability gives information about t h e number p o l y n o m i a l which a r e i n s i d e or o u t s i d e actually locating determine the continuous Consider discrete them. a real Fourier This roots the unit c i r c l e same method of a without may be u s e d t o phase. sequence transform X(co) of x ( n ) , n=0,1,...N-1 with the (DFT) = ^x(n)-exp[-j-nco] (8) to The p h a s e c a n t h e n be w r i t t e n a s for is arg[X(co)] = t a n " [ X i (co)/Xr (cu) ] tan i s the p r i n c i p a l 0<CO<VT. _ 1 1 [*] - 7 r / 2 < t a n " <ir/2, a n d L(co) y makes a r g [ » ] a c o n t i n u o u s unwrapping Figure 6. defined problem. Consider by F(co) singularities). F(co) at if goes t o - oo + oo without The singularities However, a branch singularities. point . arctangent, function Finding L(co) the arctangent the principal o f co F(co), branch, in L=0, t o have no c a n change b r a n c h e s t o L=1 o n l y i f , this can 'wrap' o n t o t h e L=1 point is change not F(co) may go t o i n f i n i t y branch Likewise, c a n change t o t h e L=-1 branches or e q u i v a l e n t l y , change which function e x h i b i t i n g any d i s c o n t i n u i t y . -co that i s t h e phase ( a n d a s s u m i n g Xi(co) I t can then arctangent of on = X i (co)/Xr (co) This F(co) goes t o branch. Consider L(CO).TT value i s an i n t e g e r function. a point + at zeros required and r e t u r n only at o f Xr(co). at these on t h e same 30 arctan(x) + Ln n L=1 J L=o J L=-1 J 0 n_ -20 r r r 20 0 X Figure The 6 - Branches function can of arctangent continuously Function change branches at InfInlt1es. branch. Thus, direction necessary Let of whether the sign to determine us or assume change that t h e r e a r e no z-plane. The product ratio Xi ( C J ) / X r (CJ) but has on z e r o s on X i ( C J ) «Xr ( C J ) no change o c c u r s . L ( C J ) the u n i t the singularities. change s i g n of X r ( c j ) . ) I f t h e p r o d u c t ' s L(CJ) is the sign and the information increased i s decreased by by is the circle in the same s i g n Consider as the a z e r o of X r ( c j ) . f o r CJ c l o s e sign This 0^CJ<7T. has as CJ i n c r e a s e s t h r o u g h Xi(cj) d o e s n ' t negative, provides X(CJ)*0 that that F ( C J ) changes L(CJ). constraint of X i ( C J ) «Xr ( C J ) not changes one. one. I f an If there positive opposite is sign (Note enough t o t h e from the no zero to sign sign 31 change by i s not changed. L(CJ) letting L ( C J ) as decrementing provides the s t a r t i n g required. has X i ( C J ) «Xr ( C J ) through the z e r o s can Pierce, become be 1963). that obtained Sturm's a polynomial the g e n e r a t i o n of polynomials imaginary it parts parts of to of,the of X ( C J ) and defined as A Sturm's in DFT the The and h e n c e of t h e s e c o n d k i n d , U ( n ) , as real in given and found terms of sufficient and phase. imaginary polynomials of the f i r s t functions of kind, COS(CJ) are p.11-24) T ( n ) = T ( n , c o s c j ) = cos(n«cj) U(n) They t h e unwrapped polynomials roots use t h e r e a l i n t e r m s of C h e b y s h e v Chebyshev and involves from two sequence The and polynomials sequence p r o v i d e s p r o c e e d s as f o l l o w s . are expressed of the sequence. Sturm Sturm real of this theory The method McGowan and Kuc to s t a r t to (Beaumont number recursively phase changes solution theorem arguments. generated ( S n y d e r , 1966, the s i g n sequence, a sequence express the second k i n d . T(n), two to determine L ( C J ) method the problem of finding or arg[0] =0 t h e number o f d i s t i n c t 1966). polynomials. information The degree, (Marden, convenient Chebyshev between of a Sturm decreasing Now theorem a r i s e s incrementing condition of Xr(cj). using p r o v i d e s a method o f f i n d i n g of The of and CJ=CJ. L(0)=0. value, unwrapping problem i s generated, in p r i n c i p l e , f r o m CJ=0 t o increase CJ L(CJ,) = U(n,coscj) = s i n [ (n+1 ) c j ] / s i n ( w ) nZO n£0 (9) (10) 32 Relevant recursion r e l a t i o n s between these polynomials are, for n£2, For n>1 there T(n) = [U(n) - U ( n - 2 ) ] (11a) T(1) = U(1) (11b) T(0) = U(0) = 1 (11c) i s the r e c u r s i o n relation U(n+1) = U ( n ) - U d ) - U(n-1) where it cos(co). set in is understood (Note that (Snyder, that U ( n ) and T ( n ) a r e f u n c t i o n s Chebyshev p o l y n o m i a l s 1966)). (I1d) are From t h e i d e n t i t i e s (8) c a n be e x p r e s s e d an (9) and Appendix Details B. on Equation ( 1 0 ) , X(co) (12) * --o where t h e a ( n ) a r e l i n e a r x(n). orthogonal as X(co) = "2,a(n)-U(n) + j • s i n (co) ^ b ( n) • U (n ) ft of combinations calculating (12) may a(n) of x(n) and be w r i t t e n and b(n) = - b(n) are given i n as X(co) = Xr(co) + j-Xi(co) (13a) X(co) = P(0,co) + j - s i n ( c o ) «P( 1 ,co) (13b) 33 X(CJ) where it is = P(0) understood p o l y n o m i a l s which are P(0) + j -sin(cj) -P( 1 ) that P(0) functions = P ( 0 , 6 j ) = Xr(cj) of (13c) and CJ. P(1) That represent is = S"a(n) - U ( n (14a) . C O S C J ) and P(1) Note = P(1,CJ) that the = 2.b(n) = Xi(cj)/sincj sign of P(1) is the -U( n , (14b) C O S C J ) same a s t h a t of Xi(cj) on 0<CJ<7T. P(0) generate and P ( 1 ) a Sturm where M ^ N - 1 . of are polynomials This sequence, sequence relationship is generated (Marden, 1966). Euclid's algorithm divisor of integers 1963). The relationship P(k-1) 1<k^M-1. The d e g r e e of P(k) or using In is the fact, for They to {P(0),P(1),...,P(M)} a polynomials function 'negative this finding polynomials a r e used a s e q u e n c e of Each polynomial merely where d e n o t e d by Sturm s e q u e n c e i s decreasing degree. The in U(n). the of CJ. remainder' relationship is g r e a t e s t common (Beaumont and Pierce, is = Q(k)-P(k) and Q(k) the p o l y n o m i a l s - P(k+1) are polynomials P and Q i s defined (15) of degree to be t h a t k. of 34 the h i g h e s t a d e g r e e Chebyshev P c o n t a i n s a t e r m U(K) sum, P i s of d e g r e e The finding recursion Q(k) and p o l y n o m i a l , U, where K in their i s the h i g h e s t sum. degree If i n the K. (15) p r o c e e d s , g i v e n P ( k - 1 ) P(k+1) s u c h and P(k), by that deg{Q(k)} = d e g { P ( k - l ) } - d e g { P ( k ) } > 0 (16) and deg{P(k+1)} < The e x i s t e n c e o f Q(k) Pierce, the 1963) long and, P(k-1) other in this algorithm of p o l y n o m i a l s . and is identical terms of l e s s e r the e l i m i n a t i o n P(k+1) a r e g u a r a n t e e d in fact, division (16) i s s a t i s f i e d and deg{P(k)} such that to that We the of Q(k)«P(k). of t h e s e terms. then P(k+l) w i l l if Q(1) and P(0) is be of d e g r e e d e g r e e N and less P(1) equivalent to I t may that term happen This i s , writing of that results (15) a s P(k-1) than P ( k - 1 ) . i s o f d e g r e e N-R, o f d e g r e e N-(N-R) = R so t h a t such degree identical. That P(k+1) = Q ( k ) - P ( k ) - and c h o o s e Q(k) highest degree are a l s o is (Beaumont Q(1)-P(1) is For example, t h e n we of take degree N 35 P(2) The highest degree but t h i s discarding divisor constant. the Sturm sequence sign Sturm fora fixed co. L(co) We These the greatest P(N-1) are given V(co), s i g n s V(co) i t i s unchanged. common = P(M) is a i n Appendix changes i n the f i x an co and c o u n t This to follows and i s p r o v e d require unwrap from the a r e p r e s e n t e d here with proofs. then .divisor of P(0) phase. and P ( 1 , c o ) sign[P(M,co) ], one, i n the L(co), the Specifically, of the up t o t h e f i n a l have no where P ( M ) and P ( 1 ) , c a n n o t Sturm change = Q(1)-P( 1 ) - P(2) common result. zeros i s the greatest for P ( M ) was g e n e r a t e d by E u c l i d ' s P(0) If below. to lead 0<co<7r, changes determine the properties two t h e o r e m s Theorem; I f P ( 0 , c o ) i s i n c r e m e n t e d by o f co, w i l l C. when a p p l i e d t o The number o f s i g n as a f u n c t i o n required Proof; P(M+1)=0. where from a d j a c e n t P ( n ) i n ( P ( 0 ) , P ( 1 ) , . . . , P ( M ) } . sequence, sequence k=M-1 t h e number o f s i g n P ( n ) have o p p o s i t e = V(co)„. will continues, is TO calculate terms The a l g o r i t h m M=N-1, If by t h e h i g h e s t other t h e o p e r a t o r V(co) w h i c h , gives function P(M) and sequence, otherwise until of t h e r e c u r s i o n define changes adjacent k P(1). and Details us some i s not guaranteed.) = Q(M)«P(M) of P ( 0 ) Let (Possibly Q, by i n c r e a s i n g P(M-1) Then - P(0) i n P(0) i s cancelled term i n Q(1)«P(l). term cancel, degree = Q(1).P(1) 0<CO<TT. algorithm: for common 36 P(1) = Q(2)-P(2) P(M-2) = Q(M-1)-P(M-1) P(M-1) = Q(M)-P(M) P(M)=0 f o r some CJ, then induction, P(M-1)=0, P(M-2)=0, hypothesis, P ( 1 ) and P ( 2 ) have If Therefore P(M)*0 - P(3) - P(M) this recursion P(1)=0, no common implies, P(0)=0. zeros f o r any CJ on 0<cj<ir a n d c a n n o t B u t , by on change by 0<CJ<7T. sign. Q.E.D. Theorem: V{P(k,cj)} zero of If and P ( i , c j ) P(0,CJ) does not change, as a f u n c t i o n P(0,CJ) Proof: Say V { P ( k , c j ) } then of zeros, CJ, except then at a . z e r o o f some P(k,a>), s a y algorithm have no common implies changes. CJ=CJ 0 . ,CJ ) 0 happen a t a i s , P(k,cj ) = 0. That Euclid's o (5)) (equation P(k-1 T h i s can only = -P(k + 1 ,CJ ) 0 or P ( k - 1 , C J ) -P(k+1 , C J ) 0 That 0 i s , members o f t h e Sturm s e q u e n c e s e p a r a t e d member have opposite signs. otherwise P(k-1,cj ) = P(k+1,cj ) P(J,CJ <0 O ) = 0 (This 0 f o r e v e r y J>k i n c l u d i n g 0 = 0 must be which by one o t h e r true because would J=M w h i c h c a n n o t be imply by the 37 previous theorem.) T h e r e f o r e V { P ( k , c j ) } c a n n o t change through a z e r o o f P ( k , w ) , 0<k<M. Q.E.D. This of t h e o r e m c a n be i l l u s t r a t e d {P(k-1) ,P(k) ,P(k+1 )} n e a r CJ=CJ . i f P(k) changes sign, By other possible inspection unchanged. {+,-,-} becoming {+,+,-} {+,+,-} b e c o m i n g {+,-,-} t o each P(n) cannot both contrary zero, the to sign adjacent 0 P(0) that individually. interest of sign 1). changes Note that 0 and P(1) a r e a l s o our h y p o t h e s i s . Note both that is t h e above P ( k ) a n d P(k+1) f o r some u> a s t h e d i v i s i o n zero zero, algorithm which is when P ( k ) i s a c t u a l l y c h a n g e s c a n n o t be c o u n t e d by c o m p a r i n g t h e s i g n s o f a n d compare the signs 1963). However, of we can ignore t h e t e r m s on e i t h e r This follows directly side from algorithm. l e t us consider the f i r s t term of t h e sequence, as w passes through a zero, CJ , of P(0,tj )=0. o t h e number t e r m s a s z e r o h a s no s i g n . Now P(0,w), becomes (i.e. t h a n one o f t h e P ( n ) a r e z e r o (Beaumont a n d P i e r c e , Euclid's i s t h e same {-,+,+} I f more imply changes {-,-,+} b e c o m i n g applies would a r e {-, + , + }. cases are we s e e be the signs f r o m + t o -, t h e s e q u e n c e {-,-,+} a n d t h e number o f s i g n The Say t h e y 0 Then by c o n s i d e r i n g 0 Since P ( 1 , O J ) * 0 by h y p o t h e s i s , can occur. o First, and V ( C J ) does n o t c h a n g e . P(0,OJ). three That i s s i t u a t i o n s of P ( 0 ) * P ( 1 ) does n o t change Second, P(0)«P(1) changes sign, from - t o 38 + and V(co) from having P(0)-P(1) is, decreases P(0) opposite changes P(1) and opposite by one. signs change By the p r e v i o u s under w h i c h V(co) equal the minus P ( 0 , co) V(co) that to keep t r a c k 2 increases it by one. t h e same s i g n theorem changes these to changes - having are the only V(co ) Therefore, from That is 0 f r o m + t o -, to +, as co 0 X i (co)/Xr(co). co Third, and from 0 t o co . of phase t h e same s i g n . number o f t i m e s P ( 0 ) « P ( 1 ) Because follows change can change. t h e number o f t i m e s increases to having from h a v i n g conditions P(1) i s , P(0) f r o m + t o - and V(co) signs. to That between = Xr(co) contains In f a c t , and points, = X i ( c o ) / s i n (co) t h e b r a n c h number L(co) = V(co). of the r e l a t i v e two P(1,co) V(co) of t h e may arctangent also be number o f b r a n c h c h a n g e s s a y co, and co , a s co g o e s 2 it used i n the f r o m co, to . Because the relation calculating V(co) the sequence Sturm (Appendix at C). does not Then evaluate t h o s e co f o r w h i c h we number of s i g n changes c a n c h o o s e a s many o r as recursive the is want methods. a purely P(n,co) the phase. and V(co) few co a s we is direct, We algebraic compute recursion f o r a l l n from 0 We i n t h e P a t e a c h co, and does n o t depend integration X(co) sampling X(co). involve which we between like. to M then c a l c u l a t e the V(co). We yielding The method on p r e v i o u s c a l c u l a t i o n s is not as d i d 39 3.5.1 Computational D e t a i l s McGowan and Kuc (1982) carrying out the e v a l u a t i o n general this compared to a s i n g l e Fast 1974, algorithm p 177). trigonometric of the of L ( C J ) will using be Fourier results sequence. of Transform from the Say we lengths N, N-1, N-2,..., must and the sequence o f p o l y n o m i a l s o u t a b o u t N /2 2 multiplications. operations. If compares w i t h about of these 2) B = N need N«log N 2 FFT. polynomials we such In f a c t , and sequence of 3 most, Thus we Each function. function need N evaluate coefficients. o f CJ we N /2 o f many evaluation a t one v a l u e trigonometric For B values (Brigham, computation i s m u l t i p l i e d by a t r i g o n o m e t r i c carry (power (FFT) 1. In slow when sequence o f , a t 2 to evaluate method. have an i n p u t N + (N-1) + (N-2) + ... + 1 o r a b o u t N /2 coefficient this program f o r computationally Then we c a l c u l a t e a Sturm polynomials a computer f u n c t i o n s and m u l t i p l i c a t i o n s i n t h e Sturm l e n g t h N. This present Thus, CJ we of evaluations B-N /2 such such o p e r a t i o n s . This operations 2 for a single we may u s e t h e F F T t o e v a l u a t e increase efficiency. Consider the polynomial P(k,cj) = a ( 0 ) - U ( 0 , c j ) where the U(n,cj) sum + a(l)-U(1,cj) = s i n [ (n+1 ) C J ) ] / s i n ( C J ) . +... Bringing sin(cj) out of yields sin (CJ) • P( k , C J ) = a(0)«sin[cj] + a ( 1 ) • s i n [ 2CJ] + . . . 40 w h i c h can be written as sin(co) 'P(k,co) = 0 « s i n [ 0 t o ] which will is reduce power of 3.6 the the 2, these will we Of Phase the accurately. can defined evaluation The As question when t h e will we naturally z e r o s of X ( z ) be well these are be closer unit circle, that zero. If amplitude becomes difference i n the (minimum circle can phase) be Poggiagliomli et is a 2 2 (Bhanu, on the if outside 1982). slowly phase a zero (maximum Clayton Thus, and the phase and any zero comes a rapid change unit is sequence speaking, vary a result. of circle, the phase introduces unit As we numerical Generally phase d e v e l o p s is That unwrapped incorrect detected. 1977; data. to which types will curve just al., as the and phase large FFT N log N. finite l e a d t o an It zero zero or if N expressions errors. easily the the However, f a r from the will near the s e q u e n c e x ( n ) , where have s e e n , determined. the input arises discontinuities to the T h e s e may a r e more s u s c e p t i b l e t o example, t h e p r o b l e m of a c c u r a t e analytic errors. Using Unwrapping only the For number a b o u t analytically. of computational will p h a s e of know x ( n ) be operations. operations consider want +... form of a F o u r i e r t r a n s f o r m . number of Difficulties We is, in + a(0)-sin[w] circle the undefined. The is just inside phase) the unit Wiggins, 1976; i f computations shift a 41 zero just result. across For zeros determination X i (co)/Xr (co) If both used The the and question Xr(co) serious inaccuracies could circle are cause the This unit small, small l a r g e changes circle. thesis A type of increases, i t s zeros zero not A l s o , the so be this c l o s e t o the able to place sequence number still of tightly probability unit i t on extra zeros sequence have sequence which will tend with are impulse train. t o become clustered that a given the c o r r e c t s i d e of i s convolved ratio. near evenly the increases that there w i l l circle a the as sequence circle. of in random v a r i a b l e s t h e n , the and ratio If are angle the a sequence. sequence v a l u e s length phase determined. i n the (1982) d i s c u s s s u c h independent of errors c o n s i s t s of a random Dickinson in i s because becomes p o o r l y the distributed problem becomes t h a t o f w h i c h t y p e s relevant in this and unit i n the a r c t a n g e n t of Xr(co) can Steiglitz circle, ill-conditioned. z e r o s c l o s e t o the be unit near is Xi(co) evaluation the short be algorithm the the a a will circle. sequence, i n t r o d u c e d and unit If finite above result holds. The circle (1977). question to For be of how incorrectly simple short 1.00001 were c o r r e c t l y about circle. 0.9 and c l o s e a zero 1.1 are must identified sequences, identified typically be was explored zeros with although to unit by Bhanu a magnitude magnitudes considered the of between c l o s e to the unit 42 3.7 C o m p a r i s o n Of Phase U n w r a p p i n g The comparison difficult. On p o o r l y conditioned data even a s h o r t sequence analysis. techniques has That 2 2 found t h a t their algorithms of initially data, padded w i t h took 2 typically Bhanu time. ran and zeros time o f a b o u t in example, was On t h e d a t a identical algorithm f o r an i n p u t used almost using for a as a this except f o r always took s e q u e n c e o f 90 p o i n t s T r i b o l e t ' s algorithm 0.02 ms on an Amdahl et a l . published results t o 2 5 6 p o i n t s , McGowan and while The Adaptation implemented 1977). power 2 adaptation. algorithm a T = 0(N ). Poggiagliolmi was does not and D i c k i n s o n algorithm yielded some i s of the order parameters. Tribolet's For ms of Using and t h e p r o g r a m Tribolet's to techniques requirements Steiglitz 2 (Tribolet, the amenable of v a r i o u s on i t s l e n g t h . i s T = 0(N log N). programs 0.048 not Sturm s e q u e n c e s , r u n t i m e Kuc's program. time. depend on e m p i r i c a l l y or a n a l y t i c a l l y . 2 the a l l work w e l l . will and time run i n T = 0 ( N « l o g N ) before computation less only algorithm McGowan and thesis complicated Tribolet, FORTRAN p r o g r a m will is i n t e r a c t i o n between t h e z e r o s o f evaluated t i m e depends on t h e d a t a FORTRAN techniques t i m e o f McGowan and Kuc's t e c h n i q u e 2 FFT t o e v a l u a t e N log N. they accuracy Computational been computational of is data their The unwrapping Thus, the comparison depend o f t h e i n p u t of data, set. must be e m p i r i c a l . The phase On w e l l c o n d i t i o n e d individual direct of Techniques Kuc's algorithm had an a v e r a g e r u n 470 V/8 c o m p u t e r . 43 3.8 Summary Phase unwrapping Iterative methods integration and provides direct phase. zeros a include polynomial In principle, near c a n be a p p r o a c h e d the unit search t h e phase between directions. algorithm, factorization. relation circle. f r o m many Number a time sequence determination adaptive theory and i t s i s unstable for 44 IV. PRINCIPAL COMPONENT ANALYSIS 4.1 I n t r o d u c t i o n We have examined t h e o r e t i c a l aspects o f t h e homomorphic transform as a p p l i e d t o c o n v o l u t i o n a l problems. The t r a n s f o r m changes c o n v o l u t i o n a l problems t o a d d i t i v e ones. developed be a p p l i e d . this for additive transform, wavelet channel common presented i n t h e next problem can be function each we o f n. of addition in anticipation have Write a suite {x,(n),x (n),...}. a We seek a way t o in an signal linear optimal between this Loeve, p. technique how case. is this I t thus application. The some further common signal extract the fashion. common in s e q u e n c e s by t h e i r of (Dhrymes, signal We may measure t h e the covariance. sequences we a r e l e d t o a t e c h n i q u e analysis that i s , i n general, different i n combination o p t i m i z e s t h e common s i g n a l component formalism Suppose 2 which seeking multi- of sequences, each sequence a component suite of i t s a shows general Using i s requested. to the principal the contains amount o f common By in sequences a The This chapter these each sequence. from problem. chapter. approached forbearance Suppose e s t i m a t i o n c a n be s t a t e d a s information provides a solution reader's p r o b l e m s may t h e n Techniques 1970, p. 5 3 ) . which known a s In f a c t , c a n be shown t o be e q u i v a l e n t t o t h e K a r h u n e n - Hotelling, or eigenvector 115). Because t h e r e e x i s t transforms a variety (Hall, 1979, of ways t o d e r i v e this 45 result, with 4.2 that in this for us b e g i n be we write development, Say although we dealt We define reflecting have a s u i t e of o f l e n g t h N. and p r i m e d e n o t e s M the will where o f an e l e m e n t be sequences, In v e c t o r the transpose the v a r i a n c e terminology the d e f i n i t i o n s {x , (n) ,x, ( n ) , . . .x.,(n)} = x ' ( n ) matrix. The literature, a f u n c t i o n of n and denotes a vector x^(n), t o the problem by d e f i n i n g s e v e r a l t e r m s . d e t e r m i n i s t i c data. sequence insight derived. of the s t a t i s t i c a l of t h i s each or t h e most Components reminiscent origins lends thesis will Principal Let is which of form underscore a vector of x ( n ) , say by var(xj) = where m^ 1/N ^[x (n)-m- i s t h e mean o f x ^ ( n ) , (1) ] 2 L defined by A/ m, We define the c o v a r i a n c e = 1/N ^x of two (n) (2) sequences x^(n), X j ( n ) as A/ c o v ( x , , x-) Equation C, (3) c a n = 1/N £[x,(n)-m:][xj(n)-fj] be e x t e n d e d to define of x ( n ) , whose e l e m e n t s a r e g i v e n the covariance by C.- = C O V ( X ; , X J ) . (3) matrix That 46 /V C = 1/N ^ [ x ( n ) - m ] [ x ( n ) - m ] ' where m i s t h e v e c t o r - of means m' = { m , , m , . . . m } 2 represents the outer product symmetric and be (Dhrymes, define this can 1970, p 3 ) . the operator (1),(2), can (Gelb, be 1974). positive compactness of random notation let (Note in probability variables.) is semidefinite E { « } t o be t h e sum 1 / N ^ J « } . operator C us that t h e o r y , as Then equations ) = E{[x- (nj-m^] } (5b) C = E{[x(n)-m][x(n)-m]'] (5c) the covariance us a t t e m p t = a'x(n), about t o maximize by t a k i n g where a linear measure variance. Thus we a t t e m t variance of y ( n ) . s e q u e n c e s t o be a where that of to f i n d is signal that a which then, of the s u i t e of the x ( n ) , say i s a vector the y(n) = a'x(n) signal a mean. combination 2 Our two t h e common s i g n a l a' = {a,,a ,...} coefficients. If of common t o b o t h , i n t e r m s of f l u c t u a t i o n s sequences (5a) 2 L m^ = E { x ; ( n ) } consider Let y(n) to and [ • 3 [ • ]' M a n d (3) become measure o f t h e s i g n a l of For with var(x- defined of v e c t o r s shown i s not the e x p e c t a t i o n we a r e n o t d e a l i n g We (4) in of constant y(n) is its maximizes from e q u a t i o n s the (5) 47 var(y) (Dhrymes, 1970) where C However, we s e e t h a t t h e r e solution, that b = sa and, structure with of a by any c o n s t a n t that If a scalar We have a p r o b l e m o f s c a l e . which i s important, s such It is not i t s magnitude. requiring constraint i s t o maximize v a r ( y ) a'a = 1. Lagrange m u l t i p l i e r s . a t o be o f u n i t = a'Ca T h i s c a n be done w i t h subject We d e f i n e t h e L a g r a n g i a n where X i s a L a g r a n g e m u l t i p l i e r . t o a and X by s e t t i n g to t h e method o f L = a'Ca + X d - a ' a ) This is a i s a'a = 1 . Now t h e p r o b l e m the of x ( n ) . s c o u l d be i n c r e a s e d t o make Thus we f i x t h e s c a l e by a r b i t r a r i l y length, matrix i s no maximum o f a ' C a . a fixed, large. (6) i s the covariance i t c o u l d be m u l t i p l i e d var(y) a r b i t r a r i l y the = a'Ca its partial (7) L i s maximized w i t h d e r i v a t i v e s equal respect to zero. yields dL/da = 0 = 2Ca - 2Xa (8a) dL/dX = 0 = 1 - a'a (8b) where we have used t h e r u l e s 1974; Hall, Premultiplying 1970) and equation 0 of v e c t o r d i f f e r e n t i a t i o n denotes (8a) by a' y i e l d s the zero (Gelb, vector. 48 a'Ca which, from multiplier (6), equates the v a r i a n c e i s the eigenvalue eigenvalues equation Thus, the Lagrange = Xa problem X and e i g e n v e c t o r s t o X, we must c h o o s e (10),-say the s o l u t i o n , (10) (Strang, 1980, p. 179-239) o f y t o be m a x i m a l , the l a r g e s t X which and t h i s satisfies X,, and i t s a s s o c i a t e d e i g e n v e c t o r , say a,. d e n o t e d by y , ( n ) , i s y , (n) = a , ' x ( n ) This i s >defined (Dhrymes, easily or finally t o be t h e f i r s t 1970, p . 3 4 ) . derived. with a. B e c a u s e we want t h e v a r i a n c e equal of y with (8a) i m p l i e s Ca is (9) X. Equation which = X The (11) principal mean and The mean, my,, i s my, = E{y,(n) } my, = my, = a,'E{x(n)} E{a,'x(n)} component variance of of x ( n ) y, are 49 my, = a1 ' m (12) Thus t h e mean o f t h e f r i r s t principal sum The v a r i a n c e o f y , i s , from o f t h e means o f x ( n ) . component i s the weighted (6) a n d (9), var(y,) Thus t h e v a r i a n c e o f t h e f i r s t the largest We other linear the signal linear properties principal first of t h e combinations they have. combination principal suite. of the Accordingly, of can y,(n). rewrite so t h a t component We x ( n ) , say By u n c o r r e l a t e d , may i s equal to t o contain the then sequences we maximal v a r i a n c e b u t i s u n c o r r e l a t e d component component eigenvalue of equation (10). expect most common (13) = X seek to the see what normalized y_(n) = a ' x ( n ) , which has with the f i r s t principal we mean c o v ( y . , y ) = 0. t h i s c o v a r i a n c e as ( u s i n g We equation (11)) cov(y,,y) E { [ y , ( n ) - m y , ] [ y (n)-my]'} cov(y,,y) E{[a,'(x(n)-m)][x(n)-m]'a} cov(y,,y) a 'E{tx(n)-m][x(n)-m]'}a cov(y,,y) a, 'Ca the c o n s t r a i n t investigate 1 may be w r i t t e n a s a , * Ca = 0 (14) 50 Proceeding a s b e f o r e we d e f i n e t h e L where X and u derivatives Lagrangian = a'-Ca + X ( l - a ' a ) + are to zero Lagrange ua/Ca multipliers. Setting partial yields dL/da = 0 = 2Ca - 2Xa + uCa, (15a) dL/dX = 0 = 1 - (15b) a'a dL/du = 0 = a , ' C a Premultiplying and equation (15c) (15a) by a' and u s i n g e q u a t i o n s (15b) (15c) y i e l d s a'Ca = X This y i s just is t h e same r e s u l t equal derivation we derived earlier: t o the Lagrange m u l t i p l i e r = X, a , by a' and u s i n g e q u a t i o n of From t h e p r e v i o u s have Ca , ' Premultiplying X. the v a r i a n c e (15c) y i e l d s 51 or 0 = a'a, which i s an e x p r e s s i o n eigenvectors. Now, are from zero." requiring equations ( 1 5 c ) and of (15a) by a' the yields = 0 (16), the f i r s t two terms i m p l i e s u = 0 and e q u a t i o n (15a) Hence from e q u a t i o n in general, = 0 (13) uX, Now, orthogonality - 2Xa,'a + ua,'Ca, ua,'Ca, and the p r e m u l t i p l y i n g equation 2a,'Ca However, (16) X, > 0 w h i c h = 0 becomes Ca = Xa and we we have returned t o the e i g e n v a l u e must c h o o s e t h e s e c o n d l a r g e s t associated principal normalized component eigenvector, i s defined as problem. eigenvalue, say a . 2 However, say X , Thus, 2 the now and i t s second 52 y (n) = a ' x(n) 2 and has the p r o p e r t i e s 2 that mean(y ) = a 'm 2 var(y ) = 2 2 X cov(y.,y ) 2 This procedure finding that variance but which components. linear be = 0 c o n t i n u e d (Dhrymes, combination i s uncorrelated Doing of = a^'x(n) i by maximal with the p r e v i o u s principal combinations o f x ( n ) , say i = 1,2,...M X^ cov(y^,y ) = 0 i * j ; J where t h e X^ The p 56) that var(yj ) = of x(n) 1970, having so g i v e s M l i n e a r y^(n) such can 2 C and y^(n) uncorrelated a r e the o r d e r e d a^ are t h e i r form a . linear set (largest associated of to smallest) eigenvalues orthonormal e i g e n v e c t o r s . sequences combinations of which x(n) are having mutually maximal 53 variance. In general y_'(n) = {y i (n) , y ( n ) , . . . ,y (n)} as 2 components of x(n) ^with we the y^(n) as define vector the of p r i n c i p a l i-th principal component. If of we w r i t e a matrix the transpose has the p r o p e r t y eigenvalues Equation y_(n). (17) t h e sum o f t h e v a r i a n c e s t h e sum of the (17) defines a linear alternately insight presented component derived problem. into o f y_(n), t h e variances transformation I t i s known a s t h e e i g e n v e c t o r , compression p. = Ax(n) the briefly. transform. as Hotelling, This between we its redundant In the s o l u t i o n This derivation transform. The For a d e t a i l e d some of x(n) elements. to utilization separate have information, order reduce This storage x(n) into d i f f e r e n t components, Karhunen- t o the optimal data provides additional ideas w i l l d e r i v a t i o n see H a l l a h i g h mutual correlation or s i g n a l , from x ( n ) can essential and e x p e c t is x(n) (1956). correlation fewer store or channel data. components, y_(n) = A x ( n ) , then of sequences. or transmission using be (1979, indicative common t o s e v e r a l space we w i s h t o r e p r e s e n t these x(n) transform 1 1 5 ) , Ahmed and Rao (1975) o r Kramer and Mathews Say of 1973). Loeve, or p r i n c i p a l be that o f C, e q u a l s (Ready and W i n t z , to a s t h e rows A we g e t y_(n) A of the e i g e n v e c t o r s transmit We discard those 54 remaining. as possible cannot are using x(n) a r e then the remaining be r e c o n s t r u c t e d components. exactly unless us d e n o t e t h e r e c o n s t r u c t e d error define E(n). total If A i s error, T, a l l of the components constructed when such discarding component by t h e i r to t h e sum o f t h e v a r i a n c e s the variances reconstruction eigenvalue contain error. principal principal component if one p r i n c i p a l to considered It follows that, o f C. obtained Thus the of the f i r s t - the proportion would have been o b t a i n e d to that using if However, regarding give the T i s equal r e p l a c e d y_(n). the r a t i o is o f y_(n) the eigenvalues will error mean v a l u e s , In p a r t i c u l a r , which also the of using the rest of components. Obviously, common We total elements information t o t h e sum o f t h e r e s t reconstruction first of these o f t h e y_(n) a r e j u s t eigenvalues the successive some o f t h e y_(n) a r e r e p l a c e d the 0 that results. the x ( n ) by x ( n ) and d e f i n e a s t h e sum o f t h e v a r i a n c e s o f transform x(n) well In g e n e r a l x(n) 0 the principal only as i n r e c o n s t r u c t i o n as E(n) = x(n) - x ( n ) . minimized the reconstructed retained. Let the The o r i g i n a l we want component, a l l of the as the r a t i o t o the u n c o r r e l a t e d the optimal this x(n). reconstruction component Thus this must using be t h e r a t i o may a l s o be o f t h e c o r r e l a t e d o r common s i g n a l signal or n o i s e . most in 55 4.3 Some P r o p e r t i e s Of P r i n c i p a l Let us behaviour examine of t h i s a 'noise' sequence. of sequences component us assume z e r o and t h a t that information for signal estimator. Let x ( n ) a s a common signal u(n), in general d i f f e r e n t about comparison Consider correlated combination the x(n) uses this u(n) a r e the case where weights. Ready a n d sum o f t h e x ( n ) and u s e i t principal other with s ( n ) . o f t h e x ( n ) w h i c h u s e s no identical the coherent t o the f i r s t w i t h each uncorrelated the component. u(n) are identically ( C O V ( U J , U J ) = f , f = c o n s t a n t ) and Then the covariance matrix elements Expanding f o r each the v a r i a n c e s of the u(n) a r e i d e n t i c a l . (1973) c a l l C^ s(n) = s ( n ) + u(n) t h e mean v a l u e s o f s ( n ) a n d A normalized linear Wintz cases t o i l l u s t r a t e the That i s , x(n) Let limiting t r a n s f o r m a s a common us w r i t e t h e s u i t e plus some Components = cov(x^Xj) = t h e sum E{[s(n)+u^(n)][s(n)+uy(n)]} yields c o v ( x ^ , x •) = v a r ( s ) + f C has 56 or, as v a r ( s ) i s a constant, cov(Xj-,Xj) = constant Thus C i s a constant zero eigenvalue. elements and coherent sum. the first all zero The f i r s t the principal rank h a v i n g eigenvector principal component. t h e same r e s u l t correlations from signal. with first of u n i t would will be This between n o i s e s makes them and w i t h the s i g n a l . will equal the into i f t h e n o i s e s were c o n s i d e r t h e c a s e where t h e n o i s e s each o t h e r identical concentrated In f a c t , follow. o n l y one non- have component A l l of the v a r i a n c e w i l l the Now matrix indicates that indistinguishable are uncorrelated Then = v a r ( s ) + dj-y »f where is d^ - a constant diagonal. the i s the Kronecker d e l t a matrix The f i r s t eigenvalue identical plus a another eigenvalue o f C,y = v a r ( s ) . and e q u a l We c a n s a y t h a t t h e principal component among a l l the component is signal and to A l l the other has added to the eigenvalues are f (Ready and W i n t z , concentrated in the 1973). first t h a t t h e n o i s e has s p r e a d uniformly components. principal the eigenvectors are orthogonal constant Thus C o f C i s X, = X + f where X i s to the constant principal equal and f i s a c o n s t a n t . coherent to the f i r s t The sum. first A l l the other and t h e r e f o r e t h e sum 57 of their signal elements i s in In these p r i n c i p a l the with each is the in is other. spread principal having this over of particular, the that much ratio of there negative Thus, the with signal and principal signal noise components. correlation of signal appear of was in the the noise in original the the signal to noise ratio. is no restriction on t h e This signal sum of component, could occur -s(n), or first eigenvalues to in a p r i n c i p a l inverted the in first may be i n v e r t e d . an no uncorrelated s p r e a d i n g out signal the had the lesser and t h e of x(n) all and elements x(n) a. is t h e most c o r r e l a t e d concentration be u s e d a s a measure of Note that c o m p o n e n t s , . w e c a n use t h e how there components. component, the p r i n c i p a l measure that component and t h e lesser principal Because of principal implies components. We c a n say uniformly successive This g e n e r a l ca-se n o i s e s - a r e c o r r e l a t e d first Signals all zero. if rest sign one or the as data. the if over some noise of a In can the more o f of led the to correlations. 4 . 4 Summary Signal linear measure this common t o combinations. leads matrix show where s e q u e n c e s can be e s t i m a t e d u s i n g Interpreting variance t o the c o v a r i a n c e m a t r i x . define the several principal information components lies. a s an information The e i g e n v e c t o r s and the of eigenvalues 58 V. 5.1 The Problem of seismic exploration The having smooth a (Tribolet, isolated wavelet 1979, Earth. The structure often called further The the u n i t sequence circle) consists of sequence i n the i s a consequence the geometry of the problem. form of a s u i t e wavelet of in (1980) general, or Earth's The time trace different. Telford et a l . data sequences i s assumed c o n s t a n t from impulses a r e , Treitel of to (See (1976) f o r information.) wavelet, problem of estimating the g i v e n the above d a t a . Solution Three L i n e s and and usually m o d e l s a wave p r o p a g a t i n g i n t h e T h i s c h a p t e r d e a l s w i t h the A impulse function, f a r from impulse the wavelet The t r a c e w h i l e the and (zeros the p h y s i c a l traces. limited w i t h an are values. available Robinson time spectrum impulse and is a p. 9 ) . non-zero Physically, 5.2 the r e c o r d e d data modeled as t h e c o n v o l u t i o n of a w a v e l e t sequence. are ESTIMATION ~ In t h e f i e l d often WAVELET methods Ulrych homomorphic wavelet e s t i m a t i o n were d i s c u s s e d by (1977): the Weiner-Levinson, methods. minimum d e l a y w a v e l e t requires of and the s e p a r a t i o n The former stationary two Wold-Kolmogorov methods r e q u i r e impulses while the o f t h e complex c e p s t r a of t h e a latter wavelet 59 and the impulse deconvolution (Wiggins, e s t i m a t e r , can It makes an the yield assumptions include as many a r b i t r a r y would like the be form suite of a necessary. We of impulse w r i t e the d a t a k or data as a wavelet and Ulrych, or simplicity, 1979). t o a v o i d making In g e n e r a l , we of these would like information while excluding For necessarily are a example, cepstrum's as i f • we length stationary, single I f the data sequences t h i s Each sequence x (n) entropy we this. segmented. w i t h an like as p o s s i b l e . assume a random, n o t If (Ooe e s t i m a t e of t h e w a v e l e t sequence can wavelet a priori to include sequence. would of minimum designed the entropy, basis. restrictions an a p r i o r i method estimate We arbitrary much have We an about sequence. on an The 1978), w h i l e not such assumption impulse to sequence. i s modeled as time are impulse sequence, already segmentation this in may the c o n v o l u t i o n the not be of a sequence. as = w(n)*i^(n) k=1,...M (1a) in vector notation x(n) where x(n), w(n) the wavelet and formed the by the and = w(n)M (n) (1b) i^(n) r e p r e s e n t , r e s p e c t i v e l y , impulse sequences. c o n v o l u t i o n of the Each same w(n) x(n) with a the data, is thus different 60 i(n). The convolutional amenable t o operation analysis^ transformation However, discussed from a c o n v o l u t i o n a l earlier, space i n (1) i s not by using i n a form most the homomorphic t h e e q u a t i o n s c a n be mapped t o an a d d i t i v e space. That i s D[w(n)*i. (n) ] = D [ w ( n ) ] + D [ U n ) ] (2a) or x ( n ) = w(n) where t h e c i r c u m f l e x now denotes i n t h e form o f a s u i t e common component, w(n), w i s h t o e s t i m a t e w(n) averaging the x(n) 1977). Averaging i g n o r e s any analysis, information it may be the a complex c e p s t r u m . of sequences, from the. x.(n). (Clayton is each a other component, the hand, improvement data. is partially, (Ulrych, i.(n). We done by Smith, technique which component data-dependent and ( U l r y c h e t a l , 1983). over a v e r a g i n g . from t h e i.(n) by of The uses Thus homomorphic separating, at low q u e f r e n c y windowing 1971). Before let w(n) a O t i s and Principal t r a n s f o r m p r o v i d e s us w i t h t h e p o s s i b i l i t y least data are been 1976; data-independent in The containing T h i s has and W i g g i n s , i n an o p t i m a l f a s h i o n an (2b) and a d i f f e r e n t information on + i_(n) discussing the p r i n c i p a l us examine a v e r a g i n g . component Assume t h a t , f o r each (P.C.) method, fixed n, the 61 elements of i_(n) can be c o n s i d e r e d a s random variables with identical probability distributions and zero mean. Then A, l / M S ^ i ( n ) w i l l tend to zero for large M (Otis and Smith, K' 1977). The e s t i m a t e o f w(n), d e n o t e d by w ( n ) , i s d e f i n e d a s k x ff M % (n) = 1/M ^x.(n) Ml w (n) e % The = 1/M equation of w (n) (4c) yields e Proceeding e + (4b) Af 2i„(n) (4c) the wavelet The inverse the wavelet as above, we + i) (n) ] i 1/M equals i s zero. w (n), K Al £_[u(n) (n) = w(n) estimated cepstrum in K (4a) cepstrum homomorphic i f the sum transform estimate. d e f i n e t h e P.C. e s t i m a t e of w(n) as = [ <£a x(n)]/2a e where the normalizing (5a) can a^ are Al M w (n) (5a) k the first t h e e s t i m a t e can be P.C. weights. shown as follows. reason for Equation be w r i t t e n as M Af M w (n) = [ ^.a w(n)+ S a i ( n ) ] / 2 a „ e or The K k k (5b) 62 w (n) = w(n) e + ^Ta which i s analogous zero, the e s t i m a t e d cepstrum It s h o u l d be wavelet's cannot If the scale noted in w (n) factor time we in this quefrency to 5.2.1 i t from The Smoothing expect in we to Again, not be w(n) (5c) we normalize cepstral have a finite i s an with they them. impulse a the number of (4c) o r ( 5 c ) , the (4c) or (5c) w (n) e i n v e r s e homomorphic w e g(n) In p r i n c i p l e the at the multiplicative set scale For the zero factor and Function equation where convolution, correlations. c e p s t r a to average with is cepstrum. arbitrarily impulse Then sum t h e complex c e p s t r u m ) . the as g ( n ) . the concerned o n l y by 6 of if original were l o s t . thesis zero (5c) the wavelet (4c) or from affecting In p r a c t i c e will property t h e examples u s e d prevent origin differ (see point equals equations will e (4c). t h a t , i n the r e c o v e r e d and sum origin, to equation s c a l e and be ^ ( n j / ^ a sum data to zero. i s not zero. and cannot Suppose t h a t , Write the sum becomes = w (n) + g (n) transform (n) = w ( n ) * g ( n ) i s t h e complex c e p s t r u m of g ( n ) . The estimate is 63 the c o n v o l u t i o n considered as estimate. result model of the t r u e wavelet a smoothing In of l i n e a r to i l l u s t r a t e the f o l l o w i n g example. of a time of l e n g t h 70 Poisson amplitude. mean The P o i s s o n and variance and t h e cepstral window o f l e n g t h length is appropriate function a a suite data spacing parameter of result i n subsequent approximates their a v e r a g e a r e shown decay indicating an to a as samples Then one w o u l d read sequences, impulses i s 11 ms, w h i c h smoothing hindsight examples. Figure function which 8. the l a c k of c e p s t r a l i s both are obtained origin). and will the 7 shows t h e windowing w i t h Note t h a t impulse. in Figure consider and G a u s s i a n , d i s t r i b u t e d i n ( c e n t r e d on t h e of the function function of s i x data are and low q u e f r e n c y 40 ms this the data the d i s t r i b u t i o n . corresponding averaging smoothing interval. in of be 1979). a 1 ms The distorts resolution consider ms. distributed this generate may i n v e r s i o n where t h e e s t i m a t e d (One may function with f o r ' p o i n t s ' . ) We i_(n), data the s i m i l a r i t y 1972; O l d e n b u r g and Samson, order g(n) which i s t h e t r u e model v i e w e d t h r o u g h a (Wiggins, 'ms' g(n). function ( I t i s worth n o t i n g from t h e f i e l d with by a Hanning The window be shown t o be the smoothing The complex c e p s t r a and Note the aliasing. exponential 64 (a) JL^_JL (b) 1 ms 60 Figure Six function 5.2.2 The Now wavelet, time random Impulse sequences (a) and their due t o c e p s t r a l a v a r a g l n g a n d w i n d o w i n g ( b ) . l e t us c o n s i d e r t h e w a v e l e t . shown w i t h 25 frequency and This using this than t h e same s c a l e . shown. i t s cepstrum t o 20 ms ms. comparisons amplitude Function smoothing Wavelet limited about 7 - Smoothing Note We i n F i g u r e 9. i t s cepstrum t h e use wavelet. Note power i t s low spectrum due spectrum pass n a t u r e . The wavelet i s quefrency justifies o f a 40 ms the the p r e v i o u s impulse The assume a mixed d e l a y of limited the higher Plotting wavelet to window f o r cepstrum's cepstra. is is i s to also Small changes to the h i g h to c o n v o l u t i o n or noise may dominate 65 0 -128 128 ms Figure 7. 8 - Complex C e p s t r a Due to Impulses These c e p s t r a c o r r e s p o n d to the impulse sequences of F i g u r e Note the e x p o n e n t i a l decay. the linear the p h a s e t o z e r o a b o v e some f r e q u e n c y and r e m o v i n g t h e l i n e a r trend to frequencies will phase c o m p u t a t i o n . that point. If T h i s c a n be a v o i d e d by the energy content setting f o r greater i s s m a l l , changes t o t h e wavelet and i t s cepstrum be s m a l l . 66 Figure 9 - An Assumed Wavelet The assumed wavelet (a) i s mixed d e l a y . I t s cepstrum (b) 1s time l i m i t e d and I t s power spectrum ( c ) Is low pass. 5.2.3 Principal Let feature Components v s us now compare t h e of the P.C. Averaging P.C. method This is input the convolution are demonstrated s e q u e n c e s o f random i m p u l s e s . the first variance cepstrum others. sequence. o f 25% o f t h a t of the noisy The P.C. the input estimate of in i s zero noise averaging. Figure 10. with The three h a s been added mean G a u s s i a n free sequence. i s noticeably different y i e l d e d weights A to discriminate our wavelet Random n o i s e The n o i s e of to is its ability between c e p s t r a . data method of to with a The from t h e [0.1,0.7,0.7] 67 Figure 10 - D i s c r i m i n a t i o n by Principal Components Input impulse sequences (a) and data derived by their convolution with a wavelet (b). The top d a t a s e q u e n c e has additive noise. C e p s t r a ( c ) show t h e n o i s y i n p u t ' s d i f f e r e n c e . The a v e r a g e and f i r s t p r i n c i p a l component of ( c ) are shown in (d). indicating wavelet shown successful discrimination of the noisy e s t i m a t e s due t o a v e r a g i n g a n d in Figure 11. the first input. The P.C. are 68 1 PC AVE. INPUT -30 0 ms 30 11 - Wavelet Estimates Figure Wavelet estimates derived from the f i r s t principal component (1 PC) and a v e r a g i n g (AVE.) of c e p s t r a of F i g u r e 10. The Input wavelet 1s shown at the bottom. To d e f i n e the o b t a i n a m i s f i t time and Thus the i t s e l f has e s t i m a t e 40. These a v e r a g i n g . for cepstrum, n u l l a had The between e s t i m a t e s h i f t e d p o i n t n u m e r i c a l e r r o r . d i f f e r e n c e wavelet 256 a are a s c a l e d i n d i c a t e m i s f i t . wavelet has an of 78 to In improvement of the we i s of t h i s versus of u n i t Since m i s f i t the sum and e s t i m a t e a of the f i r s t m i s f i t . m i s f i t i s wavelet e s t i m a t e zero m i s f i t the minimum the comparison the the e s t i m a t e s we squares of e s t i m a t e . v a r i a n c e , are 256 256 then c a l c u l a t i n g p o i n t s and example p r i n c i p a l P.C. The the the a l o n g . wavelet average component's method over 69 5.2.4 Multiple Let The us now data are wavelet 12). wavelet 15 s e q u e n c e s overlap parameter due Their The first and the second P.C. to respective P.C. bears are the first misfits by a d d i n g variance of 5% o f t h a t and misfits which first are two and The of the d a t a . P.C.s lend than the first. development resolve and this homomorphic quefrency The This average wavelet. credence noisy to our 280 than not In and paradox, and we shown Misfits 16. of the Note t h a t the f o r the n o i s e f r e e of the by our case, first examination. predicted fact, input data are e s t i m a t e i s i n the second examples. transform the to the estimates. bears c l o s e r is apparent 370. noise i s Gaussian with a The a r e 20, larger and the proper wavelet and for a noisy However, t h e b e h a v i o u r that a r e shown i n than resemblance wavelets method. generally i s unexpected P.C.s a smaller misfit random n o i s e . i s expected. P.C.s two 5.9 little (Figure The a r e 8.1, F i g u r e 14 a l o n g w i t h t h e w a v e l e t average 70 ms our i m p l y i n g an a p p r e c i a b l e the b a s i s of the p r e v i o u s r e s u l t proceed in length convolution. encouraging e x p e c t a t i o n s of t h e P.C. the c o n v o l u t i o n of of i s 6 ms the e s t i m a t e has results On sequences from a v e r a g i n g and 13. These g e n e r a t e d by impulse Poisson estimated Data a p p l y t h e s e methods t o a more g e n e r a l e x a m p l e . w i t h random The Figure Sequence It appears P.C. rather theoretical i t seems p a r a d o x i c a l . let us re-examine i t s interactions two w i t h the time To the and domains. behaviour of the homomorphic transform in the 70 (a) — w-^ (b) n- —S/W^VAA— -WyA/vv-^-— -W\Av~ —J\r—-^/v-^~->y^v— —'v-v.—J\^yv-A/^— 0 90 ms 90 ms Figure Fifteen wavelet (b). presence of random additive 12 - Multiple Sequence Data impulse sequences noise i s not w e l l t h e phase s p e c t r u m may be s t r o n g l y Clayton and Wiggins, possible that P.C. how n o i s e a f f e c t s cepstrum (a) 1976; and c o n v o l v e d understood. affected J i n and i n the second P.C. may 1983). 1975; It understanding The a p p e a r a n c e imply that However, (Buttkus, Rogers, a n a l y s i s can a i d i n our the cepstrum. with a is of of the wavelet n o i s e has induced 71 2 PC 0 -30 Figure 60 ms 13 - Wavelet Estimates Wavelet e s t i m a t e s due to f i r s t two p r i n c i p a l components (1 PC. 2 PC) and average (AVE.) of c e p s t r a of data i n F i g u r e 12. The input wavelet i s shown a t the bottom. a common first s i g n a l i n the c e p s t r a . P.C. and wavelet, appears tested noise the next most i n t h e s e c o n d P.C. by e x a m i n i n g signal appears common signal, This hypothesis the time-quefrencey i n the due t o t h e can be r e l a t i o n s h i p f o r pure inputs. In p r a c t i c e we c a n n o t noise. the This However, t i m e domain generate we c a n compare (input) with those (output). Poorly correlated matrix with s i m i l a r eigenvalues. common s i g n a l and t h e d i s t r i b u t i o n completely uncorrelated the c o r r e l a t i o n s of n o i s e i n in the sequences This quefrency have a domain covariance i s due t o t h e l a c k of u n c o r r e l a t e d noise of a over 72 (a) —\/\AA~^/\/\AIU (b) 2 PC 1 PC —sW\j\AJ^\A— AVE. —»~\J\~s\l\J\r— — a --V^AA^/^-—^— INPUT —w—y/V/-^-—'y^A*.— 60 -30 ms 90 ms Figure 14 - Noisy Data T h e d a t a o f F i g u r e 12 b u t w i t h 5% a d d i t i v e n o i s e ( a ) . (b) shows the wavelet estimates due .to t h e f i r s t two p r i n c i p a l c o m p o n e n t s (1 P C . 2 PC) a n d a v e r a g i n g of cepstra. The input w a v e l e t i s shown a t t h e b o t t o m . all t h e P.C.s. yield On t h e o t h e r (ordered) eigenvalues hand, h i g h l y c o r r e l a t e d which o c c u r s because the sequences c o n t a i n Thus, the eigenvalues changes in relative Consider Figure 15. irrelevant mostly of i n p u t s and o u t p u t s rapidly. common This signal. a r e i n d i c a t i v e of correlations. the s u i t e The n o i s e as s c a l e decrease sequences of s i x input noise s e q u e n c e s shown i s z e r o mean G a u s s i a n . f a c t o r s a r e removed in I t s variance i s i n the cepstra (Figure 73 256 ms Figure A suite 16). odd the each of Note t h e h i g h component 15 - random n o i s e of largest value. rapidly than amplitude of the cepstra due t o the phase. eigenvalues f o r the input (cepstra) (time). One must transform The c e p s t r a is normalized For be cautioned 17). T h i s against non-random decrease implies from f o r time this induced c o r r e l a t i o n dominates the f i r s t observation P.C. that implies component, the inputs that the wavelet that, although this is domain n o i s e . estimate the noise uncorrelated with may i n the It i s likely appears the inputs. may n o t be a s u i t a b l e measure o f s i m i l a r i t y domain more that the random quefrency noise comparison, inferring outputs strong by d i v i d i n g by t h e f o r the output (Figure the are plotted to a r e more h i g h l y c o r r e l a t e d t h a n has produced Correlation and two e x a m p l e s . The e i g e n v a l u e s outputs Inputs Inputs. same s c a l e a s t h e p r e v i o u s set Noise P.C. that The i n t h e second induce the wavelet a common cepstrum. 74 O •128 128 ms Figure The complex This follows We can covariance because gain matrices. cepstra 16 of Cepstra the noise of Noise Inputs of Figure t h e P.C.s a r e m u t u a l l y further insight For comparison, by 15. uncorrelated. examining the wavelet cepstral alone y i e l d s 75 a o - INPUT OUTPUT a a Figure I 5.0 1 J-o 3.0 EIGENVALUE * 17 - Eigenvalues The eigenvalues of input that those of the o u t p u t . This c o r r e l a t e d than the Input. a constant due to matrix input range of v a l u e s with alone both 0.633 0.314 0.400 0.079 -0.411 -O.130 0.314 0.359 0.241 -0.097 -0.365 -O.143 Figure 18 - (0.063). g r e a t e r than and Output 0.063. The c o v a r i a n c e i s shown i n F i g u r e signs. 0.40O 0.241 0.574 0.068 -0.219 -0.029 Several 0.079 -0.097 0.068 0.244 0.042 0.065 Covariance 18. elements than Matrix of l e t us examine due t o s i g n a l the Note t h e have a Noise the n o i s e c e p s t r a of wavelet cepstrum's the covariance those matrix -0.411 -0.130 -0.365 -0.143 -0.219 -0.029 0.042 0.065 0.496 0 . 166 0 . 166 O. 146 F i g u r e 16. correlation of the wavelet Thus c o v a r i a n c e s o f c e p s t r a due t o n o i s e be much l a r g e r Now Input ( d a r k l i n e ) d e c r e a s e more s l o w l y i m p l i e s t h a t t h e o u t p u t Is more The covariance matrix of These e l e m e n t s compare w i t h the of 0 . 0 6 3 . magnitude of the value noise of i 7.0 cepstrum alone can alone. matrices from the previous 76 example. For c l a r i t y , inputs w i l l the be c o n s i d e r e d . cepstra convolved noise, C o n l y the of Consider the impulse with the wavelet (Figure 2 0 002 0.002 -0.000 -0.001 -0.001 -0.000 0.009 0.006 0.009 0.008 0.002 0.009 19). noise, The i m p u l s e C,; 0.064 0.023 0.051 0.069 O. 104 0.061 0.066~ 0.071 0.065 0.065 0.061 0.072 0.009 0.008 0.016 -0.001 0.028 - 0 . 0 0 0 -O.OOO 0.012 -0.017 0.003 0.009 0.010 0.002 -0.005 -0.017 0.003 0.044 -0.008 19 - six Covariance from the impulses and with 5% 0 0.064 0.055 0.062 0.068 0.069 0.065 Figure first cepstra, C , generally -O.OOO -0.011 0.008 0.021 0.022 0.023 0.064 0.075 0.069 0.062 0.051 0.065 the the m a t r i c e s d e r i v e d -0.001 -0.001 -0.036 -0.046 0.014 0.007 0.039 0.035 0.035 0.092 0.02 1 0.022' 0.065 O. 127 0.075 0.055 0.023 0.071 0.006 0.016 0.016 -0.001 -0.005 0.006 from s e q u e n c e s , C „ ; from w i t h no 0 . 0 0 2 -O.OOO 0.114 -0.009 -0.009 0.017 -0.036 0.014 -0.046 0.007 -0.011 0.008 0.064 0.065 0.064 0.064 0.064 0.066 result 0.009 0.006 0.009 0.010 -0.008 0.015 (b) (C) Matrices The covariance matrices due to impulse cepstra (a), impulses convolved with a wavelet (b) and impulses convolved w i t h a w a v e l e t , w i t h 5% a d d i t i v e n o i s e (c). yield o f f d i a g o n a l terms covariance less ( 0 . 0 6 3 ) . • The with generally higher c o r r e l a t i o n s , noisy data, C 2 t the covariance convolved to the wavelet, than C,, h a s l e s s wavelet's matrix due t o structure a l l positive. The cepstral impulses than C 0 matrix has g e n e r a l l y s m a l l e r c o r r e l a t i o n s than and due C, 77 and more s t r u c t u r e , These observations Correlations the including due to matrix structure domain noise can negative can the than be interpreted impulse those dominate corrrelations. sequence due both to of as follows. contribute the wavelet. these i n the c o v a r i a n c e m a t r i x and strongly affect Thus, P.C.s amplify noise e f f e c t s in a sense, may the less P.C. to Time cepstral estimate. compared with averaging. The presence implies may that t h u s be 5.2.5 Having interactions can achieved an Clayton indicator and these This and tend to c a n c e l . converted the to can complex Wiggins, of 1976) 1977). The be of some transform windowing of the and P.C. and will, the envelope input basis of i s defined (Bracewell, i t s maxima c a n be i s smoothed s p a c e d maxima. used as the 271; an (Taner noise C e n t e r i n g windows on segments, w h i c h general, as p to reduce by the 1965, of t h e u n d e r l y i n g i m p u l s e s the t r a c e in sequence. multi-sequence done on The envelope merge c l o s e l y yields a signal the l o c a t i o n s maxima These of n o i s e dominance. homomorphic Segmentation of Sheriff, effects the weights proceed. be and will understanding ( U l r y c h e t a l , 1983). magnitude n e g a t i v e P.C. t h e c a s e where t h e d a t a a r e a s i n g l e segmentation. envelope indicator Data l e t us and cepstrum Sequence between Consider This the wavelet u s e d as an Single analysis, of b o t h p o s i t i v e destroy may the overlap. original 78 convolutional homomorphic window it model and transformation changes to a f f e c t the on into our impulse (Figure averaging and and marginally first variability the Let signal to noise segments' shown segments are 60, P.C. in were and to data window extract also 409. are convolved estimates are and 18 due Note t h a t the is dissimilar by c e p s t r a are the weights first and that P.C. impulse lack 22 of similarity cepstra to to d i d not structure. 5%, of a d d i t i v e 10% and 30% estimates are shown defined.) M i s f i t s that, of in have average in the c o r r e s p o n d i n g (Note noise noise are segments may wavelet, 25. being in Figure The the m i s f i t ' s the individual and shown to shown. estimate P.C. The Figure The maxima y i e l d used ratios.) cepstra have second with that expect before. sequence P.C.s 59, we the t h e P.C.'s examine t h e e f f e c t (Note was second us 23. impulse if data as the misfits, yielded Figure the proceed wavelet corresponding data we Envelope the c o v a r i a n c e now Once on have s i m i l a r indicate The two ms affect context. first The among t h e s e average's strongly P.C. The the weights and misfits better. the w a v e l e t . 500 The first respective average with the However, s e q u e n c e example. H a n n i n g window 21). effect the wavelet, sequence. ( F i g u r e 20). A 90 ms 1979). than t h e d e s i r e d form wavelet segments Their scale us c o n s i d e r a s i n g l e locations. are (Tribolet, c o n v o l u t i o n of a s i n g l e with a poorly understood a longer o n l y the been c o n v e r t e d Let has for shown i n different estimate Figure first 24. P.C. 30% this and The weights noise, 17 the v a r i o u s e s t i m a t e s are 79 80 (a) —-VWv— —Aw\/^— (b) 1 PC 2 PC AVE. -^/WWV INPUT ——JiA'yW— 60 •30 ms — ^ ^ - j y ^ A ^ — ^ — 90 ms Figure 21 - Segmentation and Estimation Segments of t h e d a t a i n F i g u r e 20 (a) and w a v e l e t e s t i m a t e s d u e t o f i r s t two p r i n c i p a l c o m p o n e n t s (1 P C , 2 PC) a n d a v e r a g i n g (AVE.) of c e p s t r a (b). T h e i n p u t w a v e l e t 1s shown a t t h e b o t t o m of (b). shown than i n Table the f i r s t noise, 1. The s e c o n d P.C. as the n o i s e l e v e l i t i s better than the becomes a b e t t e r e s t i m a t o r increases. average. In f a c t , Examination f o r 30% of the 81 -4— 1 pc ave. Jfb 0.212 0.248 H~ 0.242 •i^/w—— 0.227 \r$[s 0.22S |i~ 0.247 0.238 j A 0.215 — j \ f k - — 0.224 $ 0.250 0.242 A n •fr U . Jl OCA L DO 0.236 1 -f- 0.240 — 0.234 0.230 ^ I 0.243 — 0.227 . -30 Figure 22 - . i . 0 ms Cepstra . 30 of Segments The c e p s t r a o f t h e s e g m e n t s shown 1n F i g u r e 2 1 . The f i r s t p r i n c i p a l c o m p o n e n t (1 P C ) a n d a v e r a g e ( A V E . ) c e p s t r a are also shown. The weights used to form the f i r s t P . C . weights are shown n e x t t o t h e c o r r e s p o n d i n g c e p s t r a . 82 Figure A data bottom). 23 - Noisy s e q u e n c e w i t h 5%, Single Sequence Data 10% a n d 30% a d d i t i v e noise (top to 83 -30 0 60 ms Figure 24 - Wavelet Estimates The w a v e l e t e s t i m a t e s due t o t h e n o i s y d a t a of F i g u r e 23. T h e f i r s t (1 P C ) , s e c o n d (2 PC) a n d a v e r a g e ( A V E . ) e s t i m a t e s a r e shown f o r 5% ( a ) . 10% ( b ) a n d 30% ( c ) n o i s e . The I n p u t w a v e l e t i s shown a t t h e b o t t o m o f e a c h p a n e l . 84 (a) (c) (b) J\j\r- -v\/V- 0.35 — 0.21 -0.23 0.41 0.34 0.22 0.47 -0.16 —A— 0.27 0.28 0.11 0.12 0.20 0.19 -\/^A/ 0.28 —v^/v 0.26 -ju- 0.35 -0.04 0.16 0.24 —w*^- 0.20 0.26 0.22 0.29 0.24 -0.10 0.18 0.30 0.07 ~V— 0.00 0.10 0.17 0.24 0.17 0.22 0.27 0.02 0.28 -0.08 0.21 0.06 0.10 0.21 0.44 0.36 0.29 0.23 0.17 0.20 -30 0 -0.08 -0.07 —d\j^\r— 0.25 -30 0.32 30 •30 ms Figure 25 - 0 ms 30 Segment's O ms 30 Cepstra The c e p s t r a due t o the n o i s y d a t a of Figure 24. Inputs have 5% (a). 10% ( b ) , a n d 30% ( c ) n o i s e . The f i r s t principal c o m p o n e n t w e i g h t s a r e shown b e s i d e t h e c o r r e s p o n d i n g c e p s t r a . first P.C. variability from the weights increases, first to and m i s f i t s shows t h a t , the wavelet estimate the second P.C. as the weight's tends Note a l s o to that, shift f o r low 85 Wavelet Ave. 1 P.C. 2 P.C. 5% 10% 30% 39 18 75 19 74 134 190 56 32 noise (5%), the f i r s t than averaging. noise 1 - Wavelet P.C. Also, the Misfits y i e l d e d a b e t t e r wavelet f o r high a better estimate level, Misfits Noise Table yielded Estimate than P.C. noise (30%), averaging. method estimate t h e second A t an P.C. intermediate y i e l d e d a worse r e s u l t than averaging. I t may the m i s f i t there are affected as seem p a r a d o x i c a l i s lower many f o r 10% n o i s e factors setting Also, t h e phase t o z e r o relation the fact may be cautioned will, to compared against in general, the for this noise that not estimate. the estimate. the n o i s e wavelet, may a f f e c t t h e detract l e v e l , t h e P.C. adding However, Segmentation i s particular T h e s e do average inferring improve f o r 5% n o i s e . above some f r e q u e n c y for a fixed to than estimate, windowing may a f f e c t to the wavelet. that, f o r the averaged be c o n s i d e r e d . by n o i s e , and c e p s t r a l w e l l as the wavelet. noise's that, estimates One noise from must be t o an i n p u t 86 5.2.6 Low Pass In previous stabilized It Inputs may by setting seem p h a s e , we should do the cepstral averaging will adversely this the a the admitted is a logical common P.C. ignorance the magnitude. o f e n e r g y above u, this transform phase t o z e r o above some a n g l e l i k e w i s e with induce homomorphic that, having shifting affect CJ.. of This this would to zero phase. thing signal was to do. For However, i n the c e p s t r a which method of a n a l y s i s . Let us may examine further. Say a time sequence magnitude M(CJ) . equivalent to s e t t i n g T h i s may be boxcar large function, i s a sine cosine A-S(n) above CJ,. Fourier to a large u n i t y on 0<6j<co say A, f u n c t i o n and sine 1965). on The becomes Figure the log The values the original 27 magnitude larger the T h i s procedure is of A, domain call it this values wavelet above CJ, 10, 7, is a d d i t i o n of i s set to smaller our the function i s f u n c t i o n and s u b t r a c t i n g 13, Note t h e a peak v a l u e of a 0<CJ<CJ. and quefrency a s M(CJ) by on CJ,£CJ<.T) . bandpass set to v a r i o u s values a t CJ,. zero i s z e r o on shows t h e c e p s t r u m were d e r i v e d by value and 1 function with In t h e value. Fourier transform t h a t of with above CJ=CJI, i s negative CJ,<CJ<7T. to c o n v o l u t i o n with a sine A transform multiply log[M(cj)] First, i n F i g u r e 26. (Bracewell, A-S(n). a f u n c t i o n which value, modulated equivalent steps. (i.e. shown g r a p h i c a l l y boxcar log[M(cj)] a bandpass negative has S e t t i n g M(CJ) t o a s m a l l v a l u e done i n two S e c o n d , add a the obvious prevent this examples = with 0.62-7T. 3 and dominance of A « S ( n ) 0 from as A 87 F i g u r e 26 - Low Pass S i g n a l s The energy above some angle may be set to a small constant by s e t t i n g the log magnitude to a l a r g e v a l u e . This can be done by multiplying the l o g magnitude by a boxcar f u n c t i o n , then adding a bandpass f u n c t i o n . becomes l a r g e . (The z e r o quefrency point has been set to zero.) This is effect and t h e a r b i t r a r i n e s s n o t an a p p r o p r i a t e of A suggest p r o c e d u r e when P.C. analysis that is this to be u s e d on t h e c e p s t r a . 5.3 An A l t e r n a t e Solution We have seen t h a t , use of P.C.s principal The Exponential A sequence by wavelet this from w h i c h principal estimate method. estimation weighting the This of inputs estimate c a n be o b t a i n e d yields t h e most common w a v e l e t a (Buhl, from suite each of c a n be e x t r a c t e d components. i s illustrated sequence d a t a f o r wavelet o f o u t p u t s may improve 1974). This windowing. However, t h e be u t i l i z e d . quefrency de-weighting by inappropriate. allows and estimates is homomorphic t r a n s f o r m low input cepstra component method may s t i l l by et a l , on i n the presence of a d d i t i v e noise, the generated in Figure previously. 28 u s i n g the noisy An e x p o n e n t i a l multiple weighting 88 ^ • i i i -128 1 1 > 0 128 ms Figure 27 - Out of Band Energy C e p s t r a due t o t h e w a v e l e t of F i g u r e 9 but s e t t i n g t h e log magnitude c o n s t a n t above the a n g l e u i - 0 . 6 2 - i r . T h i s c o n s t a n t has b e e n s e t t o t h e o r i g i n a l v a l u e m i n u s 1 3 , 10, 7 , 3 , a n d 0 ( t o p t o bottom). factor o f 0.98 and a 5 ms c e p s t r a l wavelet estimates. These v a r i a n c e and s h i f t e d common P.C. estimate Their misfits be were u s e d t o p r o d u c e e s t i m a t e s were n o r m a l i z e d t o u n i t to align may window their defined envelope peaks. The most as the average or the f i r s t a r e 163 a n d 20 r e s p e c t i v e l y which shows a 89 Figure 28 - Time Domain P r i n c i p a l Components Wavelets estimated from the data o.f Figure 14 using c e p s t r a l windowing ( a ) . The f i r s t p r i n c i p a l component (PC) and a v e r a g e (AVE) o f t h e s e e s t i m a t e s ( b ) . I n p u t w a v e l e t Is shown a t the bottom. significant improvement T h e s e m i s f i t s compare w i t h o f t h e P.C. 20, 280 a n d 16 for average and f i r s t two c e p s t r a l yield which i s comparable t o c e p s t r a l an e s t i m a t e P.C.s. method o v e r averaging. the cepstral Thus t h i s t e c h n i q u e c a n averaging. 90 5.4 Practical The several previously window shape and l e n g t h and these exponential must be information and The length specified. length, weighting carefully, probably adjusted These cepstral factor. on the scheme has include the window shape Specification basis of a window s h o u l d be slowly of priori varying. distortion and s e q u e n c e s have g r e a t e r d i s t o r t i o n more r e l i a b l e be estimation experience. Shorter Exponential the wavelet i s a t r a d e - o f f between w a v e l e t should can done segmentation length. yield discussed p a r a m e t e r s w h i c h must be segmentation and Considerations phase not by unwrapping. have s h a r p monitoring weighting w a v e l e t ' s c e p s t r u m and cepstral from to window estimates. a t r a d e - o f f between it tend Its length corresponding separating sequence but discontinuities. the i s again The Its distorting the impulse cepstra. will and Due to tend to This technique power judicious 5.5 statistical improve can the of the scheme, more data estimate. i s not be nature user fully independent. realized only Its with flexibility careful and method for use. Summary The wavelet used the a p p l i c a t i o n of estimation i n the cepstral quefrency windowing. has the principal component been d e m o n s t r a t e d . domain or in Additive noise the may The time induce method may be domain after a common signal 91 in the will cepstra. S e t t i n g o u t o f band e n e r g y add a common t e r m these applied common signals i n the time t o the cepstra. t h e P.C. domain. method In to a small the i s more value presence of appropriately 92 VI. A I N V E R T I B I L I T Y OF THE HOMOMORPHIC TRANSFORM recently unwrapping, (Jin and published t h e homomorphic Rogers, 1983). illustrating their This should failure same examples successful within failure This i s supported to obtain with can transform However, will from t h e f a c t that an yield incorrect i f this Kuc (1983) w h i c h a d d s m u l t i p l e s to the whole inverse Fourier transforms, The first 2000 and 1999, sequence and phase operation system (except If occurs, input. which a r e separated by are this 20 points follows of transparent (Note t h a t t h e and with shown. zeros. Note in program.) aliasing. of magnitude (Figure arithmetic and e x h i b i t s some inverse exponentiation ( C a l c u l a t i o n s were p e r f o r m e d w i t h FORTRAN be This These a r e inversion. a will invertible.) t o 1024 p o i n t s inverse is integer multiples of e x p o n e n t i a t i o n . by phase the example c o n s i s t s o f two i m p u l s e s , t h e cepst.rum's delay the t h e method o f McGowan and o f it) . i s realized i s extended their i f the cepstrum the o r i g i n a l to inverse on t h e u n i t c i r c l e , p h a s e u n w r a p p i n g adds principal these transformation. 2ir the inversion. demonstrating computation. even still by e x a m p l e s be c a l c u l a t e d , s o c a n i t s i n v e r s e , t o incorrectly calculated. invertible a successful here zeros the p r e c i s i o n of the unwrapped be To s e t o u r minds a t e a s e presented f o r inputs transform note not due t o p h a s e may not occur. are suggests that, transform forward and i n v e r s e Except forward note 29). The I t s cepstrum the successful single precision The c e p s t r u m i s minimum 93 (a) 1024 (b) 512 -512 (C) 1024 (d) 1024 (e) r -512 512 (f) 1024 ms Figure 29 - I n v e r s i o n of Cepstra A n I n p u t ( a ) c o n s i s t s o f two i m p u l s e s o f a m p l i t u d e 2 0 0 0 a n d 1999. Its c e s t r u m ( b ) a n d t h e c e p s t r u m ' s I n v e r s e ( c ) show t h e r e t u r n of the o r i g i n a l Input. (d) Is t h e I n p u t s e q u e n c e o f (a) with additive noise. (e) i s the cepstum of (d) and ( f ) Is t h e i n v e r s e of (e). 94 The 5 is effect also perceptible relative nearly of This in plots. the noise delay. Again t h e s e c o n d example, The shown. with the second input, It When n o i s e having is has, the cepstrum the f i r s t and too however, such that i s reduced cepstrum a standard level magnitude of t h e impulses 2000 w h i l e 30). the shown. minimum In of a d d i t i v e n o i s e the a standard i n v e r t s t o the o r i g i n a l The Figures All thesis results of to be changed the sequence is i m p u l s e has a m a g n i t u d e by cepstrum's deviation 33% (Figure inverse are o f 20 i s a d d e d , cepstrum's appearance changes d r a m a t i c a l l y . still low inverted c o r r e c t l y . from t h i s the d e v i a t i o n of However, it input. J i n and R o g e r s (1983) a r e r e p r o d u c e d i n 31 and 32 f o r c o m p a r i s o n . of the inputs dealt with were i n v e r t i b l e . We i n the production c a n be s a t i s f i e d of J i n and R o g e r s a r e a n o m a l o u s . that of this the r e s u l t s 95 (a) 1024 (b) -512 512 (c) 1024 (d) 1024 (e) -512 512 (f) 1024 ms Figure 30 - I n v e r s i o n of Cepstra An i n p u t (a) consists of two impulses with different amplitude. I t s c e p s t r u m ( b ) a n d t h e c e p s t r u m ' s i n v e r s e ( c ) show the return of the input. (d) i s the input (a) w i t h a d d i t i v e noise. (e) Is t h e c e p s t r u m o f ( d ) a n d ( f ) Is the inverse of (e). 96 (a) 1.024 (b) -0.512 0.512 (O ' I. '1 1.024 (d) 1.024 (e) -0.512 0.512 1.024 (a) Signal with two spikes of nearly identical amplitudes, free of noise, (b) Complex cepstrum of the signal in (a), (c) Returned signal of (a), (d) Signal of (a) with noise added, (e) Complex cepstrum of the signal in (d). (f) Returned signal of (d). Figure The results Figure 29. of 31 Jin - I n v e r s i o n of and Rogers Cepstra (1983). This compares with 97 (a) 1.024 (b) -0.512 0.512 (c) 1.024 (d) 1.024 (a) Signal with two spikes whose amplitudes are quite different, free of noise, (b) Complex cepstrum of the signal in (a), (c) Returned signal of (a), (d) Signal of (a) with noise added, (e) Complex cepstrum of the signal in (d). (f) Returned signal of (d). Figure The r e s u l t s Figure 30. of 32 Jin and I n v e r s i o n of Rogers Cepstra (1983). This compares with 98 VII. The terms wavelet of DISCUSSION AND CONCLUSIONS e s t i m a t i o n problem multi-channel common has signal not require r e s t r i c t i v e The f o r m u l a t i o n i n c l u d e s a homomorphic convolutional the transform's continuous In t h e a d d i t i v e estimation. the This the wavelet's transform (cepstral) requires formulated estimation. assumptions about t o an a d d i t i v e transform been space. in does phase. t o map from a R e a l i z a t i o n of computation of a Fourier phase. space Properly the problem formulated, components y i e l d s an o p t i m a l i s one o f common the estimate method of signal principal w h i c h c a n be compared t o averaging. The successful demonstrated application f o r noise principal components averaging. This occurs signal the in estimate this cepstrum. setting signal may cepstral an i s added for solution. noisy suite the noisy poorer may instead to a small case estimate than value, a common components of i s made t o r e d u c e h a s been data induce Principal signal data, we i s estimated weighting components a r e t h e n this technique the then wavelet n o i s e e f f e c t s by another common t o the c e p s t r a . A wavelet exponential a noise space. attempt this In yield because o u t o f band e n e r g y Thus, from free data. noise-induced If of used and are from cepstral l e d to each input alternate channel windowing. t o e x t r a c t t h e most of e s t i m a t e s . an common by Principal estimate 99 APPENDIX A - THE The z-transform of a sequence x(n) i s d e f i n e d as X(z) where z i s a complex The Z-TRANSFORM ^x(n).z = (1) - h number. inverse z-transform x(n) = i s d e f i n e d as 1/2TTJJ^"X(Z)-z"" dz 2 (2) C where C i s a c o u n t e r c l o c k w i s e convergence Together, If, origin, of c l o s e d contour X ( z ) , encircling (1) and (2) f o r m the o r i g i n a z-transform i n (1) and ( 2 ) , C i s t a k e n 2[ ( >* x n r of i n the z-plane. pair. as a c i r c l e z = r • exp[ ja>], on -7r<co<jr, t h e n X ( r » e x p [ jo>]) = i n the region c e n t r e d on t h e (1) and (2) become 3*exp[-jco] (3) -tr x{n)'r~ h = ]/2nJ~X{r-exp[ jco]) «exp[ jam]da> (4) -TT These imply transform t h a t X ( r «exp[ jco]) of x(n)«r" . h to evaluate If r=1, e q u a t i o n s Setting and the discrete time Fourier Thus t h e F o u r i e r t r a n s f o r m c a n be used the z-transform on a c i r c l e c e n t r e d on t h e o r i g i n . (3) and (4) become a F o u r i e r T r a n s f o r m rrM i s r e f e r r e d Schafer, is 1975). t o as e x p o n e n t i a l weighting pair. (Oppenheim 100 APPENDIX B ~ FOURIER TRANSFORMS AS CHEBYSHEV The D i s c r e t e Time F o u r i e r Transform in terms o f Chebyshev p o l y n o m i a l s and Kuc, 1982). Consider The (DTFT) c a n be w r i t t e n of the second (McGowan and n>N. DTFT o f x ( n ) i s i n terms of r e a l F(CJ) Chebyshev cos(w), = 2 = 2 x < ) *exp[- jnu] and i m a g i n a r y x ( n ) «cos(n6j) polynomials are defined of (1) n - the as (Snyder, parts j- x(n)'Sinfnco) first kind, those of the second 1966) (3) kind as U(n,cj) = s i n [ (n+1 ) c j ] / s i n ( C J ) There i s a recursion relation betwen (2) and (3) we write (4) (3) and ( 4 ) : T(n,w) = [U(n,oj) - U(n-2,w)]/2 From e q u a t i o n s (2) as a f u n c t i o n of T(n,cj) = c o s ( n u ) and kind a s e q u e n c e x ( n ) where x(n)=0 f o r n<0 F(GJ) or, POLYNOMIALS (5) 101 A/ N F(GJ) Defining 2 < ) •T(n,u>) x = n R = R e a l { F } amd I = I m a g i n a r y { F } , we w r i t e F Consider the r e a l part of R = Dropping t h e argument R= or, using 2 CJ and u s i n g (^x(n).U(n) R = = 0, (6b) (7) x(n)-T(n.w) t e r m s and n o t i n g x(1)«sin(0) as (6): R = { ^ x ( n ) - s i n [ (n+1 ) C J ] - and (6a) = R + j•I the d e f i n i t i t i o n i n Rearranging (6a) - j« ] ^ x ( n ) *sin(na>) (5), equation (7) becomes - ^fx(n)-U(n-2)} (4), ^? x ( n ) - s i n [ (n-1 ) O J ] } / 2 s i n ( o ) that -x(0) • sin(-CJ) = x(0)«sin(cj) yields {[2.x(0)-x(2)]-sin(cj) + ^ [ x U - l ) - x ( n + l ) ] | sin(ncj) + [x(N-l).] • s i n ( NCJ ) + [x(N) sin[(N+l)cj]}/2sin(cj) , 102 or, bringing the denominator inside t h e sum and using (4) yields R = ^ a ( n ) -U(n,co) (8) «=<» where a(0) = x(0) a(n) = a(N-D a(N) Now [x(n>-x(n+2)]/2 n=1,...N-2 [X(N-1)]/2 = = x(2)/2 - x(N)/2 consider the imaginary part of (6): I = - ^ x ( n ) »sin (nco) or, as the f i r s t (9a) term i s z e r o , A/-J I = - 2 x ( n + l ) s i n [ (n+1 )o>] (9b) y> zo Setting b(n) in = -x(n+1) (8). y i e l d s A/-1 1 = 2 b < > * s i n [ (n+1 )CJ] n (10) 1 03 where we in 2 { b ( ) ' S i n [ (n + 1 ) 6 j ] / s i n ( c j ) } I = sin(oj)* I = sin(oj) • S have b ( n ) .U(n,a;) used the d e f i n i t i o n terms of Chebyshev F(CJ) = n polynomials A/ 2 a ( n ) -U(n,cj) Therefore we of t h e s e c o n d k i n d write (1) as M-7 + j-sin(cj) • 5Jb(n)-U(n,6j) where t h e a ( n ) and b ( n ) a r e g i v e n respect i v e l y . (4). by e q u a t i o n s (8) and (10) 1 04 APPENDIX C - GENERATION This sequence appendix from two C h e b y s h e v (McGowan easy and Kuc, transcription Chebyshev defined illustrates OF A STURM a method o f g e n e r a t i n g polynomials 1982). The n o t a t i o n t o a computer polynomials of i n a p p e n d i x B a n d have SEQUENCE of the a Sturm second kind i s designed to allow program. the second the r e c u r s i o n kind, U(n), are relation (Snyder, 1966) U ( n ) - U d ) = U(n+1) + U(n-1) We start where with t h e two p o l y n o m i a l s a(k,n) polynomial polynomials P(0) = ^a(0,n)-U(n) (2) P(1) = 2 (3) refers P(k). of to The a { 1 the Sturm decreasing P ( 0 ) and P ( 1 ) by E u c l i d ' s .n).U(n) n-th coefficient sequence of degree, generated r e c u r s i v e l y from (Marden, general a of the k-th sequence algorithm is P(k+1) = Q ( k ) - P ( k ) - P ( k - 1 ) The (1) 1966) (4) p o l y n o m i a l c a n be w r i t t e n a s P(k) = ^Ta(k,n)-U(n) (5) 1 05 and the final polynomial P(N) In the degree = a(N,0)-U(0) r e c u r s i o n ( 4 ) , l e t us define where t h e For d are Q(k) as the first = { let us begin ( 2 ) , (3) and (7) with into an example by (4) generating yields 2a(l,n).U(n)}.{d(l,0).U(0)+d(l,l).U(l)]} a(0,n).U(n)} expanding, P(2) relation (7) constants. -{ S The the = d(k,0)-U(0) + d(k,l).U(1) clarity Inserting P(2) or, (6) polynomial Q(k) P(2). as (1) = { f a ( l , n ) - d ( l ,0)-U(n)}.U(0) + { S - { ^ a ( O . n ) .U(n)} and the a( 1 , n ) - d d , 1 ) - U ( n ) } - U ( 1) fact t h a t U(0)=1 yield 106 P(2) = ^>(1 ,n)-d( 1 ,0)-U(n) + {a(1,1)-d(1,1)}-U(0) + 2 {a(1,n-1)-d(1,1)+a(1,n+l)-d(1,1)}-U(n) + {a(1,N-2)-d(1,1)}-U(N-1) + {a(1 N-1)-d(1,1)}-U(N) - f S a(0,n)-U(n) h-o Collecting like terms P(2) = { a ( 1 , 0 ) - d ( 1 , 0 ) yields + a(1,l)-d(1,1) - a(0,0)} A/-7 ^{a(1,n)-d(1,0) * + [ a ( 1 , n - 1 ) + a ( 1 , n + 1 ) ]-d(1,1) - { a ( 1 , N - 1 ) - d ( 1 , 0 ) + a ( 1 , N - 2 ) - d ( 1 , 1 ) - a ( 0 , N - 1 )} -U(N-1 ) + + {a(,1 ,N-1)-d( 1 ,1 )-a(0,N)}-U(N) We c h o o s e Q by in U(N) By a(0,n)}-U(n) and finding d ( 1 , 0 ) and d ( l , l ) (N-1) c a n c e l . a(0,N)/a(1,N-1) d(1,0) = {a(0,N-O polynomial as we can the terms That i s , d( 1 , 1) = induction such t h a t - write a(1,N-2)-d(1,1)}/a(1,N-1) the general term f o r the k-th 1 07 P(k) {a(k-1,0)-d(k-1,0)+a(k-1,1)-d(k-1,1)-a(k-2,0)} = + U(k-1,n).d(k-1,0)+[a(k-1,n-1)+a(k-1,n+1)] S - and the f i n a l P(N) p o l y n o m i a l as {a(N-1,0)-d(N-1,0) + = - The a(k-2,n)}-U(n) a(N-1,1)-d(N-1,1) a(N-2,0)}.U(0) k - t h p o l y n o m i a l c a n be w r i t t e n a s P(k) = £ a(k,n)-U(n) k=2,...N where a(k,0) and, = a ( k - 1 , 0 ) - d ( k - 1 , 0 ) + a ( k - 1 , 1 ) - d ( k - 1 , 1 ) - a ( k - 2 , 0 ) (8a) f o r n>0 a(k,n) = a(k-1,n)-d(k-1,0)+[a(k-1,n-l)+a(k-1,n+1)]-d(k-1,1) - a(k-2,n) (8b) where, f o r k=0,...N-2, d(k+1,l) = a(k,N-k)/a(k+1,N-k-1) (9a) 108 and d(k+1,0) = [a(k,N-k-1)-a(k+1,N-k-2)>d(k+1,1) ] /a(k+1,N-k-1) (9b) The all recursion the prime. relations (8) a n d (9) a l l o w t h e If this case i s the case, the h i g h e s t degree t h e next by more t h a n continues of P ( k ) , k>0, p r o v i d e d P ( 0 ) and P ( 1 ) a r e n o t r e l a t i v e l y a P ( k ) i s z e r o and a d i v i s i o n this computation one. by z e r o i n (9) w i l l polynomial generated The r e c u r s i o n f o r subsequent coefficient then calculations. occur. has a d e g r e e of For reduced s k i p s one i t e r a t i o n and 109 BIBLIOGRAPHY 1. Ahmed, N. and Rao, K.R., 1975, O r t h o g o n a l T r a n s f o r m s f o r D i g i t a l S i g n a l P r o c e s s i n g : S p r i n g e r - V e r l a g , New Y o r k . 2. Beaumont, R.A. and P i e r c e , R.S., 1963, The A l g e b r a i c F o u n d a t i o n s o f M a t h e m a t i c s : A d d i s o n - W e s l e y , U.S.A. 3. Bhanu, B., 1977, C o m p u t a t i o n o f Complex C e p s t r u m : and E . E . t h e s i s , M.I.T., M a s s a c h u s e t t s . 4. B o g e r t , B.P., H e a l y , M.J.R. and T u k e y , J.W., 1962, The Q u e f r e n c y A l a n y s i s o f Time S e r i e s f o r E c h o e s : C e p s t r u m , P s e u d o - A u t o c o v a r i a n c e , C r o s s - C e p s t r u m and Saphe C r a c k i n g : P r o c . Symp. on Time S e r i e s A n a l y s i s , E d . M. R o s e n b l a t t , J o h n W i l e y and s o n s , I n c . , N.Y. 5. B r a c e w e l l , R., 1965, The F o u r i e r T r a n s f o r m and I t s A p p l i c a t i o n s : M c G r a w - H i l l , U.S.A. 6. B r i g h a m , E.O., Prentice-Hall, 7. B u h l , P., S t o f f a , P.L. and B r y a n , G.M., 1974, The A p p l i c a t i o n o f Homomorphic D e c o n v o l u t i o n t o S h a l l o w Water M a r i n e S e i s m o l o g y - P a r t I I : R e a l D a t a : G e o p h y s i c s , v . 39, p. 417-426. 8. B u t t k u s , B., 1975, Homomorphic F i l t e r i n g - T h e o r y and P r a c t i c e : G e o p h y s i c a l P r o s p e c t i n g , v . 23, p . 712-748. 9. C l a y t o n , R.W. and W i g g i n s , R.A., 1976, S o u r c e Shape E s t i m a t i o n and D e c o n v o l u t i o n o f T e l e s e i s m i c Bodywaves: G e o p h y s . J.R. a s t r . S o c , v . 47, p. 151-177. 10. 1974, The F a s t F o u r i e r I n c . , New J e r s e y . M.Sc. Transform: Dhrymes, P . J . , 1970, E c o n o m e t r i c s , S t a t i s t i c a l F o u n d a t i o n s and A p p l i c a t i o n s : S p r i n g e r - V e r l a g , New York. 11. G e l b , A. e t a l . , 1974, A p p l i e d O p t i m a l A n a l y t i c S c i e n c e s C o r p o r a t i o n , M.I.T., E s t i m a t i o n : The Massachusetts. 12. H a l l , E . L . , 1979, Computer Image P r o c e s s i n g and R e c o g n i t i o n : Academic P r e s s , I n c . , New Y o r k . 13. J i n , D . J . and R o g e r s , J.R., 1983, Homomorphic D e c o n v o l u t i o n : G e o p h y s i c s , v . 48, p. 1014-1016. 14. Kramer, H.P. and Mathews, M.V., 1956, L i n e a r ' C o d i n g F o r T r a n s m i t t i n g a S e t o f C o r r e l a t e d S i g n a l s : I.R.E. T r a n s , on I n f o r m a t i o n T h e o r y , v . I T - 2 , p. 41-46. 15. L i n e s , L.R. and U l r y c h , T . J . , 1977, The O l d and t h e New i n S e i s m i c D e c o n v o l u t i o n and W a v e l e t E s t i m a t i o n : 110 Geophysical Prospecting, v . 25, p . 512-540. 16. L i p s c h u t z , S., 1974, L i n e a r A l g e b r a : S e r i e s , M c G r a w - H i l l , London. Schaum's O u t l i n e 17. Marden, M., 1966, Geometry o f P o l y n o m i a l s : M a t h e m a t i c a l S o c i e t y , Rhode I s l a n d . 18. McGowan, R. a n d Kuc, R., 1982, A D i r e c t R e l a t i o n Between a S i g n a l Time S e r i e s a n d I t s Unwrapped P h a s e : I . E . E . E . T r a n s , on A.S.S.P., v . ASSP-30, p . 719-726. 19. O l d e n b u r g , D.W. a n d Samson, J . C , 1979, I n v e r s i o n o f I n t e r f e r o m e t r i c D a t a From C y l i n d r i c a l l y Symmetric, R e f r a c t i o n l e s s P l a s m a s : J . O p t . S o c . Am., v . 69, p . 927942. 20. Ooe, M. a n d U l r y c h , T . J . , 1979, Minimum E n t r o p y D e c o n v o l u t i o n w i t h an E x p o n e n t i a l Transformation: G e o p h y s i c a l P r o s p e c t i n g , v . 27, p . 458-473. 21. Oppenheim, A.V., 1965, S u p e r p o s i t i o n i n a C l a s s o f N o n l i n e a r S y s t e m s : T e c h . R e p t . 432, R e s e a r c h L a b o r a t o r y o f E l e c t r o n i c s , M.I.T., M a s s a c h u s e t t s . 22. Oppenheim, A.V. a n d S c h a f e r , R.W., 1975, D i g i t a l P r o c e s s i n g : P r e n t i c e - H a l l , I n c . New J e r s e y . 23. Oppenheim, A.V., S c h a f e r , R.W. a n d Stockham, T.G., 1968, N o n l i n e a r F i l t e r i n g o f M u l t i p l i e d and C o n v o l v e d S i g n a l s : P r o c I . E . E . E . , v . 56, p . 1264-1291. 24. O t i s , R.M. a n d S m i t h , R.B., 1977, Homomorphic D e c o n v o l u t i o n by L o g S p e c t r a l A v e r a g i n g : G e o p h y s i c s , v . 42, p . 1146-1157. 25. P o g g i a g l i o l m i , E . , B e r k h o u t , A . J . a n d Boone, M.M., Phase U n w r a p p i n g , P o s s i b i l i t i e s a n d L i m i t a t i o n s : G e o p h y s i c a l P r o s p e c t i n g , v . 30, p . 281-291. 26. Ready, P . J . a n d W i n t z , P.A., 1973, I n f o r m a t i o n E x t r a c t i o n , S.N.R. Improvement, and D a t a C o m p r e s s i o n i n M u l t i s p e c t r a l Imagery: I . E . E . E . T r a n s , on C o m m u n i c a t i o n s , v . COM-21, p 1123-1130. 27. R o b i n s o n , E.A. and T r e i t e l , S., 1980, G e o p h y s i c a l A n a l y s i s : P r e n t i c e - H a l l , I n c . , New J e r s e y . 28. S c h a f e r , R.W., 1969, Echo Removal by D i s c r e t e G e n e r a l i z e d L i n e a r F i l t e r i n g : T e c h . R e p t . 466, MIT R e s e a r c h L a b o r a t o r y o f E l e c t r o n i c s , MIT, M a s s a c h u s e t t s . 29. S n y d e r , M.A., 1966, C h e b y s h e v Methods i n N u m e r i c a l A p p r o x i m a t i o n : P r e n t i c e - H a l l , I n c . , New J e r s e y . American Signal 1982, Signal 111 30. S t e i g l i t z , K. and by F a c t o r i z a t i o n : 30, p. 984-991. 31. S t o f f a , P.L., B u h l . , P . and B r y a n , G.M., 1974, The A p p i c a t i o n s of Homomrphic D e c o n v o l u t i o n t o S h a l l o w - W a t e r M a r i n e S e i s m o l o g y - P a r t I : M o d e l s : G e o p h y s i c s , v . 39, p. 401-416. 32. S t r a n g , G., 1980, L i n e a r A l g e b r a and e d . ) : Academic P r e s s . New Y o r k . 33. T a n e r , M.T. and S h e r i f f , R.E., 1977, A p p l i c a t i o n o f A m p l i t u d e , F r e q u e n c y , and O t h e r A t t r i b u t e s t o S t r a t i g r a p h i c and H y d r o c a r b o n D e t e r m i n a t i o n : i n S e i s m i c S t r a t i g r a p h y - A p p l i c a t i o n s to Hydrocarbon E x p l o r a t i o n , E d . C.E. P a y t o n , Memoir 26 o f A.A.P.G, T u l s a , p. 301328. 34. T e l f o r d , W.M., G e l d a r t , L.P., S h e r i f f , R.E. and K e y s , D.A., 1976, A p p l i e d G e o p h y s i c s : Cambridge U n i v e r s i t y P r e s s , New Y o r k . 35. T r i b o l e t , J.M., 1977, A New Phase Unwrapping A l g o r i t h m : I . E . E . E . T r a n s , on A.S.S.P., v . ASSP-25, p. 170-177. 36. T r i b o l e t , J.M., 1979, S e i s m i c A p p l i c a t i o n s o f Homomorphic S i g n a l P r o c e s s i n g : P r e n t i c e - H a l l , I n c , Jersey. New U l r y c h , T . J . , 1971, A p p l i c a t i o n o f Homomorphic D e c o n v o l u t i o n t o S e i s m o l o g y : G e o p h y s i c s , v . 36, p. 660. 650- 37. D i c k i n s o n , B., 1982, Phase Unwrapping I . E . E . E . T r a n s , on A.S.S.P., v ASSP- Its Applications (2 38. U l r y c h , T . J . , L e v y , S., O l d e n b u r g , D.W., and J o n e s , I . F . , 1983, A p p l i c a t i o n s o f t h e K a r h u n e n - L o e v e Transformation i n R e f l e c t i o n Seismology: presented at t h e 5 3 r d A n n u a l M e e t i n g o f t h e S.E.G., L a s V e g a s . 39. W i g g i n s , R.A., 1972, The G e n e r a l L i n e a r I n v e r s e P r o b l e m : I m p l i c a t i o n s o f S u r f a c e Waves and F r e e O s c i l l a t i o n s f o r E a r t h S t r u c t u r e : Reviews o f G e o p h y s i c s and Space P h y s i c s , v . 10, p . 251-285. 40. W i g g i n s , R.A., 1978, Minimum E n t r o p y D e c o n v o l u t i o n : G e o e x p l o r a t i o n , v . 16, p. 21-35.
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Multi-channel homomorphic wavelet estimation
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Multi-channel homomorphic wavelet estimation Lane, Mark Christopher 1983
pdf
Page Metadata
Item Metadata
Title | Multi-channel homomorphic wavelet estimation |
Creator |
Lane, Mark Christopher |
Publisher | University of British Columbia |
Date Issued | 1983 |
Description | Wavelet estimation can be posed as a multi-channel common information problem. Each channel of data is modeled as the convolution of a wavelet with an impulse sequence. A homomorphic transform maps the data from a convolutional to an additive space. The mapping may also effect partial separation of wavelet and impulses. In the additive space the wavelet can be estimated using averaging. This is termed cepstral averaging. This thesis reviews the homomorphic transform and provides a synthesis and comparison of the techniques available for its realization. The method of principal components for wavelet estimation is proposed as an alternative to cepstral averaging. The effect of noise on this method is investigated. The investigation shows that noise may cause principal components to produce estimates which are inferior to cepstral averaging. For these cases an alternate solution is proposed in which principal components are used in the original convolutional space. A wavelet is estimated by homomorphic separation for each data channel. Principal components may then be used to define a best estimate from this suite of estimates. |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2010-05-15 |
Provider | Vancouver : University of British Columbia Library |
Rights | For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use. |
DOI | 10.14288/1.0052934 |
URI | http://hdl.handle.net/2429/24714 |
Degree |
Master of Science - MSc |
Program |
Geophysics |
Affiliation |
Science, Faculty of Earth, Ocean and Atmospheric Sciences, Department of |
Degree Grantor | University of British Columbia |
Campus |
UBCV |
Scholarly Level | Graduate |
Aggregated Source Repository | DSpace |
Download
- Media
- 831-UBC_1984_A6_7 L35.pdf [ 4MB ]
- Metadata
- JSON: 831-1.0052934.json
- JSON-LD: 831-1.0052934-ld.json
- RDF/XML (Pretty): 831-1.0052934-rdf.xml
- RDF/JSON: 831-1.0052934-rdf.json
- Turtle: 831-1.0052934-turtle.txt
- N-Triples: 831-1.0052934-rdf-ntriples.txt
- Original Record: 831-1.0052934-source.json
- Full Text
- 831-1.0052934-fulltext.txt
- Citation
- 831-1.0052934.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.831.1-0052934/manifest