UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Multi-channel homomorphic wavelet estimation Lane, Mark Christopher 1983

Your browser doesn't seem to have a PDF viewer, please download the PDF to view this item.

Item Metadata

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

MULTI-CHANNEL HOMOMORPHIC WAVELET ESTIMATION by MARK CHRISTOPHER LANE B . S c , The U n i v e r s i t y Of V i c t o r i a , 1980 A THESIS SUBMITTED IN PARTIAL FULFILMENT OF THE REQUIREMENTS FOR THE DEGREE OF MASTER OF SCIENCE in THE FACULTY OF GRADUATE STUDIES Department Of Geophysics And Astronomy We accept t h i s t h e s i s as conforming to the r e q u i r e d standard THE UNIVERSITY OF BRITISH COLUMBIA November 1983 © Mark Ch r i s t o p h e r Lane, 1983 In p r e s e n t i n g t h i s t h e s i s i n p a r t i a l f u l f i l m e n t of the requirements f o r an advanced degree at the U n i v e r s i t y o f B r i t i s h Columbia, I agree t h a t the L i b r a r y s h a l l make i t f r e e l y a v a i l a b l e f o r r e f e r e n c e and study. I f u r t h e r agree t h a t p e r m i s s i o n f o r e x t e n s i v e copying o f t h i s t h e s i s f o r s c h o l a r l y purposes may be granted by the head of my department or by h i s o r her r e p r e s e n t a t i v e s . I t i s understood t h a t copying or p u b l i c a t i o n o f t h i s t h e s i s f o r f i n a n c i a l g a i n s h a l l not be allowed without my w r i t t e n p e r m i s s i o n . Department O f G e o p h y s i c s and A s t r o n o m y The U n i v e r s i t y of B r i t i s h Columbia 1956 Main M a l l Vancouver, Canada V6T 1Y3 Date J/^. /?, /403 DE-6 (3/81) i i A b s t r a c t Wavelet e s t i m a t i o n can be posed as a mu l t i - c h a n n e l common in f o r m a t i o n problem. Each channel of data i s modeled as the c o n v o l u t i o n of a wavelet with an impulse sequence. A homomorphic transform maps the data from a c o n v o l u t i o n a l to an a d d i t i v e space. The mapping may a l s o e f f e c t p a r t i a l s e p a r a t i o n of wavelet and impulses. In the a d d i t i v e space the wavelet can be estimated using a v e r a g i n g . T h i s i s termed c e p s t r a l a veraging. T h i s t h e s i s reviews the homomorphic transform and pro v i d e s a s y n t h e s i s and comparison of the techniques a v a i l a b l e f o r i t s r e a l i z a t i o n . The method of p r i n c i p a l components f o r wavelet e s t i m a t i o n i s proposed as an a l t e r n a t i v e to c e p s t r a l a v e r a g i n g . The e f f e c t of noise on t h i s 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 shows that noise may cause p r i n c i p a l components to produce estimates which are i n f e r i o r to c e p s t r a l a v e r a g i n g . For these cases an a l t e r n a t e s o l u t i o n i s proposed i n which p r i n c i p a l components are used i n the o r i g i n a l c o n v o l u t i o n a l space. A wavelet i s estimated by homomorphic s e p a r a t i o n f o r each data channel. P r i n c i p a l components may then be used to d e f i n e a best estimate from t h i s s u i t e of estimates. Table of Contents A b s t r a c t i i L i s t of Tables v L i s t of F i g u r e s v i Acknowledgements v i i i I. INTRODUCTION 1 I I . HOMOMORPHIC TRANSFORMS 3 2.1 I n t r o d u c t i o n 3 2.2 G e n e r a l i z e d S u p e r p o s i t i o n 4 2.3 The C h a r a c t e r i s t i c System 7 2.4 The Inverse 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 System 12 2.6 P r o p e r t i e s Of The Complex Cepstrum 13 2.7 Computational 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 18 I I I . PHASE UNWRAPPING 19 3.1 I n t r o d u c t i o n 19 3.2 P r i n c i p a l Value 23 3.3 I n t e g r a t i o n 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 Theory ....28 3.5.1 Computational D e t a i l s 39 3.6 D i f f i c u l t i e s Of Phase Unwrapping 40 3.7 Comparison Of Phase Unwrapping Techniques 42 3.8 Summary 43 i v 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 Components 45 4.3 Some P r o p e r t i e s Of P r i n c i p a l Components 55 4.4 Summary 57 V. WAVELET ESTIMATION 58 5.1 The Problem 58 5.2 A S o l u t i o n 58 5.2.1 The Smoothing Fu n c t i o n 62 5.2.2 The Wavelet 64 5.2.3 P r i n c i p a l Components vs Averaging 66 5.2.4 M u l t i p l e Sequence Data 69 5.2.5 S i n g l e Sequence Data 77 5.2.6 Low Pass Inputs 86 5.3 An A l t e r n a t e S o l u t i o n 87 5.4 P r a c t i c a l C o n s i d e r a t i o n s 90 5.5 Summary 90 VI. INVERTIBILITY OF THE HOMOMORPHIC TRANSFORM 92 V I I . DISCUSSION AND CONCLUSIONS 98 APPENDIX A - THE Z-TRANSFORM 99 APPENDIX B - FOURIER TRANSFORMS AS CHEBYSHEV POLYNOMIALS 100 APPENDIX C - GENERATION OF A STURM SEQUENCE . . 104 BIBLIOGRAPHY 109 V L i s t of T a b l e s 1 . W a v e l e t M i s f i t s 85 v i L i s t of F i g u r e s 1. Canonic Representation of a Homomorphic System 6 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 8 3. Representation of the System D?? 12 4. Complex Fu n c t i o n i n the Complex Plane 20 5. P r i n c i p a l Value and Unwrapped Phase 22 6. Branches of arctangent F u n c t i o n 30 7. Smoothing F u n c t i o n 64 8. Complex Cepstra Due to Impulses 65 9. An Assumed Wavelet 66 10. D i s c r i m i n a t i o n by P r i n c i p a l Components 67 11. Wavelet Estimates 68 12. M u l t i p l e Sequence Data 70 13. Wavelet Estimates ..71 14. Noisy Data 72 15. Noise Inputs 73 16. Cepstra of Noise 74 17. Eigenvalues of Input and Output 75 18. Covariance Matrix of Noise ..75 19. Covariance M a t r i c e s 76 20. S i n g l e Sequence Data 79 21. Segmentation and E s t i m a t i o n 80 22. Cepstra of Segments 81 23. Noisy S i n g l e Sequence Data 82 24. Wavelet Estimates 83 25. Segment's Cepstra 84 v i i 26. Low Pass S i g n a l s 87 27. Out of Band Energy 88 28. Time Domain P r i n c i p a l Components 89 29. I n v e r s i o n of Cepstra 93 30. I n v e r s i o n of Cepstra 95 31. I n v e r s i o n of Cepstra 96 32. I n v e r s i o n of Cepstra 97 v i i i Acknowledgement At t h i s p o i n t i n a t h e s i s , space i s a l l o c a t e d f o r p a n e g y r i z a t i o n . T h e r e f o r e , f o l l o w i n g the course of an u n t o l d m u l t i t u d e , I w i l l do so. My s u p e r v i s o r , P r o f e s s o r Tadeusz U l r y c h , has managed to guide me through the l a b y r i n t h of academic rese a r c h to the completion of t h i s t h e s i s . T h i s he has done by s h a r i n g f r e e l y h i s wisdom, e x p e r t i s e and p h i l o s o p h i e s . I am g r a t e f u l f o r the o p p o r t u n i t i e s given me, by Tad, to meet many people who make resea r c h what i s i s . I thank Brad Prager, p a r t i c u l a r l y f o r s h a r i n g h i s v ast knowledge of computers and h i s programs. Shlomo Levy, I a i n Jones and Kerry S t i n s o n p r o v i d e d me with many s t i m u l a t i n g d i s c u s s i o n s about the t o p i c s d e a l t with i n t h i s t h e s i s , as w e l l as myriad o t h e r s . The comraderie which accompanies students p r o g r e s s i n g through the system together i s important. A c c o r d i n g l y , and g l a d l y , I thank J u l i a n (Bandito) Cabrera, Don White, Lynda F i s k and Bob J e f f r e y . Without them l i f e would have been much d u l l e r . C h r i s t i a n Sutherland has s u r v i v e d my many l u c u b r a t i o n s and pro v i d e d much needed balance and ix frame. I am deeply g r a t e f u l . I thank the U n i v e r s i t y of B r i t i s h Columbia and the N a t u r a l E n g i n e e r i n g and Research C o u n c i l f o r t h e i r f i n a n c i a l support i n the form of f e l l o w s h i p s . 1 I. INTRODUCTION In any i n v e s t i g a t i o n , i t i s the understanding and i n s i g h t gained which i s paramount. "Let a l l t h i n g s be done unto e d i f y i n g . " 1 The o r i g i n a l purpose of t h i s t h e s i s was to i n v e s t i g a t e wavelet e s t i m a t i o n by combining a n o n - l i n e a r mapping and an optimal l i n e a r system. However, i t has become p r i m a r i l y an i n v e s t i g a t i o n of the mapping, the system and t h e i r i n t e r a c t i o n . Wavelet e s t i m a t i o n can be posed as a m u l t i - c h a n n e l common in f o r m a t i o n problem. Each channel of data i s modeled as the c o n v o l u t i o n of a wavelet with an impulse sequence. The wavelet remains constant while the impulse sequence v a r i e s from channel to channel. A n o n - l i n e a r homomorphic transform can be used to map the data from a c o n v o l u t i o n a l to an a d d i t i v e space. T h i s mapping may a l s o e f f e c t p a r t i a l or t o t a l s e p a r a t i o n of the wavelet and impulses. In the a d d i t i v e space the problem can be analyzed with l i n e a r techniques. We propose the use of p r i n c i p a l components to f u r t h e r separate the wavelet and impulses. In p r a c t i c e t h i s approach y i e l d e d some unexpected r e s u l t s . R e s o l u t i o n of these n e c e s s i t a t e d re-examination of the homomorphic transform and p r i n c i p a l components. T h i s l e d to a c l e a r e r understanding of t h e i r i n t e r a c t i o n . 1 St. Paul, 1 C o r i n t h i a n s 14,26, King James V e r s i o n of the Holy B i b l e . 2 C h a p t e r I I o f t h i s t h e s i s e x a m i n e s t h e h o m o m o r p h i c t r a n s f o r m : i t s u n d e r l y i n g t h e o r y , d e v e l o p m e n t , p r o p e r t i e s a n d p r a c t i c a l r e a l i z a t i o n . R e a l i z a t i o n i n c l u d e s t h e c a l c u l a t i o n o f t h e c o n t i n o u s p h a s e o f a F o u r i e r t r a n s f o r m . T h i s i s t e r m e d p h a s e u n w r a p p i n g a n d i s d i s c u s s e d i n C h a p t e r I I I . P h a s e u n w r a p p i n g i s e s s e n t i a l l y a n u m e r i c a l a n a l y s i s p r o b l e m a l t h o u g h i t m a y i n c o r p o r a t e d i v e r s e m a t h e m a t i c a l r e l a t i o n s . A n e w a p p r o a c h , o f t h e o r e t i c a l i n t e r e s t , i s e x a m i n e d i n d e t a i l . H a v i n g m a p p e d f r o m a c o n v o l u t i o n a l t o a n a d d i t i v e s p a c e , a v e r a g i n g c a n b e u s e d f o r w a v e l e t e s t i m a t i o n . A n a l t e r n a t e m e t h o d i s p r o v i d e d b y p r i n c i p a l c o m p o n e n t a n a l y s i s . T h i s i s d e v e l o p e d i n C h a p t e r I V i n t e r m s o f o p t i m a l i n f o r m a t i o n e x t r a c t i o n . T h e a p p l i c a t i o n o f p r i n c i p a l c o m p o n e n t s t o w a v e l e t e s t i m a t i o n i s d e a l t w i t h i n C h a p t e r V . F o r t h e c a s e o f d a t a c o r r u p t e d b y a d d i t i v e n o i s e , p r i n c i p a l c o m p o n e n t s m a y y i e l d p o o r e r e s t i m a t e s t h a n a v e r a g i n g . T o a v o i d t h i s p r o b l e m , p r i n c i p a l c o m p o n e n t s m a y b e u s e d i n t h e o r i g i n a l c o n v o l u t i o n a l s p a c e . F o r e a c h c h a n n e l o f d a t a a w a v e l e t e s t i m a t e m a y b e f o u n d b y h o m o m o r p h i c s e p a r a t i o n . P r i n c i p a l c o m p o n e n t s m a y t h e n b e u s e d t o d e f i n e a b e s t e s t i m a t e f r o m t h e s e i n d i v i d u a l e s t i m a t e s . T h i s p r o c e d u r e i s a l s o d i s c u s s e d i n C h a p t e r V . R e c e n t l i t e r a t u r e h a s c a s t d o u b t s o n t h e h o m o m o r p h i c t r a n s f o r m ' s i n v e r t i b i l i t y . C h a p t e r V I b r i e f l y a d d r e s s e s t h i s c o n c e r n . 3 I I . HOMOMORPHIC TRANSFORMS 2.1 I n t r o d u c t i o n 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 i s that of se p a r a t i n g s i g n a l s which have been combined i n a known way. The theory of mappings has proven to be u s e f u l i n i t s s o l u t i o n . Consider a t r a n s f o r m a t i o n which maps combined s i g n a l s i n t o d i f f e r e n t r e g i o ns of a space. These regions may be separated and i n d i v i d u a l l y mapped back to the o r i g i n a l space, y i e l d i n g the separate input s i g n a l s . Even i f the s i g n a l s are not separated, v a r i o u s techniques may be more amenable to t h i s s e p a r a t i o n i n the transform space than i n the o r i g i n a l space. The problem then becomes t h a t of f i n d i n g an ap p r o p r i a t e t r a n s f o r m a t i o n . We do not know, a p r i o r i , t h a t such a transform e x i s t s and must c o n s i d e r each case i n d i v i d u a l l y . One approach to f i n d i n g an a p p r o p r i a t e transform i s provided by the theory of homomorphic systems. The theory of homomorphic systems was f o r m a l i z e d by Oppenheim (1965). He c o n s i d e r e d c l a s s e s of n o n l i n e a r systems which obey a p r i n c i p l e of g e n e r a l i z e d s u p e r p o s i t i o n . These systems may be represented by a l g e b r a i c a l l y l i n e a r t r a n s f o r m a t i o n s between input and output v e c t o r spaces and are thus c a l l e d homomorphic systems ( L i p s c h u t z , 1974, p. 123). The g e n e r a l i z e d s u p e r p o s i t i o n d e f i n e s the r e l a t i o n s h i p between a r u l e f o r combining s i g n a l s i n the input space and that f o r combining them i n the output space. L i n e a r systems are a 4 s p e c i a l case of homomorphic systems. 2.2 G e n e r a l i z e d S u p e r p o s i t i o n For the purposes of t h i s d i s c u s s i o n , a l l input s i g n a l s are c o n s i d e r e d to be d i s c r e t e and thus may be c o n s i d e r e d as sequences or as m u l t i d i m e n s i o n a l v e c t o r s . The n o t a t i o n x(n) i n d i c a t e s that x i s a f u n c t i o n of the d i s c r e t e v a r i a b l e n. Consider a system d e f i n e d by the l i n e a r t r a n s f o r m T. Let x,(n) and x 2 ( n ) be input s i g n a l s and a and b be constant s c a l a r s . Then the system T, by d e f i n i t i o n , s a t i s f i e s the s u p e r p o s i t i o n r e l a t i o n T[a-x,<n) + b - x 2 ( n ) ] = a'T[x,(n)] + b - T [ x 2 ( n ) ] (1) Note that the order of s c a l a r m u l t i p l i c a t i o n , s i g n a l a d d i t i o n and system t r a n s f o r m a t i o n of s i g n a l s i s i r r e l e v a n t . T h i s s u p e r p o s i t i o n r e l a t i o n shows the s u i t a b i l i t y of a l i n e a r t r ansform to the s e p a r a t i o n of s i g n a l s combined by a d d i t i o n . I t i s p o s s i b l e to d e f i n e more general r u l e s f o r combining s i g n a l s and f o r combining s i g n a l s with s c a l a r s . In t h i s case the c orresponding r u l e s f o r combining the s i g n a l transforms w i l l be, i n g e n e r a l , d i f f e r e n t from those used to combine the input s i g n a l s . Let * represent a r u l e f o r combining input s i g n a l s and : represent a r u l e f o r combining a s c a l a r with an input s i g n a l . S i m i l a r l y , l e t + represent a r u l e f o r combining output s i g n a l s and j_ represent a r u l e f o r combining output s i g n a l s with s c a l a r s . Then a system H s a t i s f i e s a g e n e r a l i z e d 5 p r i n c i p l e of superposition i f H[a:x, (n)*b:x 2 (n) ] = a^Htx^n)] + bj_H[x 2(n)] (2) Comparison of equations (1) and (2) shows that linear systems s a t i s f y a generalized superposition with the operations : and as m u l t i p l i c a t i o n and the operations * and + as addition. There are a variety of systems which obey a generalized p r i n c i p l e of superposition. The mathematical r e s t r i c t i o n s and the formalism of such systems were developed by Oppenheim (1975) and the d e t a i l s w i l l not be given here. The class of systems sp e c i f i e d by equation (2) can be interpreted as a l g e b r a i c a l l y l i n e a r transformations from an input vector space to an output vector space. The rules for combining signals correspond to vector addition and those for combining scalars with signals correspond to scalar m u l t i p l i c a t i o n . A l l systems of t h i s class can be represented as a cascade of three systems known as the canonic representation. In t h i s thesis we w i l l consider that homomorphic system in which * represents discrete convolution and + represents addition. The operation : represents a rule for combining scalars with input signals and i s best defined through the output operation _^  which represents scalar m u l t i p l i c a t i o n . In the special case that the scalar i s an integer m, : corresponds to the convolution of the signal with i t s e l f m times. In the interest of r e a d a b i l i t y , l e t us denote the output o p e r a t i o n s + and j_ by + and • r e s p e c t i v e l y . Then equation (2) becomes H [ a : x , ( n ) * b : x 2 ( n ) ] = a«H[x,(n)] + b - H [ x 2 ( n ) ] (3) The canonic r e p r e s e n t a t i o n of t h i s homomorphic system i s shown in F i g u r e 1. x(n) Figure 1 - Canonic Representation of a Homomorphic System 0 maps from a convolutional space to an additive space. P i s a homomorphic system which transforms from a c o n v o l u t i o n a l space to an a d d i t i v e space. I t i s d e f i n e d by the r e l a t i o n s h i p D [ a : x , ( n ) * b : x 2 ( n ) ] = a«D[x,(n)] + b«D[x 2(n)] (4) L i s a l i n e a r system, d e f i n e d by the r e l a t i o n i n equation ( 1 ) . The system D~1 i s the i n v e r s e of D and i s d e f i n e d by the r e l a t i o n s h i p 7 D" 1{a-D[x,(n)] + b«D[x 2(n)]} = a.x, ( n ) * b : x 2 ( n ) (5) In g e n e r a l the system D s p e c i f i e s the canonic r e p r e s e n t a t i o n and i s c a l l e d the c h a r a c t e r i s t i c system. The canonic r e p r e s e n t a t i o n shows t h a t , once D i s f i x e d , the problem of s e p a r a t i n g s i g n a l s which have been combined by c o n v o l u t i o n i s reduced to that of l i n e a r f i l t e r i n g . The s o l u t i o n , t h e r e f o r e , l i e s i n the s p e c i f i c a t i o n of the l i n e a r system L. 2.3 The C h a r a c t e r i s t i c System The c h a r a c t e r i s t i c system D serves to transform from a c o n v o l u t i o n a l space to an a d d i t i v e space. I t can be decomposed i n t o i t s c o n s t i t u e n t transforms. R e c a l l that the z-transform of the c o n v o l u t i o n of two s i g n a l s i s equal to the product of the z-transforms of the i n d i v i d u a l s i g n a l s (Appendix A ) . That i s Z [ x , ( n ) * x 2 ( n ) ] = z[x,(n) ]-Z[x 2(n) ] (6) I t i s a homomorphic system which maps from a c o n v o l u t i o n a l to a m u l t i p l i c a t i v e space. The system D has a canonic r e p r e s e n t a t i o n shown i n F i g u r e 2 and may be r e a l i z e d i n three s t e p s . 8 x(n) Figure 2 - Representation of the C h a r a c t e r i s t i c System D F i r s t , the z-transform maps an input sequence, x ( n ) , from a c o n v o l u t i o n a l space i n t o a continuous f u n c t i o n , X ( z ) , i n a m u l t i p l i c a t i v e space. Next, the complex l o g a r i t h m maps t h i s f u n c t i o n from the m u l t i p l i c a t i v e space i n t o another continuous f u n c t i o n , X ( z ) , i n an a d d i t i v e space. F i n a l l y , the i n v e r s e z-transform maps t h i s continuous f u n c t i o n from the a d d i t i v e space i n t o a sequence, x ( n ) , i n another a d d i t i v e space. T h i s output, x ( n ) , i s c a l l e d the complex cepstrum. The word cepstrum was proposed by Bogert, Healy and Tukey (1963) as a paraphrase of the word spectrum. T h i s , along with fo u r t e e n other words, was proposed f o r use i n t h i s a d d i t i v e space because "although strange at f i r s t s i g h t , [they] c o n s i d e r a b l y reduce c o n f u s i o n on b a l a n c e . " 1 The cepstrum i s d e f i n e d as the power spectrum of the l o g a r i t h m of the power spectrum of a s i g n a l . Oppenheim, Schafer and 1 Bogert et a l , 1963, "The Quefrency A l a n y s i s of Time S e r i e s f o r Echoes: Cepstrum, Pseudo-Autocovariance, Cross-Cepstrum and Saphe C r a c k i n g , " ,in Proc. Symp. on Time S e r i e s A n a l y s i s (N.Y.: John Wiley and Sons), p. 209. 9 Stockham (1968) added the word complex to emphasize the use of the complex transforms and the complex l o g a r i t h m i n the system D. We have i m p l i c i t l y assumed that the complex l o g a r i t h m maps from a m u l t i p l i c a t i v e space to an a d d i t i v e space. T h i s means that i s must be d e f i n e d such that l o g [ X , ( z ) - X 2 ( z ) ] = l o g [ X , ( z ) ] + l o g [ X 2 ( z ) ] (7) We have a l s o assumed that X(z) has an i n v e r s e z-transform. Thus X(z) must be the v a l i d z-transform of some sequence, x ( n ) . The unique d e f i n i t i o n of x(n) r e q u i r e s the s p e c i f i c a t i o n of a region of convergence of X ( z ) . We may r e s t r i c t x(n) and x(n) to be r e a l and such that the regions of convergence of X(z) and X(z) i n c l u d e the u n i t c i r c l e . That i s , X(z) and X(z) are a n a l y t i c i n a region i n c l u d i n g the u n i t c i r c l e . In the case that X(z) or X(z) i s not a n a l y t i c on the u n i t c i r c l e , an a l t e r n a t e contour of i n t e g r a t i o n may be used (Appendix A ) . Let us denote the e v a l u a t i o n of X(z) on the u n i t c i r c l e , z=exp(jw), by X ( j u ) . In terms of r e a l and imaginary p a r t s X ( J C J ) = Xr(jw) + j - X i ( j w ) (8) where j i s the imaginary u n i t , j 2 = -1. The requirement that x(n) be r e a l i m p l i e s that Xr(jo>) i s an even f u n c t i o n of w and Xi(jw) i s an odd f u n c t i o n of I t f u r t h e r i m p l i e s that X ( J C J ) 10 i s p e r i o d i c i n co with p e r i o d 2ir. A n a l y t i c i t y of X(z) on the u n i t c i r c l e i m p l i e s that X(jco) must be a continous f u n c t i o n of a). We have X(jw) = log[X(jo>)] (9a) X(jw) = log|X(j<j)| + j-arg[X( j<j) ] (9b) Comparison of equations (9a,b) with equation (8) i m p l i e s that Xr(jw) = log|X( jo>)| (10a) Xi (jw) - arg[x( jco) ] (10b) T h e r e f o r e , l o g | X ( j a ) | and arg[X(jcj)] must be continuous f u n c t i o n s of a. C o n t i n u i t y of Xr(jco) = log|X(jw)| i s assured, provided X(z) has no zeros on the u n i t c i r c l e , by the a n a l y t i c i t y of X(jco) on the u n i t c i r c l e . Note that zeros of X(z) become pole s of X ( z ) . C o n t i n u i t y of Xi(jw) = arg[X(jco)] depends on the d e f i n i t i o n of the complex l o g a r i t h m . Note that X(jco), the z-transform e v a l u a t e d on the u n i t c i r c l e , i s j u s t the d i s c r e t e time F o u r i e r transform. The e s s e n t i a l problem with uniqueness and a n a l y t i c i t y of the complex lo g a r i t h m r e s u l t from the ambiguity i n d e f i n i n g the argument (phase) of a complex f u n c t i o n . Any i n t e g e r m u l t i p l e of 2n may be added to the p r i n c i p a l phase without a f f e c t i n g the phase's being r e p r e s e n t a t i v e of that f u n c t i o n . 11 However, o n l y one c o n t i n u o u s phase can be r e p r e s e n t a t i v e of the f u n c t i o n . A l s o , i n g e n e r a l , on l y t h i s c o n t i n u o u s phase w i l l s a t i s f y the p r o p e r t y of a d d i t i o n s p e c i f i e d i n e q u a t i o n ( 7 ) . The procedure of f i n d i n g t h i s unique phase i s commonly termed phase unwrapping. T h i s t e r m i n o l o g y r e s u l t s from the f a c t t h a t , i n the n u m e r i c a l computat ion of a phase , the p r i n c i p a l v a l u e i s o b t a i n e d . Thus the phase i s s a i d to have wrapped back i n t o the p r i n c i p a l v a l u e ' s range of (-n,n). Phase unwrapping has been , and c o n t i n u e s to b e , a major f a c t o r in the c a l c u l a t i o n of complex c e p s t r a . From a t h e o r e t i c a l p o i n t of v iew, the c o n t i n u o u s phase can be u n i q u e l y d e f i n e d . The d e f i n i t i o n and computat ion of the unwrapped phase a re d e a l t w i th e lsewhere in t h i s t h e s i s . The requ i rements tha t X i ( jw ) be c o n t i n u o u s , odd and p e r i o d i c i m p l i c i t l y c o n s t r a i n the a l l o w a b l e input sequences x ( n ) . In p a r t i c u l a r , x(n) must be such that X ( jo ) has a phase of z e r o at u=it. I t must a l s o have a p o s i t i v e mean v a l u e . The former f o l l o w s d i r e c t l y from the above r e q u i r e m e n t s , n o t i n g that X i ( jw) i s the phase of X ( jw ) . The l a t t e r f o l l o w s from the oddness and c o n t i n u i t y r e q u i r e m e n t s . These imply tha t X i ( j w ) , and hence the phase of X(ja>), i s ze ro at the o r i g i n . Thus the v a l u e of X(joj) i s equa l to i t s magnitude and cannot be n e g a t i v e . The p o s s i b i l i t y of a ze ro magnitude i s p r e c l u d e d because X(z ) has no ze ros on the u n i t c i r c l e . Now, the v a l u e of X(jw) at the o r i g i n i s equa l to the mean v a l u e of x(n) which t h e r e f o r e must be p o s i t i v e . I t i s worth n o t i n g t h a t , fo r a sequence x(n) hav ing a n e g a t i v e mean v a l u e , the phase of 12 X(jo>) at CJ=0 i s an odd m u l t i p l e of ir and i s d i s c o n t i n u o u s . 2.4 The Inverse C h a r a c t e r i s t i c System The i n v e r s e c h a r a c t e r i s t i c system D"1 i s d e f i n e d i n equation ( 5 ) . I t transforms from an a d d i t i v e space to a c o n v o l u t i o n a l space and has the canonic r e p r e s e n t a t i o n shown i n F i g u r e 3. The transform i s r e a l i z e d by c a s c a d i n g a z-transform, complex e x p o n e n t i a l and i n v e r s e z-transform. The complex e x p o n e n t i a l i s unique and, i f X(z) i s a n a l y t i c on the u n i t c i r c l e , so i s e x p [ X ( z ) ] . Figure 3 - Representation of the System D"1 The inverse characterst1c system D"' maps from an additive space to a convolutional space. 2.5 The L i n e a r System The l i n e a r system L i s the only system i n H which i s u n s p e c i f i e d once D i s f i x e d . I t s s p e c i f i c a t i o n w i l l depend on the c h a r a c t e r i s t i c s of the i n p u t s , the mapping D and the intended outcome of the procedure. 13 2.6 P r o p e r t i e s Of The Complex Cepstrum We w i l l c o n s i d e r here only f i n i t e l e n g t h input sequences x(n) with p o s i t i v e index n. Such sequences have z-transforms with no s i n g u l a r i t i e s i n the z plane and can be represented as (Schafer, 1969) where |a(k)|<1 and |b(k)|<1. The a(k) are the m. zeros i n s i d e the u n i t c i r c l e and the b(k) are the m0 zeros o u t s i d e the u n i t c i r c l e . The f i r s t f a c t o r , A, i s a constant and the second f a c t o r r e p r e s e n t s a s h i f t of the input sequence by m0 p o s i t i o n s . Let us c o n s i d e r the f i r s t two f a c t o r s i n equation (11). For r e a l sequences A i s r e a l , and i f A i s p o s i t i v e i t c o n t r i b u t e s only to x(0) (Oppenheim and Schafer, 1975). I f A i s n e g a t i v e , i t s c o n t r i b u t i o n to x(n) i s more complicated. I t can be shown ( T r i b o l e t , 1979) that the s i g n of A i s equal to the s i g n of the mean value of x ( n ) . Thus, n o r m a l i z i n g the input by m u l t i p l i c a t i o n by +1 or -1 to make A p o s i t i v e i s c o n s i s t e n t with s a t i s f y i n g the c o n t i n u i t y requirement of Xi(jcj) at the o r i g i n . In the second f a c t o r , the exponent of z, m0, g i v e s r i s e to a term x,(n) in the complex cepstrum of X(z) (11) 14 X i ( n ) = -[m 0 «cos(irn) ]/n (12) ( U l r y c h , 1971). I f X(z) has many zeros o u t s i d e the u n i t c i r c l e , m0 can be l a r g e and ic, (n) may dominate the complex cepstrum. Thus i t i s important to s h i f t the input sequence to remove t h i s term. T h i s i s e q u i v a l e n t to removing the l i n e a r phase component of X(jco) and i s c o n s i s t e n t with s a t i s f y i n g the requirements that Xi(jw) be continuous and p e r i o d i c . S e v e r a l p r o p e r t i e s of the complex cepstrum have been examined (Oppenheim and Scha f e r , 1975; T r i b o l e t , 1979). Those most r e l e v a n t are o u t l i n e d below. 1. The complex cepstrum decays at l e a s t as f a s t as 1/n. That i s where C i s a constant and d i s the maximum of |a(k)| and 2. I f x(n) i s minimum phase (no zeros o u t s i d e the u n i t c i r c l e ) then |x(n)| < C-|d n/n| - oo <n< ao |b(k)|. x(n) = 0 n<0 and x(n) can be c a l c u l a t e d r e c u r s i v e l y , d i r e c t l y from x ( n ) . 3. If x(n) i s maximum phase (no zeros i n s i d e the u n i t c i r c l e ) then 15 x(n) = 0 , n>0 and x(n) can be c a l c u l a t e d r e c u r s i v e l y , d i r e c t l y from x ( n ) . 4. I f x(n) i s of f i n i t e d u r a t i o n , x(n) w i l l n e v e r t h e l e s s be of i n f i n i t e d u r a t i o n . 5. The complex cepstrum of a sequence x(n) whose spectrum i s smooth tends to concentrate near n=0. T h i s p r o p e r t y i s a r e s u l t of the f a c t t h at a smooth spectrum r e s u l t s from a sequence whose zeros are f a r from the u n i t c i r c l e . Thus |a(k)| and |b(k)| are sma l l and p r o p e r t y 1 i m p l i e s that x(n) decays r a p i d l y with n. 6. Adding an impulse at the o r i g i n to a cepstrum i s e q u i v a l e n t to s c a l i n g the time sequence. Consider a complex cepstrum which i s non-zero only at the o r i g i n and has magnitude A. Then a p p l y i n g the i n v e r s e t r a n s f o r m y i e l d s a time sequence which i s exp[A] at the o r i g i n and zero elsewhere. Since a d d i t i o n maps to c o n v o l u t i o n and t h i s c o n v o l u t i o n merely s c a l e s the sequence by exp[A], the pr o p e r t y f o l l o w s . Note that the s c a l e f a c t o r i s always p o s i t i v e . R e c a l l that we have excluded input sequences having zeros on the u n i t c i r c l e . In the case that x(n) has f i n i t e d u r a t i o n , X(z) has a r e g i o n of convergence which i n c l u d e s the u n i t c i r c l e . In the general case the contour of i n t e g r a t i o n of X(z) can be changed by e x p o n e n t i a l l y weighting the input sequence (Appendix A). That i s , x(n) becomes w(n) = x(n)-expfa-n] 16 where a i s a r e a l c o n s t a n t . The sequence w(n) has the z-transform W(z) = X(z-exp[-a]) or W(z) = X(z'b) where b=exp[-a]. Thus, the zeros and poles of X(z) are moved r a d i a l l y by the f a c t o r b. I f we have a c o n v o l u t i o n a l input, x ( n ) = x t ( n ) * x 2 ( n ) , then x(n)*exp[a*n] = x,(n)•exp[a-n]*x 2(n)«exp[a«n] Thus, e x p o n e n t i a l weighting of the c o n v o l u t i o n a l of sequences i s e q u i v a l e n t to the c o n v o l u t i o n of e x p o n e n t i a l l y weighted sequences. E x p o n e n t i a l weighting can be used to make a mixed phase sequence e i t h e r minimum phase or maximum phase by moving a l l i t s zeros and poles e i t h e r i n s i d e the u n i t c i r c l e or out s i d e the u n i t c i r c l e . Then the s p e c i a l p r o p e r t i e s of the minimum and maximum phase sequences can be e x p l o i t e d . 1 7 2.7 Computational 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 We have d i s c u s s e d the c h a r a c t e r i s t i c system, D, i n terms of continuous transforms. In f a c t , the z-transform, when eva l u a t e d on the u n i t c i r c l e , i s j u s t the F o u r i e r transform. In p r a c t i c e we cannot compute continous transforms and must use a d i s c r e t e r e p r e s e n t a t i o n . We use the d i s c r e t e F o u r i e r transform to evaluate samples of X(jw). T h i s computation lead s to an a l i a s e d cepstrum (Oppenheim and S c h a f e r , 1975). That i s , i f x(n) i s the true cepstrum, the c a l c u l a t e d cepstrum, x 0 ( n ) , w i l l be x 0 ( n ) = "2.x(n+k«N) where N i s the number of p o i n t s i n the o r i g i n a l sequence x ( n ) . S i n c e , i n g e n e r a l , x(n) i s of i n f i n i t e d u r a t i o n , x 0 ( n ) w i l l be an a l i a s e d v e r s i o n of x ( n ) . However, as x(n) decays e x p o n e n t i a l l y , appending zeros to the input sequence w i l l i n c r e a s e N and reduce the a l i a s i n g e r r o r . Schafer (1969) has shown that one may reduce c e p s t r a l a l i a s i n g by e x p o n e n t i a l l y weighting the input sequence. T h i s has the e f f e c t of e x p o n e n t i a l l y weighting the complex cepstrum, causing i t to decay more r a p i d l y . Examples of t h i s e f f e c t may be found i n S t o f f a et a l . ( l 9 7 4 ) . The most d i f f i c u l t p a r t , i n p r a c t i c e , i n c a l c u l a t i n g the complex cepstrum l i e s i n the computation of the argument or phase of X(ja>). We must compute samples of the continous phase of X(jo>) to d e f i n e Xi(jcj) as in equation (10b). If we 18 use the r e l a t i o n ARG[X(jco)] = t a n " 1 [Xi (jw)/Xr (jco) ], where Xi(j w ) and Xr(jco) represent the imaginary and r e a l p a r t s of X(jcu), we o b t a i n the p r i n c i p a l value of arg[X(ja>)]. In g e n e r a l , the p r i n c i p a l value of the phase i s d i s c o n t i n u o u s and thus w i l l not s u f f i c e . Methods of computing samples of the continuous phase are d e a l t with i n the next c h a p t e r . i For the s p e c i a l case that the input sequence i s e i t h e r minimum or maximum phase, the complex cepstrum can be c a l c u l a t e d , without a l i a s i n g , d i r e c t l y from x ( n ) . Oppenheim and Schafer (1975) pr o v i d e d e t a i l s of the method used. 2.8 Summary We have d i s c u s s e d a n o n - l i n e a r homomorphic transform which maps c o n v o l u t i o n s i n t o a d d i t i o n s . T h i s transform may e f f e c t the s e p a r a t i o n of convolved s i g n a l s . S e v e r a l p r o p e r t i e s of t h i s transform were presented. There are r e s t r i c t i o n s on the a l l o w a b l e input s i g n a l s due to the use of complex logarithms and i n v e r s e z-transforms. There are computational problems i n computing the complex l o g a r i t h m which w i l l be d e a l t with i n the next chapter. 19 I I I . P H A S E UNWRAPPING 3.1 I n t r o d u c t i o n The problem of phase unwrapping can be approached from many d i r e c t i o n s . Let us use a g e o m e t r i c a l approach to e l u c i d a t e the concepts. Say we have a r e a l sequence, x ( n ) , with F o u r i e r t r a n s f o r m X(w). X(co) i s a complex valued f u n c t i o n of co and may be represented by a path i n the complex plane. Let Xr(cp) and Xi (co) represent the r e a l and imaginary p a r t s of X(co) and l e t them d e f i n e the axes of a complex plane as i n F i g u r e 4. As to i n c r e a s e s from zero, X(co) t r a c e s out a path i n the complex plane as shown. The p r i n c i p a l phase of X(co) , A R G [ X ( C J ) ] , i s d e f i n e d as the angle a s t r a i g h t l i n e from the o r i g i n to X(a>) makes with the p o s i t i v e r e a l a x i s . T r a d i t i o n a l l y t h i s angle i s measured p o s i t i v e upward and negative downward. The angles +ir and -ir c o i n c i d e on the n e g a t i v e r e a l a x i s . I f the path of X(co) c r o s s e s the negative r e a l a x i s , A R G[X ( c o ) ] undergoes a jump d i s c o n t i n u i t y of 27r. I f X(co) i s moving down, the jump i s from +IT to -it. I t i t i s moving up, the jump i s from -it to +ir. We wish to d e f i n e a unique phase, a r g [ ( X ( w ) ] , which changes c o n t i n u o u s l y everywhere i n c l u d i n g the negative r e a l a x i s . T h i s phase w i l l g e n e r a l l y be unbounded and w i l l i n c r e a s e or decrease c o n t i n u o u s l y a c r o s s t h i s a x i s . arg[X ( co ) ] can be d e f i n e d by adding to ARG[X(co)] an i n t e g e r number of m u l t i p l e s of 2TT. That i s 20 Xi(uj) Figure 4 - Complex Function in the Complex Plane As u increases, F(u) traces out a path in the complex plane. arg[X(o>)] = ARG[X(a>)] + L(G>)-2TT where L(a>) i s an i n t e g e r which makes a r g [ « ] c o n t i n u o u s . L ( u ) can change on l y a t d i s c o n t i n u i t i e s in ARG[•]. These occur on l y when Xr(w)=0 and Xi(cj) changes s i g n . The c o n t i n u o u s phase , a rg tX (c j ) ] , i s c a l l e d the unwrapped phase of X ( C J ) and the a c t of f i n d i n g L(o>) i s c a l l e d phase unwrapping. Note t h a t i f X(o>) = 0 f o r some u> t h e r e can be no unique p h a s e . T h i s problem i s i n t i m a t e l y r e l a t e d to the unique d e f i n i t i o n of the complex l o g a r i t h m . T h i s f o l l o w s as arg[X(w)] i s the imaginary p a r t of the l o g a r i t h m of X ( C J ) . I t i s a l s o r e l a t e d to the requi rement t h a t , i f two complex f u n c t i o n s a r e . m u l t i p l i e d , t h e i r phases w i l l be added. In 21 gen e r a l ARGtX.tw)] + A R G [ X 2 ( G > ) ] * A R G[X.(w) + X 2(o>)] but arg[X,(w)] + arg[X 2(o>)] = arg[X, + X 2(w)] Note that i f one uses the d e f i n i t i o n ARG[X(w)] = tan" 1 [Xi (w)/Xr (OJ) ] then ARG[ • ] i s d e f i n e d on (-7r/2,7r/2) and the d i s c o n t i n u i t i e s occur at p o i n t s on the imaginary a x i s . T h i s f o l l o w s because the r a t i o X i (w)/Xr(CJ) does not have the s i g n i n f o r m a t i o n of i t s c o n s t i t u e n t s and thus tan _ 1 [ ' 3 i s p r o j e c t e d onto the r i g h t h a l f plane. Then L(o>) can change when X(w) c r o s s e s the imaginary a x i s (Xr(w)=0). An example of the p r i n c i p a l value of a phase and the corresponding unwrapped phase are shown i n Fi g u r e 5. Va r i o u s approaches have been put f o r t h to determine the unwrapped phase. Schafer (1969) used an a l g o r i t h m which searched f o r d i s c o n t i n u i t i e s i n the p r i n c i p a l phase, A R G ( C J ) , and removed them by a d d i t i o n of n-27r using the a p p r o p r i a t e n. A technique to e x p l o i t the a n a l y t i c d e f i n i t i o n of the phase as the i n t e g r a l of i t s d e r i v a t i v e was developed by T r i b o l e t (1977). Bhanu (1977) used a d i f f e r e n t technique to e x p l o i t 22 ARG[X(CJU)] —n _ -2n _ F i g u r e 5 - P r i n c i p a l V a l u e 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 i n t o ( - T T ) , T T ) w h i l e the unwrapped phase ( a r g ) Is c o n t i n u o u s . t h i s same d e f i n i t i o n . S t e i g l i t z and Di c k i n s o n (1982) proposed polynomial f a c t o r i z a t i o n as a means of computing the phase. A novel approach, of some t h e o r e t i c a l i n t e r e s t , has been developed by. McGowan and Kuc (1982) using number theory 23 a p p l i e d to polynomials. These methods w i l l be presented b r i e f l y . The method of McGowan and Kuc w i l l be d i s c u s s e d i n more d e t a i l to e l u c i d a t e the concepts i n v o l v e d . Some d i s c u s s i o n of the p r o p e r t i e s of f u n c t i o n s which a f f e c t phase unwrapping w i l l be presented. 3.2 P r i n c i p a l Value The method of phase unwrapping by p r o c e s s i n g the p r i n c i p a l value i s due to Schafer (1969). The r e l a t i o n between the p r i n c i p a l value, ARG(w), and the continuous phase, arg(OJ) , i s A R G ( C J ) = {arg(w) }M0D2TT That i s , the continuous phase, modulo 2TT, i s the p r i n c i p a l phase. Schafer's a l g o r i t h m searches f o r d i s c o n t i n u i t i e s i n ARG(a>) caused by the modulo o p e r a t i o n and removes them by the a p p r o p r i a t e a d d i t i o n or s u b t r a c t i o n of 2ir. The a l g o r i t h m s t a r t s at arg(O), which i s 0 f o r r e a l x ( n ) , and searches f o r d i s c o n t i n u i t i e s at i n c r e a s i n g CJ. The user must s p e c i f y a t h r e s h o l d w i t h i n which adjacent samples must l i e f o r the phase to be c o n s i d e r e d continuous. In g e n e r a l , the a l g o r i t h m becomes, more accurate as the spacing between adjacent samples i s reduced. Any m i s - i d e n t i f i e d d i s c o n t i n u i t i e s w i l l a f f e c t the unwrapping at l a r g e r w. More r e c e n t l y , P o g g i a g l i o l m i et a l . (1982) proposed another 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. I t a l s o 24 searches f o r d i s c o n t i n u i t i e s , but a d a p t i v e l y r e c a l c u l a t e s phase val u e s between the o r i g i n a l values of the F o u r i e r transform. They suggest s h i f t i n g the input sequence such that i t s f i r s t moment i s zero before c a l c u l a t i n g the transform. T h i s i s an attempt to remove part of the l i n e a r phase component and may reduce the number of d i s c o n t i n u i t i e s . 3.3 I n t e g r a t i o n Of The D e r i v a t i v e V a r i o u s methods to compute the unwrapped phase use a numerical s o l u t i o n of an a n a l y t i c e x p r e s s i o n . E s s e n t i a l l y , these methods are attempts to do numerical i n t e g r a t i o n with c o n s t r a i n t s . Let x(n) be a sequence and X(co) be i t s F o u r i e r transform. Then, as i n the previous chapter, we d e f i n e where the complex l o g a r i t h m i s d e f i n e d as i n the p r e v i o u s chapter and equation (2b) represents a proper F o u r i e r transform. The d e r i v a t i v e of X(co) i s w e l l d e f i n e d and i s given by X(co) = log[X(co)] (2a) or X(co) = log|X(cu)| + j »arg[X(co) ] (2b) X* (co) = d{log[X(co) ]}/dco (3a) 25 X' ( C J ) = log'|X(cj)| + j«arg'[X(«)] (3b) X ' (w) = X ' (C J ) / X ( C J ) (3c) where the prime denotes d i f f e r e n t i a t i o n with respect to C J . Comparison of equations (3b) and (3c) shows that the d e r i v a t i v e of arg[X(cj)] i s equal to the imaginary p a r t of X' ( C J ) / X ( C J ) . I t can be shown (Oppenheim and Schafer, 1975) that arg'[X(cj)] = {Xr (CJ)-Xi ' (CJ ) - X i ( C J )-Xr ' (CJ) }/| X ( C J ) | 2 (4) The phase, arg[X(cj)], can be unambiguously d e f i n e d as with i n i t i a l c o n d i t i o n arg[X(0)]=0. As mentioned i n the prev i o u s chapter, a r g [ X ( t o ) ] must be continuous and odd. Equation (5) e s t a b l i s h e s c o n t i n u i t y . Oddness w i l l r e s u l t i f the mean phase d e r i v a t i v e i s zero. That i s If t h i s i s not t r u e , removal of the l i n e a r phase component w i l l make i t so. I m p l i c i t y , i t has been assumed that we have independent knowledge of ar g ' [ X ( c j ) ] . In f a c t , t h i s can be computed d i r e c t l y from x(n) using the r e l a t i o n ( T r i b o l e t , (5) o ( 6 ) o 26 1977) Xr'(co) + j-Xi'(co) = -j.F{n-x(n)} (7) where F{«} denotes the F o u r i e r transform. Thus, computing arg'[X(co)] by using equations (7) and (4) and i n t e g r a t i n g as i n equation (5) y i e l d s the unwrapped phase. The computational problem i s the numerical i n t e g r a t i o n of equation (5) . The most s t r a i g h t - f o r w a r d s o l u t i o n i s to use a t r a p e z o i d a l r u l e f o r i n t e g r a t i o n . Again, t h i s method improves as the spacing between adjacent samples of arg'[X(co)] decreases. T h i s approach a l s o s u f f e r s from e r r o r propagation i n that an e r r o r i n i n t e g r a t i o n at u>=u, w i l l a f f e c t the r e s u l t at co>co, . A more s o p h i s t i c a t e d method of c a r r y i n g out the i n t e g r a t i o n was developed by T r i b o l e t (1977) and i s c a l l e d a daptive i n t e g r a t i o n . The p r i n c i p a l phase ARG(co) as w e l l as the phase d e r i v a t i v e , arg' (co) are c a l c u l a t e d . The phase d e r i v a t i v e i s then i n t e g r a t e d as i n (5) u s i n g a t r a p e z o i d a l r u l e between, say, co, and co 2. T h i s computed value of the phase i s then compared to the p r i n c i p a l v a l u e . The c o n s t r a i n t that these v a l u e s must be c l o s e to each other to w i t h i n an i n t e g e r m u l t i p l e of 2ir i s added. That i s | (ARG[X(co) ] + 2 T T - L ( C J ) } - arg[ (X(co) ] | < E for some i n t e g e r L(co), where E i s some small p o s i t i v e number 27 s u p p l i e d by the user. The a l g o r i t h m checks va l u e s of L(co) to s a t i s f y the i n e q u a l i t y . I f no L(co) s u f f i c e s , another p o i n t of the phase d e r i v a t i v e i s c a l c u l a t e d midway between co, and co 2. The t r a p e z o i d a l r u l e i s a p p l i e d between co, and (co2-co,)/2 and the procedure i t e r a t e d u n t i l the i n e q u a l i t y i s s a t i s f i e d . The procedure i s i n i t i a t e d by t a k i n g a DFT of x(n) and n«x(n) to ob t a i n ARG[X(co)] and arg'[X(co)]. Then the i n t e g r a t i o n s t a r t s at co=0 and proceeds to higher co, r e c a l c u l a t i n g p o i n t s as needed. T h i s method a l s o s u f f e r s from e r r o r propagation. An e r r o r i n determining L(co,) w i l l c a r r y through f o r co>co, . Bhanu (1977) i n v e s t i g a t e d T r i b o l e t ' s a l g o r i t h m and compared i t to a l t e r n a t e t e c hniques. He a l s o proposed and i n v e s t i g a t e d the use of piecewise polynomial i n t e r p o l a t i o n of the phase d e r i v a t i v e as a b e t t e r r u l e of i n t e g r a t i o n . T h i s r e q u i r e s the e v a l u a t i o n of the second d e r i v a t i v e of the phase, computed from F{n 2«x(n)} as w e l l as the f i r s t d e r i v a t i v e and the p r i n c i p a l v a l u e . There i s a t r a d e - o f f of more i n i t i a l computation f o r in c r e a s e d accuracy i n the i n t e g r a t i o n . Note t h a t , although the p r i n c i p a l value of the arctangent i s c o n tained i n ( - T T/2 , TT/2 ) , by r e t a i n i n g the s i g n i n f o r m a t i o n of the terms i n tan" 1 [Xi (co)/Xr (co) ] a value unique w i t h i n ( - 7 T , < T ) can be assigned. 28 3.4 F a c t o r i z a t i o n S t e i g l i t z and D i c k i n s o n (1982) proposed the use of polynomial f a c t o r i z a t i o n of the z-transform as a " r e l i a b l e and a c u r a t e " 1 method of phase computation. Once the r o o t s of a z-transform are known, i t can be f a c t o r e d i n t o the product of polynomials of degree one ( d i p o l e s ) . These d i p o l e s may then be e v a l u a t e d on the u n i t c i r c l e to form the F o u r i e r transform. The phase of the product i s equal to the sum of the phases of each of the d i p o l e s . The continuous phase of each d i p o l e , and thus the product, can be c a l c u l a t e d unambiguously. With t h i s method, the problem becomes that of f i n d i n g r o o t s of a complex polynomial. S t e i g l i t z and D i c k i n s o n used a Newton-Raphson r o o t - f i n d i n g a l g o r i t h m . They ignored the problem of repeated r o o t s and saddle p o i n t s on the grounds that these r a r e l y occur i n p r a c t i c e . 3 . 5 Number Theory McGowan and Kuc (1982) proposed a method of d e f i n i n g the continuous phase based on number theory. The theory w i l l be d i s c u s s e d i n d e t a i l and r e l e v a n t p r o o f s presented. In l i n e a r system c o n t r o l theory the problem of s t a b i l i t y i s important. There are methods which w i l l determine the s t a b i l i t y of a system without a c t u a l l y f i n d i n g which p a r t of the system i s u n s t a b l e . A p a r t i c u l a r method of determining 1 S t e i g l i t z and D i c k i n s o n , 1982, "Phase Unwrapping by F a c t o r i z a t i o n , " I.E.E.E. Trans. A.S.S.P., v. ASSP-30, No. 6, p. 984. 29 s t a b i l i t y g i v e s i n f o r m a t i o n about the number of ro o t s of a polynomial which are i n s i d e or o u t s i d e the u n i t c i r c l e without a c t u a l l y l o c a t i n g them. T h i s same method may be used to determine the continuous phase. Consider a r e a l sequence x ( n ) , n=0,1,...N-1 with the d i s c r e t e F o u r i e r transform (DFT) X(co) = ^ x ( n ) - e x p [ - j - n c o ] (8) to The phase can then be w r i t t e n as arg[X(co)] = tan" 1 [Xi (co)/Xr (cu) ] + L ( C O ) . T T f o r 0<CO<VT. t a n _ 1 [ * ] i s the p r i n c i p a l value arctangent, that i s -7r/2<tan" y<ir/2, and L(co) i s an i n t e g e r f u n c t i o n of co which makes arg[»] a continuous f u n c t i o n . F i n d i n g L(co) i s the phase unwrapping problem. Consider the arctangent f u n c t i o n i n Fig u r e 6. Consider a p o i n t on the p r i n c i p a l branch, L=0, d e f i n e d by F(co) = Xi (co)/Xr (co) (and assuming Xi(co) to have no s i n g u l a r i t i e s ) . T h i s p o i n t can change branches to L=1 only i f F(co) goes to + oo . I t can then 'wrap' onto the L=1 branch at - oo without e x h i b i t i n g any d i s c o n t i n u i t y . L i k e w i s e , i f F(co) goes to - c o , t h i s p o i n t can change to the L=-1 branch. The arctangent can change branches only at s i n g u l a r i t i e s of F(co), or e q u i v a l e n t l y , at zeros of Xr(co). However, a branch change i s not r e q u i r e d at these s i n g u l a r i t i e s . F(co) may go to i n f i n i t y and r e t u r n on the same 3 0 arctan(x) + Ln n 0 L=1 J r r n _ L=o J L=-1 J r - 2 0 0 20 X F i g u r e 6 - Branches of a r c t a n g e n t F u n c t i o n The f u n c t i o n can c o n t i n u o u s l y change branches a t I n f I n l t 1 e s . branch. Thus, whether or not F ( C J ) changes s i g n and the d i r e c t i o n of the sig n change p r o v i d e s the i n f o r m a t i o n necessary to determine L ( C J ) . Let us assume that X ( C J ) * 0 on 0^CJ<7T. T h i s i s the c o n s t r a i n t t h a t there are no zeros on the u n i t c i r c l e i n the z-plane. The product X i ( C J ) «Xr ( C J ) has the same s i g n as the r a t i o X i ( C J ) / X r ( C J ) but has no s i n g u l a r i t i e s . Consider the s i g n of X i ( C J ) «Xr ( C J ) as CJ i n c r e a s e s through a zero of Xr(cj). (Note that Xi(cj) doesn't change s i g n f o r CJ c l o s e enough to the zero of Xr(cj).) I f the product's s i g n changes from p o s i t i v e to negative, L ( C J ) i s i n c r e a s e d by one. If an opposite sig n change occurs. L ( C J ) i s decreased by one. If there i s no sig n 31 change L ( C J ) i s not changed. L ( C J , ) i s generated, i n p r i n c i p l e , by l e t t i n g CJ i n c r e a s e from CJ=0 to CJ=CJ. and incrementing or decrementing L ( C J ) as r e q u i r e d . The c o n d i t i o n arg[0] = 0 p r o v i d e s the s t a r t i n g v a l u e , L(0)=0. Now the problem of phase unwrapping has become that of f i n d i n g the s i g n changes of X i ( C J ) «Xr ( C J ) through the zeros of Xr(cj). A s o l u t i o n to t h i s problem can be obtained using Sturm's theorem (Beaumont and P i e r c e , 1963). Sturm's theorem a r i s e s i n number theory and p r o v i d e s a method of f i n d i n g the number of d i s t i n c t r e a l r o o t s of a polynomial between two arguments. The method i n v o l v e s the g e n e r a t i o n of a Sturm sequence, a sequence of polynomials of d e c r e a s i n g degree, generated r e c u r s i v e l y from two given polynomials (Marden, 1966). McGowan and Kuc use the r e a l and imaginary p a r t s o f , t h e DFT to s t a r t the sequence. They found i t convenient to express the Sturm sequence i n terms of Chebyshev polynomials. The Sturm sequence p r o v i d e s s u f f i c i e n t i n f o r m a t i o n to determine L ( C J ) and hence the unwrapped phase. The method proceeds as f o l l o w s . The r e a l and imaginary p a r t s of X ( C J ) are expressed i n terms of Chebyshev polynomials of the second ki n d . Chebyshev polynomials of the f i r s t k ind, T ( n ) , and of the second kind, U(n), as f u n c t i o n s of C O S ( C J ) are d e f i n e d as (Snyder, 1966, p.11-24) T(n) = T(n,coscj) = cos(n«cj) nZO (9) U(n) = U(n,coscj) = s i n [ (n+1 )cj]/sin(w) n£0 (10) 32 Relevant r e c u r s i o n r e l a t i o n s between these polynomials a r e , f o r n£2, T(n) = [U(n) - U(n-2)] (11a) T(1) = U(1) (11b) T(0) = U(0) =1 (11c) For n>1 there i s the r e c u r s i o n r e l a t i o n U(n+1) = U ( n ) - U d ) - U(n-1) (I1d) where i t i s understood that U(n) and T(n) are f u n c t i o n s of cos(co). (Note that Chebyshev polynomials are an orthogonal set (Snyder, 1966)). From the i d e n t i t i e s (9) and (10), X(co) i n (8) can be expressed as X(co) = "2,a(n)-U(n) + j • s i n (co) ^ b ( n) • U (n ) (12) ft --o * where the a(n) are l i n e a r combinations of x(n) and b(n) = -x ( n ) . D e t a i l s on c a l c u l a t i n g a(n) and b(n) are given i n Appendix B. Equation (12) may be w r i t t e n as X(co) = Xr(co) + j-Xi(co) (13a) X(co) = P(0,co) + j-sin(co) «P( 1 ,co) (13b) 33 X (CJ ) = P(0) + j -s in(c j ) -P( 1 ) (13c) where i t i s unders tood t h a t P(0) and P(1) r e p r e s e n t p o l y n o m i a l s which are f u n c t i o n s of CJ . That i s P(0) = P (0,6j) = Xr(cj) = S"a(n) -U(n . C O S C J ) (14a) Note that the s i g n of P(1) i s the same as tha t of Xi(cj) on 0<CJ<7T. P(0) and P(1) a re p o l y n o m i a l s in U ( n ) . They a re used to generate a Sturm sequence , denoted by { P ( 0 ) , P ( 1 ) , . . . , P ( M ) } where M^N-1. T h i s Sturm sequence i s a sequence of p o l y n o m i a l s of d e c r e a s i n g d e g r e e . Each p o l y n o m i a l i s a f u n c t i o n of C J . The sequence i s genera ted u s i n g the ' n e g a t i v e remainder ' r e l a t i o n s h i p (Marden, 1966) . In f a c t , t h i s r e l a t i o n s h i p i s merely E u c l i d ' s a l g o r i t h m f o r f i n d i n g the g r e a t e s t common d i v i s o r of i n t e g e r s or p o l y n o m i a l s (Beaumont and P i e r c e , 1963) . The r e l a t i o n s h i p i s and P ( 1 ) = P ( 1 , C J ) = Xi(cj)/sincj = 2.b(n) -U( n, C O S C J ) (14b) P ( k - 1 ) = Q ( k ) - P ( k ) - P(k+1) (15) where 1<k^M-1. P(k) and Q(k) are p o l y n o m i a l s of degree k. The degree of the p o l y n o m i a l s P and Q i s d e f i n e d to be that of 34 the highest degree Chebyshev poly n o m i a l , U, i n t h e i r sum. I f a P c o n t a i n s a term U(K) where K i s the h i g h e s t degree i n the sum, P i s of degree K. The r e c u r s i o n (15) proceeds, given P(k-1) and P(k), by f i n d i n g Q(k) and P(k+1) such that deg{Q(k)} = deg{P(k-l)} - deg{P(k)} > 0 (16) and deg{P(k+1)} < deg{P(k)} The e x i s t e n c e of Q(k) and P(k+1) are guaranteed (Beaumont and P i e r c e , 1963) and, i n f a c t , t h i s a l g o r i t h m i s e q u i v a l e n t to the long d i v i s i o n of p o l y n o m i a l s . We choose Q(k) such that (16) i s s a t i s f i e d and such that the h i g h e s t degree term of P(k-1) i s i d e n t i c a l to that of Q(k)«P(k). I t may happen that other terms of l e s s e r degree are a l s o i d e n t i c a l . T h i s r e s u l t s i n the e l i m i n a t i o n of these terms. That i s , w r i t i n g (15) as P(k+1) = Q(k)-P(k) - P(k-1) then P(k+l) w i l l be of degree l e s s than P(k-1). For example, i f P(0) i s degree N and P(1) i s of degree N-R, then we take Q(1) of degree N-(N-R) = R so that Q(1)-P(1) i s of degree N and 35 P (2 ) = Q ( 1).P ( 1 ) - P (0 ) The h i g h e s t degree term i n P (0 ) i s c a n c e l l e d by the hi g h e s t degree term i n Q ( 1)«P(l). ( P o s s i b l y some other terms w i l l c a n c e l , but t h i s i s not guaranteed.) The a l g o r i t h m c o n t i n u e s , d i s c a r d i n g Q, by i n c r e a s i n g k u n t i l k=M-1 where P (M+1)=0. Then P ( M -1 ) = Q (M)«P (M) and P (M) i s the g r e a t e s t common d i v i s o r of P (0 ) and P ( 1 ) . I f M=N -1, P(N -1 ) = P (M) i s a con s t a n t . D e t a i l s of the r e c u r s i o n are given i n Appendix C. Let us d e f i n e the operator V(co) which, when a p p l i e d to the Sturm sequence, g i v e s the number of sign changes i n the sequence f o r a f i x e d co. TO c a l c u l a t e V ( c o ) , f i x an co and count s i g n changes from adjacent P(n) i n (P ( 0 ),P(1),...,P ( M )}. I f adjacent P(n) have opposite s i g n s V(co) i s incremented by one, otherwise i t i s unchanged. The number of s i g n changes i n the Sturm sequence, as a f u n c t i o n of co, w i l l determine L ( c o ) , the f u n c t i o n r e q u i r e d to unwrap the phase. S p e c i f i c a l l y , L(co) = V(co)„ . T h i s f o l l o w s from the p r o p e r t i e s of the Sturm sequence and i s proved below. We r e q u i r e two theorems to l e a d up to the f i n a l r e s u l t . These are presented here with p r o o f s . Theorem; I f P(0 ,co) and P(1 ,co) have no common zeros f o r 0<co<7r, then sign[P ( M , c o ) ] , where P (M) i s the g r e a t e s t common .divisor of P (0 ) and P ( 1 ) , cannot change f o r 0<CO<TT. Proof; P (M) was generated by E u c l i d ' s a l g o r i t h m : P (0) = Q(1)-P( 1 ) - P (2) 36 P(1) = Q (2 )-P (2 ) - P (3 ) P(M-2) = Q(M -1)-P(M -1) - P(M) P(M-1) = Q(M)-P(M) If P(M)=0 f o r some CJ, then t h i s r e c u r s i o n i m p l i e s , by i n d u c t i o n , P(M -1)=0, P(M -2)=0, P (1 )=0 , P ( 0 ) = 0 . But, by hyp o t h e s i s , P (1 ) and P (2 ) have no common zeros on 0<CJ<7T. Therefore P(M )*0 f o r any CJ on 0<cj<ir and cannot change s i g n . Q.E.D. Theorem: If P ( 0 , C J ) and P ( i , c j ) have no common zeros, then V{P(k,cj)} does not change, as a f u n c t i o n of CJ, except at a zero of P ( 0 , C J ) . Proof: Say V{P(k ,cj)} changes. T h i s can only happen at a zero of some P(k,a>), say C J = C J 0 . That i s , P(k , c j o ) = 0. E u c l i d ' s a l g o r i t h m then i m p l i e s (equation ( 5 ) ) P(k -1 , C J 0 ) = -P(k + 1 , C J 0 ) or P(k -1 , C J 0 ) -P(k+1 , C J 0 ) < 0 That i s , members of the Sturm sequence separated by one other member have opposite s i g n s . (This must be tr u e because otherwise P ( k - 1 , c j 0 ) = P(k+1,c j 0 ) = 0 which would imply P ( J , C J O ) = 0 f o r every J>k i n c l u d i n g J=M which cannot be by the 37 pre v i o u s theorem.) The r e f o r e V{P(k,cj)} cannot change through a zero of P(k,w), 0<k<M. Q.E.D. Th i s theorem can be i l l u s t r a t e d by c o n s i d e r i n g the s i g n s of {P(k-1) ,P(k) ,P(k+1 )} near CJ=CJ 0 . Say they are {-, + , + }. Then i f P(k) changes s i g n , from + to -, the sequence becomes {-,-,+} and the number of si g n changes i s the same ( i . e . 1). The other p o s s i b l e cases are {-,-,+} becoming {-,+,+} {+,-,-} becoming {+,+,-} {+,+,-} becoming {+,-,-} By i n s p e c t i o n we see that the number of s i g n changes i s unchanged. I f more than one of the P(n) are zero the above a p p l i e s to each P(n) i n d i v i d u a l l y . Note that P(k) and P(k+1) cannot both be zero f o r some u>0 as the d i v i s i o n a l g o r i t h m would imply P(0) and P(1) are a l s o both zero, which i s co n t r a r y to our h y p o t h e s i s . Note that when P(k) i s a c t u a l l y zero, s i g n changes cannot be counted by comparing the signs of adjacent terms as zero has no s i g n . However, we can ignore the 0 and compare the si g n s of the terms on e i t h e r s i d e (Beaumont and P i e r c e , 1963). T h i s f o l l o w s d i r e c t l y from E u c l i d ' s a l g o r i t h m . Now l e t us co n s i d e r the f i r s t term of the sequence, P(0,w), as w passes through a zero, C J 0 , of P(0,O J). That i s P(0,tj o)=0. Since P(1,O J o)*0 by hyp o t h e s i s , three s i t u a t i o n s of i n t e r e s t can occur. F i r s t , P(0)*P(1) does not change s i g n , and V(C J) does not change. Second, P(0)«P(1) changes from - to 38 + and V(co) decreases by one. That i s , P(0) and P (1 ) change from having o p p o s i t e s i g n s to having the same s i g n . T h i r d , P ( 0)-P (1 ) changes from + to - and V(co) i n c r e a s e s by one. That i s , P(0) and P (1 ) change from having the same s i g n to having o p p o s i t e s i g n s . By the p r e v i o u s theorem these are the only c o n d i t i o n s under which V(co) can change. T h e r e f o r e , V (co 0 ) i s equal to the number of times P ( 0)«P (1 ) changes from + to -, minus the number of times i t changes from - to +, as co i n c r e a s e s from 0 to co 0 . Because P ( 0 , co) = Xr(co) and P(1 ,co) = X i ( c o)/sin (co) i t f o l l o w s that V(co) c o n t a i n s the branch number of the arctangent of X i (co)/Xr(co). In f a c t , L(co) = V(co) . V(co) may a l s o be used to keep t r a c k of the r e l a t i v e number of branch changes i n the phase between two p o i n t s , say co, and co 2, as co goes from co, to co 2 . Because the r e l a t i o n between X(co) and V(co) i s d i r e c t , c a l c u l a t i n g V(co) does not i n v o l v e sampling X (co) . We compute the Sturm sequence which i s a p u r e l y a l g e b r a i c r e c u r s i o n (Appendix C). Then we ev a l u a t e P(n ,co) f o r a l l n from 0 to M at those co f o r which we want the phase. We then c a l c u l a t e the number of s i g n changes i n the P at each co, y i e l d i n g V(co). We can choose as many or as few co as we l i k e . The method i s not r e c u r s i v e and does not depend on pr e v i o u s c a l c u l a t i o n s as d i d the i n t e g r a t i o n methods. 39 3.5.1 Computational D e t a i l s McGowan and Kuc (1982) present a computer program f o r c a r r y i n g out the e v a l u a t i o n of L ( C J ) using t h i s method. In gene r a l t h i s a l g o r i t h m w i l l be c o m p u t a t i o n a l l y slow when compared to a s i n g l e Fast F o u r i e r Transform (FFT) (Brigham, 1974, p 177). T h i s r e s u l t s from the computation of many t r i g o n o m e t r i c 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 the e v a l u a t i o n of the Sturm sequence. Say we have an input sequence of le n g t h N. Then we c a l c u l a t e a Sturm sequence o f , at most, N polynomials of lengths N, N-1, N-2,..., 1. Thus we evaluate N + (N-1) + (N-2) + ... + 1 or about N 2/2 c o e f f i c i e n t s . Each c o e f f i c i e n t 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 f u n c t i o n . Thus, to e v a l u a t e the sequence of polynomials at one value of CJ we must c a r r y out about N 2/2 t r i g o n o m e t r i c f u n c t i o n e v a l u a t i o n s and m u l t i p l i c a t i o n s . For B v a l u e s of CJ we need B-N2/2 such o p e r a t i o n s . If B = N we need N 3/2 such o p e r a t i o n s . T h i s compares with about N«log 2N such o p e r a t i o n s f o r a s i n g l e (power of 2) FFT. In f a c t , we may use the FFT to evaluate these polynomials and i n c r e a s e e f f i c i e n c y . Consider the polynomial P(k,cj) = a(0)-U(0,cj) + a(l)-U(1,cj) +... where U(n,cj) = s i n [ (n+1 ) C J ) ] / s i n ( C J ) . B r i n g i n g sin(cj) out of the sum y i e l d s s i n ( C J ) • P( k ,C J ) = a(0)«sin[cj] + a ( 1 ) • s i n [ 2CJ] + . . . 40 which can be w r i t t e n as sin(co) 'P(k,co) = 0«sin [0to] + a ( 0 ) - s i n [ w ] +... which i s i n the form of a F o u r i e r transform. Using the FFT w i l l reduce the number of o p e r a t i o n s . For example, i f N i s a power of 2, these o p e r a t i o n s w i l l number about N 2 l o g 2 N . 3.6 D i f f i c u l t i e s Of Phase Unwrapping We w i l l c o n s i d e r only the problem of a c c u r a t e data. That i s , we want the phase of the input sequence x ( n ) , where we know x(n) a c c u r a t e l y . As we have seen, the unwrapped phase can be d e f i n e d a n a l y t i c a l l y . However, f i n i t e numerical e v a l u a t i o n of the a n a l y t i c e x p r e s s i o n s i n t r o d u c e s computational e r r o r s . These may l e a d to an i n c o r r e c t r e s u l t . The q u e s t i o n n a t u r a l l y a r i s e s as to which types of sequence are more s u s c e p t i b l e to these e r r o r s . G e n e r a l l y speaking, when the zeros of X(z) are f a r from the u n i t c i r c l e , the phase w i l l be w e l l determined. I t w i l l vary slowly and any d i s c o n t i n u i t i e s w i l l be e a s i l y d e t e c t e d . As a zero comes c l o s e r to the u n i t c i r c l e , the phase develops a r a p i d change near that zero. If the zero i s on the u n i t c i r c l e the amplitude becomes zero and the phase i s undefined. The d i f f e r e n c e i n the phase curve i f a zero i s j u s t i n s i d e (minimum phase) or j u s t o u t s i d e (maximum phase) the u n i t c i r c l e can be l a r g e (Bhanu, 1977; Clayton and Wiggins, 1976; P o g g i a g l i o m l i et a l . , 1982). Thus, i f computations s h i f t a 41 zero j u s t a c r o s s the u n i t c i r c l e , s e r i o u s i n a c c u r a c i e s c o u l d r e s u l t . For zeros near the u n i t c i r c l e the problem of phase d e t e r m i n a t i o n i s i l l - c o n d i t i o n e d . T h i s i s because the r a t i o X i (co)/Xr (co) used i n the arctangent becomes p o o r l y determined. If both Xi(co) and Xr(co) are s m a l l , small e r r o r s i n the e v a l u a t i o n of Xr(co) can cause l a r g e changes i n the r a t i o . The q u e s t i o n becomes that of which types of sequence have zeros c l o s e to the u n i t c i r c l e . A type of sequence which w i l l be r e l e v a n t i n t h i s t h e s i s c o n s i s t s of a random impulse t r a i n . S t e i g l i t z and D i c k i n s o n (1982) d i s c u s s such a sequence. If the sequence values are independent random v a r i a b l e s then, as the sequence l e n g t h i n c r e a s e s , i t s zeros tend to become evenly d i s t r i b u t e d i n angle and t i g h t l y c l u s t e r e d near the u n i t c i r c l e . A l s o , the p r o b a b i l i t y i n c r e a s e s that there w i l l be a zero so c l o s e to the u n i t c i r c l e that a given a l g o r i t h m w i l l not be able to p l a c e i t on the c o r r e c t side of the c i r c l e . If t h i s sequence i s convolved with a short sequence, a f i n i t e number of e x t r a zeros are i n t r o d u c e d and the above r e s u l t s t i l l h o l d s . The q u e s t i o n of how c l o s e a zero must be to the u n i t c i r c l e to be i n c o r r e c t l y i d e n t i f i e d was e x p l o r e d by Bhanu (1977). For simple short sequences, zeros with a magnitude of 1.00001 were c o r r e c t l y i d e n t i f i e d although magnitudes between about 0.9 and 1.1 are t y p i c a l l y c o n s i d e r e d c l o s e to the u n i t c i r c l e . 42 3.7 Comparison Of Phase Unwrapping Techniques The comparison of phase unwrapping techniques i s d i f f i c u l t . On w e l l c o n d i t i o n e d data they w i l l a l l work w e l l . On p o o r l y c o n d i t i o n e d data, t h e i r accuracy w i l l depend on the i n d i v i d u a l data s e t . The i n t e r a c t i o n between the zeros of even a short sequence i s complicated and not amenable to d i r e c t a n a l y s i s . Thus, the comparison of v a r i o u s techniques must be e m p i r i c a l . Computational time requirements of some techniques has been eva l u a t e d e m p i r i c a l l y or a n a l y t i c a l l y . The computational time of McGowan and Kuc's technique does not depend of the input data, only on i t s l e n g t h . Using a power of 2 FFT to evaluate Sturm sequences, run time i s of the order of N 2 l o g 2 N . That i s T = 0 ( N 2 l o g 2 N ) . S t e i g l i t z and D i c k i n s o n found that t h e i r a l g o r i t h m t y p i c a l l y ran i n T = 0 ( N 2 ) . The alg o r i t h m s of T r i b o l e t , Bhanu and P o g g i a g l i o l m i et a l . i n i t i a l l y run i n T = 0(N«log 2N) before a d a p t a t i o n . Adaptation time depends on the data and the program parameters. McGowan and Kuc's a l g o r i t h m was implemented using a FORTRAN program. T r i b o l e t ' s a l g o r i t h m was p u b l i s h e d as a FORTRAN program ( T r i b o l e t , 1977). On the data used f o r t h i s t h e s i s the programs y i e l d e d i d e n t i c a l r e s u l t s except f o r computation time. T r i b o l e t ' s a l g o r i t h m almost always took l e s s time. For example, f o r an input sequence of 90 p o i n t s padded with zeros to 256 p o i n t s , McGowan and Kuc's a l g o r i t h m took 0.048 ms while T r i b o l e t ' s a l g o r i t h m had an average run time of about 0.02 ms on an Amdahl 470 V/8 computer. 43 3.8 Summary P h a s e u n w r a p p i n g c a n be a p p r o a c h e d f r o m many d i r e c t i o n s . I t e r a t i v e m e t h o d s i n c l u d e s e a r c h a l g o r i t h m , a d a p t i v e i n t e g r a t i o n a n d p o l y n o m i a l f a c t o r i z a t i o n . Number t h e o r y p r o v i d e s a d i r e c t r e l a t i o n b e t w e e n a t i m e s e q u e n c e a nd i t s p h a s e . I n p r i n c i p l e , t h e p h a s e d e t e r m i n a t i o n i s u n s t a b l e f o r z e r o s n e a r t h e u n i t c i r c l e . 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 of the homomorphic transform as a p p l i e d to c o n v o l u t i o n a l problems. The transform changes c o n v o l u t i o n a l problems to a d d i t i v e ones. Techniques developed f o r a d d i t i v e problems may then be a p p l i e d . Using t h i s transform, wavelet e s t i m a t i o n can be s t a t e d as a m u l t i -channel common i n f o r m a t i o n problem. The formalism i s presented i n the next c h a p t e r . T h i s chapter shows how t h i s problem can be approached i n the general case. I t thus p r o v i d e s a s o l u t i o n i n a n t i c i p a t i o n of i t s a p p l i c a t i o n . The reader's forbearance i s requested. Suppose we have a s u i t e of sequences, each sequence a f u n c t i o n of n. Write { x , ( n ) , x 2 ( n ) , . . . } . Suppose f u r t h e r that each of these sequences c o n t a i n s some common s i g n a l i n a d d i t i o n to a component which i s , i n g e n e r a l , d i f f e r e n t i n each sequence. We seek a way to e x t r a c t the common s i g n a l from the s u i t e i n an optimal f a s h i o n . We may measure the amount of common s i g n a l between sequences by t h e i r c o v a r i a n c e . By seeking a l i n e a r combination of the sequences which optimi z e s the common s i g n a l we are l e d to a technique known as p r i n c i p a l component a n a l y s i s (Dhrymes, 1970, p. 53). In f a c t , t h i s technique can be shown to be e q u i v a l e n t to the Karhunen-Loeve, H o t e l l i n g , or ei g e n v e c t o r transforms ( H a l l , 1979, p. 115). Because there e x i s t a v a r i e t y of ways to d e r i v e t h i s 45 r e s u l t , t h a t which lends the most i n s i g h t to the problem d e a l t with i n t h i s t h e s i s w i l l be d e r i v e d . 4.2 P r i n c i p a l Components Let us begin by d e f i n i n g s e v e r a l terms. The terminology i s r e m i n i s c e n t of the s t a t i s t i c a l l i t e r a t u r e , r e f l e c t i n g the o r i g i n s of t h i s development, although the d e f i n i t i o n s w i l l be f o r d e t e r m i n i s t i c data. Say we have a s u i t e of M sequences, each sequence a f u n c t i o n of n and of l e n g t h N. In v e c t o r form we w r i t e {x , (n) ,x, ( n ) , . . .x.,(n)} = x'(n) where underscore denotes a v e c t o r and prime denotes the transpose of a v e c t o r or matrix. We d e f i n e the v a r i a n c e of an element of x ( n ) , say x ^ ( n ) , by var(xj) = 1/N ^[xL(n)-m- ]2 ( 1 ) where m^  i s the mean of x ^ ( n ) , d e f i n e d by A/ m, = 1/N ^ x (n) (2) We d e f i n e the c o v a r i a n c e of two sequences x ^ ( n ) , X j ( n ) as A/ cov(x, , x-) = 1/N £[x,(n)-m:][xj(n)-fj] (3) Equation (3) can be extended to d e f i n e the c o v a r i a n c e matrix C, of x ( n ) , whose elements are given by C.- = C O V ( X ; , X J ) . That 46 /V C = 1/N ^ [ x ( n ) - m ] [ x ( n ) - m ] ' (4) where m i s the vector - of means m' = { m , , m 2 , . . . m M } and [ • 3 [ • ]' r e p r e s e n t s the outer p roduct of v e c t o r s ( G e l b , 1974) . C i s symmetric and can be shown to be p o s i t i v e s e m i d e f i n i t e (Dhrymes, 1970, p 3 ) . For compactness of n o t a t i o n l e t us d e f i n e the o p e r a t o r E { « } to be the sum 1 / N ^ J « } . (Note that t h i s i s not the e x p e c t a t i o n o p e r a t o r i n p r o b a b i l i t y t h e o r y , as we a re not d e a l i n g w i th random v a r i a b l e s . ) Then e q u a t i o n s ( 1 ) , ( 2 ) , and (3) become var(x-L ) = E{[x- (n j -m^] 2 } (5a) m^ = E{x ; (n) } (5b) C = E { [ x ( n ) - m ] [ x ( n ) - m ] ' ] (5c) We can c o n s i d e r the c o v a r i a n c e of two sequences to be a measure of the s i g n a l common to b o t h , where t h a t s i g n a l i s d e f i n e d in terms of f l u c t u a t i o n s about a mean. Let us attempt to maximize the common s i g n a l of the s u i t e of sequences by t a k i n g a l i n e a r c o m b i n a t i o n of the x ( n ) , say y(n) = a ' x ( n ) , where a ' = { a , , a 2 , . . . } i s a v e c t o r of c o n s t a n t c o e f f i c i e n t s . Our measure of the s i g n a l i n y(n) i s i t s v a r i a n c e . Thus we attemt to f i n d that a which maximizes the v a r i a n c e of y ( n ) . If y (n) = a ' x ( n ) t h e n , from e q u a t i o n s (5) 47 v a r ( y ) = a'Ca (6) (Dhrymes, 1970) where C i s the cov a r i a n c e matrix of x ( n ) . However, we see that there i s no maximum of a'Ca. I f a i s a s o l u t i o n , i t c o u l d be m u l t i p l i e d by any constant s c a l a r s such that b = sa and, with a f i x e d , s c o u l d be i n c r e a s e d to make var(y ) a r b i t r a r i l y l a r g e . We have a problem of s c a l e . I t i s the s t r u c t u r e of a which i s important, not i t s magnitude. Thus we f i x the s c a l e by a r b i t r a r i l y r e q u i r i n g a to be of u n i t l e n g t h , that i s a'a = 1 . Now the problem i s to maximize var ( y ) = a'Ca subj e c t to the c o n s t r a i n t a'a = 1. T h i s can be done with the method of Lagrange m u l t i p l i e r s . We d e f i n e the Lagrangian where X i s a Lagrange m u l t i p l i e r . L i s maximized with respect to a and X by s e t t i n g i t s p a r t i a l d e r i v a t i v e s equal to zero. T h i s y i e l d s L = a'Ca + Xd-a'a) (7) dL/da = 0 = 2Ca - 2Xa ( 8 a ) dL/dX = 0 = 1 - a'a ( 8 b ) where we have used the r u l e s of v e c t o r d i f f e r e n t i a t i o n (Gelb, 1974; H a l l , 1970) and 0 denotes the zero v e c t o r . P r e m u l t i p l y i n g equation (8a) by a' y i e l d s 48 a'Ca = X (9) which, from (6), equates the v a r i a n c e of y with the Lagrange m u l t i p l i e r X. Equation (8a) i m p l i e s which i s the eigenvalue problem (Strang, 1980, p. 179-239) with eigenvalues X and e i g e n v e c t o r s a. Because we want the v a r i a n c e of y to be maximal, and t h i s i s equal to X, we must choose the l a r g e s t X which s a t i s f i e s equation (10),-say 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,. Thus, the s o l u t i o n , denoted by y , ( n ) , i s T h i s i s >defined to be the f i r s t p r i n c i p a l component of x(n) (Dhrymes, 1970, p. 34). The mean and v a r i a n c e of y, are e a s i l y d e r i v e d . The mean, my,, i s Ca = Xa (10) y , (n) = a,'x(n) (11) my, = E{y,(n) } my, = E{a,'x(n)} my, = a,'E{x(n)} or f i n a l l y 49 my, = a 1 ' m (12) Thus the mean of the f r i r s t p r i n c i p a l component i s the weighted sum of the means of x ( n ) . The v a r i a n c e of y, i s , from (6) and Thus the v a r i a n c e of the f i r s t p r i n c i p a l component i s equal to the l a r g e s t eigenvalue of equation (10). We expect the f i r s t p r i n c i p a l component to c o n t a i n the most common s i g n a l of the s u i t e . We may then i n v e s t i g a t e other l i n e a r combinations of the sequences to see what p r o p e r t i e s they have. A c c o r d i n g l y , we seek the normalized l i n e a r combination of x ( n ) , say y_(n) = a'x(n), which has maximal v a r i a n c e but i s u n c o r r e l a t e d with the f i r s t p r i n c i p a l component y , ( n ) . By u n c o r r e l a t e d , we mean cov(y.,y) = 0. We can r e w r i t e t h i s c o v a r i a n c e as (usi n g equation (11)) (9), v a r ( y , ) = X (13) cov(y,,y) E{[y,(n)-my, ] [ y (n)-my]'} cov(y,,y) E{[a,'(x(n)-m)][x(n)-m]'a} a 1'E{tx(n)-m][x(n)-m]'}a cov(y,,y) cov(y,,y) a, 'Ca so that the c o n s t r a i n t may be w r i t t e n as a , * Ca = 0 (14) 50 Proceeding as before we d e f i n e the Lagrangian L = a'-Ca + X ( l-a'a) + u a / C a where X and u are Lagrange m u l t i p l i e r s . S e t t i n g p a r t i a l d e r i v a t i v e s to zero y i e l d s dL/da = 0 = 2Ca - 2Xa + uCa, (15a) dL/dX = 0 = 1 - a'a (15b) dL/du = 0 = a ,'Ca (15c) P r e m u l t i p l y i n g equation (15a) by a' and using equations (15b) and (15c) y i e l d s a'Ca = X T h i s i s j u s t the same r e s u l t d e r i v e d e a r l i e r : the v a r i a n c e of y i s equal to the Lagrange m u l t i p l i e r X. From the pr e v i o u s d e r i v a t i o n we have Ca , ' = X, a , P r e m u l t i p l y i n g by a' and using equation (15c) y i e l d s 51 or 0 = a'a, (16) which i s an e x p r e s s i o n r e q u i r i n g the o r t h o g o n a l i t y of the e i g e n v e c t o r s . Now, p r e m u l t i p l y i n g equation (15a) by a' y i e l d s 2a,'Ca - 2Xa,'a + ua,'Ca, = 0 However, from equations (15c) and (16), the f i r s t two terms are zero." Hence ua,'Ca, = 0 and from equation (13) uX, = 0 Now, i n g e n e r a l , X, > 0 which i m p l i e s u = 0 and equation (15a) becomes Ca = Xa and we have returned to the eigenvalue problem. However, now we must choose the second l a r g e s t e i g e n v a l u e , say X 2, and i t s a s s o c i a t e d normalized e i g e n v e c t o r , say a 2 . Thus, the second p r i n c i p a l component i s d e f i n e d as 52 y 2 ( n ) = a 2 ' x(n) and has the p r o p e r t i e s that mean(y 2) = a2'm v a r ( y 2 ) = X 2 c o v ( y . , y 2 ) = 0 T h i s procedure can be continued (Dhrymes, 1970, p 56) by f i n d i n g t h a t l i n e a r combination of x(n) having maximal v a r i a n c e but which i s u n c o r r e l a t e d with the p r e v i o u s p r i n c i p a l components. Doing so g i v e s M l i n e a r combinations of x ( n ) , say y^(n) = a^'x(n) i = 1,2,...M such that var(yj i) = X^ cov(y^,y ; J) = 0 i * j where the X^ are the ordered ( l a r g e s t to s m a l l e s t ) eigenvalues of C and a^ are t h e i r a s s o c i a t e d orthonormal e i g e n v e c t o r s . The y^(n) form a set of sequences which are mutually u n c o r r e l a t e d . l i n e a r combinations of x(n) having maximal 53 v a r i a n c e . In gen e r a l we d e f i n e y_'(n) = {y i (n) ,y 2 ( n ) , . . . ,y (n)} as the ve c t o r of p r i n c i p a l components of x(n) ^with y^(n) as the i - t h p r i n c i p a l component. If we w r i t e the transpose of the e i g e n v e c t o r s as the rows of a matrix A we get y_(n) = Ax(n) (17) A has the prop e r t y that the sum of the v a r i a n c e s of y_(n), the eigenvalues of C, equals the sum of the v a r i a n c e s of x(n) (Ready and Wintz, 1973). Equation (17) d e f i n e s a l i n e a r t r a n s f o r m a t i o n from x(n) to y_(n). I t i s known as the e i g e n v e c t o r , H o t e l l i n g , Karhunen-Loeve, or p r i n c i p a l component tr a n s f o r m . T h i s transform can be a l t e r n a t e l y d e r i v e d as the s o l u t i o n to the optimal data compression problem. T h i s d e r i v a t i o n p r o v i d e s a d d i t i o n a l i n s i g h t i n t o the transform. The e s s e n t i a l ideas w i l l be presented b r i e f l y . For a d e t a i l e d d e r i v a t i o n see H a l l (1979, p. 115), Ahmed and Rao (1975) or Kramer and Mathews (1956). Say we have x(n) and expect a high mutual c o r r e l a t i o n between i t s elements. T h i s c o r r e l a t i o n i s i n d i c a t i v e of redundant i n f o r m a t i o n , or s i g n a l , common to s e v e r a l sequences. In order to reduce storage space or t r a n s m i s s i o n channel u t i l i z a t i o n we wish to represent x(n) using fewer data. We separate x(n) i n t o d i f f e r e n t components, y_(n) = Ax(n), d i s c a r d some of these components, then s t o r e or transmit those 54 remaining. The o r i g i n a l x(n) are then r e c o n s t r u c t e d as w e l l as p o s s i b l e using the remaining components. In general x(n) cannot be r e c o n s t r u c t e d e x a c t l y unless a l l of the components are r e t a i n e d . Let us denote the r e c o n s t r u c t e d x(n) by x 0 ( n ) and d e f i n e the e r r o r i n r e c o n s t r u c t i o n as E(n) = x(n) - x 0 ( n ) . We a l s o d e f i n e the t o t a l e r r o r , T, as the sum of the v a r i a n c e s of E ( n ) . If A i s c o n s t r u c t e d such that the t o t a l e r r o r i s minimized when d i s c a r d i n g s u c c e s s i v e elements of y_(n) the p r i n c i p a l component transform r e s u l t s . I t f o l l o w s t h a t , i f some of the y_(n) are r e p l a c e d by t h e i r mean v a l u e s , T i s equal to the sum of the v a r i a n c e s of these r e p l a c e d y_(n). However, the v a r i a n c e s of the y_(n) are j u s t the eigenvalues of C. Thus the eigenvalues c o n t a i n i n f o r m a t i o n r e g a r d i n g the r e c o n s t r u c t i o n e r r o r . In p a r t i c u l a r , the r a t i o of the f i r s t -e igenvalue to the sum of the r e s t w i l l g i v e the p r o p o r t i o n of the r e c o n s t r u c t i o n which would have been obtained using the f i r s t p r i n c i p a l component to that obtained using the r e s t of the p r i n c i p a l components. Obviously, i f we want the optimal r e c o n s t r u c t i o n u s i n g only one p r i n c i p a l component, t h i s component must be the most common to a l l of the x ( n ) . Thus t h i s r a t i o may a l s o be cons i d e r e d as the r a t i o of the c o r r e l a t e d or common s i g n a l i n x(n) to the u n c o r r e l a t e d s i g n a l or n o i s e . 55 4.3 Some P r o p e r t i e s Of P r i n c i p a l Components Let us examine some l i m i t i n g cases to i l l u s t r a t e the behaviour of t h i s transform as a common s i g n a l e s t i m a t o r . Let us w r i t e the s u i t e of sequences x(n) as a common s i g n a l s(n) p l u s a 'noise' component u(n), i n general d i f f e r e n t f o r each sequence. That i s , x(n) = s(n) + u(n) Let us assume that the mean valu e s of s(n) and u(n) are zero and that the v a r i a n c e s of the u(n) are i d e n t i c a l . A normalized l i n e a r combination of the x(n) which uses no in f o r m a t i o n about the x(n) uses i d e n t i c a l weights. Ready and Wintz (1973) c a l l t h i s the coherent sum of the x(n) and use i t fo r comparison to the f i r s t p r i n c i p a l component. Consider the case where the u(n) are i d e n t i c a l l y c o r r e l a t e d with each other ( C O V ( U J , U J ) = f, f = constant) and u n c o r r e l a t e d with s ( n ) . Then the c o v a r i a n c e matrix C has elements C^ = c o v ( x ^ X j ) = E { [ s ( n ) + u ^ ( n ) ] [ s ( n ) + u y ( n ) ] } Expanding the sum y i e l d s cov(x^,x •) = v a r ( s ) + f 56 or, as v a r ( s ) i s a constant, cov(Xj-,Xj) = constant Thus C i s a constant matrix of u n i t rank having only one non-zero e i g e n v a l u e . The f i r s t e i g e n v e c t o r w i l l have i d e n t i c a l elements and the f i r s t p r i n c i p a l component w i l l equal the coherent sum. A l l of the v a r i a n c e w i l l be c o n c e n t r a t e d i n t o the f i r s t p r i n c i p a l component. In f a c t , i f the no i s e s were a l l zero the same r e s u l t would f o l l o w . T h i s i n d i c a t e s t h at the c o r r e l a t i o n s between n o i s e s makes them i n d i s t i n g u i s h a b l e from s i g n a l . Now c o n s i d e r the case where the no i s e s are u n c o r r e l a t e d with each other and with the s i g n a l . Then = var (s) + dj-y »f where d^ - i s the Kronecker d e l t a and f i s a con s t a n t . Thus C i s a constant matrix p l u s a another constant added to the d i a g o n a l . The f i r s t e i genvalue of C i s X, = X + f where X i s the eigenvalue of C,y = v a r ( s ) . A l l the other eigenvalues are i d e n t i c a l and equal to the constant f (Ready and Wintz, 1973). We can say that the s i g n a l has concentrated i n the f i r s t p r i n c i p a l component and that the noise has spread u n i f o r m l y among a l l the p r i n c i p a l components. The f i r s t p r i n c i p a l component i s equal to the coherent sum. A l l the other e i g e n v e c t o r s are orthogonal to the f i r s t and t h e r e f o r e the sum 57 of t h e i r e lements i s z e r o . T h i s i m p l i e s t h a t there i s no s i g n a l in these p r i n c i p a l components . In the g e n e r a l ca-se n o i s e s - a re c o r r e l a t e d w i th s i g n a l and w i th each o t h e r . We can say tha t the most c o r r e l a t e d s i g n a l i s in the f i r s t p r i n c i p a l component and the u n c o r r e l a t e d n o i s e i s spread u n i f o r m l y over a l l the p r i n c i p a l components. S i g n a l s h a v i n g l e s s e r and l e s s e r c o r r e l a t i o n appear in s u c c e s s i v e p r i n c i p a l components . Because of t h i s c o n c e n t r a t i o n of s i g n a l in the f i r s t p r i n c i p a l component, and the s p r e a d i n g out of the n o i s e over a l l the p r i n c i p a l components , .we can use the e i g e n v a l u e s as a measure of how much s i g n a l was i n the o r i g i n a l d a t a . In p a r t i c u l a r , the r a t i o of the f i r s t to the sum of the r e s t can be used as a measure of the s i g n a l to n o i s e r a t i o . Note t h a t there i s no r e s t r i c t i o n on the s i g n of the e lements of a . Thus , in a p r i n c i p a l component, one or more of the x(n) may be i n v e r t e d . T h i s c o u l d occur i f some of the x(n) had an i n v e r t e d s i g n a l - s ( n ) , or i f the n o i s e l e d to n e g a t i v e c o r r e l a t i o n s . 4 .4 Summary S i g n a l common to s e v e r a l sequences can be e s t i m a t e d u s i n g l i n e a r c o m b i n a t i o n s . I n t e r p r e t i n g v a r i a n c e as an i n f o r m a t i o n measure l e a d s to the c o v a r i a n c e m a t r i x . The e i g e n v e c t o r s of t h i s mat r i x d e f i n e p r i n c i p a l components and the e i g e n v a l u e s show where the i n f o r m a t i o n l i e s . 58 V. WAVELET ESTIMATION 5.1 The Problem ~ In the f i e l d of seismic e x p l o r a t i o n the recorded data are o f t e n modeled as the c o n v o l u t i o n of a wavelet with an impulse sequence. The wavelet i s a time l i m i t e d f u n c t i o n , u s u a l l y having a smooth spectrum (zeros f a r from the u n i t c i r c l e ) ( T r i b o l e t , 1979, p. 9). The impulse sequence c o n s i s t s of i s o l a t e d non-zero v a l u e s . P h y s i c a l l y , the wavelet models a wave propagating i n the E a r t h . The impulse sequence i s a consequence of the Earth's s t r u c t u r e and the p h y s i c a l geometry of the problem. The data are o f t e n a v a i l a b l e i n the form of a s u i t e of time sequences c a l l e d t r a c e s . The wavelet i s assumed constant from t r a c e to t r a c e while the impulses a r e , i n g e n e r a l , d i f f e r e n t . (See Robinson and T r e i t e l (1980) or T e l f o r d et a l . (1976) f o r f u r t h e r information.) T h i s chapter d e a l s with the problem of e s t i m a t i n g the wavelet, given the above data. 5.2 A S o l u t i o n Three methods of wavelet e s t i m a t i o n were d i s c u s s e d by L i n e s and U l r y c h (1977): the Weiner-Levinson, Wold-Kolmogorov and homomorphic methods. The former two methods r e q u i r e a minimum delay wavelet and s t a t i o n a r y impulses while the l a t t e r r e q u i r e s the s e p a r a t i o n of the complex c e p s t r a of the wavelet 59 and the impulse sequence. The method of minimum entropy d e c o n v o l u t i o n (Wiggins, 1978), while not designed as a wavelet e s t i m a t e r , can y i e l d such an estimate (Ooe and U l r y c h , 1979). I t makes an assumption about the entropy, or s i m p l i c i t y , of the impulse sequence. We would l i k e t o a v o i d making these assumptions on an a r b i t r a r y b a s i s . In g e n e r a l , we would l i k e to i n c l u d e as much a p r i o r i i n f o r m a t i o n while e x c l u d i n g as many a r b i t r a r y r e s t r i c t i o n s as p o s s i b l e . For example, i f • we have an a p r i o r i estimate of the wavelet cepstrum's l e n g t h we would l i k e to i n c l u d e t h i s . We assume a random, not n e c e s s a r i l y s t a t i o n a r y , impulse sequence. If the data are a s i n g l e time sequence, t h i s sequence can be segmented. If the data are a l r e a d y i n the form of a s u i t e of sequences t h i s segmentation may not be necessary. Each sequence i s modeled as the c o n v o l u t i o n of a wavelet with an impulse sequence. We w r i t e the data as x k ( n ) = w(n)*i^(n) k=1,...M (1a) or i n v e c t o r n o t a t i o n x(n) = w(n)M (n) (1b) where x ( n ) , w(n) and i^(n) r e p r e s e n t , r e s p e c t i v e l y , the data, the wavelet and the impulse sequences. Each x(n) i s thus formed by the c o n v o l u t i o n of the same w(n) with a d i f f e r e n t 60 i ( n ) . The c o n v o l u t i o n a l o p e r a t i o n i n (1) i s not i n a form most amenable to a n a l y s i s ^ However, by using the homomorphic t r a n s f o r m a t i o n d i s c u s s e d e a r l i e r , the equations can be mapped from a c o n v o l u t i o n a l space to 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) + i_(n) (2b) where the c i r c u m f l e x denotes a complex cepstrum. The data are now i n the form of a s u i t e of sequences, each c o n t a i n i n g a common component, w(n), and a d i f f e r e n t component, i.(n). We wish to estimate w(n) from the. x.(n). T h i s has been done by averaging the x(n) (Clayton and Wiggins, 1976; O t i s and Smith, 1977). Averaging i s a data-independent technique which ignores any i n f o r m a t i o n i n the data. P r i n c i p a l component a n a l y s i s , on the other hand, i s data-dependent and uses i n f o r m a t i o n i n an optimal f a s h i o n ( U l r y c h et a l , 1983). Thus i t may be an improvement over averaging. The homomorphic transform p r o v i d e s us with the p o s s i b i l i t y of s e p a r a t i n g , at l e a s t p a r t i a l l y , w(n) from the i.(n) by low quefrency windowing (U l r y c h , 1971). Before d i s c u s s i n g the p r i n c i p a l component (P.C.) method, l e t us examine averaging. Assume t h a t , f o r each f i x e d n, the 61 elements of i_(n) can be c o n s i d e r e d as random v a r i a b l e s with i d e n t i c a l p r o b a b i l i t y d i s t r i b u t i o n s and zero mean. Then A, l / M S ^ i k ( n ) w i l l t e n d to zero f o r l a r g e M ( O t i s and Smith, K x' 1977). The estimate of w(n), denoted by w f f(n), i s d e f i n e d as M % (n) = 1/M ^ x . ( n ) (4a) M l K Al we (n) = 1/M K£_[u(n) + i)i (n) ] (4b) Af % (n) = w(n) + 1/M 2 i „ ( n ) (4c) The estimated cepstrum equals the wavelet cepstrum i f the sum in equation (4c) i s zero. The i n v e r s e homomorphic transform of w e(n) y i e l d s w e(n), the wavelet e s t i m a t e . Proceeding as above, we d e f i n e the P.C. estimate of w(n) as M Al w e(n) = [ <£a kx(n)]/2a (5a) where the a^ are the f i r s t P.C. weights. The reason f o r n o r m a l i z i n g the estimate can be shown as f o l l o w s . Equation (5a) can be w r i t t e n as M Af M w e(n) = [ ^ .a Kw(n)+ S a k i k ( n ) ] / 2 a „ (5b) or 62 w e(n) = w(n) + ^ T a ^ ( n j / ^ a (5c) which i s analogous to equation ( 4 c ) . Again, i f the sum i s zero, the estimated cepstrum equals the wavelet cepstrum. I t should be noted t h a t , i n the o r i g i n a l c o n v o l u t i o n , the wavelet's s c a l e and time o r i g i n were l o s t . In p r i n c i p l e they cannot be recovered and we w i l l not be concerned with them. If the sum i n equations (4c) or (5c) i s an impulse a t the o r i g i n , w e(n) w i l l d i f f e r from w(n) only by a m u l t i p l i c a t i v e s c a l e f a c t o r (see p r o p e r t y 6 of the complex cepstrum). For the examples used i n t h i s t h e s i s we a r b i t r a r i l y set the zero quefrency p o i n t to zero to normalize the s c a l e f a c t o r and prevent i t from a f f e c t i n g c e p s t r a l c o r r e l a t i o n s . 5.2.1 The Smoothing F u n c t i o n In p r a c t i c e we have a f i n i t e number of data and cannot expect the impulse c e p s t r a to average to zero. Suppose t h a t , i n equation (4c) or (5c), the sum i s not z e r o . Write the sum as g ( n ) . Then (4c) or (5c) becomes w e(n) = w (n) + g (n) with i n v e r s e homomorphic tra n s f o r m we (n) = w(n)*g(n) where g(n) i s the complex cepstrum of g ( n ) . The estimate i s 63 the c o n v o l u t i o n of the true wavelet with g ( n ) . g(n) may be co n s i d e r e d as a smoothing f u n c t i o n which d i s t o r t s the estimate. ( I t i s worth n o t i n g the s i m i l a r i t y of t h i s t o a r e s u l t from the f i e l d of l i n e a r i n v e r s i o n where the estimated model i s the true model viewed through a r e s o l u t i o n f u n c t i o n (Wiggins, 1972; Oldenburg and Samson, 1979). In order to i l l u s t r a t e t h i s smoothing f u n c t i o n c o n s i d e r the f o l l o w i n g example. (One may c o n s i d e r the data as samples of a time f u n c t i o n with a 1 ms i n t e r v a l . Then one would read 'ms' f o r 'p o i n t s ' . ) We generate a s u i t e of s i x data sequences, i_(n), of l e n g t h 70 ms. The data are impulses which are Poisson d i s t r i b u t e d i n spacing and G a u s s i a n , d i s t r i b u t e d i n amplitude. The Poisson parameter i s 11 ms, which i s both the mean and v a r i a n c e of the d i s t r i b u t i o n . F i g u r e 7 shows the data and the corresponding smoothing f u n c t i o n obtained by c e p s t r a l averaging and low quefrency windowing with a Hanning window of l e n g t h 40 ms (centred on the o r i g i n ) . The window len g t h i s a r e s u l t of h i n d s i g h t and w i l l be shown to be a p p r o p r i a t e i n subsequent examples. Note that the smoothing f u n c t i o n approximates an impulse. The complex c e p s t r a and t h e i r average are shown i n F i g u r e 8. Note the e x p o n e n t i a l decay i n d i c a t i n g the lack of c e p s t r a l a l i a s i n g . 64 JL^ _JL 1 (a) (b) ms 60 F i g u r e 7 - Smoothing F u n c t i o n S i x random Impulse sequences (a) and t h e i r smoothing f u n c t i o n due t o c e p s t r a l a v a r a g l n g and windowing ( b ) . 5.2.2 The Wavelet Now l e t us c o n s i d e r the wavelet. We assume a mixed delay wavelet, shown with i t s cepstrum i n F i g u r e 9. The wavelet i s time l i m i t e d to 20 ms and i t s cepstrum i s quefrency l i m i t e d to about 25 ms. T h i s j u s t i f i e s the use of a 40 ms window f o r comparisons using t h i s wavelet. Note the cepstrum's higher amplitude than the p r e v i o u s impulse c e p s t r a . P l o t t i n g i s to the same s c a l e . The power spectrum of the wavelet i s a l s o shown. Note i t s low pass nature. Small changes to the high frequency spectrum due to c o n v o l u t i o n or noise may dominate 65 -128 0 128 ms Figure 8 - Complex Cepstra Due to Impulses These cepstra correspond to the impulse sequences of Figure 7. Note the exponential decay. the l i n e a r phase c o m p u t a t i o n . T h i s can be a v o i d e d by s e t t i n g the phase t o z e r o above some fr e q u e n c y and removing the l i n e a r t r e n d t o t h a t p o i n t . I f the energy c o n t e n t f o r g r e a t e r f r e q u e n c i e s i s s m a l l , changes t o the wavelet and i t s cepstrum w i l l be s m a l l . 66 Figure 9 - An Assumed Wavelet The assumed wavelet (a) is mixed delay. Its cepstrum (b) 1s time lim i t e d and Its power spectrum (c) Is low pass. 5.2.3 P r i n c i p a l Components vs Averaging Let us now compare the P.C. method to av e r a g i n g . A fe a t u r e of the P.C. method i s i t s a b i l i t y to d i s c r i m i n a t e between c e p s t r a . T h i s i s demonstrated i n F i g u r e 10. The input data are the c o n v o l u t i o n of our wavelet with three sequences of random impulses. Random noise has been added to the f i r s t sequence. The noise i s zero mean Gaussian with a va r i a n c e of 25% of that of the noise f r e e sequence. The cepstrum of the noisy input i s n o t i c e a b l y d i f f e r e n t from the ot h e r s . The P.C. estimate y i e l d e d weights of [0.1,0.7,0.7] 67 F i g u r e 10 - D i s c r i m i n a t i o n by P r i n c i p a l Components Input i m p u l s e s e q u e n c e s (a ) and d a t a d e r i v e d by t h e i r c o n v o l u t i o n w i t h a w a v e l e t ( b ) . The top d a t a sequence has a d d i t i v e n o i s e . C e p s t r a ( c ) show the 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 ) a r e shown i n ( d ) . i n d i c a t i n g s u c c e s s f u l d i s c r i m i n a t i o n of the nois y i n p u t . The wavelet estimates due to averaging and the f i r s t P.C. are shown in F i g u r e 1 1 . 6 8 1 PC AVE. INPUT -30 0 30 ms Figure 11 - Wavelet Estimates Wavelet estimates derived from the f i r s t p r i n c i p a l component (1 PC) and averaging (AVE.) of cepstra of Figure 10. The Input wavelet 1s shown at the bottom. T o o b t a i n a n u m e r i c a l c o m p a r i s o n o f t h e e s t i m a t e s w e d e f i n e a m i s f i t e r r o r . T h e m i s f i t i s t h e s u m o f s q u a r e s o f t h e d i f f e r e n c e b e t w e e n t h e w a v e l e t a n d t h e e s t i m a t e . T h e w a v e l e t a n d e s t i m a t e a r e f i r s t s c a l e d t o u n i t v a r i a n c e , t h e n t i m e s h i f t e d f o r m i n i m u m m i s f i t . S i n c e w e a r e c a l c u l a t i n g a 2 5 6 p o i n t c e p s t r u m , t h e w a v e l e t e s t i m a t e i s 2 5 6 p o i n t s l o n g . T h u s t h e n u l l e s t i m a t e h a s a m i s f i t o f 2 5 6 a n d t h e w a v e l e t i t s e l f h a s a z e r o m i s f i t . I n t h i s e x a m p l e t h e a v e r a g e e s t i m a t e h a d a m i s f i t o f 7 8 v e r s u s t h e p r i n c i p a l c o m p o n e n t ' s 4 0 . T h e s e i n d i c a t e a n i m p r o v e m e n t o f t h e P . C . m e t h o d o v e r a v e r a g i n g . 69 5.2.4 M u l t i p l e Sequence Data Let us now apply these methods to a more general example. The data are 15 sequences generated by the c o n v o l u t i o n of our wavelet with random impulse sequences of l e n g t h 70 ms ( F i g u r e 12). The Poisson parameter i s 6 ms implying an a p p r e c i a b l e wavelet o v e r l a p due to the c o n v o l u t i o n . The wavelets estimated from averaging and the f i r s t two P.C.s are shown i n F i g u r e 13. T h e i r r e s p e c t i v e m i s f i t s are 8.1, 5.9 and 370. The f i r s t P.C. estimate has a s m a l l e r m i s f i t than the average and the second P.C. bears l i t t l e resemblance to the wavelet. These r e s u l t s are encouraging and lend credence to our e x p e c t a t i o n s of the P.C. method. On the b a s i s of the p r e v i o u s r e s u l t f o r a n o i s y input we proceed by adding random n o i s e . The noise i s Gaussian with a v a r i a n c e of 5% of that of the d a t a . The n o i s y data are shown in F i g u r e 14 along with the wavelet e s t i m a t e s . M i s f i t s of the average and f i r s t two P.C.s are 20, 280 and 16. Note that the m i s f i t s are g e n e r a l l y l a r g e r than f o r the noise f r e e case, which i s expected. However, the behaviour of the f i r s t two P.C.s i s unexpected and bears c l o s e r examination. I t appears that the proper wavelet estimate i s i n the second P.C. r a t h e r than the f i r s t . T h i s i s not p r e d i c t e d by our t h e o r e t i c a l development and examples. In f a c t , i t seems p a r a d o x i c a l . To r e s o l v e t h i s apparent paradox, l e t us re-examine the homomorphic transform and i t s i n t e r a c t i o n s with the time and quefrency domains. The behaviour of the homomorphic transform i n the 70 — w - ^ n-(a) 0 90 ms —S/W^ VAA— -WyA/vv-^ -— -W\Av~ —J\r-—-^/v-^~->y^v— —'v-v.—J\^ yv-A/^ — (b) ms 90 F i g u r e 12 - M u l t i p l e Sequence Data F i f t e e n random i m p u l s e s e q u e n c e s (a ) and c o n v o l v e d w i t h a w a v e l e t ( b ) . presence of a d d i t i v e noise i s not w e l l understood. However, the phase spectrum may be s t r o n g l y a f f e c t e d (Buttkus, 1975; Clayton and Wiggins, 1976; J i n and Rogers, 1983). I t i s p o s s i b l e that P.C. a n a l y s i s can a i d i n our understanding of how noise a f f e c t s the cepstrum. The appearance of the wavelet cepstrum i n the second P.C. may imply that noise has induced 71 2 PC -30 0 60 ms Figure 13 - Wavelet Estimates Wavelet estimates due to f i r s t two pr i n c i p a l components (1 PC. 2 PC) and average (AVE.) of cepstra of data in Figure 12. The input wavelet is shown at the bottom. a common s i g n a l i n the c e p s t r a . T h i s s i g n a l appears in the f i r s t P.C. and the next most common s i g n a l , due to the wavelet, appears i n the second P.C. T h i s hypothesis can be t e s t e d by examining the time-quefrencey r e l a t i o n s h i p f o r pure noise i n p u t s . In p r a c t i c e we cannot generate completely u n c o r r e l a t e d n o i s e . However, we can compare the c o r r e l a t i o n s of noise in the time domain (input) with those i n the quefrency domain (o u t p u t ) . Poorly c o r r e l a t e d sequences have a covariance matrix with s i m i l a r e i g e n v a l u e s . T h i s i s due to the lack of a common s i g n a l and the d i s t r i b u t i o n of u n c o r r e l a t e d noise over 72 —\/\AA~^ /\/\AIU —sW\j\AJ^\A— —»~\J\~s\l\J\r—a— --V^ AA^ /^ -—^ — —w—y/V/-^-—'y^A*.— (a) ms 90 2 PC 1 PC AVE. INPUT (b) -30 60 ms F i g u r e 14 - N o i s y D a t a The d a t a o f F i g u r e 12 but w i t h 5% a d d i t i v e n o i s e ( a ) . (b) shows t h e w a v e l e t e s t i m a t e s due .to the f i r s t two p r i n c i p a l components (1 PC . 2 PC) and a v e r a g i n g o f c e p s t r a . The input w a v e l e t i s shown a t t h e b o t t o m . a l l the P.C.s. On the other hand, h i g h l y c o r r e l a t e d sequences y i e l d (ordered) eigenvalues which decrease r a p i d l y . T h i s occurs because the sequences c o n t a i n mostly common s i g n a l . Thus, the eigenvalues of inputs and outputs are i n d i c a t i v e of changes i n r e l a t i v e c o r r e l a t i o n s . Consider the s u i t e of s i x input noise sequences shown i n Fig u r e 15. The noise i s zero mean Gaussian. I t s v a r i a n c e i s i r r e l e v a n t as s c a l e f a c t o r s are removed i n the c e p s t r a (Figure 73 256 ms F i g u r e 15 - N o i s e I n p u t s A s u i t e o f random n o i s e I n p u t s . 16). Note the high amplitude of the c e p s t r a and the strong odd component due to the phase. The c e p s t r a are p l o t t e d to the same s c a l e as the pr e v i o u s two examples. For comparison, each set of eigenvalues i s normalized by d i v i d i n g by the l a r g e s t v a l u e . The eigenvalues f o r the output decrease more r a p i d l y than f o r the input ( F i g u r e 17). T h i s i m p l i e s that the outputs ( c e p s t r a ) are more h i g h l y c o r r e l a t e d than the inputs ( t i m e ) . One must be cautioned a g a i n s t i n f e r r i n g that the transform has produced non-random outputs from random i n p u t s . C o r r e l a t i o n may not be a s u i t a b l e measure of s i m i l a r i t y i n the quefrency domain f o r time domain n o i s e . I t i s l i k e l y that t h i s noise induced c o r r e l a t i o n dominates the f i r s t P.C. The obs e r v a t i o n that the wavelet estimate appears i n the second P.C. i m p l i e s t h a t , although the noise may induce a common component, t h i s i s u n c o r r e l a t e d with the wavelet cepstrum. 74 •128 O ms 128 F i g u r e 16 - C e p s t r a of N o i s e The complex c e p s t r a o f the n o i s e Inputs of F i g u r e 15. T h i s f o l l o w s because the P.C.s are mutually u n c o r r e l a t e d . We can gain f u r t h e r i n s i g h t by examining c e p s t r a l covariance m a t r i c e s . For comparison, the wavelet alone y i e l d s 75 a o - INPUT OUTPUT a a 1 I i J-o 3.0 5.0 7.0 EIGENVALUE * F i g u r e 17 - E i g e n v a l u e s o f Input and Output The e i g e n v a l u e s o f i n p u t ( d a r k l i n e ) d e c r e a s e more s l o w l y t h a t t h o s e o f t h e o u t p u t . T h i s i m p l i e s t h a t t h e o u t p u t Is more c o r r e l a t e d t h a n t h e I n p u t . a constant matrix with the value 0.063. The co v a r i a n c e matrix due to input noise alone i s shown i n F i g u r e 18. Note the range of values of both s i g n s . S e v e r a l elements have a 0 . 6 3 3 0 . 3 1 4 0 . 4 0 0 0 . 0 7 9 - 0 . 4 1 1 - O . 1 3 0 0 . 3 1 4 0 . 3 5 9 0 . 2 4 1 - 0 . 0 9 7 - 0 . 3 6 5 - O . 1 4 3 0 . 4 0 O 0 . 2 4 1 0 . 5 7 4 0 . 0 6 8 - 0 . 2 1 9 - 0 . 0 2 9 0 . 0 7 9 - 0 . 0 9 7 0 . 0 6 8 0 . 2 4 4 0 . 0 4 2 0 . 0 6 5 - 0 . 4 1 1 - 0 . 3 6 5 - 0 . 2 1 9 0 . 0 4 2 0 . 4 9 6 0 . 166 - 0 . 1 3 0 - 0 . 1 4 3 - 0 . 0 2 9 0 . 0 6 5 0 . 166 O. 146 F i g u r e 18 - C o v a r i a n c e M a t r i x of N o i s e The c o v a r i a n c e m a t r i x o f t h e n o i s e c e p s t r a of F i g u r e 16. T h e s e e l e m e n t s compare w i t h t h e w a v e l e t c e p s t r u m ' s c o r r e l a t i o n o f 0 . 0 6 3 . magnitude gre a t e r than the co v a r i a n c e of the wavelet cepstrum ( 0 . 0 6 3 ) . Thus co v a r i a n c e s of c e p s t r a due to noise alone can be much l a r g e r than those due to s i g n a l alone. Now l e t us examine the matrices from the previous 76 example. For c l a r i t y , only the r e s u l t from the f i r s t s i x inputs w i l l be c o n s i d e r e d . Consider the ma t r i c e s d e r i v e d from the c e p s t r a of the impulse sequences, C„; from the impulses convolved with the wavelet with no n o i s e , C,; and with 5% n o i s e , C 2 ( F i g u r e 19). The impulse c e p s t r a , C 0, g e n e r a l l y 0 002 0 . 0 0 2 0 . 0 0 2 0 . 1 1 4 - 0 . 0 0 0 - 0 . 0 0 9 - 0 . 0 0 1 - 0 . 0 3 6 - 0 . 0 0 1 - 0 . 0 4 6 - 0 . 0 0 0 - 0 . 0 1 1 -O.OOO - 0 . 0 0 1 - 0 . 0 0 9 - 0 . 0 3 6 0 . 0 1 7 0 . 0 1 4 0 . 0 1 4 0 . 0 3 9 0 . 0 0 7 0 . 0 3 5 0 . 0 0 8 0 . 0 2 1 - 0 . 0 0 1 -O.OOO - 0 . 0 4 6 - 0 . 0 1 1 0 . 0 0 7 0 . 0 0 8 0 . 0 3 5 0 . 0 2 1 0 . 0 9 2 0 . 0 2 2 0 . 0 2 2 ' 0 . 0 2 3 0 . 0 6 4 0 . 0 6 5 0 . 0 6 4 0 . 0 6 4 0 . 0 6 4 0 .066~ 0 . 0 6 5 0 . 0 6 4 0 . 0 6 4 0 . 0 6 4 0 . 0 6 6 O. 127 0 . 0 7 5 0 . 0 5 5 0 . 0 2 3 0 . 0 7 1 0 . 0 7 5 0 . 0 6 9 0 . 0 6 2 0 . 0 5 1 0 . 0 6 5 0 . 0 5 5 0 . 0 6 2 0 . 0 6 8 0 . 0 6 9 0 . 0 6 5 0 . 0 2 3 0 . 0 5 1 0 . 0 6 9 O. 104 0 . 0 6 1 0 . 0 7 1 0 . 0 6 5 0 . 0 6 5 0 . 0 6 1 0 . 0 7 2 (b) 0 . 0 0 9 0 . 0 0 6 0 . 0 0 9 0 . 0 0 8 0 . 0 0 2 0 . 0 0 9 0 . 0 0 6 0 . 0 0 9 0 . 0 0 8 0 . 0 0 2 0 . 0 0 9 0 . 0 1 6 0 . 0 1 6 - 0 . 0 0 1 - 0 . 0 0 5 0 . 0 0 6 0 . 0 1 6 0 . 0 2 8 -O.OOO - 0 . 0 1 7 0 . 0 0 9 - 0 . 0 0 1 - 0 . 0 0 0 0 . 0 1 2 0 . 0 0 3 0 . 0 1 0 - 0 . 0 0 5 - 0 . 0 1 7 0 . 0 0 3 0 . 0 4 4 - 0 . 0 0 8 0 . 0 0 6 0 . 0 0 9 0 . 0 1 0 - 0 . 0 0 8 0 . 0 1 5 (C) F i g u r e 19 - C o v a r i a n c e M a t r i c e s The c o v a r i a n c e m a t r i c e s due t o i m p u l s e c e p s t r a ( a ) , i m p u l s e s c o n v o l v e d w i t h a w a v e l e t (b) and i m p u l s e s c o n v o l v e d 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 ) . y i e l d o f f d i a g o n a l terms l e s s than the wavelet's c e p s t r a l c ovariance (0.063). • The c o v a r i a n c e matrix due to impulses convolved with the wavelet, C,, has l e s s s t r u c t u r e than C 0 and g e n e r a l l y higher c o r r e l a t i o n s , a l l p o s i t i v e . The matrix due to noisy data, C 2 t has g e n e r a l l y smaller c o r r e l a t i o n s than C, 77 and more s t r u c t u r e , i n c l u d i n g negative c o r r r e l a t i o n s . These o b s e r v a t i o n s can be i n t e r p r e t e d as f o l l o w s . C o r r e l a t i o n s due to the impulse sequence c o n t r i b u t e l e s s t o the matrix s t r u c t u r e than those due to the wavelet. Time domain noise can dominate both of these i n the c e p s t r a l c o v a r i a n c e matrix and s t r o n g l y a f f e c t the P.C. e s t i m a t e . Thus, i n a sense, P.C.s may a m p l i f y noise e f f e c t s compared with a v e r a g i n g . The presence of both p o s i t i v e and negative P.C. weights i m p l i e s that the wavelet cepstrum w i l l tend to c a n c e l . These may thus be used as an i n d i c a t o r of noise dominance. 5.2.5 S i n g l e Sequence Data Having achieved an understanding of some of the i n t e r a c t i o n s between the homomorphic tra n s f o r m and P.C. a n a l y s i s , l e t us proceed. Consider the case where the data are a s i n g l e sequence. T h i s can be converted to a multi-sequence input by segmentation. Segmentation can be done on the b a s i s of the envelope ( U l r y c h et a l , 1983). The envelope i s d e f i n e d as the magnitude of the complex s i g n a l ( B r a c e w e l l , 1965, p 271; Cl a y t o n and Wiggins, 1976) and i t s maxima can be used as an i n d i c a t o r of the l o c a t i o n s of the u n d e r l y i n g impulses (Taner and S h e r i f f , 1977). The envelope i s smoothed to reduce noise e f f e c t s and merge c l o s e l y spaced maxima. C e n t e r i n g windows on these maxima y i e l d s the t r a c e segments, which may o v e r l a p . T h i s windowing w i l l , i n g e n e r a l , destroy the o r i g i n a l 78 c o n v o l u t i o n a l model and has a p o o r l y understood e f f e c t on the homomorphic t r a n s f o r m a t i o n ( T r i b o l e t , 1979). However, i f the window changes on a longer s c a l e than the wavelet, we expect i t to a f f e c t only the impulse sequence. Once the data have been converted i n t o the d e s i r e d form we proceed as b e f o r e . Let us c o n s i d e r a s i n g l e sequence example. The data are the c o n v o l u t i o n of a s i n g l e 500 ms impulse sequence convolved with our wavelet ( F i g u r e 20). Envelope maxima y i e l d window l o c a t i o n s . A 90 ms Hanning window was used to e x t r a c t 18 segments ( F i g u r e 21). The wavelet and estimates due to averaging and the f i r s t and second P.C.s are a l s o shown. T h e i r r e s p e c t i v e m i s f i t s are 60, 59, and 409. Note that the average and f i r s t P.C. have s i m i l a r m i s f i t s , the P.C.'s being m a r g i n a l l y b e t t e r . The second P.C. estimate i s d i s s i m i l a r to the wavelet. The corresponding c e p s t r a are shown in F i g u r e 22 with the weights y i e l d e d by the f i r s t P.C. The lack of v a r i a b i l i t y among these weights and the m i s f i t ' s s i m i l a r i t y to the average's i n d i c a t e t h a t the impulse c e p s t r a d i d not s t r o n g l y a f f e c t the c o v a r i a n c e s t r u c t u r e . Let us now examine the e f f e c t of a d d i t i v e noise i n t h i s c o n t e x t . The data with 5%, 10% and 30% noise are shown in F i g u r e 23. (Note that i n d i v i d u a l segments may have d i f f e r e n t s i g n a l to noise r a t i o s . ) The wavelet, average estimate and f i r s t two P.C. estimates are shown in F i g u r e 24. The segments' c e p s t r a and the c o r r e s p o n d i n g f i r s t P.C. weights are shown i n F i g u r e 25. (Note t h a t , f o r 30% n o i s e , 17 segments were defined.) M i s f i t s of the v a r i o u s estimates are 79 80 —-VWv— —Aw\/^ — -^ /WWV ——JiA'yW— — ^ ^ - j y ^ A ^ — ^ — (a) 90 ms 1 PC 2 PC A V E . INPUT (b) •30 60 ms F i g u r e 21 - S e g m e n t a t i o n and E s t i m a t i o n 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 due t o f i r s t two p r i n c i p a l components (1 PC, 2 PC) and a v e r a g i n g ( A V E . ) o f c e p s t r a ( b ) . The i n p u t w a v e l e t 1s shown a t the bot tom o f ( b ) . shown i n Table 1. The second P.C. becomes a b e t t e r estimator than the f i r s t as the noise l e v e l i n c r e a s e s . In f a c t , f o r 30% no i s e , i t i s b e t t e r than the average. Examination 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 j 0.238 A 0.215 — j \ f k - — 0.224 $ 0.250 A 0.242 n OCA •fr Jl U . L D O 1 0.236 - f - 0.240 0.230 — ^ I 0.234 — 0.243 . . i . . 0.227 -30 0 30 ms F i g u r e 22 - C e p s t r a of Segments The c e p s t r a o f t h e segments 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 component (1 PC) and a v e r a g e ( A V E . ) c e p s t r a a r e a l s o shown. The w e i g h t s u s e d t o f o r m the f i r s t P . C . w e i g h t s a r e shown nex 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 F i g u r e 23 - N o i s y S i n g l e Sequence Data A d a t a s e q u e n c e w i t h 5%, 10% and 30% a d d i t i v e n o i s e ( t o p t o b o t t o m ) . 83 -30 0 60 ms F i g u r e 24 - Wave le t E s t i m a t e s The w a v e l e t e s t i m a t e s due t o the n o i s y d a t a of F i g u r e 2 3 . The f i r s t (1 P C ) , s e c o n d (2 PC) and 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 ) and 30% ( c ) n o i s e . The Input w a v e l e t i s shown a t t h e b o t t o m of e a c h p a n e l . 84 J\j\r-(a) (b) 0.21 0.34 0.22 0.27 0.20 0.19 0.28 0.18 0.16 0.30 0.20 0.26 0.22 0.29 0.21 0.21 0.25 0.17 -30 0 30 ms — —A— - \ /^A/ — v ^ / v —w*^-—d\j^\r— 0.35 0.41 - 0 . 1 6 0.11 0.12 - 0 . 1 0 0.35 - 0 . 0 4 0.24 0.10 0.24 0.22 0.02 - 0 . 0 7 0.06 0.44 0.29 0.20 •30 0 30 ms (c) -v\/V-- j u -~ V — 0.32 - 0 . 2 3 0.47 0.28 0.24 0.26 - 0 . 0 8 0.07 0.00 0.17 0.17 0.27 0.28 - 0 .08 0.10 0.36 0.23 -30 O 30 ms F i g u r e 25 - S e g m e n t ' s C e p s t r a The c e p s t r a due t o t h e n o i s y d a t a of F i g u r e 2 4 . Inputs have 5% ( a ) . 10% ( b ) , and 30% ( c ) n o i s e . The f i r s t p r i n c i p a l component w e i g h t s a r e shown b e s i d e the c o r r e s p o n d i n g c e p s t r a . f i r s t P.C. weights and m i s f i t s shows t h a t , as the weight's v a r i a b i l i t y i n c r e a s e s , the wavelet estimate tends to s h i f t from the f i r s t to the second P.C. Note a l s o t h a t , f o r low 85 Wavelet Estimate M i s f i t s Noise Ave. 1 P.C. 2 P.C. 5% 39 19 190 10% 18 74 56 30% 75 134 32 Table 1 - Wavelet M i s f i t s noise ( 5 % ) , the f i r s t P.C. y i e l d e d a b e t t e r wavelet estimate than a v e r a g i n g . A l s o , f o r high noise (30%), the second P.C. y i e l d e d a b e t t e r estimate than a v e r a g i n g . At an intermediate noise l e v e l , the P.C. method y i e l d e d a worse r e s u l t than averaging. I t may seem p a r a d o x i c a l t h a t , f o r the averaged estimate, the m i s f i t i s lower f o r 10% noise than f o r 5% n o i s e . However, there are many f a c t o r s to be c o n s i d e r e d . Segmentation i s a f f e c t e d by n o i s e , and c e p s t r a l windowing may a f f e c t the nois e as w e l l as the wavelet. A l s o , f o r t h i s p a r t i c u l a r wavelet, s e t t i n g the phase to zero above some frequency may a f f e c t the no i s e ' s r e l a t i o n to the wavelet. These do not d e t r a c t from the f a c t t h a t , f o r a f i x e d n o i s e l e v e l , the P.C. estimates may be compared to the average estimate. One must be caut i o n e d a g a i n s t i n f e r r i n g that adding n o i s e to an input w i l l , i n g e n e r a l , improve the esti m a t e . 86 5 . 2 . 6 Low Pass Inputs In p r e v i o u s examples the homomorphic transform was s t a b i l i z e d by s e t t i n g the phase to zero above some angle CJ.. I t may seem obvious t h a t , having admitted ignorance of t h i s phase, we should do l i k e w i s e with the magnitude. T h i s would prevent the s h i f t i n g of energy above u, to zero phase. For c e p s t r a l averaging t h i s i s a l o g i c a l t h i n g to do. However, t h i s w i l l induce a common s i g n a l i n the c e p s t r a which may a d v e r s e l y a f f e c t the P.C. method of a n a l y s i s . Let us examine t h i s f u r t h e r . Say a time sequence has a F o u r i e r transform with magnitude M(CJ) . S e t t i n g M(CJ) to a small value above CJ=CJI, i s e q u i v a l e n t to s e t t i n g log[M ( c j ) ] to a l a r g e negative v a l u e . T h i s may be done i n two s t e p s . F i r s t , m u l t i p l y log[M ( c j ) ] by a boxcar f u n c t i o n , ( i . e . u n i t y on 0<6j<co 1 and zero on CJ,£CJ<.T) . Second, add a bandpass f u n c t i o n which i s zero on 0<CJ<CJ. and a l a r g e negative v a l u e , say A, on CJ,<CJ<7T. T h i s procedure i s shown g r a p h i c a l l y i n F i g u r e 2 6 . The F o u r i e r transform of the boxcar i s a sine f u n c t i o n and that of the bandpass f u n c t i o n i s a cosine modulated sine f u n c t i o n with peak value A, c a l l i t A-S(n) (Bra c e w e l l , 1 9 6 5 ) . In the quefrency domain t h i s i s e q u i v a l e n t to c o n v o l u t i o n with a s i n e f u n c t i o n and a d d i t i o n of A-S(n). A becomes l a r g e r as M(CJ) i s set to s m a l l e r values above CJ,. F i g u r e 27 shows the cepstrum of our wavelet with the l o g magnitude set to v a r i o u s values above CJ, = 0 .62-7T. The values were d e r i v e d by s u b t r a c t i n g 13, 10, 7, 3 and 0 from the o r i g i n a l value at CJ,. Note the dominance of A«S(n) as A 87 Figure 26 - Low Pass Signals 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 large value. This can be done by m u l t i p l y i n g the log magnitude by a boxcar funct ion , then adding a bandpass func t ion . becomes l a r g e . (The zero quefrency p o i n t has been set to zero.) T h i s e f f e c t and the a r b i t r a r i n e s s of A suggest that t h i s i s not an a p p r o p r i a t e procedure when P.C. a n a l y s i s i s to be used on the c e p s t r a . 5.3 An A l t e r n a t e S o l u t i o n We have seen t h a t , i n the presence of a d d i t i v e n o i s e , the use of P.C.s on c e p s t r a i s i n a p p r o p r i a t e . However, the p r i n c i p a l component method may s t i l l be u t i l i z e d . The homomorphic transform a l l o w s f o r wavelet e s t i m a t i o n by low quefrency windowing. E x p o n e n t i a l weighting of inputs and de-weighting of outputs may improve the estimate (Buhl, et a l , 1974). A wavelet estimate can be obtained from each input sequence by t h i s method. T h i s y i e l d s a s u i t e of estimates from which the most common wavelet can be e x t r a c t e d by p r i n c i p a l components. Th i s i s i l l u s t r a t e d i n F i g u r e 28 using the noisy m u l t i p l e sequence data generated p r e v i o u s l y . An e x p o n e n t i a l weighting 88 ^ • i i i 1 1 > -128 0 128 ms F i g u r e 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 the l o g m a g n i t u d e c o n s t a n t above t h e 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 the o r i g i n a l v a l u e minus 13, 10, 7 , 3 , and 0 ( t o p to b o t t o m ) . f a c t o r of 0.98 and a 5 ms c e p s t r a l window were used to produce wavelet e s t i m a t e s . These estimates were normalized to u n i t v a r i a n c e and s h i f t e d to a l i g n t h e i r envelope peaks. The most common estimate may be d e f i n e d as the average or the f i r s t P.C. T h e i r m i s f i t s are 163 and 20 r e s p e c t i v e l y which shows a 89 F i g u r e 28 - T ime Domain P r i n c i p a l Components W a v e l e t s e s t i m a t e d f rom t h e d a t a o.f F i g u r e 14 u s i n g c e p s t r a l w indowing ( 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 ) . Input w a v e l e t Is shown a t t h e b o t t o m . s i g n i f i c a n t improvement of the P.C. method over averaging. These m i s f i t s compare with 20, 280 and 16 f o r the c e p s t r a l average and f i r s t two c e p s t r a l P.C.s. Thus t h i s technique can y i e l d an estimate which i s comparable to c e p s t r a l a v e r a g i n g . 90 5.4 P r a c t i c a l C o n s i d e r a t i o n s The p r e v i o u s l y d i s c u s s e d wavelet e s t i m a t i o n scheme has s e v e r a l parameters which must be s p e c i f i e d . These i n c l u d e the segmentation window shape and l e n g t h , c e p s t r a l window shape and length and e x p o n e n t i a l weighting f a c t o r . S p e c i f i c a t i o n of these must be done c a r e f u l l y , on the b a s i s of a p r i o r i i n f o r m a t i o n and experience. The segmentation window should be slowly v a r y i n g . I t s l e n g t h i s a t r a d e - o f f between wavelet d i s t o r t i o n and sequence l e n g t h . S h o r t e r sequences have g r e a t e r d i s t o r t i o n but tend to y i e l d more r e l i a b l e phase unwrapping. The c e p s t r a l window should probably not have sharp d i s c o n t i n u i t i e s . I t s l e n g t h can be a d j u s t e d by monitoring the corresponding e s t i m a t e s . E x p o n e n t i a l weighting i s again a t r a d e - o f f between d i s t o r t i n g the wavelet's cepstrum and s e p a r a t i n g i t from the impulse c e p s t r a . Due to the s t a t i s t i c a l nature of the scheme, more data w i l l tend to improve the estimate. T h i s technique i s not user independent. I t s f l e x i b i l i t y and power can be f u l l y r e a l i z e d only with c a r e f u l and j u d i c i o u s use. 5.5 Summary The a p p l i c a t i o n of the p r i n c i p a l component method for wavelet e s t i m a t i o n has been demonstrated. The method may be used in the quefrency domain or i n the time domain a f t e r c e p s t r a l windowing. A d d i t i v e noise may induce a common s i g n a l 91 i n the c e p s t r a . S e t t i n g out of band energy to a small value w i l l add a common term to the c e p s t r a . In the presence of these common s i g n a l s the P.C. method i s more a p p r o p r i a t e l y a p p l i e d i n the time domain. 92 VI. INVERTIBILITY OF THE HOMOMORPHIC TRANSFORM A r e c e n t l y p u b l i s h e d note suggests t h a t , due to phase unwrapping, the homomorphic transform may not be i n v e r t i b l e ( J i n and Rogers, 1983). T h i s note i s supported by examples i l l u s t r a t i n g t h e i r f a i l u r e to o b t a i n a s u c c e s s f u l i n v e r s i o n . T h i s f a i l u r e should not occur. To set our minds at ease these same examples are presented here demonstrating t h e i r s u c c e s s f u l forward and in v e r s e t r a n s f o r m a t i o n . Except f o r inputs with zeros on the u n i t c i r c l e , i f the forward transform can be c a l c u l a t e d , so can i t s i n v e r s e , to w i t h i n the p r e c i s i o n of the computation. If the phase i s unwrapped i n c o r r e c t l y an i n c o r r e c t cepstrum w i l l be c a l c u l a t e d . However, even i f t h i s occurs, the i n v e r s e transform w i l l s t i l l y i e l d the o r i g i n a l i n p u t . T h i s f o l l o w s from the f a c t that phase unwrapping adds i n t e g e r m u l t i p l e s of 2ir to the p r i n c i p a l phase (except the method of McGowan and Kuc (1983) which adds m u l t i p l e s of it) . These are tran s p a r e n t to the in v e r s e o p e r a t i o n of e x p o n e n t i a t i o n . (Note that the whole i n v e r s e system i s r e a l i z e d by t h i s e x p o n e n t i a t i o n and F o u r i e r transforms, which are i n v e r t i b l e . ) The f i r s t example c o n s i s t s of two impulses, of magnitude 2000 and 1999, separated by 20 p o i n t s (Figure 29). The sequence i s extended to 1024 p o i n t s with z e r o s . I t s cepstrum and the cepst.rum's inverse are shown. Note the s u c c e s s f u l i n v e r s i o n . ( C a l c u l a t i o n s were performed with s i n g l e p r e c i s i o n a r i t h m e t i c i n a FORTRAN program.) The cepstrum i s minimum delay and e x h i b i t s some a l i a s i n g . 93 (a) 1024 (b) -512 512 (C) 1024 (d) 1024 (e) r -512 512 (f) ms 1024 F i g u r e 29 - I n v e r s i o n of C e p s t r a An Input (a ) c o n s i s t s of two i m p u l s e s o f a m p l i t u d e 2000 and 1999. I t s c e s t r u m (b) and t h e c e p s t r u m ' s I n v e r s e ( c ) show the r e t u r n o f t h e o r i g i n a l I n p u t . (d ) Is the Input sequence of (a ) w i t h a d d i t i v e n o i s e . (e ) i s t h e cepstum o f (d) and ( f ) Is t h e i n v e r s e o f ( e ) . 94 The e f f e c t of a d d i t i v e n o i s e with a standard d e v i a t i o n of 5 i s a l s o shown. T h i s n o i s e l e v e l i s too low to be p e r c e p t i b l e i n the p l o t s . I t has, however, changed the r e l a t i v e magnitude of the impulses such that the sequence i s n e a r l y minimum d e l a y . Again the cepstrum i n v e r t e d c o r r e c t l y . In the second example, the f i r s t impulse has a magnitude of 2000 while the second i s reduced from t h i s by 33% (Figure 30). The input, cepstrum and the cepstrum's i n v e r s e are shown. When noise having a standard d e v i a t i o n of 20 i s added, the cepstrum's appearance changes d r a m a t i c a l l y . However, i t s t i l l i n v e r t s to the o r i g i n a l i n p u t . The r e s u l t s of J i n and Rogers (1983) are reproduced i n F i g u r e s 31 and 32 f o r comparison. A l l of the inputs d e a l t with i n the p r o d u c t i o n of t h i s t h e s i s were i n v e r t i b l e . We can be s a t i s f i e d that the r e s u l t s of J i n and Rogers are anomalous. 95 (a) 1024 (b) -512 512 (c) 1024 (d) 1024 (e) -512 512 (f) ms 1024 F i g u r e 30 - I n v e r s i o n o f C e p s t r a An 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 w i t h d i f f e r e n t a m p l i t u d e . I t s c e p s t r u m (b) and the 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 o f the i n p u t . (d) i s the input (a) w i t h a d d i t i v e n o i s e . (e ) Is t h e c e p s t r u m of (d ) and ( f ) Is the i n v e r s e 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). F i g u r e 31 - I n v e r s i o n of C e p s t r a The r e s u l t s of J i n and R o g e r s ( 1 9 8 3 ) . T h i s compares w i t h F i g u r e 2 9 . 97 (a) (b) -0.512 (c) (d) 1.024 0.512 1.024 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). F i g u r e 32 - I n v e r s i o n o f C e p s t r a The r e s u l t s o f J i n and Rogers ( 1 9 8 3 ) . T h i s compares w i t h F i g u r e 3 0 . 98 V I I . DISCUSSION AND CONCLUSIONS The wavelet e s t i m a t i o n problem has been formulated i n terms of m u l t i - c h a n n e l common s i g n a l e s t i m a t i o n . T h i s does not r e q u i r e r e s t r i c t i v e assumptions about the wavelet's phase. The f o r m u l a t i o n i n c l u d e s a homomorphic transform to map from a c o n v o l u t i o n a l to an a d d i t i v e ( c e p s t r a l ) space. R e a l i z a t i o n of the transform r e q u i r e s the computation of a F o u r i e r transform's continuous phase. In the a d d i t i v e space the problem i s one of common s i g n a l e s t i m a t i o n . P r o p e r l y formulated, the method of p r i n c i p a l components y i e l d s an optimal estimate which can be compared to averag i n g . The s u c c e s s f u l a p p l i c a t i o n of t h i s technique has been demonstrated f o r noise f r e e d ata. In the n o i s y data case p r i n c i p a l components may y i e l d a poorer estimate than averaging. T h i s occurs because noise may induce a common s i g n a l i n the c e p s t r a l space. P r i n c i p a l components then estimate t h i s n o i s e - i n d u c e d s i g n a l i n s t e a d of the wavelet cepstrum. I f an attempt i s made to reduce n o i s e e f f e c t s by s e t t i n g out of band energy to a small value, another common s i g n a l i s added to the c e p s t r a . Thus, f o r noisy data, we are l e d to an a l t e r n a t e s o l u t i o n . A wavelet i s estimated from each input channel by exp o n e n t i a l weighting and c e p s t r a l windowing. P r i n c i p a l components are then used to e x t r a c t the most common estimate from t h i s s u i t e of e s t i m a t e s . 99 APPENDIX A - THE Z-TRANSFORM The z-transform of a sequence x(n) i s d e f i n e d as X(z) = ^ x ( n ) . z - h (1) where z i s a complex number. The i n v e r s e z-transform i s d e f i n e d as x(n) = 1 / 2 T T J J ^ " X ( Z ) - z " " 2 dz (2) C where C i s a co u n t e r c l o c k w i s e c l o s e d contour i n the region of convergence of X ( z ) , e n c i r c l i n g the o r i g i n i n the z-plane. Together, (1) and (2) form a z-transform p a i r . I f , i n (1) and (2), C i s taken as a c i r c l e c e n t r e d on the o r i g i n , z = r • exp[ ja>], on -7r<co<jr, then (1) and (2) become X(r»exp[ jo>]) = 2[x(n>*r 3*exp[-jco] (3) -tr x{n)'r~h = ]/2nJ~X{r-exp[ jco]) «exp[ jam]da> (4) -TT These imply that X ( r «exp[ jco]) i s the d i s c r e t e time F o u r i e r transform of x(n)«r" h. Thus the F o u r i e r transform can be used to e v a l u a t e the z-transform on a c i r c l e c e n t r e d on the o r i g i n . If r=1, equations (3) and (4) become a F o u r i e r Transform p a i r . S e t t i n g rrM i s r e f e r r e d to as e x p o n e n t i a l weighting (Oppenheim and S c h a f e r , 1975). 100 APPENDIX B ~ FOURIER TRANSFORMS AS CHEBYSHEV POLYNOMIALS The D i s c r e t e Time F o u r i e r Transform (DTFT) can be w r i t t e n i n terms of Chebyshev polynomials of the second kind (McGowan and Kuc, 1982). Consider a sequence x(n) where x(n)=0 f o r n<0 and n>N. The DTFT of x(n) i s F ( G J ) = 2 x< n) *exp[- jnu] (1) or, i n terms of r e a l and imaginary p a r t s F ( C J ) = 2 x(n) «cos(n6j) - j - x (n) 'S infnco) (2) Chebyshev polynomials of the f i r s t kind, as a f u n c t i o n of cos(w), are d e f i n e d as (Snyder, 1966) T(n,cj) = cos(nu) (3) and those of the second kind as U(n,cj) = s i n [ (n+1 )cj]/sin (CJ ) (4) There i s a r e c u r s i o n r e l a t i o n betwen (3) and ( 4 ) : T(n,w) = [U(n,oj) - U(n-2,w)]/2 (5) From equations (2) and (3) we wr i t e 101 N A/ F ( G J ) = 2 x < n ) •T(n,u>) - j« ]^x(n) *sin(na>) (6a) D e f i n i n g R = Real{F} amd I = Imaginary{F}, we w r i t e (6a) as F = R + j • I (6b) Consider the r e a l part of ( 6 ) : R = 2 x(n)-T(n.w) (7) Dropping the argument CJ and using equation ( 5 ) , (7) becomes R = ( ^ x ( n ) . U ( n ) - ^fx(n)-U(n-2)} or, using the d e f i n i t i t i o n i n ( 4 ) , R = { ^ x ( n ) - s i n [ (n+1 ) C J] - ^? x ( n ) - s i n [ (n-1 ) O J ] }/2sin(o) Rearranging terms and n o t i n g that -x(0) • s i n ( - C J ) = x(0)«sin(cj) and x(1)«sin(0) = 0, y i e l d s R = {[ 2 . x ( 0 ) - x ( 2 ) ] - s i n ( c j ) + ^ [ x U - l ) - x(n+l) ] | sin(ncj) + [x(N-l).] • s i n ( NCJ ) + [ x ( N ) , s i n [ ( N + l ) c j ] } / 2 s i n ( c j ) 102 or, b r i n g i n g the denominator i n s i d e the sum and using (4) y i e l d s where R = ^ a(n) -U(n,co) «=<» (8) a ( 0 ) = x ( 0 ) - x ( 2 ) / 2 a(n) = [x(n>-x(n+2 ) ] /2 a ( N - D = [ X ( N - 1 ) ] / 2 a ( N ) = x(N ) / 2 n=1,...N-2 Now c o n s i d e r the imaginary p a r t of ( 6 ) : I = - ^ x ( n ) »sin (nco) (9a) or, as the f i r s t term i s zero, A/-J I = - 2 x ( n + l ) s i n [ (n+1 )o>] y> zo S e t t i n g (9b) b(n) = -x(n+1) (10) i n (8). y i e l d s A/-1 1 = 2 b< n> *s i n [ (n+1 )CJ] 1 03 I = sin(oj)* 2 { b ( n ) 'Sin[ (n + 1 )6j]/sin(cj)} I = sin(oj) • S b(n) .U(n,a;) where we have used the d e f i n i t i o n ( 4 ) . T h e r e f o r e we w r i t e (1) i n terms of Chebyshev polynomials of the second kind as A/ M-7 F (C J ) = 2 a ( n ) - U ( n , c j ) + j-sin(cj) • 5 J b(n)-U(n,6j) where the a(n) and b(n) are given by equations (8) and (10) resp e c t i v e l y . 1 04 APPENDIX C - GENERATION OF A STURM SEQUENCE T h i s appendix i l l u s t r a t e s a method of gen e r a t i n g a Sturm sequence from two Chebyshev polynomials of the second kind (McGowan and Kuc, 1982). The n o t a t i o n i s designed to allow easy t r a n s c r i p t i o n to a computer program. Chebyshev polynomials of the second kin d , U(n), are d e f i n e d i n appendix B and have the r e c u r s i o n r e l a t i o n (Snyder, 1966) U ( n ) - U d ) = U(n+1) + U(n-1) (1) We s t a r t with the two polynomials P(0) = ^ a ( 0 , n ) - U ( n ) (2) P(1) = 2 a { 1 .n).U(n) (3) where a(k,n) r e f e r s to the n-th c o e f f i c i e n t of the k-th polynomial P ( k ) . The Sturm sequence i s a sequence of polynomials of d e c r e a s i n g degree, generated r e c u r s i v e l y from P(0) and P(1) by E u c l i d ' s a l g o r i t h m (Marden, 1966) P(k+1) = Q(k)-P(k) - P(k-1) (4) The ge n e r a l polynomial can be w r i t t e n as P(k) = ^Ta(k,n)-U(n) (5) 1 0 5 and the f i n a l polynomial as P(N) = a(N,0)-U(0) (6) In the r e c u r s i o n (4), l e t us d e f i n e the Q(k) as the f i r s t degree polynomial Q(k) = d(k,0)-U(0) + d ( k , l ) . U ( 1 ) (7) where the d are c o n s t a n t s . For c l a r i t y l e t us begin with an example by generating P(2). I n s e r t i n g (2), (3) and (7) i n t o (4) y i e l d s P(2) = { 2 a ( l , n ) . U ( n ) } . { d ( l , 0 ) . U ( 0 ) + d ( l , l ) . U ( l ) ] } -{ S a(0,n).U(n)} or, expanding, P(2) = { f a ( l , n ) - d ( l ,0)-U(n)}.U(0) + { S a( 1 , n ) - d d , 1)-U(n)}-U( 1) - { ^ a ( O . n ) .U(n)} The r e l a t i o n (1) and the f a c t t h a t U(0)=1 y i e l d 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 fN-1)-d(1,1)}-U(N) - S a(0,n)-U(n) h-o C o l l e c t i n g l i k e terms y i e l d s P(2) = {a(1 , 0 )-d(1 , 0 ) + 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 (0,n)}-U(n) + {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 choose Q by f i n d i n g d(1,0) and d ( l , l ) such that the terms i n U(N) and (N-1) c a n c e l . That i s , d( 1 , 1) = a(0,N)/a(1,N-1) d(1,0) = { a ( 0 , N - O - a (1,N-2)-d(1,1 ) } / a (1,N-1) By i n d u c t i o n we can w r i t e the general term f o r the k-th polynomial as 1 07 P(k) = { a ( k - 1 , 0 ) - d ( k - 1 , 0 ) + a ( k - 1 , 1 ) - d ( k - 1 , 1 ) - a ( k - 2 , 0 ) } + S U(k - 1,n).d(k - 1,0)+[a(k - 1,n - 1)+a(k - 1,n+ 1 ) ] - a(k-2,n)}-U(n) and the f i n a l polynomial as P(N) = {a(N -1,0)-d(N -1,0) + a ( N - 1 , 1 ) - d ( N - 1 , 1 ) - a(N-2,0)}.U(0) The k-th polynomial can be w r i t t e n as P(k) = £ a(k,n)-U(n) k=2,...N where a(k,0) = a ( k - 1 , 0 ) - d ( k - 1 , 0 ) + a ( k - 1 , 1 ) - d ( k - 1 , 1 ) - a ( k - 2 , 0 ) (8a) and, 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 r e c u r s i o n r e l a t i o n s (8) and (9) allow the computation of a l l the P ( k ) , k>0, pr o v i d e d P(0) and P(1) are not r e l a t i v e l y prime. I f t h i s i s the case, the h i g h e s t degree c o e f f i c i e n t of a P(k) i s zero and a d i v i s i o n by zero i n (9) w i l l o ccur. For t h i s case the next polynomial generated has a degree reduced by more than one. The r e c u r s i o n then s k i p s one i t e r a t i o n and contin u e s f o r subsequent c a l c u l a t i o n s . 109 BIBLIOGRAPHY 1. Ahmed, N. and Rao, K.R., 1975, Orthogonal Transforms 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 York. 2. Beaumont, R.A. and P i e r c e , R.S., 1963, The A l g e b r a i c Foundations of Mathematics: Addison-Wesley, U.S.A. 3. Bhanu, B., 1977, Computation of Complex Cepstrum: M.Sc. and E.E. t h e s i s , M.I.T., Massachusetts. 4. Bogert, B.P., Healy, M.J.R. and Tukey, J.W., 1962, The Quefrency A l a n y s i s of Time S e r i e s f o r Echoes: Cepstrum, Pseudo-Autocovariance, Cross-Cepstrum and Saphe C r a c k i n g : Proc. Symp. on Time S e r i e s A n a l y s i s , Ed. M. Ros e n b l a t t , John Wiley and sons, Inc., N.Y. 5. B r a c e w e l l , R., 1965, The F o u r i e r Transform and I t s A p p l i c a t i o n s : McGraw-Hill, U.S.A. 6. Brigham, E.O., 1974, The F a s t F o u r i e r Transform: P r e n t i c e - H a l l , Inc., New J e r s e y . 7. Buhl, P., S t o f f a , P.L. and Bryan, G.M., 1974, The A p p l i c a t i o n of Homomorphic Deconvolution t o Shallow-Water Marine Seismology-Part I I : Real Data: Geophysics, v. 39, p. 417-426. 8. Buttkus, B., 1975, Homomorphic F i l t e r i n g - Theory and P r a c t i c e : G eophysical 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 Wiggins, R.A., 1976, Source Shape E s t i m a t i o n and Deconvolution of T e l e s e i s m i c Bodywaves: Geophys. J.R. a s t r . S o c , v. 47, p. 151-177. 10. Dhrymes, P.J., 1970, Econometrics, S t a t i s t i c a l Foundations 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. Gelb, A. et a l . , 1974, A p p l i e d Optimal E s t i m a t i o n : The A n a l y t i c Sciences C o r p o r a t i o n , M.I.T., 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 Press, Inc., New York. 13. J i n , D.J. and Rogers, J.R., 1983, Homomorphic Deconvolution: Geophysics, 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 For T r a n s m i t t i n g a Set of C o r r e l a t e d S i g n a l s : I.R.E. Trans, on Information Theory, v. IT-2, p. 41-46. 15. L i n e s , L.R. and U l r y c h , T .J., 1977, The Old and the New in Seismic Deconvolution and Wavelet E s t i m a t i o n : 110 G e o p h y s i c a l P r o s p e c t i n g , 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 : Schaum's O u t l i n e S e r i e s , McGraw-Hill, London. 17. Marden, M., 1966, Geometry of Polynomials: American Mathematical S o c i e t y , Rhode I s l a n d . 18. McGowan, R. and 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 and I t s Unwrapped Phase: I.E.E.E. Trans, on A.S.S.P., v. ASSP-30, p. 719-726. 19. Oldenburg, D.W. and Samson, J . C , 1979, I n v e r s i o n of I n t e r f e r o m e t r i c Data 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 Plasmas: J . Opt. Soc. Am., v. 69, p. 927-942. 20. Ooe, M. and U l r y c h , T.J., 1979, Minimum Entropy Deconvolution with an E x p o n e n t i a l T r a n s f o r m a t i o n : 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 of Nonlinear Systems: Tech. Rept. 432, Research Laboratory of E l e c t r o n i c s , M.I.T., Massachusetts. 22. Oppenheim, A.V. and Schafer, R.W., 1975, D i g i t a l 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 , Inc. New J e r s e y . 23. Oppenheim, A.V., Schafer, R.W. and Stockham, T.G., 1968, Nonlinear F i l t e r i n g of M u l t i p l i e d and Convolved S i g n a l s : Proc I.E.E.E., v. 56, p. 1264-1291. 24. O t i s , R.M. and Smith, R.B., 1977, Homomorphic Deconvolution by Log S p e c t r a l Averaging: Geophysics, v. 42, p. 1146-1157. 25. P o g g i a g l i o l m i , E., Berkhout, A.J. and Boone, M.M., 1982, Phase Unwrapping, P o s s i b i l i t i e s and L i m i t a t i o n s : Geophysical P r o s p e c t i n g , v. 30, p. 281-291. 26. Ready, P.J. and Wintz, P.A., 1973, Information E x t r a c t i o n , S.N.R. Improvement, and Data Compression i n M u l t i s p e c t r a l Imagery: I.E.E.E. Trans, on Communications, v. COM-21, p 1123-1130. 27. Robinson, E.A. and T r e i t e l , S., 1980, Geo p h y s i c a l S i g n a l A n a l y s i s : P r e n t i c e - H a l l , Inc., New J e r s e y . 28. Schafer, 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 : Tech. Rept. 466, MIT Research Laboratory of E l e c t r o n i c s , MIT, Massachusetts. 29. Snyder, M.A., 1966, Chebyshev Methods i n Numerical Approximation: P r e n t i c e - H a l l , Inc., New J e r s e y . 111 30. S t e i g l i t z , K. and D i c k i n s o n , B., 1982, Phase Unwrapping by F a c t o r i z a t i o n : I.E.E.E. Trans, on A.S.S.P., v ASSP-30, p. 984-991. 31. S t o f f a , P.L., Buhl.,P. and Bryan, G.M., 1974, The A p p i c a t i o n s of Homomrphic Deconvolution to Shallow-Water Marine Seismology-Part I:Models: Geophysics, v. 39, p. 401-416. 32. Strang, G., 1980, L i n e a r Algebra and I t s A p p l i c a t i o n s (2 ed.): Academic P r e s s . New York. 33. Taner, M.T. and S h e r i f f , R.E., 1977, A p p l i c a t i o n of Amplitude, Frequency, and Other A t t r i b u t e s to S t r a t i g r a p h i c and Hydrocarbon Determination: i n Seismic 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 , Ed. C.E. Payton, Memoir 26 of A.A.P.G, T u l s a , p. 301-328. 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 Keys, D.A., 1976, A p p l i e d Geophysics: Cambridge U n i v e r s i t y Press, New York. 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. Trans, on A.S.S.P., v. ASSP-25, p. 170-177. 36. T r i b o l e t , J.M., 1979, Seismic A p p l i c a t i o n s of 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 , Inc, New J e r s e y . 37. U l r y c h , T.J., 1971, A p p l i c a t i o n of Homomorphic Deconvolution to Seismology: Geophysics, v. 36, p. 650-660. 38. U l r y c h , T.J., Levy, S., Oldenburg, D.W., and Jones, I.F., 1983, A p p l i c a t i o n s of the Karhunen-Loeve Transformation i n R e f l e c t i o n Seismology: presented at the 53rd Annual Meeting of the S.E.G., Las Vegas. 39. Wiggins, R.A., 1972, The General L i n e a r Inverse Problem: I m p l i c a t i o n s of Surface Waves and Free 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 of Geophysics and Space P h y s i c s , v. 10, p. 251-285. 40. Wiggins, R.A., 1978, Minimum Entropy Deconvolution: Ge o e x p l o r a t i o n , v. 16, p. 21-35. 

Cite

Citation Scheme:

        

Citations by CSL (citeproc-js)

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>
                        
                    
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:
https://iiif.library.ubc.ca/presentation/dsp.831.1-0052934/manifest

Comment

Related Items