Open Collections

UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Maximum entropy spectral analysis of free oscillations of the earth : the 1964 Alaska event 1976

You don't seem to have a PDF reader installed, try download the pdf

Item Metadata

Download

Media
UBC_1976_A6_7 D38.pdf [ 4.37MB ]
UBC_1976_A6_7 D38.pdf
Metadata
JSON: 1.0052997.json
JSON-LD: 1.0052997+ld.json
RDF/XML (Pretty): 1.0052997.xml
RDF/JSON: 1.0052997+rdf.json
Turtle: 1.0052997+rdf-turtle.txt
N-Triples: 1.0052997+rdf-ntriples.txt
Citation
1.0052997.ris

Full Text

MAXIMUM ENTROPY SPECTRAL ANALYSIS OF FREE OSCILLATIONS OF THE EARTH: THE 1964 ALASKA EVENT by John C. Davies fi.Sc, University of B r i t i s h Columbia, 1973 A THESIS SUBMITTED IN PARTIAL FULFILMENT OF THE REQUIREMENTS FOR THE DEGREE OF MASTER OF SCIENCE in the Department of Geophysics and Astronomy We accept t h i s thesis as conforming to the reguired standard The University Of B r i t i s h Columbia March, 1976 (c) John C. Davies In presenting- th i s thes is in pa r t i a l fu l f i lment of the requirements for an advanced degree at the Un ivers i ty of B r i t i s h Columbia, I agree that the L ibrary sha l l make it f ree ly ava i l ab le for reference and study. I fur ther agree that permission for extensive copying of th is thesis for scho lar ly purposes may be granted by the Head of my Department or by his representat ives. It is understood that copying or pub l i ca t ion of th i s thes i s fo r f i nanc ia l gain sha l l not be allowed without my wr i t ten permiss ion. Department of /̂ {̂ X̂ / t \\ The Univers i ty of B r i t i s h Columbia 2075 Wesbrook Place Vancouver, Canada V6T 1W5 Date VflfajL-K.Jlf. i ABSTRACT The DCLA gravimeter recording obtained after the 1964 Alaska earthguake, has been subjected to various, recently developed data processing techniques. In p a r t i c u l a r , the maximum entropy method (MEM) of spectral estimation and f i l t e r design was u t i l i s e d . Prediction f i l t e r s were used to extend the o r i g i n a l t i d a l gravimeter recording, so as to avoid information loss caused by tapering the record prior to f i l t e r i n g . Time adaptive prediction error f i l t e r s were then used to locate noise bursts, or 'gl i t c h e s ' , which are present i n the f i l t e r e d record. The record was then deglitched using two methods, one a predictive approach, and the other involving a d i v i s i o n of the record by a weighted envelope function. Power spectra using the c l a s s i c a l periodogram approach, as well as MEM, were calculated for the f i l t e r e d records, both before and af t e r deglitching. This analysis resulted i n new values f o r the s p l i t t i n g parameter '^', and the centre frequency, for various free o s c i l l a t i o n modes. A dramatic increase i n the signal to noise r a t i o was also observed after the f i l t e r e d records were deglitched. The presence of core mode o s c i l l a t i o n s was also investigated, but no evidence for these undertones was fcund. Instead, numerous peaks attributed to either instrument n o n - l i n e a r i t i e s , or barometric pressure e f f e c t s , were found i n the freguency range 0.12 to 1.20 cycles per hour. i i i TABLE OF CONTENTS Page ABSTRACT i TABLE OF CONTENTS i i i LIST OF TABLES V LIST OF FIGURES v i ACKNOWLEDGEMENTS v i i i CHAPTER I . INTRODUCTION 1-1. , Purpose 1 1-2. Data D e s c r i p t i o n 2 1-3. Free O s c i l l a t i o n s 7 1-4. S p e c t r a l S p l i t t i n g 10 I - 5. Core Undertones 12 CHAPTER I I . THEORETICAL BACKGROUND I I - 1. C l a s s i c a l S p e c t r a l E s t i m a t i o n A. , Periodogram 14 B. A u t o c o r r e l a t i o n F u n c t i o n 15 I I - 2 . Maximum Entropy Method A. Information and Entropy 16 B. Maximum Entropy S o l u t i o n 18 C. A n a l y t i c a l I n t e g r a t i o n 19 D. Z-transform Approach 22 E. H e u r i s t i c S o l u t i o n 24 .1.1-3. Levinson Recursion 26 I I - 4 . A u t o r e g r e s s i v e Models and MEM 28 I I - 5. Time Adaptive MEM 31 CHAPTER I I I . DATA ANALYSIS PROCEDURES I I I - 1. N e c e s s i t y of Tapering 35 .111-2. P r e d i c t i o n 36 I I I - 3 . F i l t e r i n g 36 i v I I I - 4 . D e g l i t c h i n g A. I d e n t i f i c a t i o n 46 B. P r e d i c t i o n i n t o Gaps 47 C. D e g l i t c h i n g by Envelope 52 I I I - 5. Power S p e c t r a l A n a l y s i s A. Periodgram and MEM Spectra 61 B. Complex Demodulation 64 CHAPTER IV. RESULTS IV- 1. Free O s c i l l a t i o n Mode Frequencies A. oS^ mode 68 B. 0 S 3 mode 78 C. 0S^. mode 86 D. 0S« mode 93 E. 0S f e mode 99 IV-2. Core Undertones 100 CHAPTER V SUMMARY 110 REFERENCES 112 APPENDIX 1 115 V LIST OF TABLES PAGE Table 1. 0S^ Frequencies 69 Table 2. 0S 3 Frequencies 79 Table 3. 0S ¥ Frequencies 87 Table 4. 0SS Frequencies 9U Table 5. „S6 Frequencies 99 v i LIST OF FIGURES Page Figure 1. Unprocessed gravimeter data 4 Figure 2. Periodogram power of unprocessed data 6 Figure 3. Predicted gravimeter data 38 Figure 4. High pass f i l t e r e d predicted data 41 Figure 5. Band pass f i l t e r e d predicted data 43 Figure 6. High pass f i l t e r e d unpredicted data 45 Figure 7. Prediction error trace of high passed record 49 Figure 8. Prediction error trace of band passed record 51 Figure 9. Gap prediction test case 54 Figure 10. Deglitched high passed record 56 Figure 11. Deglitched band passed record 57 Figure 12. Prediction error trace of deglitched high passed record 59 Figure 13. Prediction error trace of deglitched band passed record 60 Figure 14. Deglitched high passed record using envelope sguared 63 Figure 15. Power spectra of 0S a mode 73 Figure 16. Power spectra of 0Sj mode 82 Figure 17. Power spectra of aS+ mode 89 Figure 18. Power spectra of 0SS mode 96 F i g u r e 19. Power s p e c t r a of 0S& mode F i g u r e 20. Periodogram of u n d e g l i t c h e d band passed r e c o r d F i g u r e 21. Periodogram of d e g l i t c h e d band passed r e c o r d v i i i ACKNOWLEDGEMENTS I would l i k e to thank Dr. T. J. Ulrych for many valuable discussions that took place during the course of t h i s research project. Dr. Ulrych also devoted time and e f f o r t i n the c r i t i c a l appraisal of t h i s t h e s i s , as well as provided the necessary f i n a n c i a l support to see t h i s project to an end. My thanks also go to Dr. R. M. Clowes who gave a second opinion of t h i s thesis and to Rob Clayton and George Spence who dropped several useful programming hints my way. Much use was made of the University of B r i t i s h Columbia computing centre and the computing costs were born by a National Research Council of Canada grant made to Dr. T. J. Dlrych (#67-1804). 1 I. INTBODUCTION 1. PUBPOSE T h i s r e s e a r c h p r o j e c t was undertaken i n order t o i n v e s t i g a t e e a r t h f r e e o s c i l l a t i o n data u t i l i s i n g r e l a t i v e l y new technigues i n time s e r i e s a n a l y s i s . The i n v e s t i g a t i o n focused on the w e l l known UCI& gravimeter data obtained by Dr. L. S l i c h t e r a f t e r the Alaska earthquake of March 27, 1964 { S l i c h t e r 1967b). Maximum entropy technigues were used f o r two purposes i n t h i s study: 1. To compute power s p e c t r a of fundamental s p h e r o i d a l o s c i l l a t i o n modes, and t o use the improved r e s o l u t i o n f e a t u r e s of the maximum entropy method to study the a s s o c i a t e d s p e c t r a l s p l i t t i n g . 2. To design f i l t e r o p e r a t o r s necessary t o i d e n t i f y and p r e d i c t s e c t i o n s of the re c o r d which were made untenable by noise g l i t c h e s . The d e g l i t c h i n g c f the data s e t was an important s t e p i n redu c i n g the nois e power of t h i s r e c o r d , so t h a t the presence and d e t e c t a b i l t i y of core undertones c o u l d be i n v e s t i g a t e d . 2 2. DATA. DESCRIPTION The data used i n t h i s study originated from UCLA gravimeter #4, one of two Lacoste-Romberg t i d a l gravimeters that were i n operation i n Los Angeles at the time of the 1964 Alaska earthquake. This seismic event i s the largest to be extensively recorded with long-period instruments, and i s the primary source for much of the present information on normal modes. The data has a sampling rate of 1/12 of an hour and extends for 18.5 days. Free o s c i l l a t i o n information i s superimposed on high amplitude earth tid e data,*and, i n the unprocessed form, the normal mode data i s not observable when the record i s plotted (Figure #1). The data set was recorded on a Lacoste-Romberg t i d a l gravimeter, which in simple terms, i s a mass on the end of a spring-supported beam. A photo-electric c e l l detects motions in the mass,• relays the information to a servomechanism, which i n turn raises or lowers the upper end of the spring by means of a screw. The angular position of the screw i s recorded d i g i t a l l y i n units of 0.1 microgals. Lacoste-Romberg gravimeters generally have d r i f t rates of between 1 and 20 microgals per day. The noise power for t h i s data set increases 28 db. from 0.4 cycles per hour to 0.0 cycles per hour (Figure #2). The f i r s t step i n the data processing was to eliminate the t i d a l frequencies from the record using a f i l t e r i n g scheme. However, since the free o s c i l l a t i o n information s t a r t s at the beginning of the t i d a l record, any frequency O r i g i n a l unprocessed gravimeter data recorded at Los Angeles during the 1964 Alaska earthguake. The record consists of 5320 points with a sampling rate of 1/12th of an hour.  FIGURE 2. Unsmcothed periodogram .power spectrum of the unprocessed gravimeter data. Ihe Nyguist freguency i s at 6.0 cycles per hour. The location of seme of the fundamental spheroidal o s c i l l a t i o n modes i s indicated, filso shewn are the basic diurnal and semi-diurnal t i d e s . CM I I ' I 1 I I 1 I 0.0 0.857 1.714 2.571 3.429 4.2B6 5.143 6.0 cr> FREQUENCY CY/HR. FIGURE 2 7 domain f i l t e r i n g method w i l l r e s u l t i n t a p e r i n g of the f r o n t and t a i l end of the data s e t . I t i s p a r t i c u l a r l y important not t o taper the r e c o r d i n t h i s manner, as high amplitude normal mode i n f o r m a t i o n at the onset of the r e c o r d , s i l l be l o s t . To circumvent t h i s problem, a scheme based on p r e d i c t i v e f i l t e r i n g was used (Olrych e t a l . 1973). Another problem with the data becomes immediately e v i d e n t when low-pass t i d a l f i l t e r i n g i s accomplished. The f i l t e r e d r e c o r d i s contaminated by n o i s e b u r s t s or g l i t c h e s which are caused by d i g i t i z a t i o n e r r o r s and a f t e r s h o c k f e a t u r e s (Wiggins and M i l l e r 1972). To e l i m i n a t e these unwanted noi s e f e a t u r e s , a method of data p r e d i c t i o n was used to f i l l i n reasonable data values (Olrych and C l a y t o n 1975) . 3. FEEJ OSCILLATIONS Free o s c i l l a t i o n s have t h e i r t h e o r e t i c a l beginning with Lamb (1882), who f i r s t c o n s i d e r e d the problem of normal modes of a uniform sphere. Love (1911) then extended the development to i n c l u d e a s e l f - g r a v i t a t i n g uniform sphere. B e n i o f f e t a l . (1954) f i r s t r e p o r t e d e a r t h f r e e o s c i l l a t i o n s a f t e r observing a normal mode with a p e r i o d of 57 minutes. In 1961 v a r i o u s groups r e p o r t e d observing a number of normal modes e x c i t e d by the C h i l e a n earthquake of 1960 (Alsop e t a l . 1961; B e n i o f f et a l . 1961; Bogert 1961; and Ness et a l . 1961). T h e o r e t i c a l i n t e r e s t i n c r e a s e d g r e a t l y at t h i s p o i n t , and the mathematical theory f o r a heterogeneous e a r t h model was developed, a l l o w i n g a comparison between 8 experimental and t h e o r e t i c a l r e s u l t s . I n v e r s i o n schemes were developed, which when u t i l i s e d , y i e l d e d v e l o c i t y - d e p t h models f o r the e a r t h , which i n turn enabled s e p a r a t i o n of Lame's parameters as a f u n c t i o n of depth {Press 1968; Haddon and B u l l e n 1967; Backus and G i l b e r t 1967; and Dziewonski and G i l b e r t 1973) . There are two fundamental types of f r e e o s c i l l a t i o n s of a sphere: 1. S p h e r o i d a l : r a d i a l displacements always e x i s t . 2. T o r o i d a l ; displacements are e n t i r e l y t a n g e n t i a l . S p h e r i c a l geometry can be broken down i n t o three components, l a t i t u d e { 0) , l o n g i t u d e ( f) , and r a d i u s ( r ) . S l i c h t e r (1967a) g i v e s the displacements f o r these t h r e e components f o r t o r o i d a l modes as: 1. 1^=0 2. O 0 = (r,£,m) c o s ( e ) - » P™(sin © ) X "dexp f - i ojUt+mcar-*/! \ 3. oy =•-v„ ( r , £ , a ) 3 vTisxn e) X exp (-i^ft+may-Y} ) 9 For spheroidal modes the displacment components are; 1. 0^ = 0^ (r,^,m) P™(sin©) X exp (-i ct>s {t+mittj-1/}) 2. 0 e = W^Cr,Hfm) 9Posing ) X exp (-idOgCt+mt^-1^}) 3. 0? = W^r^m) (cos©)- 1- pJVin© ) X J i ^ p l - l ^ f t - S r a ^ - i / } ) P̂ 1: spherical surface harmonic of degree *J0* and order . «m« (m=0,1,2...l) . For the r a d i a l functions V^O^and , • n• denotes the overtone number, where the fundamental overtone i s ' n=0*. Spheroidal free o s c i l l a t i o n s are represented by the symbols ^S^ and t o r o i d a l o s c i l l a t i o n s are represented by . The 0 S 0 mode i s e n t i r e l y a r a d i a l motion of the earth, expanding and contracting much as i f a balloon were being i n f l a t e d and deflated p e r i o d i c a l l y . Another spheroidal mode, eS^ , i s nicknamed the f o o t b a l l mode and represents a change of shape of the earth from prolate to oblate and back again. Nodal points do ex i s t and i t i s possible some modes w i l l not be observed when the recording station i s located at one of these points of minimum amplitude. Some free o s c i l l a t i o n modes are not possible for the earth. . 0S ( i s a 10 r i g i d body o s c i l l a t i o n and i s not o b s e r v a b l e . i s a r i g i d body t r a n s l a t i o n about the p o l a r a x i s and has not yet been observed, and f o r the mode n T ^ a l l the displacements are e g u a l t o zero. 4. SPECTRAL SPLITTING• A f t e r the r e c o r d s from the 1960 C h i l e a n earthquake were examined, s p l i t t i n g of some of the fundamental s p h e r o i d a l modes was observed. An analogous phenomenon i n n u c l e a r p h y s i c s i s the Zeeman e f f e c t where atomic s p e c t r a l l i n e s show s p l i t t i n g i n the presence of a magnetic f i e l d . R o t a t i o n and e l l i p t i c i t y of the e a r t h were b e l i e v e d to be two p o s s i b l e causes f o r t h i s behaviour. S u i t a b l e theory was developed and i t was found t h a t r o t a t i o n would cause the s p l i t t i n g of the c e n t r a l peak i n t o 2^+1 peaks, and e l l i p t i c i t y would cause a s p l i t t i n g i n t o ^+1 peaks. Only s p h e r o i d a l modes have been observed t o show s p l i t t i n g , as t o r s i o n a l modes attenuate too g u i c k l y , r e s u l t i n g i n d i s s i p a t i v e broadening of the s p e c t r a l l i n e s . The f o l l o w i n g f o r m u l a t i o n taken from Alterman e t a l . (1974) p r e d i c t s s p l i t t i n g e f f e c t s t o f i r s t order i n e l l i p t i c i t y and second order i n r o t a t i o n . (1 + ^ + «4 + m 2^) \ = 0 ( ^ ) 2 + E £ LU^ = s p l i t peak freguency 11 m = number of the s p l i t peak •A- = angular freguency of the earth's rotation 1X3 e, = eigenfreguency of S or T £ = e f f e c t of earth's e l l i p t i c i t y A,B,C,D,E = constants The centre frequency of a group of s p l i t peaks w i l l d i f f e r by from the frequency of an unsplit peak. Also, because of the e f f e c t of 'm2', s p l i t t i n g w i l l be asymmetric around the c e n t r a l peak (n=0). For lower order modes, the rotation e f f e c t i s much qreater than the e l l i p t i c i t y e f f e c t , so that the f i r s t order s p l i t t i n g e f f e c t in e l l i p t i c i t y i s of the same magnitude as the second order s p l i t t i n g e f f e c t in ro t a t i o n . For the higher order modes, t h i s equivalence does not hold. Lateral heterogeneities can also give r i s e to s p l i t t i n g or widening of spectral l i n e s , but the e f f e c t s cannot be well predicted or quantized. Up to the present, only conventional power spectral techniques have been used i n the study of free o s c i l l a t i o n s p l i t t i n g . MEM seems to be a l o q i c a l choice i n studying these e f f e c t s because of the higher resolution and optimally smooth nature of the HEM power spectrum. Bolt and Currie (1975) and Olrych and Bishop (1975) have shown the improvement in frequency resolution when using MEM, Bolt and Currie (1975) studied the t o r s i o n a l eigenfreguencies recorded at Trieste a f t e r the 1960 Chilean earthquake, and reported an increase i n precision and number of peaks detected. 12 5. CORE UNDERTONES One o f t h e aims o f t h i s p r o j e c t was t o l o o k once a g a i n f o r c o r e mode o s c i l l a t i o n s i n t h e f r e g u e n c y r a n g e 0.0 t o 1.0 c y c l e s p e r h o u r . S l i c h t e r (1961) r e p o r t e d a s p e c t r a l peak a t .698 c y c l e s p e r h o u r i n t h e r e c o r d o f t h e 1960 C h i l e e v e n t . He a l s o p r e d i c t e d t h a t t h i s p a r t i c u l a r mode s h o u l d show s p e c t r a l s p l i t t i n g peaks a t .741-and .659 c y c l e s p e r h o u r . The 1964 A l a s k a e a r t h q u a k e y i e l d e d no s u c h p e a k s , and i t i s now g e n e r a l l y b e l i e v e d t h a t t h e peak r e p o r t e d i n 1961 was n o t v a l i d . I t has been shown ( C r c s s l e y , p e r s o n a l c o m m u n i c a t i o n 1976) t h a t t h e p e r i o d o f t h e S l i c h t e r mode i s v e r y model d e p e n d e n t , and t h e c o r r e l a t i o n between S l i c h t e r * s o b s e r v a t i o n and p r e d i c t e d p e r i o d s i s v e r y low. S p h e r o i d a l u n d e r t o n e s c a n o n l y e x i s t i n a s u b - a d i a b a t i c c o r e , as any r a d i a l p a r t i c l e m o t i o n i s p o s s i b l e o n l y when t h e medium i s g r a v i t a t i o n a l l y s t a b l e . R o t a t i o n a l c o u p l i n g between t o r s i o n a l and s p h e r o i d a l modes r e s t r i c t s t h e t h e o r e t i c a l r a n g e o f t h e u n d e r t o n e p e r i o d s t o between 0 and 12 h o u r s . At l o n g e r p e r i o d s , s t a r t i n g a t a b o u t 72 h o u r s , R o s s b y waves c a n t h e o r e t i c a l l y e x i s t ( C r o s s l e y 1975). The most i n t e r e s t i n g o f c o r e u n d e r t o n e s i s t h e S l i c h t e r mode, which i s a t r a n s l a t i o n a l m o t i o n o f t h e s o l i d i n n e r c o r e d r i v i n g a r e t u r n m o t i o n i n t h e f l u i d o u t e r c o r e . J a c k s o n and S l i c h t e r (1974) b e l i e v e t h e y c o u l d d e t e c t s u c h a m o t i o n i f i t p r o d u c e d a s i g n a l o f 0.08 m i c r o g a l s on a L a c o s t e - R o m b e r g g r a v i m e t e r a t t h e S o u t h P o l e . T h i s s i g n a l a m p l i t u d e c o r r e s p o n d s t o a c o r e m o t i o n o f 6 c e n t i m e t e r s w i t h an i n n e r - 13 outer core d e n s i t y c o n t r a s t of 0.3 g/cm- 3. T h i s minimum d e t e c t a b l e s i g n a l would be c o r r e s p o n d i n g l y l a r g e r f o r s t a t i o n s at other l o c a t i o n s on the e a r t h . The o b s e r v a t i o n of a core undertone such as the S l i c h t e r mode would-be-extremely u s e f u l i n a c c u r a t e l y determining the d e n s i t y c o n t r a s t between the i n n e r and outer core. At t h i s time i t i s unknown whether core undertones do e x i s t , and whether they are observable phenomena. T h i s l a c k of knowledge i s due to t h r e e f a c t o r s : 1. The l e v e l of e x c i t a t i o n of core modes i s b e l i e v e d to be q u i t e low (Smith 1974). 2. Undertones are non-degenerate with r e s p e c t to the azimuthal order number 'm*. The number of undertones i s t h e r e f o r e '2n+1* g r e a t e r than the corresponding overtones, making i d e n t i f i c a t i o n o f i n d i v i d u a l modes d i f f i c u l t ( C r o s s l e y 1975). 3. I t i s p o s s i b l e t h a t many of the undertones w i l l have p e r i o d s packed around an upper l i m i t of 12 hours. This l i m i t i s based upon c a l c u l a t i o n s done by C r o s s l e y (1975) who used a formulated computing scheme t o a r r i v e at t h i s r e s u l t . T h i s 12 hour p e r i o d corresponds to t h a t of the s e m i - d i u r n a l t i d e , and so power from t h i s t i d a l peak w i l l tend to obscure the s p e c t r a l r e c o r d . 14 I I . TJJOBJTICJL BACKGROUND 1. CLASSICAL SPECTBAI ESTIMATION Orthodox methods of spectral estimation often make u n r e a l i s t i c assumptions about the extension of the sampled process beyond the windowed region. The maximum entropy method of spectral estimation was developed s p e c i f i c a l l y to obviate t h i s r e s t r i c t i v e supposition. The two ' c l a s s i c a l * methods of obtaining a power spectrum are outlined below. A. Periodggram A function f ( t ) can be expanded in a Fourier series i f f ( t ) i s piecewise continuous i n the i n t e r v a l -1/2 < t < T/2 and i s perio d i c with period T , Then f ( t ) can be expressed as a sum of sine and cosine functions. f (t) = 1/2aG+ ̂  (ancos«;nt + b^sinuj^t) -n=i n=0,1,2... n=1,2,3... The amplitude spectrum for the process f (t) i s given as 15 ! f n | =Va2 • b* Due to-its'Conditions of existence, the periodogram method assumes a periodic extension of the data beyond the sampled region. B » M i s s °££§1§ l i s s fjyssliSJQ Seiner (1930) f i r s t shoved the r e l a t i o n s h i p between the autocorrelation function in the time domain and the power spectrum i n the freguency domain. The f i r s t step i n obtaining a power spectrum using t h i s method i s to obtain a suitable estimate of the autocorrelation funtion 8 ( t ) , where fi (t) — 0 f o r | t | > 8. This i s eguivalent to saying that an i n f i n i t e non-zero autocorrelation funtion H (t) i s multiplied by some weighting function H (t) where I (t) •- 0 for | t j > N, The product E(t)W(t) i s then transformed to the freguency domain. The estimated spectrum i s therefore the convolution of the Fourier transform of the weighting function with the true spectrum. Considerable work has been put into trying to design the best possible weighting function and i t s transform. There are two d i f f i c u l t i e s involved i n the design, the f i r s t being the tradeoff that exists between resolution and variance i n the freguency domain, and the second being the requirement that the spectrum be non-negative. As Burg(1975) points out, although window theory has a certain neatness to i t , i t i s an a r t i f i c e imposed by the assumption that B(t) - 0-for-It| > N. 16 " I f one were not blinded by the mathematical elegance of the conventional approach, making unfounded assumptions as to the value of unmeasured data, and changing data values that one knows would be t o t a l l y unacceptable from a common sense, and hopefully a s c i e n t i f i c point of view" 1. 2. MAXIMUM ENTBOPY METHOD A. Information and Entropy Burg's r a t i o n a l in developing the theory of the maximum entropy method (MEM), was, that since there i s an i n f i n i t e set of non-negative power spectra that correspond to a given autocorrelation function, the best spectrum would be the one that corresponds to the most random time s e r i e s . The concept of maximum entropy spectral analysis is'based upon t h e o r e t i c a l developments made i n s t a t i s t i c a l mechanics and information theory. One of the great advantages of the MEM technigue, i s i t s a b i l i t y to resolve cl o s e l y spaced spectral l i n e s i n a much better manner than conventional methods. From a s t a t i s t i c a l viewpoint, MEM can be shown to correspond to a s i t u a t i o n with the most random behaviour. Following Ulrych and Bishop (1975) i t can be seen that i f there are 1M* d i f f e r e n t things •m,-' which could possibly occur, a l l with p r o b a b i l i t i e s ' p L ' f then i f a l l the *p^*s' 1 J. P. Burg, Maximum Entropy Spectral Analysis, Ph.D. Thesis, Stanford University 1975. 17 a r e t h e same, t h e n no i n f o r m a t i o n about the system has been g a i n e d . I f one of the 'p^'s* i s d i f f e r e n t , t h e n some i n f o r m a t i o n about the p r o c e s s has been a g u i r e d . A r e l a t i o n s h i p between i n f o r m a t i o n and p r o b a b i l i t y can be e x p r e s s e d as I = k l o g < V P c ) k = 1 when l o g i s base 2. I f t h e o c c u r r e n c e s a r e summed over a l o n g p e r i o d T then Iror/ic = MP,T l o g d / P , ) + P aT log(1/P^) + ...) The average e n t r o p y g i v e n by Shannon (1948) i s = k f P . l o g ( P L ) II-1 l-1 I f a l l t h e • P,;* but one a r e e g u a l t o z e r o , and t h i s e x c e p t i o n i s e g u a l t o one, then H = 0, A g a i n , a l l knowledge about the p r o c e s s i s a v a i l a b l e , and t h e r e i s no u n c e r t a i n t y i n t h e system. I f H i s g r e a t e r than z e r o , a measure of u n c e r t a i n t y i s a v a i l a b l e . The e n t r o p y r a t e o r e n t r o p y per u n i t sample i s d e f i n e d as l i m H/ (N + 1) 18 and can be shown to be proportional to In P (f) df P (f) = power spectrum II-2 f., = Nyguist freguency 8. Maximum Entropy Solution To solve the maximum entropy v a r i a t i o n a l problem, the following argument i s taken from Burg(1975). The function P (f) that maximizes eguation II-2 must be obtained so that the constraint eguations are s a t i s f i e d . B (n) i s the autocorrelation function. Another condition i s the Fourier transform relationship between the autocorrelation and the power spectrum. The constraint eguation II-3 w i l l be s a t i s f i e d i f the B (n) i n equation II-3 are equivalent to the B(n) i n eguation 11-4. To solve t h i s problem, lagrangian m u l t i p l i e r s can be used, but what follows i s a simpler derivation. Substituting II-4 into II-2 gives (-N < n < N) II-3 P(f) = 1/2f„ 2 B(n) exp {-i2 7T f t) II-4 19 h << \ ln[ 1/24 2. B («) -«*P ( - i 2 n tn a t) ] df II-5 Maximizing II-5 with respect to E(n) r e s u l t s i n h,(.1/2f N y P-» (f) exp<-i2 n fsn t) df = 0 for |s| > N II-6 where P - 1 (f) i s the inverse of the power spectrum P ( f ) . Expanding P"~1 (f) i n terms of a Fourier series gives P-* (f) = ^ \exp (-i2 n fn a t) .'. P(f) = 1/j£ A 5exp (-i2 ft ts a t) II-7 II-7 Is the same eguation that would r e s u l t from the v a r i a t i o n a l s o l u t i o n , where the ? v*s are the Lagrangian m u l t i p l i e r s . From II-6 i t can be seen that A1N=0 for In j > N. The next task i s to express the /\ 's i n a suitable form. Two methods are presented, one an a n a l y t i c a l integration scheme, and the other a z-transform argument. C. M a l y i i c a l Integration Substituting eguation II-7 into eguation II-3 gives ft exp(i27r fnA t) df = fi<n) (-N < n < N) J ^ x p ( - i 2 7T fsat) II-8 Letting z = exp(-i2ff f / l t ) which i s the normal z-transform equivalence i t can be shown that 20 dz = - i 2 T i a t exp(-i2Hf At) df = -±2Tf£it z df and df = -dz/i2n£t z. Equation II-8 can now be written as 1 / 2 7 T i 4 t ^ z - ^ - V J : A s z s dz = R(n) I I - 9 This integration curve i s counterclockwise around the unit c i r c l e i n the complex z plane. From II-7 and the p r i n c i p l e of polynomial f a c t o r i z a t i o n i t i s possible to say that 2 \zs = [PNZJt]-» [1 • a, z +...+ a^z^ [1 + a%-»- +...• a*z~ w3 - [P NAt]-» Z a s z S £ a * 11-10 Substituting 11-10 into II-9 re s u l t s i n an expression f o r B (n). fi(n) = P»/ L z~^ - 1 dz - N < n < N 2-iTi ™ i a, z s £ as z " 5 11-11 Since £ (n) = \ P ( f ) z - ^ and i f 11-11 i s summed over a so that i\--o 2rri J 2 a 5 2 5 Z a s z " dz 21 i t i s possible to write ,Z.a* R(n-r) = P/o C z^-* 5. a t ^ d z ^ 2~~T3 i a zs ^ a r Z - 5 = P P C ds r>0 2/ri J 1" a, z s .11-12 Knowing that 2 a z i s analytic for z < 1, and for r > 1, s«e> z ^ - 1 i s analytic for z < 1, i t i s evident that l z ^ - i / - s a z2] dz = 0 for r > 1 S - O and also that j"f a^R (n-r) = 0 for r > 1 11-13 T\-0 Now that the case for r > 1 has been solved a l l that need be done i s obtain a solution for the r = 0 case. Since [1/27ri] ^ [ f ( z ) / z ] dz = f(0) i t can be shown that 1/2 7Ti ^ ( Z a 5 z s ) - * / z dz = f (0) where f (z) - ( 2 a^z a chas been already defined as M", so that f (0) = 1. Therefore 11-12 becomes 22 B(n-r) = E w for r = 0 11-14 T l v - o Taking the complex conjugates of 11-13 and 11-14 we get ^ B ( n - r ) a„= Pw r = 0 II-15 Z B(n-r) a~ = 0 r > 1 11-16 An expression for the power P(f) can also be written by substituting 11-10 into II-7 to give P(f) = P^At/jr a s z 5 £a*z- 5 11-17 S-.o which can be rewritten in the f a m i l i a r form N P(f) = P N / ( 2 ^ 11 • 2 V e x P (-i2 77 fs o. t) | 2) 11-18 s-o D. ^tra n s f o r m approach An alternative method to obtain equations 11-15 and 11-16 i s through a z-transform approach. Substituting 11-14 and 11-10 i n t o II-6 gives S*o s-o = (V2f N)2 E<r) zA~ 11-19 A.s-»0 KJ ^ Multiplying by ^ a ^ z - ^ gives 23 S=<5 T K O 1 A.---0O = £ a*z-*£ S (n-r) z*-^ -n-.o = 2 [|<S(n-r) ] z-^ A* ̂  V\_ . 11-20 When s = 0 then 11*0 11-21 and when s > 0, powers of z do not match so that 0 =2 a^B (n-r) 11-22 Taking the complex conjugates of 11-21 and 11-22, equations 11-15 and 11-16 are again obtained. Writing 11-15 and 11-16 in matrix form gives I B (0) B(1) I B(1) H(0) I I . . I B(N) R(N-1) . B(H) . B(N-1) B(0) i 1 J I a, | I . I I • I I 0 | I • I ! • I I 0 I 11-23 which i s the eguation for obtaining the N-t'h point prediction f i l t e r . From 11-18 i t can be seen that to obtain a spectral estimate P (£), a l l that i s required i s a prediction error f i l t e r and the error power cf prediction. These twc values are qiven by equation 11-23. 2H E. Heuristic Solution From Ulrych et a l . (1973) a h e u r i s t i c solution to the maximum entropy problem can be rep l i c a t e d . I f there exists a time series x (t) which has a Fourier transform X ( t c ) , there i s some f i l t e r i (t) that whitens x ( t ) . The maximally white time series that corresponds to a par t i c u l a r power spectrum has already been shown to be the one that exhibits maximum entropy. If the transfer function of i (t) i s I (to), then | X ( L U ) I(u>)| equals a constant 'k». Also | X ( « ; ) |2 = k V I K w ) I2 i s a power estimate. Any stochastic,non-deterministic, stationary time series can be represented as a moving-average process. x(t) = Z*>T e^_ s b Q > 0 b 0z + b,2 + b^2 •+...< oo e^. = a white noise s e r i e s . Therefore x(t) i s the convolution of some f i l t e r or wavelet b (t) with a white noise ser i e s e ( t ) . b ( t ) - 1 i s the f i l t e r that whitens x ( t ) . In other words, i f a(t) * b(t) = £ (t) , then x (t) * a(t) - e ( t ) . a(t) i s a deconvolution or spiking f i l t e r and can be determined from 25 .8 a •= 1 11-24 B i s an N by N autocorrelation matrix a = a, >a^ ,...a n i s the prediction error f i l t e r 1 = (1,0,0,.. .0) The prediction f i l t e r (t) = 8(t) - a (t) so 11-24 can be written as = ft" — 1 # -a | , — a ̂ , •. • — a PN= (P^,0,0,...0) f»-l P w - r Q + 2 R L = e r r o r power of prediction. L-I From the power estimate IX (w) 12 = k V I K w i l 2 i t can be seen that u P(f) = P w /H + 2 rKexp (-12 77fk) |2 which i s equivalent to equation 11-28. From the above discussion, the fundamental difference between HEM and conventional power spectral analysis can be seen. Maximum entropy spectral analysis gives the spectrum of a model that f i t s or approximates a f i t to the process, a sample of which i s represented by the data. Conventional 26 f o u r i e r transform methods, on the other hand, obtain the spectrum of the sample. 3. LEVINSON J J C O B S I C N Burg developed a solution to eguation 11-23 based upon the Levinson recursion algorithm. The solution i s also consistent with the concept of maximum entropy. The accuracy of Burg's method i s greater that that of the normal Levinson method, and i t i s also executed about 1/3 f a s t e r . The development i s shown below. A prediction f i l t e r i s run over the data set i n both the forward and reverse d i r e c t i o n s , and the power output from t h i s operation i s minimized. The f i l t e r i s not run off the ends of the data, and so no assumptions about the continuation of the data outside the sampling i n t e r v a l are made. This i s the major feature that i s related to maximizing the entropy of the process. For the two point f i l t e r ( 1,Y), the power output from running the f i l t e r over the data i s = I/2(N - D i2 (x d + ( + y x - ) 2 • t - l The value of P^ can be shown to be a minimum when = 2 | ' ( x i x U ( ) / | ' ( x . ) 2 11-25 27 For the two point case, eguation 11-23 i s written as | B(0) B{1) | | 1 | I Pol I 1 1 1 = 1 1 I Rd) E{0) I | y I I 0 I ^ i s obtained from eguation 11-25. The f i r s t autocorrelation A; value B(0) i s estimated by B(0) = 1/H ~Z. (x / ) 2 , and sc L-i B(1) can be calculated from B(1) = - ^ ( B ( 0 ) ) , and P^ can be obtained from P^= B(0) (1 - ^ 2 ) . For the three point f i l t e r (1 , V,, ^ ) the power output after running the f i l t e r over the data i s *-a P~ = (1/2 (N-2)) 2 { ( X l + a + X l + lV a+X (: j ^ ) 2 When ^ i s set egual to ]f(1 +^ 3), where Y i s the f i l t e r c o e f f i c i e n t obtained for the two-point case, then the value of Xjthat minimizes P 3 can be obtained. The f i l t e r { 1# )fn *Xi)'^3^ w i l l also be minimum phase. Eguation 11-23 for the three point case i s written as 28 1 B(0) B<1) B (2) J J 1 I | 0 | 1 B(1) B(0) B(1) | | y | + Yjl Y I = | B(2) B(1) E (0) | \ 0 I 1 1 | I P ^ l 1 4,1 1 o | + y 3i 0 | ! 4*1 I F A I Zl a= B{2) + B(1))f & - - v p * The advantage of t h i s recursive method i s that only one variable, i n t h i s case }fs * n a s t o b e adjusted so as to obtain the minimum power output. The recursion i s continued u n t i l the desired f i l t e r length i s obtained. The values for the prediction f i l t e r { 1, x* , o" a#» # o^ j ) and the error poser of prediction P N, are then plugged into equation 11-18 tc obtain an estimate of the power P ( f ) . As Lacoss(1971) and Burg(1972) pointed out; the power estimate P(f) i s actually a spectral density estimate. If a more r e l i a b l e estimate of the spectrum (in terms of r e l a t i v e peak amplitudes) i s required, the maximum entropy estimate should be integrated. I P 3 I = 1 0 | 1 0 | 4. ASJ'CREGRESSIVE MODELS AND MEM One of the c r i t i c a l factors i n determining the maximum entropy spectrum i s deciding upon which length f i l t e r 29 o p e r a t o r should be used i n the s p e c t r a l estimate c a l c u l a t i o n . Too high an order l e a d s t o i n s t a b i l i t y f a i r l y q u i c k l y , whereas too low an order gives a poorly r e s o l v e d e s t i m a t e . In some cases prior'knowledge of the spectrum (throuqh c o n v e n t i o n a l methods) can be very u s e f u l i n p i c k i n g a f i l t e r o r d e r . Due to the work of VandenBos (1971), who showed the equ i v a l e n c e of a u t o r e g r e s s i v e time s e r i e s modelling and the p r i n c i p l e of MEM, the l a r g e g u a n t i t y of l i t e r a t u r e on a u t o r e g r e s s i v e modelling can be a p p l i e d to the f i l t e r l e n g t h e s t i m a t i o n problem. In p a r t i c u l a r the works of Akaike (1969* 1970) and Gersch(1970) are worth n o t i n g . The a u t o r e g r e s s i v e model f o r an order M, i s w r i t t e n and i t can be seen t h a t the present value x(t) of the time s e r i e s i s being produced from the known past values x̂ ., ,x^.a -,.. x^.w. a^ i s an a d d i t i v e white no i s e s e r i e s c a l l e d the i n n o v a t i o n and has zero mean and a v a r i a n c e of t Y ^ 2 . I f the z-transform o f 11-26 i s taken, the r e s u l t can be w r i t t e n as X(z) - X(z) [ ^ z +\z* +...+ = A(z) and t h e r e f o r e | X ( z ) | 2 = |A(z )|V I 1 -(k,z - ^ z 2 < ^ z V H-27 30 Equation 11-27 can be rewritten as P ( f ) . = 2 ^ / 1 1 -£^.exp(-i2^fj) I 2 which i s equivalent to 11-18, the expression for the maximum entropy spectrum, with the equivalence of autoregressive modellinq and MEM shown, the f i n a l prediction error c r i t e r i o n of Akaike can be used to help determine the f i l t e r Order. The FPE i s defined as the mean square prediction error. 11-28 of xf The FPE can be shown to be (Ulrych and Bishop 1975) FPE^ » p + (M~1)/N - (M + 1)} S2 M i s the order of f i l t e r S 2 i s the minimum residual sum of squares M In many cases the FPE c r i t e r i o n qives an excellent idea cf the f i l t e r length, but for the data being studied i n t h i s project, i t was found not to be as useful. In numerous examples a guess of the optimum order was made, in some cases afte r several test runs were performed. M was never taken to be greater than N/2 where N i s the length of the data set. FPE = E [ (xt - x t ) 32 A x = estimated prediction 31 5. TIME ADAPTIVE MEM As s e l l as handling stationary time s e r i e s , maximum entropy can be expanded into the realm of non-stationarity i n space and time. The use of time-adaptive MEM was considered here as a possible means of i d e n t i f i c a t i o n of the noise bursts i n the data. These gl i t c h e s can be thought of as reverberant peaks superimposed on a non-stationary time s e r i e s . Time adaptive MEM i s a deconvolution scheme that designs a prediction error f i l t e r i n order to locate areas of u n p r e d i c t a b i l i t y . Seismic a r r i v a l s i n r e f l e c t i o n records, and the noise g l i t c h e s i n this data set, can both be thought of as regions of un p r e d i c t a b i l i t y . The deconvolved record should show spikes corresponding to the location of the g l i t c h e s , and a f t e r t h i s i d e n t i f i c a t i o n has been made, the noise regions can be removed by some predictive method. Normal time adaptive schemes are handled by time gates, and two possible approaches can be used. 1. The autocorrelation lags for each gate can be calculated, the prediction error f i l t e r s can be determined from these lags using the Levinson recursion algorithm, and the f i l t e r s can then be applied to the data. 2. The Burg approach can be used to calculate the prediction error operators i n the time gates, and the f i l t e r can be applied to the data. 32 Burg (1971) has developed an alternative to the time gating approach, where the prediction error c o e f f i c i e n t s are continuously updated as the f i l t e r i s run across the data. This eliminates the problem of picking time gates to cover only stationary data segments. A brief description of the processing scheme follows. The parameters discussed are described i n the comment section of the program l i s t e d i n APPENDIX 1. A variable NWABM controls the length of a warmup segment of data. A prediction error operator i s calculated on th i s 'stationary' data segment in the same manner as for the normal Burg technigue where the f i l t e r length i s set egual to 'LCN +1' points. This f i l t e r i s then run over the f i r s t 'NWARM' points to fi n d the prediction error of the * NHARM • 1* point. The f i l t e r i s new recomputed using only the 'LCN* pcints of the segment from 'N8ARM - LCN +1* to * NWARM + 1'. The prediction error f o r the 'NWARM + 2* point i s then found. The process i s continued u n t i l the end of the data segment i s reached. If the parameter 'IFIGAP' i s greater than 1, then only every *IFLGAP* point i s used i n the c a l c u l a t i o n of the f i l t e r and the determination of the prediction errors. After the prediction error trace for every •IFLGAP* point has been calculated, intermediate prediction error values are computed;,There i s also a parameter *NZC* which controls the number of r e f l e c t i o n c o e f f i c i e n t s set to zero. This has the e f f e c t of creating a f i l t e r that predicts further than unit distance, while s t i l l maintaining minimum phase requirements. 33 As discussed e a r l i e r , the optimum prediction error f i l t e r i s the-one with the minimum output power as i t i s run over the data. For the time adaptive case, the minimum average squared power must be set by the constant updating of the f i l t e r c o e f f i c i e n t s . A recursive scheme as outlined for the stationary MEM i s u t i l i s e d to obtain the N + 1 point f i l t e r from the N point f i l t e r . Only one parameter C ̂  i s adjusted to minimize the average sguared error. corresponds to a r e f l e c t i o n c o e f f i c i e n t i n an i d e a l l y layered medium. Burg (1972) has shown that the average power i s P̂  = CV(2N)} f [F. (k) + C d B ^ ( k ) ] 2 + [Cj Frf (k) + Bj (k) ] 2 which i s the sum of the sguared forward and backward prediction errors. The value cf which minimizes the power j> ± i s given as which i s the negative average cross-power of the forward and backward prediction error, divided by the averaged auto-powers of the forward and backward prediction errors. To handle nonstationarity the cross-powers and auto-powers are weighted so that the s t a t i s t i c s close to the point predicted 34 w i l l be more i m p o r t a n t t h a n s t a t i s t i c s f u r t h e r away. The w e i g h t i n g i s e x p o n e n t i a l and i s b o t h t i m e and s p a c e v a r i a b l e H e i g h t (k,m) = e x p ( - k d T / r r - m/^ ) % = r e l a x a t i o n t i m e (seconds) % - r e l a x a t i o n d i s t a n c e ( t r a c e s ) The r e l a x a t i o n t i m e i s 1/e t h e d e c a y t i m e . S i n c e t h e r e i s a t r a d e o f f between spee d y a d a p t i o n and a c c u r a t e c o m p u t a t i o n o f C: , t h e use o f MEM i s p a r t i c u l a r y a d v a n t a g e o u s i n t h a t i t a l l o w s c a l c u l a t i o n o f c o e f f i c i e n t s on s h o r t d a t a segments. A " F o r t r a n IV l i s t i n g o f t h e program i s i n c l u d e d i n A p p e n d i x #1. 3 5 I I I . DATA ANALYSIS PROCEDURES As discussed in Chapter I I , two major problems exist with the time series under investigation. One i s the fact i that the onset of the free o s c i l l a t i o n information i n t h i s data corresponds to the.beginning of the record. The second d i f f i c u l t y i s the presence of high amplitude noise g l i t c h e s which occur i n the f i l t e r e d records. Both of these problems can be overcome by use of prediction f i l t e r s designed with the Burg algorithm. 1. NECESSITY OF TAPERING In order to accomplish freguency domain f i l t e r i n g , i t i s es s e n t i a l that the data set be zero mean and tapered in some manner before transformation to the freguency domain. Otherwise unnecessary noise power due to step d i s c o n t i n u i t i e s i n the data w i l l be added to the spectrum. If a data set i s not zero mean and tapered, i t can be broken down into two parts, one being the signal information, and the other being a D.C. component. This D.C. sign a l i s i n e f f e c t a step function which represents unwanted information. The Fourier transform of a step function u (t) = 1/2 + 1/2 sgn (t) i s <p [u(t) ] = ifSlcu) + 1/itu 36 Therefore a step function w i l l add a D.C. component, as well as a freguency dependent component that i s higher at low freguencies. Therefore the necessity of removing the mean and tapering the data set before transforming to the freguency domain, i s clear. 2. J?lJDICTION Since the free o s c i l l a t i o n information s t a r t s at the onset of the record, any tapering w i l l remove desired information. To eliminate t h i s problem the data set was predicted forward and backwards i n time so that necessary information was preserved. The 500 point prediction f i l t e r was designed using the Burg maximum entropy algorithm. The o r i g i n a l 5320 points were predicted up to 8092 points, which represents an extension of the data 1386 points both forward and backwards. This record was then tapered with a 5% cosine b e l l taper (Figure #3). 3. FILTERING Two f i l t e r e d data sets were obtained from the predicted o r i g i n a l data. One was a high pass f i l t e r e d record with a cutoff freguency of 0.66 cycles per hour. This record was used to study the fundamental free o s c i l l a t i o n modes and the associated spectral s p l i t t i n g . The second data set was a band pass f i l t e r e d record with cutoffs of 0.-12 and 1.20 cycles per hour. This band passed data was used for the investigation of Predicted gravimeter data extended forward and backward i n time. Prediction f i l t e r length was 500 points and 1386 points were added in each d i r e c t i o n . The position of the o r i g i n a l record i s indicated. 38 < Pred ic t ion FIGURE 3 4 0 0 0 yqals 3 9 the core undertones. Both f i l t e r i n g operations were performed i n the freguency domain by a smoothed boxcar m u l t i p l i c a t i o n operator. The seguence of f i l t e r design i s as follows: 1. The boxcar function with the desired cutoffs i s constructed in the freguency domain. 2. This boxcar i s then transformed to the time domain and the time function i s truncated and tapered to some desired length. 3 . This truncated time function i s retransformed back to the freguency domain, re s u l t i n g i n a smoothed boxcar. The amount of smoothing i n the freguency domain depends upon the amount the time function has been truncated i n the time domain. The longer the time domain function, the sharper the cutoff i n the freguency domain. Figure #4 i s the high pass f i l t e r e d predicted record and Figure #5 i s the band passed predicted record. The noise g l i t c h e s are immediately evident i n both data sets, and can be seen to have very high amplitudes in r e l a t i o n to the s i g n a l . A comparison should be made between Figure #4 and Figure #6, which i s the high pass f i l t e r e d not predicted uo High pass f i l t e r e d predicted record. The f i l t e r cutoff i s at 0.66 cycles per hour. The noise glitches are immediately evident, and have a s e l l defined character. 12 h r s FIGURE 4 1 0 0 /QaiS Band pass f i l t e r e d predicted record. The f i l t e r cutoffs are at 0.12 and 1.20 cycles per hour. Glitches are observable, but not as well defined as for the high passed record. 43 High pass f i l t e r e d unpredicted o r i g i n a l record. The tapering, which i s necessary to accomplish the f i l t e r i n g , i s immediately evident at the onset of the record, where there i s a large reduction i n the si g n a l amplitude. 45 | 12 hrs FIGURE 6 100 )/gals 46 Alaska record. This i s the record that would r e s u l t from tapering and f i l t e r i n g the t i d a l record before prediction forward and backward in time. The difference i n the information available at the beginning of the record i s obvious, and much better s i g n a l to noise r a t i o s would be expected i n the power spectrum. 4. DEGLITCHING As seen on both Figures #4 and #5, high amplitude g l i t c h e s are present i n the f i l t e r e d records. These gli t c h e s contribute a s i g n i f i c a n t amount cf noise in the frequency range 0.0 to 2.0 cycles per hour (Wiggins and M i l l e r 1972). A. I d e n t i f i c a t i o n The f i r s t step in elimination of these g l i t c h e s i s i d e n t i f i c a t i o n . For t h i s record v i s u a l i d e n t i f i c a t i o n can i s o l a t e most of the noise areas at l e a s t for the high-passed record (Figure #4). For the band passed record, the g l i t c h e s do not have the same well defined character and are harder to spot, p a r t i c u l a r l y i n the e a r l i e r part of the record. To aid i n t h i s i d e n t i f i c a t i o n , i t was thought that i t would be useful to run the time-adaptive MEM routine over the f i l t e r e d records to see i f regions of u n p r e d i c t a b i l i t y were located. As mentioned i n Chapter I I , the g l i t c h e s should show up as areas of r e l a t i v e l y large prediction error. After experimenting with the length of the f i l t e r 47 o p e r a t o r , and the time r e l a x a t i o n f a c t o r (ZTAU), r e s p e c t i v e v a l u e s of 100 and 5.0 were chosen. There i s no hard and f a s t r u l e f o r determining these v a l u e s , although Burg (1972) does make some recommendations f o r ZTAU. F i g u r e s #7 and #8 show the p r e d i c t i o n e r r o r t r a c e s f o r the high passed and band passed r e c o r d s . The f i r s t 100 p o i n t s of each r e c o r d r e p r e s e n t a warm-up p e r i o d where the f i l t e r c o e f f i c i e n t s are s e t to some i n i t i a l v a l u e . Areas B and C on F i g u r e #7 are a l s o warm- up s e c t i o n s , as the long l e n g t h of t h i s r e c o r d r e g u i r e d i t s d i v i s i o n i n t o t h r e e segments. As hoped f o r , the g l i t c h e s show up very w e l l i n both examples. B. P r e d i e t i o n Into Gaps In order to remove these g l i t c h e s from the r e c o r d s , a p r e d i c t i o n scheme based on one demonstrated by U l r y c h and C l a y t o n (1975) was used. I t i n v o l v e s the c a l c u l a t i o n of p r e d i c t i o n f i l t e r c o e f f i c i e n t s u s i n g the Burg a l g o r i t h m . However i n s t e a d of d e s i g n i n g the f i l t e r on data only i n one segment, t h i s gap p r e d i c t i o n method c r e a t e s a f i l t e r frcm segments of data t h a t are not connected. An averaging e f f e c t i s thus achieved. Time adaptive prediction error trace of high pass f i l t e r e d predicted data. Glitches show up as areas of u n p r e d i c t a b i l i t y . The time relaxation f a c t o r , ZTAO, equals 5.0 and the f i l t e r length equals 100 points. Areas A, B, and C are f i l t e r warmup sections, where no prediction occurs. 49 A k ^ T B 4- 12 hrs FIGURE 7 100 ygal 50 Time adaptive prediction error trace of band pass f i l t e r e d predicted data. The time relaxation f a c t o r , ZTAU, equals 5,0 and the f i l t e r length equals 100 points. 51 J L . 1 FIGURE 8 7 5 KS3'5 52 Once the f i l t e r has been designed, i t i s run over the data such that the prediction i s done into the gap frcm both sides. The predictions are then tapered with a ramp function and combined i n a f i n a l r e s u l t . In Figure #9, a harmonic si g n a l (0.05 and 0.06 Hz.) with 20% noise i s shown bj the s o l i d l i n e . The dashed l i n e shows the prediction of the middle 120 points of data using only information from the f i r s t and l a s t 40 points of data. The r e s u l t i s guite good and studies by Ulrych and Clayton (1975) have shown that the improvement i n the spectral resolution i s excellent. The r e s u l t s for deglitching both the high passed and band passed data sets can be seen i n Figures #10 and #11. The prediction error traces for these deglitched records are shown i n Figures #12 and #13. There i s a remarkable improvement i n the records, and t h i s can be seen v i s u a l l y i n both the deglitched and prediction error traces. c« Deglitching Ey Envelope An alternative method for deglitching was suggested by Wiggins (personal communication), and t h i s technigue was applied to the f i r s t 3000 points of the high passed record. Dividing the data by the sguare of i t s envelope was thought of as-a possible means of removing the high amplitudes of the noise bursts, even though the freguency content would not be altered. A modification of t h i s scheme was applied, i n that a function E 1 was created where Test case of the.gap prediction routine. A mixed harmonic of 0.05 Hz. and 0.06 Hz. with 20S added noise. The s o l i d l i n e represents the o r i g i n a l data, and the dotted l i n e represents the predicted data. The f i r s t 40 and l a s t 40 points were used to predict the middle 120 points. The f i l t e r length was chosen to be 35. 54 a oi FIGURE 9 FIGURES JO AND 1J. Deglitched high pass and band pass f i l t e r e d predicted record. Glitches were i d e n t i f i e d , then predicted out using a f i l t e r operator determined from the Burg algorithm. 56 I 12 h r s FIGURE 10 100 ygals 57 1IG0HES 12 AND _13_. Prediction error traces of both deglitched high passed and band passed predicted records. The areas of unpr e d i c t a b i l i t y have been larg e l y eliminated by the deglitching prccess. 59 -wW *<—r> vr>̂ *~ y B c -,/-v>v-V- 1 2 h r s 1 FIGURE 12 100 y g a l s FIGURE 13 7 5 ygals 61 E' = E/W.I. E = envelope of the process H.L. = some predetermined water l e v e l parameter. In t h i s case, the water l e v e l was not constant but was weighted e x p o n e n t i a l l y a t the beginning of the r e c o r d so as to c o r r e s p o n d 1 t o the high s i g n a l amplitude i n t h i s s e c t i o n of the data. I f any data f e l l above the water l e v e l , the f u n c t i o n E 1 was s e t egual to E/W.L. I f any data value f e l l below the water l e v e l , the f u n c t i o n E* was s e t egual to 1.0. i . e . E* = 1.0 i f E» < 1.0 E« = E/I.L. i f E» > 1.0 The r e c o r d was then d i v i d e d by (E*) 2 and t h e r e d u c t i o n i n g l i t c h amplitude was very s u c c e s s f u l (Figure #14). There i s a l s o a n o t i c e a b l e drop i n the n o i s e content i n the frequency domain (see Chapter I V ) . 5. POHEJ SPECIEJL M i l l S I S A» Peridogram and MEM Spectra In the c a l c u l a t i o n of power s p e c t r a f o r t h i s p r o j e c t , both c o n v e n t i o n a l and MEM s p e c t r a were used. In a l l cases, the periodograms were were not smoothed. The periodogram was Deglitched high passed predicted record. Deglitching was accomplished by dividing the record by the sguare of i t s envelope. The amplitudes of the gl i t c h e s have been greatly reduced, but the freguency content of the gl i t c h e s i s s t i l l evident. 63 1 FIGURE 14 1 0 0 /9 a l s 64 extremely useful i n helping to determine the optimum order of the prediction error f i l t e r needed, i n the c a l c u l a t i o n of the MEM spectra, as i t gave prior knowledge of the r e s u l t . In almost a l l cases, the HIM spectra f o r each free o s c i l l a t i o n mode were calculated several times using d i f f e r e n t f i l t e r orders and data lengths. In most cases the FPE gave no r e a l i s t i c i n d i c a t i o n of the f i l t e r order to be used. B. Complex Demodulation In the analysis of the spectra of i n d i v i d u a l free o s c i l l a t i o n modes, a technigue known as complex demodulation was used. This method involves a s h i f t i n g of the freguency of i n t e r e s t down to 1 u> = 0", then low pass f i l t e r i n g the data, and f i n a l l y resampling the data at a coarser i n t e r v a l . This method i s p a r t i c u l a r l y useful when calcul a t i n g MEM spectra, as fewer points are used in the c a l c u l a t i o n of the f i l t e r c o e f f i c i e n t s . Two sources of error are thus eliminated: 1. Roundoff errors leading to inaccuracies p a r t i c u l a r l y with harmonic terms. 2. Order of magnitude differences between the power of d i f f e r e n t freguency components which re s u l t s in numerical i n s t a b i l i t i e s (Olrych and Bishop 1975). The technigue of complex demodulation i s expressed 65 schematically below. Suppose a time series x(t) i s comprised of two freguency bands, and i t i s necessary that peak 'A* be i s o l a t e d . 0 8 ft - co - u><i-uiA ou-.o w 8 The complex time series C(t) •= (x(t),iO) i s multiplied by the factor exp ( - i ^ n at) , where cuA i s the angular freguency of peak *A* and a t i s the sampling i n t e r v a l . This m u l t i p l i c a t i o n s h i f t s peak »A' to <w=0, and the negative peak •A» to ou - -2«v». CO -LO " A The data can then be low pass f i l t e r e d so as to only retain the freguency band around ^ = 0. -co ui-.o Resampling i n the time domain can now be performed with no a l i a s i n g , and a stable MEM spectrum can be calculated for the 66 freguency component of interest. This above operation can be expressed mathematically. Childers and Pao (1972) show that i f a time series consists of a set of transient harmonics f(t ) = [e-^-TVV sinco 0(t^T)] e-<^~TVc = transient decay factor then multiplying f (t) by the demodulation factor Q-cto0-k gives f (t) = { e-<*"7>/r/2i [edu,°<*-T> - e - ^ o < £-r> ] e-iovij = e-H-T)/t/2i [ e - t a , e T - e-^to 0-t e C a 3° r], This new series f ( t ) has a 2 ooG freguency component which can be f i l t e r e d out to leave f ( t ) = ( e - t ^ V r /2i) | e - ^ r ) 67 Iv. RESQITS 1. FREE 0SCI1L&TICN MODE F R E Q U E N C I E S The deglitched and undeglitched high pass f i l t e r e d data sets were used in the examination of free o s c i l l a t i o n multiplets and the associated s p l i t t i n g parameters. Periodogram and MEM spectra were calculated using various data lengths and f i l t e r orders. The values f o r the s p l i t t i n g parameter * ft * were calculated using the f i r s t order s p l i t t i n g approximation given by Backus and Gilbe r t (1961), where 0^= cu0+ mfiJX "V s s p l i t peak freguency «3o= centre peak freguency m = s p l i t t i n g order number /S = s p l i t t i n g parameter rotation rate of the earth The following sections d e t a i l the measured freguency values and calculated s p l i t t i n g parameters for the fundamental spheroidal modes 0 S ^ , Q S 3 ' , Q S ^ , E S S , and E S 6 . Comparisons are also made between the spectra f o r undeglitched, deglitched, undemodulated, and unenhanced data sets. Remarks are alsc made about the e f f e c t of varying the f i l t e r order f o r 68 c a l c u l a t i o n of HEM spectra, and the higher resolution c a p a b i l i t y of MEM spectra as compared to periodogram spectra. In a l l cases,•the power spectra shown were calculated from the complex demodulated data sets. A. mode The periodogram spectra for the undeglitched and deglitched records are shown in Figure 15a and Figure 15b. The noise l e v e l s are indicated, and the decrease i n noise power was calculated to be 5 db. For the undeglitched spectrum (Figure 15a) i t i s only possible to i d e n t i f y the ±1 peaks, whereas the deglitched spectrum (Figure 15b) shows the ±1 peaks as well as the +2 peak. Figure 15c i s the periodogram power of the unenhanced (not predicted) high pass f i l t e r e d record (Figure 1). This spectrum shows lower sig n a l tc noise l e v e l s than 15a, as a r e s u l t of information loss at the beginning of the record. The advantage of predictive f i l t e r i n g i n the analysis of the record can be seen. Figure 15d i s shown just as a check of the complex demodulation technigue. For a periodogram there should be no difference between the spectrum calculated for the demodulated data and the one calculated for the undemodulated data. I f Figures 15a and 15d are compared, t h i s equivalence can be seen. Figures 15e and 15f show MEM spectra for the deglitched record with orders of 63 and 150 respectively. The 69 greater d e t a i l in 15f i s the r e s u l t of choosing a higher f i l t e r order. A possible -2 peak has also been i d e n t i f i e d on t h i s plot. Figures 15e and 15f show considerably l e s s noise, and better resolved peaks than Figure 15g, which i s the HEM spectrum (order=110) for the undeglitched data set. Interestingly enough, the MEM spectra, Figure 15h (order=45) , and Figure 15i (order=55) , f or the f i r s t 1200 points of the record, seem to show both the ±1 and ±2 peaks, whereas the corresponding periodogram (Figure 15j) shows only the ±1 peaks. This i s a good example of the increased resolution c a p a b i l i t y of the MEM spectral estimation method. For the c S a mode, the •,ar=01 peak cannot be observed, as i t corresponds to a nodal point located at los Angeles. Therefore, the ±1, and ±2 peaks are a l l that w i l l be observed in t h i s record. The freguency of the '0' peak i s calculated as the average of the ±1 peaks. Table 1 d e t a i l s the frequencies observed f o r the CS^ mode. TABLE JM Periodogram 5320 Points Ondeglitched 15a freg,, Jcy^hrJ. & +1 1.130 ± .001 0 1. 11*4 ± .001 -1 1.097 ± .001 .360 ± .007 .432 ± .007 70 Periodogram 5320 P o i n t s D e g l i t c h e d 15b f r e g . l£XZhLk & +2 1.159 ± .002 .528 ± .008 +1 1.129 ± .001 .336 ± .007 0 1.114 ± .001 -1 1.099 ± .001 .384 ± .007 MEM order=63 5320 P o i n t s D e g l i t c h e d 15e f r e g . J c y / h r l & +2 1.158 ± .002 .516 ± .008 +1 1.130 ± .001 .360 ± .007 0 1.115 ± .001 -1 1.099 ± .001 .384 ± .007 MEM order=J[50 5320 P o i n t s D e g l i t c h e d J 5 f freg_. Jcy/hr). g +2 1.157 ± .002 .504 ± .008 +1 1.130 ± .001 .360 ± .007 0 1.115 ± .001 -1 1.099 ± .001 .384 ± .007 -2 1.074 ± .002 .492 ± .008 MEM order=110 5320 P o i n t s Undeglitched 15g f r e g A Jcy2kEl ~ ~& +2 1.159 ± .001 .528 ± .005 +1 1.130 ± .001 .360 ± .007 0 1. 114 ± .001 -1 1.098 ± .001 .408 ± .007 MEM order=45 1200 P o i n t s Undeglitched JJ5h f reg.. SciZhLi. & +2 1.164 ± .002 .588 ± .008 +1 1.129 ± .002 .336 ± .011 0 1.115 ± .001 -1 1.100 ± .001 .360 ± .007 -2 1.073 ± .002 .504 ± .008 71 order=55 J200 P o i n t s Undeglitched J 5 i Ireg A J c y ^ h r l " +2 1.162 ± .002 .564 ± .008 + 1 1. 134 ± .002 .455 ± .011 0 1.117 ± ^001 -1 1.100 ± .001 .360 ± .007 -2 1.078 ± .002 .444 ± .008 Periodogram 1200 P o i n t s Undeglitched 15j "~ l £ S g A i c y / h r l _ £ _ +1 1.131 ± .001 .384 ± .007 0 1.115 ± .001 -1 1.099 ± .002 .384 ± .011 4Z§£§gg Values from MEM s p e c t r a freg j . i S l i h r l - 0 - + 2 1.160 ± .003 .540 ± .035 + 1 1.131 ± .002 .374 ± .046 0 1. 115 ± .001 -1 1.099 ± .001 .379 ± .020 -2 1.075 ± .003 .480 ± .032 Average @ o v e r a l l = .439 ± .081 The e r r o r s given f o r the average val u e s are standard d e v i a t i o n s . G i l b e r t and Backus(1965) gave a value of 1.116 c y c l e s per hour f o r the freguency of the 0 m o d e , and a value of 0.397 f o r the s p l i t t i n g parameter. Derr (1969) gave valu e s of 1.114 c y c l e s per hour and 0.397. These compare t o the average values of 1,115 c y c l e s per hour and 0.439 obtained i n t h i s p r o j e c t . However, the values f o r * 01 t h a t were c a l c u l a t e d f o r each peak were net c o n s t a n t . As a i g g i n s and M i l l e r (1972) noted, t h e r e appears to be a dependence of the s p l i t t i n g parameter « ^ ' on the s p l i t t i n g order number 'm*. The l a t t e r a. Periodogram 5320 points b. Periodogram 5320 points c. Periodogram 5320 points d. Periodogram 5320 points €. HEM order=63 5320 points f. MEM order=150 5320 points g. MEM crder=110 5320 points h. MEM crder=15 1200 points i . MEM order=55 1200 points j . Periodogram 1200 points Undeglitched D e g l i t c h e d Unenhanced Undemodulated D e g l i t c h e d D e g l i t c h e d Undeglitched Undeglitched Undeglitched Undeglitched 73 •x. ZD CC f— CJ". lUc.-Q_ LD a LLJ LO •—ITJ-5 ^ >z ce o r I.oi T r 1.05 1.09 1.13 FREQUENCY IN CY/HR FIGURE 15a 1.17 ZD OC 3̂<o LJo" Cu to O UJ in I—tit <£° o r 1.01 -i 1— 1 r 1.05 1.09 1.13 1.17 FREQUENCY IN CY/HR FIGURE 15b 1.21 I — 1.01 T 1 r 1.05 1.09 1.13 FREQUENCY IN CY/HR i FIGURE 15d 1.21 75 03 3= ZD CC CJ" 3 UJo" D_ CO Q UJ CO >— d° z: cc o z: o o" 1.0) -i r—— 1 — i — 1.D5 1.D9 ] . ) 3 I . 1 7 FREQUENCY IN CY/HR FIGURE 15e 1 . 2 1 U J o " Q_ to a UJ co l—IT z: cc: o CM a a" -2 -1 - 1 1.05 — ! — 1 .09 FREQUENCY FIGURE 15f IN CY/HR 1.21 76 77 78 authors gave values for » of 0.408, 0.410, and 0.522 for the -1, +1, and +2 peaks, which compare reasonably well to the values of 0.379, 0.374, and 0.540 obtained here. B. „ S„ mode o - 3 The periodogram power spectra, for the undeglitched and deglitched high passed records are shown i n figures 16a and 16b. As observed for the 0 m o d e , there i s a very noticeable reduction in the noise content for the deglitched record. Figure 16c i s the MEM spectrum (order=130) for the undeglitched data set, and should be compared to Figures 16d, 16e, and 16f, which are the MEM spectra for the corresponding deglitched data (orders 70, 130, and 150). The e f f e c t of varying the length of the calculated prediction operator can be seen i n these l a s t three fig u r e s , where the spectral resolution increases as the order of the f i l t e r increases. However, going to too high an order w i l l r e s u l t i n an unstable spectrum. This i s starting to happen with the 150 point MEM spectrum, as the noise i s beginning to take on a very spiky appearance. Deglitching the record by dividing by the sguare of the envelope i s demonstrated in the next two f i g u r e s , 16g and 16h. Figure 16g shows the periodogram power of 3000 points of undeglitched data, and 16h shows the periodogram power f c r the deglitched 3000 point segment. There i s a noticeable reduction i n the noise content using t h i s method of removing 79 t h e g l i t c h e s , b u t i t i s n o t as e f f e c t i v e a s t h e p r e d i c t i v e method o f d e g l i t c h i n g . Ey f a r , t h e b e s t s p e c t r u m o f t h e C S 3 mode i s d i s p l a y e d i n f i g u r e 1 6 f , which shows t h e s p l i t t i n g o f t h e major peak i n t o ± 1 , ± 2 , and a p o s s i b l e -3 peak. T a b l e #2 d e t a i l s some o f t h e o b s e r v e d f r e g u e n c y m u l t i p l e t s . T f l E L J #2 Periodogram 5320 Points Deglitched J6b freg. J c y / h r l & +2 1.703 ± .001 .204 ± .005 •1 1.695 ± .001 .216 ± .007 0 1.686 ± .001 -1 1.680 ± .001 .144 ± .007 -2 1.670 ± .001 .192 ± .005 order=V30 5320 P o i n t s D e g l i t c h e d J6c f r e g ^ SZlZhEL ~& +2 1.704 ± .001 .216 ± .005 +1 1.694 ± .001 .192 ± .007 0 1.685 ± .001 -1 1.680 ± .001 .144 ± .007 -2 1.667 ± .001 .228 ± .005 MEM order=70 5320 Points Deglitched J6d Jreg.. J c y / h r l & +2 1.705 ± .001 .228 ± .005 •1 1.691 ± .002 .120 ± .011 0 1.687 ± .001 -1 1.682 ± .001 .096 ± .007 -2 1.669 ± .001 .204 ± .005 80 MEM order=J30 5320 P o i n t s D e g l i t c h e d J6e +2 1.704 ± .001 .216 ±' .005 +1 1.693 ± .001 . 168 ±- .007 0 1.686 ± .001 -1 1.680 ± .001 . 144 ± .007 -2 1.670 ± .001 .192 ± .005 -3 1.665 ± .001 . 168 ± .004 MEM o r d e r ^ l 5 0 5320 P o i n t s D e g l i t c h e d Ireg. J c j r / h r l +2 1.704 ± .001 .216 ± .005 + 1 1.694 ± .001 .192 ± .007 0 1.686 ± .001 -1 1.681 ± .001 . 120 ± .007 -2 1.671 ± .001 .180 ± .005 -3 1.665 ± .001 .168 ± .004 Periodogram 5320 P o i n t s Undeglitched freg^. J c ^ h r l 0 +2 1.705 ± .001 .228 ± .005 +1 1.694 ± .001 .168 ± .007 0 1.687 ± .001 -1 1.682 ± .001 .096 ± .007 -2 1.673 ± .001 . 156 ±- .005 -3 1.667 ± .002 . 152 ± .006 Average Values fron i MEM s p e c t r a f£§3i JSlZferl +2 1.704 ± .001 .220 ±- .007 •1 1.693 i .001 .168 ± .033 0 1.686 ± .001 -1 1.681 ± .001 .126 ± .024 -2 1.669 ± .002 .201 ± .020 -3 1.665 ± .001 . 168 ±- .009 Average 0 o v e r a l l = . 177 ± .040 G i l b e r t and Backus (1965) gave a value of 1.688 c y c l e s per hour f o r the freguency of D S 3 and a value cf .1839 f o r 81 IIGUEE 16 a. Periodogram 5320 points Undeglitched b. Periodogram 5320 points Deglitched c MEM order=130 5320 points Undeglitched d. MEM order=70 5320 points Deglitched e. MEM order=130 5320 points Deglitched f. MEM order=150 5320 points Deglitched 9* Periodogram 3000 points Undeglitched h. Periodogram 3000 points Deglitched  83 03 ZZi CC UJa" CL. CO O UJ CO >r. cc o o" 1 . 5 8 1 . 6 2 1 . 6 6 1 . 7 N 1 . 7 4 FREQUENCY IN CY/HR FIGURE 16c 1 . 7 8 TI zz> cc I— O < o U J o ' Q_ CO O UJ CO »—irt 5?° CC o V -2 -1 -1 +2 I • 1 r - J — i 1 . 5 B 1 . 6 2 1 . 6 6 1 . 7 FREQUENCY IN CY/HR FIGURE 16d 1 . 7 4 1 . 7 8 84 or L3<°. UJo" cu cn a UJ in e x a : o 1.58 r - 1 1 62 1.66 1 -7 FREQUENCY IN CY/HR FIGURE 16e 1.74 1.78 11 111,1 , , FIGURE 16f l i s 1 62 1 .66 1.7 1.74 1.78 1 FREQUENCY IN CY/HR  86 the s p l i t t i n g parameter. Derr (1969) gave values of 1.687 cycles per hour and .187. These compare with the average values of 1.686 and .177 obtained i n t h i s project, k dependence of * ^ * on »m' i s again observed. C. S„ mode O-Cf. ; The periodogram power spectra for the undeglitched and deglitched high passed record are shown in Figures 17a and 17b. S p l i t t i n g can be observed in both, and although i d e n t i f i c a t i o n of the number of the s p l i t peaks i s d i f f i c u l t , a f a i r l y good explanation of the peaks has been made. Figure 17c i s the periodogram spectrum of the unenhanced data, and i f t h i s i s compared with Figure 17a, the advantage of extending the data set before tapering can be seen. Figures 17d and 17e are MEM spectra of the undeglitched (order=130) and the deglitched (order=150) records. Deglitching by dividing by the sguare of the envelope of the record i s i l l u s t r a t e d in Figures 17f and 17g, where 17f i s the periodogram power cf 3000 points of undeglitched data, and 17g i s the power of the corresponding deglitched data. There i s however, very l i t t l e difference between these two spectra, i n d i c a t i n g that t h i s method of removing the noise g l i t c h e s i s freguency dependent. Table #3 gives values for some of the observed freguency s p l i t s . + 1 0 5320 g o i n t s D e g l i t c h e d f r e g . J c y / h r l ~& 17b + 4 2.359 +2 2.343 2.335 2.328 -1 2.321 -3 2.309 .002 .002 .001 .001 .001 .002 186 ± .005 180 ± .008 168 ± .007 168 ± 152 ± .007 .006 MEM order=130 5320 P o i n t s ~ Ireg. J c y ^ h r l Ondeglitched 17d +4 2.359 ± .001 . 186 ± .003 •2 2.341 ± .001 .156 ± .005 +1 2.335 ± .001 .168 ± .007 0 2.328 ± .001 -1 2.321 ± .001 . 168 ± .007 MEM order5l50 5320 P o i n t s D e g l i t c h e d f reg.. j c y / h r j +4 2.359 ± .001 .186 ± .003 +2 2.341 ± .001 .156 ± .005 +1 2.333 ± .001 . 120 .007 0 2.328 ± .001 -1 2.321 ± .001 .168 ±- .007 -3 2.310 ± .001 .144 ± .004 Average Values from MEM s p e c t r a f reg.. SQlZhLl +4 2.359 ± .001 .186 ± .003 +2 2.341 ± .001 .156 ± .005 +1 2.334 ± .001 .144 ± .034 0 2.328 ± .001 -1 2.321 ± .001 .168 ± .007 -3 2.310 ± .001 . 144 ± .004 Average 0 o v e r a l l = .161 ± 0.021 • 17e The value f o r the Q S„ ce n t r e freguency given by 88 a. Periodogram 5320 points Undeglitched b. Periodogram 5320 points Deglitched c * Periodogram 5320 points Unenhanced d. HEM order»130 5320 points Undeglitched e. MEM order=150 5320 points Deglitched £. Periodogram 3000 points Undeglitched g. Periodogram 3000 points Deglitched 89 90 91  93 G i l b e r t and Backus (1965) i s 2.335 cycles per hour and Derr (1969) gave a value of 2.327 cycles per hour. Derr also gave a value of .103 for the s p l i t t i n g parameter which i s smaller than the value of .161 obtained here. D. S„ mode o-$ The periodogram poser spectra for the undeglitched and deglitched high passed data are shown i n Figures 18a and 18b. There i s a noticeable improvement i n the signal to noise r a t i o , but i d e n t i f i c a t i o n of in d i v i d u a l peaks remains a d i f f i c u l t matter. In a l l the spectra observed, two major symmetrical s p l i t s sere observed, and i t Is believed that these are blended +1, +2, +3 and -1, -2,,and -3 spectral l i n e s , which are not properly resolved. Therefore, for the 0 § G mode, only the *m=0' freguency was determined, and i t was taken as the midpoint of the two major peaks on each spectrum examined. Figure 18c i s the MEM (order=130) spectrum of the undeglitched data set and Figure 18d shows the corresponding MEM (crder=130) spectrum for the deglitched data. The noise reduction due to deglitching i s evident in the MEM spectra as well as the periodograms. Figure 18e i s the periodogram power for a 760 point undeglitched data segment. The observed frequencies f o r the 0Sg. mode are l i s t e d in Table #4. 94 H S U i i i £sriodogram 5320 Points Undeglitched 18a f reg.. lUllhLl 0 3.028 ± .002 £S£io^£3ilS 5320 Points Deglitched 18b jE£Sai l£lZhLl 0 3.027 ± .002 I I S order=J30 5320 Points Undeglitched 18c f£§g. i c y ^ h r l 0 3.029 ± .002 MEM order=130 5320 Points Deglitched 1.8d freg.. Jcvyhr), 0 3.028 ± .002 Periodogram 760 Points Undeglitched 18e ££§3.2. J c y / h r l 0 3.028 ± .002 Average Values from MEM spectra f r eg.. J c y ^ h r l 0 3.029 ± .002 This average value for the centre freguency of the 0 m o d e compares with that of 3.034 cycles per hour given by G i l b e r t and Backus (1965) and 3.027 cycles per hour given by Derr (1969). 95 IJGUEE 18 a. Periodogram 5320 po i n t s Undeglitched b. Periodogram 5320 p o i n t s D e g l i t c h e d c. HEH order=130 5320 p o i n t s Undeglitched a . HEH. order=130 5320 p o i n t s D e g l i t c h e d e. Perioaogram 760 po i n t s Undeglitched 9 6 97 98 99 E. Q S 6 mode figures 19a and 19b show the periodogram power spectra of the undeglitched and deglitched high passed data. The noise reduction due to removing the glitches i s not as noticeable as for other free o s c i l l a t i o n modes. There i s some evidence f o r s p l i t t i n g i n the deglitched spectrum, but the s p l i t t i n g number cannot be determined. Figures 19c and 19d are MEM power spectra for the undeglitched (order=110) and deglitched (order=55) records. The high resolution c a p a b i l i t y of MEM i s i l l u s t r a t e d in Figures 19e and 19f. The f i r s t i s a plot of the periodogram power of a 760 point undeglitched data segment, and the second i s an MEM (order=10) spectrum for the same data segment. The increased resolution i s guite obvious. Table #5 gives the observed values for the QS f e center freguency. TABLE #5 Periodogram 5320 Points Ondeglitched J9a fre g ^ J c y / h r i 0 3.739 ± .004 £e£i°dojjr.am. 532.0 Points Deglitched 19b ~ treg."jc3i/hrl ~ 0 3.739 ± .004 100 MEM orderfjHO 5320 Points OndeoJ-itched 19c fregT J c y / h r l 0 3.74 0 ± .001 MEM ord«r=55 5320 Points Deglitched 19d 0 3.740 ± .001 Periodogram 760 Points Undeglitched J9e fregj, J c y / h r l 0 3.738 ± .003 MEM ordexfJHD 760 Points Undeglitched 19f fl§g. J c y ^ h r l 0 3.743 ± .001 MSiaaS Values from JEM spectra freg.. J c y ^ h r l 0 3.741 ± .002 This average value compares with the value of 3.749 cycles per hour given by Gilbert and Backus (1965) and 3.735 cycles per hour given by Derr (1969). 2 i CORE UNDERTONES A further goal of t h i s thesis was to see i f the deglitching of the band passed Alaska record could r e s u l t i n the i d e n t i f i c a t i o n of any free o s c i l l a t i o n ' s of the earth's inner core. The period range studied was r e s t r i c t e d to between 0.833 and 8.333 hours. 8.333 hours was chosen as the lower cutoff period, as i t was f e l t that the overpowering 101 I2GURJ 19 a. Periodograo 5320 points Undeglitched b. Perioaogram 5320 points Deglitched c. HEM oraer=110 5320 points Unaeglitchea a. MEM order=55 5320 points Deglitched Perioaogram 760 points Undeglitched f . MEM order=10 760 points Undeglitched 102 FIGURE 19b FREQUENCY IN CY/HR 103  105 presence of the d i u r n a l and the s e m i d i u r n a l t i d e s would contaminate the r e c o r d at l o n g e r p e r i o d s . The p r e d i c t e d i n t e r p o l a t e d r e c o r d was bandpassed and c o n v e n t i o n a l , unsmoothed periodogram power s p e c t r a were c a l c u l a t e d both before and a f t e r removing the g l i t c h e s . F i g u r e s 20 and 21 show the s p e c t r a f o r the u n d e g l i t c h e d and d e g l i t c h e d r e c o r d s r e s p e c t i v e l y . The l o c a t i o n of the S mode i s marked at the upper end of the freguency s c a l e . Removing the n o i s e b u r s t s i n the time domain r e s u l t e d i n an amazing r e d u c t i o n i n the noise power i n the frequency domain e s p e c i a l l y , at the lower f r e g u e n c i e s . There are n o t i c e a b l e peaks i n both s p e c t r a , but they become much more obvious i n the spectrum of the d e g l i t c h e d r e c o r d . The p e r i o d s of these unexpected peaks are exact i n t e g e r m u l t i p l e s of the b a s i c l u n a r and s o l a r t i d e , and extend up to about 0.65 c y c l e s per hour. An i n i t i a l r e a c t i o n was t h a t these peaks were some t i d a l overtone e f f e c t , such as shallow water t i d e s found i n harbours and b a s i n s . However, i t i s very d o u b t f u l t h a t these high order ocean t i d e ? e f f e c t s w i l l be d e t e c t e d by a l a n d based gravimeter. Because of the high amplitudes of these peaks r e l a t i v e to 0S% , i t i s even more u n l i k e l y t h a t they r e p r e s e n t core undertone o s c i l l a t i o n s . The most reasonable e x p l a n a t i o n i s t h a t they are due to e i t h e r instrument non- l i n e a r i t i e s , or barometric pressure e f f e c t s , such as those r e p o r t e d by Warburton et a l . (1975). The noise power f o r the u n d e g l i t c h e d spectrum goes from 106 FIGURES 20 AND 21 Periodogram power s p e c t r a of both u n d e g l i t c h e d (Figure 20) and d e g l i t c h e d (Figure 21) band pass f i l t e r e d p r e d i c t e d r e c o r d . The f i l t e r c u t o f f s are at 0.12 and 1.20 c y c l e s per hour. The l o c a t i o n s of the 0S^ f r e e o s c i l l a t i o n mode and the t e r - d i u r n a l t i d e (M3) are marked.  o o lO O M o Q. a cc z: cc o ID o oo 0.125 -1 1 — 0.225 0.325 —1 0.425 0.525 0.625 FREQUENCY 0.725 (CY/HR) 0.825 0.925 -1 1.025 1.125 1.225 109 40.0 x 10*-at 1.125 Hz. to 5 x 10 6 at .725 Hz. The noise power for the deglitched spectrum goes from 12 x 10* at .125 Hz. to 2 x 10» at .725 Hz. This drop i n noise power i s 5.2 db. at .125 Hz. and 4.0 db. at .725 Hz. 110 V. SUMHAJY A t i d a l gravimeter recording from the 1964 Alaska earthquake was subjected to various data processing techniques.in:an attempt to improve the quality of the record. An e f f o r t was also made to try and locate core undertone o s c i l l a t i o n s . The findings of t h i s project can be summarized as follows: 1. A time adaptive prediction error trace provides a useful t o o l for i d e n t i f y i n g and locating noise bursts i n time series 2. The method of removing the unwanted noise by predictive deglitching, r e s u l t s i n an excellent improvement of the signal to noise r a t i o of t h i s data set. 3. The maximum entropy method of spectal estimation gives some excellent results for the power spectra of various free o s c i l l a t i o n modes. Unfortunately, estimation of the length of the f i l t e r operator reguired i n the spectral c a l c u l a t i o n remains a d i f f i c u l t task. Prior knowledge of the spectrum through a conventional spectral technique, i s very useful in aidinq i n the determination of the optimum f i l t e r order. 111 In several instances, the maximum entropy method of spectral estimation gave good r e s u l t s f o r the s p l i t t i n g of various free o s c i l l a t i o n modes. Certainly more s p l i t peaks were observed using MEM as opposed to the periodogram approach. The central freguencies for the 0 S a , 0S 3; oS^ , 0Sg-, and the QSfe modes are 1. 115, 1.686, 2.328, 3.029, and 3.741 cycles per hour respectively. The s p l i t t i n g parameters for the 0 S a , CS^, and Q S f modes are .439, .177, and .161. No evidence f o r core undertones was discovered during the analysis of the band pass f i l t e r e d record. The freguency band below 0,12 cycles per hour i s strongly dominated by t i d a l freguencies, and the band from 0,12 to 0.65 cycles per hour i s dominated by either instrument n o n - l i n e a r i t i e s or barometric pressure e f f e c t s . 112 REFERENCES Akaike, H. 1 9 6 9 . Eower spectrum estimation through autoregressive model f i t t i n g . Ann. Instjs, Statist.. Vol. 2 1 , p. 4 0 7 . Akaike, H. 1 9 7 0 . S t a t i s t i c a l predictor i d e n t i f i c a t i o n . JLSSJL l i s t i Statist,. J a t h A f Vol. 2 2 , p. 2 0 3 . Alsop, A. L i , G. H. Sutton, and M. Ewing. 1 9 6 1 . Measurement of Q for very long period free o s c i l l a t i o n s . JAGJ_RJ_, Vol. 6 6 , p. 6 3 1 . Alterman, Z. S., Y, Eyal, and A. fl. Merzer. 1 9 7 4 . Cn free o s c i l l a t i o n s of the earth. Geophysical Surveys, Vol. 1, p. 4 0 9 . Backus, G. and F . Gilbert. 1 9 6 1 . The rota t i o n a l s p l i t t i n g of the free o s c i l l a t i o n s of the earth. Froc A Nat A Acad.. Sci.. OJ.SJ.A_. , Vol. 4 7 , p. 3 6 2 . Backus, G. and F. Gilb e r t . 1 9 6 7 . Numerical applications of a formalism f o r geophysical inverse problems. Geophysic.. Jour. Roy. astr. Soc.., Vol. 1 3 , p. 2 7 7 7 ~ ~ ~ ~ Benioff, H. F., B. Gutenberg, and C. F. Richter, 1 9 5 4 . Progress report, Seismological Laboratory, C a l i f o r n i a I n s t i t u t e of Technology. Transj. Anu Geophys,. On., Vol. 3 5 , p. 9 7 9 . , Benioff, H. F., F. Press, and S. Smith. 1 9 6 1 . Excitation of the free o s c i l l a t i o n s of the earth. Jjj.Gj.Ri, Vol. 6 6 , p. 6 0 5 . Bogert, B. P. 1 9 6 1 . An observation of free o s c i l l a t i o n s of the earth. J.G.R., Vol. 6 6 , p. 6 4 3 . Bolt, B. A., and R. G. Currie. 1 9 7 5 . Maximum entropy estimates of earth to r s i o n a l eigenperiods from 1 9 6 0 Trieste data. Geophysic. Jour. Roy A astr. Soc.,, Vol. 4 0 , p. 1077. " Burg, J . P. and D. C. Riley. 1 9 7 1 . Time and space adaptive deconvolution f i l t e r . S.E.G., Preprint. Burg, J. P. 1 9 7 2 . The relat i o n s h i p between maximum entropy and maximum l i k e l i h o o d spectra. Ge oj_ hj_si cs , Vol. 3 7 , p. 3 7 5 . Burg, J . P. 1 9 7 5 . Maximum entropy spectral analysis. Unpublished Ph.D. Thesis, Stanford. 113 Childers, B. G. and Min-Tai Pao. 1972. Complex demodulation fox transient wavelet detection and extraction. I A E X E A I A on Audio and Ii§£troacoustics, Vol. AU-20, p. 295 Crossley, D. J . 1975. Core undertones with rotation. GeoDhysic. Jour.. Boy. ast r . Soc,, Vol. 42, p. 200." Derr, J . S. 1969. free o s c i l l a t i o n observations through 1968. B u l l A Seism A Soc. Am., Vol. 59, p. 2079. Dziewonski, A. M. and F. G i l b e r t . 1973. Observations of normal modes from 84 recordings of the Alaskan earthguake of 1964 March 28. Geo£hysie A Jour^ Eoy x a s t r . Soc,, Vol. 35, p. 401. Gersch, H. 1970. Estimation of the autoregressive parameters of a mixed autoregressive moving average time s e r i e s . l i J i l i J . Trans A on Automatic Control, Vol. AC-18, ~p~ ~267. G i l b e r t , F. and G. Backus. 1965. The r o t a t i o n a l s p l i t t i n g of the free o s c i l l a t i o n s of the earth. Beyiews of Geg£hysics, Vol. 3, p. 1. Haddon, B. A. and K. E. Bullen. 1967., Derivation of an earth model from free o s c i l l a t i o n data. Proc x Hat. Acad. S c i . JJ AS.A A, Vol. 58, p. 846." Jackson, B. V. and L. B. S l i c h t e r . 1974. The residual daily earth t i d e s . JLs.G.B., Vol. 79, p. 1711. Lacoss, B. T. 1971. Data adaptive spectral analysis methods. GeofhjjsJ.cs, Vol. 36, p. 661. Lamb, H. 1882. Proc. Math. Soc.., Vol. 13, p. 189. Love, A. E. H. 1911. Some problgms of ggodynamics. Cambridge University Press, Cambridge. Ness, N. F;, J . C. Harrison, and L. B. S l i c h t e r . 1961. Observations of the free o s c i l l a t i o n s of the earth. J.iGAB., Vol. 66, p. 62. Press, F. 1968. Earth models obtained by Monte Carlo inversion. j3.G AB., Vol. 73, p. 5223. Shannon, C. E. 1948. A mathematical theory of communication. B e l l . Syst A Tech* Jour., Vol. 27, p. 379. 114 S l i c h t e r , I. B. 1961. The fundamental free mode of the earth's inner core. Prcc A Nat.. Acad.. Sci.. Û Ŝ A.., Vol. 47, p. 186. ~~ S l i c h t e r , 1. B. 1967a. Free o s c i l l a t i o n s of the earth. Dictionary of Geophysics, ed. S. K. Runcorn, Pergamon Press, Oxford. S l i c h t e r , L. B. 1967b. Spherical o s c i l l a t i o n s of the earth. Geophysic... Jour. Roy., a s t r x Soc., Vol. 14, p. ~ 1 7 l T Smith, M. L. 1974. The normal modes of a rotat i n g e l l i p t i c a l earth. Ph.D. Thesis, Princeton. Ulrych, T. J . , D. E. Smylie, 0. Jensen> and G. K. C. Clarke. 1973. Predictive f i l t e r i n g and smoothing of short records by using maximum entropy. JJLGJUIA# Vol. 78, p. 4959. Ulrych, T. J . , and 1. N. , Bishop. 1975. Maximum entropy spectral analysis and autoregressive decomposition. Reviews of Geophysics and Space Physics, Vol. 13, p. 183. Ulrych, T. J . , and B. 8. Clayton. 1975. Time series modelling and maximum entropy. Preprint. VandenBos, A. 1971. Alternative interpretation of maximum entropy spectral analysis. IJJSAJSAIIA Trans A Infcrm.. Theory JCorresplj., Vol. I T - l 7, p. 493. Harburton, B. J . , C. Beaumont, and J. M. Goodkind. 1975. The eff e c t of ocean tide leading on tides of the s o l i d earth observed with the superconducting gravimeter. Geophysics. Jourj. Bgy x astr.. Soc.__, Vol. 43, pT~ 707." Weiner, N. 1930. Generalized harmonic analysis. Acta A Math.. Vol. 55, p. 117. Wiggins, R. A., and S. P. M i l l e r . 1972. New noise reduction technigue applied to long period o s c i l l a t i o n s from the Alaskan earthquake. B u l l . Seism. Soc. Am.., Vol. 62, p. 471. 115 APPENDIX 1 Fortran source l i s t i n g of the time adaptive maximum entropy deconvolution scheme. i 116 SUBROUTINE B A F L ( L O U T , C X , I F L G A P , L C N • N W A R M , I S T R T , Z T A U , X T A U , N Z E R O ) C C SUBROUTINE BAFL C THE BURG ADAPT IVE F I L T E R LADDER:AN ADAPT IVE OR T IME C VARYING F IXED LEAD PRED ICT ION ERROR PROCESSOR C ADJUSTMENT OF EACH REFLECT ION C O E F F I C I E N T IS MADE EVERY C JUMP STATE ATTEMPTING TO M IN IM IZE THE STAGE OUTPUT POWER C INPUTS C X ( L X ) = I N I T I A L DATA C LCN=LAST NON-ZERO REFLECT ION COEFF . C IFLGAP=NUMBER OF GAPS BETWEEN F I L T E R C O E F F I C I E N T S C SETT ING IFLGAP=1 DOES NOT GAP THE SPECTRUM AND T R I E S C TO OPERATE ON THE ENT IRE SPECTRUM. C NWARM= CURAT ION OF STATIONARY CYCLE C ISTRT=START CF STATIONARY GATS C ZTAU=TEMPORAL RELAXAT ION FACTOR COMMON / B L K 3 / X U 7 7 4 ) C XTAU=SPAT IAL RELAXAT ION FACTOR C OUTPUTS C C XtLOlT)=FORWARD. ERROR PRED ICT ION TRACE C C(LCN)=FORWARD STATE VECTOR C B(LCN)=BACKWARD STATE VECTOR C C ( L C N ) = R E F L E C T I O N C O E F F . AT EACH STAGE C CX=R EFLECT ION C O E F F . INTEGRATED IN SPACE AND TIME C DEN(LCN) = STAGE AUTOPOWER C NUM(LCN) = STAGE CROSSPOWER C C C A L L I N G BAFL F IRST SETS UP THE LOOPING AND PASS ING C ARRAYS FOR THE PART ICULAR PROBLEM AS S P E C I F I E D BY C L C N G I F L G A P . THEN IT COMPUTES A SHORT!LENGTH=NWARM) EST IMATE C OF THE REFLECT ION C O E F F . S ER I E S IN ORDER TO START THE C ADAPTION OUT WITH SOME REASONABLE NUMBERS C THEN IT LOADS UP THE CX ARRAY WITH THE I N I T I A L VALUES C AND PASSES INTO ENTRY ' B A F L G O * C THE USUAL ENTRY IS BAFL GO WHICH F I R ST I N I T I A L I Z E S THE C BACKWARD ARRAY THEN PASSES TO THE MAIN ROUTINE C THE CX SER IES IS UPDATED EVERY I F L G A P DATA POINTS AND IN THE C INTERMEDIATE STEPS THE OUTPUT ARE INTERPOLATED OF PROCESSED C AS THOUGH THE R .C . WERE STATIONARY C WRITTEN BY J . P . BURG REAL N U M ( 2 0 O ) , D E N { 2 0 0 ) , B { 2 0 0 ) , F ( 2 0 0 ) t C ( 2 0 0 ) , E M ( 2 0 0 0 ) , * E P ( 2 0 0 C ) , C X ( L C N , L O U T ) , A ( 2 0 0 ) OAT A A , B , C / 6 0 0 * 0 . / FTEST=1 . LCNP1=LCN+1 LCNP2=LCN+2 I F G M 1 = I F L G A P - 1 L B S P = L O U T - I F L G A P NZERP1=NZER0+1 NEND=LCN-NZERP1 DO 80 K=LBSP,LOUT 80 E P ( K ) = 0 . 117 10 c 12 11 33 e c c c- c c 1 0 0 0 c c 1 0 1 0 BEGIN STATIONARY WARMUP DO 10 I=1,NWARM E M U ) = X ( I * I F L G A P + I S T R T ) E P ( I ) = X ( I * I F L G A P + I S T R T ) DO 11 J =2 , LCNP1 D E N ( J ) = 0 . NUM(J ) = Oc DO 12 I=J,NWARM D E N ( J ) = D E N ( J ) + E P ( I ) * 5 P { I ) + E M ( I - J + i ) * E M ( I - J + l ) N U M U ) = N U M U ) + EP ( I )*EM( I-J + l ) DIV=NWARM-J+1 NUM{J) = N L M ( J ) / D I V D E N ( J ) = D E N(J ) / D I V C U ) = - 2 . * N U M l J ) / D E N ( J ) DO 11 I=J,NWARM E P I = £ P ( I > E P ( I ) = E P I + C U J * E M { I - J + l ) KS= I FLGAP*LCN+1 EM( I-J+1 ) = E M ( I - J + l ) + C ( J ) * E P I DO 8 J=1 ,LCN DO 8 K=1,KS C X ( J , K ) = C(J+1 ) DO 33 J=1 ,LCN DEN { J ) = D E N U + l ) N U M U ) = N U M U + l ) END WARM UP CYCLE SET RELAXAT ION TIMES D L X = E X P ( - 1 . / X T A U ) D L = E X P ( - 1 . / Z T A U ) DRX=1. -DLX DR= l o -DL ENTRY BAFLGO -USUAL ENTRY- I N I T I A L I Z E BACKWARD VECTOR B ( l ) = X ( K S - I F L G A P ) A ( l ) = l . DO 1000 J = 2 , L C N B(J) = X ( K S - J * I F I G A P > A(J) = 0 . 0 DO 1000 I=2,J A ( I ) = A ( I ) + C X ( J , K S ) * A { J - I + l ) 3(J) = B U ) + A C I ) * X ( K S + { I-J + 1 ) * I F L G A P ) BEGIN MAIN LOOP DO 5000 K = K S t l O U T , I F L G A P Z=X(K) DO 1010 J = l t N Z E R P l F I J ) = Z DO 2000 J=NZERP1 , LCN DEN { J ) = ( F ( J ) * * 2 + B ( J ) * * 2 ) * D R + D £ N ( J ) * D L N U M ( J ) = F{J ) * B (J)*DR+NUM{J)*DL I F ( F T E S T . L E . l . l ) CX(J,K)=-2.*NUM{J)/DEN(J) C X ( J t K ) = -2o*DRX*NUM{JJ/DEN(J}+CX{J,K)*DLX 200C F(J+1) = F { J ) + C X ( J t K ) * B { J ) X ( K ) = F ( L C N P l ) DO 3000 JR = 11NEND J=LCN - JR 3000 B l J + 1 J = B ( J ) + C X ( J t K ) * F { J) IF(NZERO.EG.O) GO TQ 5000 DO 4000 JR=1,NZER0 J=NZERP1~JR 4000 B C J + l ) = B(J) 5000 B(1)=Z C C END OF MAIN LOOP C C NOW GO BACK TO FILL IN THE GAPS IF( I F L G A P . E Q . l ) GO TO 9050 DO 9000 L = l t I F G M l KSTART=KS+L DO 9000 K=KSTART»LOUT,IFLGAP KC=K-L Z=X(K) DO 6000 J=l,NZERPl 6000 F ( J ) = Z DO 7000 J=NZERP1,LCN 70C0 F(J+1) = F< J)+CX( J,KC)*B('J) XlK) = F ( L C N P l ) DO 8000 JR=1»NEND J=LCN-JR 8000 B ( J+1) = B ( J )+CX(J,KC)*F(J) IF(NZERO.EC.O) GC. TO 9000 DO 8050 JR=1,NZER0 J=NZERP1 -JR 8050 B ( J + l ) = B ( J ) 9000 B ( l ) = Z 9050 FTEST=FTEST+1. RETURN END

Cite

Citation Scheme:

    

Usage Statistics

Country Views Downloads
China 19 18
United States 7 0
City Views Downloads
Beijing 19 2
Seattle 6 0
Ashburn 1 0

{[{ mDataHeader[type] }]} {[{ month[type] }]} {[{ tData[type] }]}

Share

Share to:

Comment

Related Items