Alias-Free Spectral Estimation of Stochastic Processes Citation Adegbola, Mashood Olayide (1971) Alias-Free Spectral Estimation of Stochastic Processes. Dissertation (Ph.D.), California Institute of Technology. doi:10.7907/ymat-1n64. https://resolver.caltech.edu/CaltechTHESIS:10112017-132121017 Abstract A scheme for the practical estimation of power spectrum from randomly-timed samples is proposed and investigated for wide-sense stationary point processes. The sampling process {t n } is assumed to be stationary point process statistically independent of the sampled process X(t). Stationarity of {t n } admits that joint statistics of t k , t k+n do not depend on k. Closed form analytical formulae are derived for the spectral window Q m (f) and for cov{S ^ (f r ), S ^ (f q )}, var{S ^ (f r )} for the particular case of independent identically distributed sampling intervals. Results confirm the alias-free character of the Poisson sampling scheme even for non-bandlimited spectra. It is shown further that for Gaussian processes with very smooth spectra Poisson sampling process can yield more reliable estimates (i.e., with a smaller variance) than the well known method of periodic sampling. Item Type: Thesis (Dissertation (Ph.D.)) Subject Keywords: (Electrical Engineering and Applied Mathematics) Degree Grantor: California Institute of Technology Division: Engineering and Applied Science Major Option: Electrical Engineering Minor Option: Applied Mathematics Thesis Availability: Public (worldwide access) Research Advisor(s): Martel, Hardy Cross (advisor) Caughey, Thomas Kirk (advisor) Thesis Committee: Caughey, Thomas Kirk (chair) George, Nicholas A. Martel, Hardy Cross Todd, John Whitham, Gerald Beresford Defense Date: 15 April 1971 Funders: Funding Agency Grant Number Caltech UNSPECIFIED Record Number: CaltechTHESIS:10112017-132121017 Persistent URL: https://resolver.caltech.edu/CaltechTHESIS:10112017-132121017 DOI: 10.7907/ymat-1n64 Default Usage Policy: No commercial reproduction, distribution, display or performance rights in this work are provided. ID Code: 10504 Collection: CaltechTHESIS Deposited By: Benjamin Perez Deposited On: 11 Oct 2017 21:01 Last Modified: 29 May 2024 21:50 Thesis Files Preview PDF - Final Version See Usage Policy. 13MB Repository Staff Only: item control page

ALIAS-FREE SPECTRAL ESTI~iATION OF STOCHASTIC PROCESSES Thesis by M. Olayide Adegbola In Partial Fulfillment of the Requirements For the Degree of · Doctor of Philosophy Ca~ifornia Institute of Technology Pasadena, California 1971 (Submitted April 15, 1971) ii @)1971 Mashood Olayide Adegbola ALL RIGHTS RESERVED -iii- ACKNOWLEDGMENTS I would like to thank my former thesis advisor, Professor H. C. Martel, for his valued help during the early phase of this investigation. I am deeply indebted to my present advisor, Professor T. K. Caughey, whose guiding inspiration and technical assistance made this work possible despite some difficulties. My sincere thanks to Professor F. J. Beutler of the University of Michigan for making his cited works available to me. I am also grateful to Professor Nicholas George for his interest and support throughout my stay at the Institute. The financial support provided by the Institute through numerous assistantships and fellowships is gratefully acknowledged. Further, I wish to thank Mrs. Ruth Stratton very much for her skillful efforts in preparing this work for presentation. Others in the Electrical Engineering Staff to whom I am grateful for assistance are Mrs. Aldre'Khoury, Mrs. Carol Norris, Mrs. Karen Current and Mrs. Paula Samazan. -iv- This work is dedicated to my cousin Dr. H. M. A. Onitiri who provided the motivation. -v- ABSTRACT A scheme for the practical estimation of power spectrum from randomly-timed samples is proposed and investigated for widesense stationary stochastic processes. The sampling process {t } n is assumed to be a stationary point process statistically indepenX(t) • Stationarity of {t } admits tk , tk+n do not depend on k . Closed dent of the sampled process that joint statistics of n form analytical formulae are derived for the spectral window ~(f) and for cov{S(f ),S(f )} , var{S(f )} r q r for the particular case of independent identically distributed sampling intervals. Results confirm the alias-free character of the Poisson sampling scheme even for non-bandlimited spectra. It is shown further that for gaussian processes with very smooth spectra Poisson sampling process can yield more reliable estimates (i.e., with a smaller variance) than the well-known method of periodic sampling. -viTable of Contents INTRODUCTION 1 I. EMPIRICAL SPECTRAL ESTIMATION 4 1.1 Statement of the Problem 4 1.2 Designing Estimation Algorithms for Practical Use 6 1.2.1 Review of Numerical Integration Schemes 6 1.2.2 Algorithms for Estimating from Randomly Sampled Data 8 1.3 II. The Sampling Process {t } n 10 FIRST ORDER STATISTICS 12 2.1 Introduction 12 2.2 Analytical Derivation of the Spectral Window 12 2.3 Desired Properties of the Spectral Window 14 2.4 Windows Resulting from Some Common Sampling Densities 16 2.5 2.4.1 Periodic Sampling: The Blackman-Tukey Window 17 2.4.2 Poisson Sampling Process: Window 17 2.4.3 Rectangular Sampling Process: Rectangular Window 19 Details on the Periodic and Poisson Windows 20 2.5.1 21 The Poisson Integrability and Absolute 1ntegrability u p Qm(w), ~(w) for Large 2.5.2 Behavior of 2.5.3 Computation of Side Lobes m 24 25 2.6 Aliasing of Spectral Estimates 30 2.7 Modification of Spectral Windows 32 2.7.1 Application of the Bartlett Function to Nonperiodic Sampling Scheme 2.7.2 Application of Hanning, Hamming Functions 39 to Nonperiodic Sampling Scheme 2.8 · Aliasing: Extension to Nonuniform Sampling and Criteria for Alias-Free Estimation 33 41 . -vii- III. VARIABILITY AND COVARIABILITY 46 3.1 Introduction 46 3.2 Analytical De rivation of cov[S(f ) , S(f )], var{ S(f )} r q r Derivation of the Blackman-Tukey Result for Periodic Sampling 3.3 3.4 Results for the Poisson Sampling Scheme 3.4.1 An Exact Result for var{S(f )} r 3.4.2 lim jf-fj-+-0) r 3.4.3 3.5 IV. q cov[S(f ),S(f )] r q cov[S(f ) ,S(f )] r q An Approximation for var{S(f r )} 47 54 57 59 61 62 Special Cases of Practical Interest 73 3.5.1 Estimation of Spectral Spikes 73 3.5.2 Estimation of Slow-Varying Spectra 76 SUMMARY AND CONCLUSIONS 80 4.1 Summary and Conclusions 80 4.2 Poisson vs. Periodic Sampling 81 4.3 Practical Implications of the Random Sampling Scheme 83 4.3.1 Communications Technology 84 4.3.2 Seismic Data Processing 84 · 4.3.3 Oceanography 84 4.3.4 Structural Design 85 APPENDIX A 86 APPENDIX B 90 REFERENCES 94 -1- INTRODUCTION Power spectral techniques ha ve found many useful applications not only in communication engineering but also in such diverse fields as Astronomy, Meteorology, Structural Dynamics, Oil Exploration and Economics. From measurements of power spectra, it has been possible in the case of linear systems to obtain useful estimates of the stochastic r elations generating such time series (see, for example, reference [19)). In turn, identification of linear systems via spec- tral analysis provides a means of constructive fault-finding and subsequent modification of the design of radio receivers, aircraft and other linear (or approximately linear) systems which are subject to stochastic excitation. In areas of application involving natural sys- tems (which cannot be modified) scientists and engineers have b een able by extrapolation to predict responses of those systems to well- defined random excitations. Symbols from a communication source [10], radar echoes from distant planets [12],[13] , swell recordings from distant storms [11], reflections from seismic explosions [19], wind velocities in atmospheric turbulence, and even day-to-day price fluctuation in the stock market are all examples of random signals with which engineers and scientists have to cope. In general, amplitudes of these signals are random and in some cases significant readings of these amplitudes arrive at random times. For many of the applications mentioned above analog methods of spectral estimation have proved inadequate due, for the most part, to the degree of frequency resolution demanded by modern noise studies. Since the advent of digital computers, discrete-time estimation of power -2- spectra has been widely carried out using the standard procedure of R. B. Blackman and J. W. Tukey [1]. For a certain class of signals, the standard procedure utilizing ampli tude readings taken periodically has produced useful results, while for others it has led to erroneous estimates. Of particular concern is the error due to aliasing or unwanted contribution (by "folding over" replicas of the true spectrum) f rom frequencies which are even multiples of the Nyquist frequency (fN = l/2~t). The attractive possibility of using randomly-timed samples as a means of eliminating aliasing was suggested by Blackman and Tukey [1] and investigated by H. S. Shapiro and R. A. Silverman [9], F. J. Beutler and 0. A. Z. Leneman [2]-[5], all of whom have shown that Poisson sampling is alias-free even for non-bandlimited spectra. Whereas the works cited above have provided the groundwork for this investigation, much-desired processing algorithms utilizing this particular technique together with qualitative analysis are still lacking. In this :thesis, a scheme for the practical estimation of power spectra from randomly-timed samples is outlined and fully investigated for a real wide-sense stationary stochastic process requires the sampling process {t } n admits that joint statistics of The scheme to be a stationary point process whose statistics are independent of those of {tn} X(t) • X(t) • tk+n' tk Stationarity of are independent of k. Chapter I provides the basic formulation of the problem together with relevant assumptions on the sampling process • An algorithm for practical estimation of power spectra from non-uniform sampling is outlined briefly. In Chapter II the first-order statistics of the estimator are fully investigated. An analytical expression for the -3characteristic spectral window is developed in its most general form. Methods of modifying spectral windows are discussed and windows arising from Hanning, Hamming and Barlett modifications are derived analytically. The notion of aliasing is generalized and testing criteria are discussed. Results on covariability and variability of estimates are presented in Chapter III. Closed form expressions of the Blackman-Tukey type are derived for the Poissson sampling scheme. It is further shown that for very smooth spectra, Poisson sampling scheme can achieve a better variance than the method of periodic sampling. For wildly fluctuating spectra, on the other hand, results show that Poisson sampling is equally unreliable. -4Chapter I EMPIRICAL SPECTRAL ESTIMATION 1.1 Statement of th e Problem In the empirical estimation of power spectral density S (f) of a X real stationary stochastic process X(t) , the following procedure is customary: (i) X(t) (ii) A processing algorithm is set up to utilize the sampled data is sampled at prescribed times X(t ),X(t )···X(tN) 1 2 along with appropriate information on the sample intervals to obtain an estimate at some frequency (iii) ately. of f t ,t •••tN . 1 2 s" x (fr ) of the noise spectrum r " The estimator is evaluated and modified appropri- S (f ) x r Evaluation consists in an analysis of the mean and variance Sx(fr) . The need to modify the estimate usually arises during the analysis of the mean of " Sx(fr) . Sampling may be done periodically as in the Blackman and Tukey algori thm or in a non-uniform manner as has been suggested and studied in the literature (1]-(5],(9]. In practice, a statistical description can be imposed on these intervals by either (a) sampling the process X(t) according to some well known distribution like Poisson, periodic, rectangular; or (b) observing, in the case of natural phenomena like ocean waves, the arrival times of records and approximating their statistical behavior as best as possible. The goodness of estimating spectra by methods described above depends to some extent on how accurate the sampling interval statistics -5have been approximated. In general, one does not need to know the actual sampling times for practical estimation, but a statistical description is required for a qualitative analysis. The scheme to be developed later in this thesis admits of data whose sampling intervals are additively random. Thus, it will be less sensitive to sampling jitter tha n the procedure of R. B. Blackman and J. W. Tukey in which the sampling intervals are fixed and equal. The equations fundamental to spectral estimation procedures for a real stationary (wide-sense) stochastic process X(t) are: 00 sX (f) J R (T) cos WT dT (1.11) X -00 R (T) X ... (1.12) E {X(t +T) X(t)} where SX (f) , is the power spectral density R (T) X is the autocorrelation function of X(t) For ergodic processes ensemble averaging implied in equation 1.12 may be replaced with time averaging defined by T R (T) X = lim T-+oo 1 T I X(t+T) X(t) dt (1.13) 0 Equations 1.11 and 1.13 underlie theoretical spectral estimation in continuous time or in discrete time with infinite data. Since empiri- cal spectral estimation is carried out in discrete time with finite data, a great deal of the efforts in the area is directed to developing -6improved sampling schemes. In what follows, we will develop some processing algorithms for randomly sampled data, modify one of them and evaluate its first and second order statistics. 1.2 Designing Estimation Algorithms for Practical Use The first step in designing a spectral estimator for a stationary process X(t) which is assumed ergodic is to obtain a suitable numerical approximation to the defining equation: T-T oo S(f ) • 2 J lim T:Tf X(t+T) X(t) cos 2TI f "C dt d"C r r OT-+oo (1. 21) 0 Numerical integration schemes for doing this abound in the literature, varying in complexity from the rectangular and trapezoidal approximations to quadrature formulas of the Lagrangian and Gaussian types (cf. Todd [22]). In this section some estimation algorithms are developed from first and second order integration schemes. The most suitable of these is selected and modified for subsequent evaluation. 1.2.1 Review of Numerical Integration Schemes Let f(t) be Riemann integrable in the interval [a,b] and let f( t ), n • 0,1,2· • •m n be samples of f(t) at the times t t ••• t • o' 1 m The area defined by the integral b I -J !:. f(t) dt a can be approximated by either of the following sums of rectangles -7- IF -I IB = I indicate that the base of corresponding rectangle m-1 t:, The subscripts on n=O f(t n )(tn+1- tn) (1. 22a) f(t n )(t n - t n- 1) (1.22b) m t:, I n=l is obtained by taking the forward step (F) or the backward step (B). A better approximation to the integral is obtained from an average of the last two sums, namely: (1. 23) which is none other than the trapezoidal approximation of an integral. To see this, we write out the sum (letting f n ~n+l ~ tn+l - tn , ~ f(t n )) and rearrange to . have m-1 • ~ {foal + L n=l f (a + a +l) + f . a } n n n m m Higher order integration schemes result from fitting the sampled data to polynomials of degree less than the number of sample points. Details leading to well known quadrature furmulas can be found in texts on numerical analysis (see, for example, pp. 59-61). Todd [22], As a rule polynomial fitting can be regarded as weighting the rectangles in the sums defined by 1. 22a, 1 .• 22b. -81.2.2 Algorithms for Estimating from Randomly Sampled Data Let X(t . ), i= 0,1,2,""",N 1 ary stochastic (real) process be samples from a wide-sense station- {X ( t)} and i=0,1,2,···,N t.' 1 be points of the stationary point process {tk} whose properties we will define later. Numerical approximations to the double integral in equation 1.21 will proceed as follows. approximate the inner integral as For simplicity, we will (letting T-T 1 T-T I '\.. = X(t+T) X(t) dt T - Na 0 and apply equations 1 . 22a, 1.22b and 1.23 to the outer integral to (w obtain respectively r = 2nf r ): m-1 N-n SF(fr) = 2 n~O (N=n)a k~l Xk+n ~cos wr(tk+n- tk)a~ ak+n+l (1.24a) (1. 24b) and m-1 + \ N-n 1 L (N- n)a n=l . k~l Xk+n~cos wr(tk+n- tk)(ak~+ ak+n+l)ak (1.25) Each of the above is a valid estimator as far as the integral in 1 . 21 is concerned and , in fact, for equi-spaced samples . §T(fr) is Blackman and Tukey's estimator (To see this, put tk+n - tk = n~t , '. ' -9ak+n = ak+l = ~t , tk = k~t in equation 1.25.) We remark in passing that there are infinitely many such valid estimators corresponding to the various integration schemes already mentioned. However, in order to obtain useful averaging filters or spectral windows, as they are commonly referred to in the literature, it is sometimes necessary to modify these estimators somewhat. sis showed that using either For example, an earlier pilot analy- " SF(fr) for general non- uniform s ampling schemes is equivalent to averaging with windows which are not integrable. " SB(fr) , on the other hand, averages with a Poisson window whose area is zero and whose bandwidth is therefor e infinite . Detail discussion on the averaging capability of the spec- tral windows is presented in Chapter 2. In general, modification of estimators can be viewed as conveniently weighting the rectangles in the sum of 1.22a, 1.22b. A systematic way of choosing these weights is that provided by higher order schemes of polynomial fitting. 61). (See, for example, Todd, pp. ' 59- Since higher order schemes will usually present analytical dif- ficulties, any choice of weights which does not hinder analysis considerably will be quite acceptable. Further accuracy in the inte- gration scheme can then be achieved by taking samples closer on the average. A simple such modification which we present here for further investigation assigns an approximate weight of 2 to the first rectangle and a weight of l to the other rectangles in the sum of equation 1.22b. Thus we have m· S" ( f ) = 2 r L\ n=l N-n 1 -----:-("' --:-) k=L\_ N-n 1 J: ) Xk+nx.. cos wr ( tk+n- tk ) ( ~+n+ aunl ~ k (l. 26) -10where n = 1 ={1 0 0 nl n f 1 The Blackman and Tukey estimator is easily derived from equation 1.26 above, by excluding the term m m defined by with the sum Z n=O n-1 a6nl and replacing the sum I m z {a0 + 2 n=O 1.3 The Sampling Process m-1 I a n=l n + a m} {t0 } employed in the scheme outlined {t } The sampling sequence n above is assumed stationary with independently distributed intervals. Furthermore, the processes cally independent. {t } n and {X(t )} n are assumed statisti- It is implied by the stationarity of the {t } n that the joint statistics of the respective number of points in any set of intervals are invariant under a translation of these intervals . Additional details on the theory of stationary point processes are provided by Beutler and Leneman, Ref. [3]. relevant (i) Summarized below are some properties and assumptions on the sampling sequence ~ t~+l- t~ =a~,~ =1,2,·· · ributed with mean {t } . n are independent and identically dis- a , common probability density characteristic function ~(iw) ~ E[exp(-iwa~)] p ( ~) , and -11- (ii) t~+V ~ t~ , v ~ 0 with equality iff V = 0 • + ) X(t )]} (iii) E[X(t~+v) X(t~)] - Et {EX (X(t ~V ~ with E , E t X inter- changeable. (iv) ~ * (iw) =/: , ~(-iw) where * implies complex conjugate. Furthermore, ~(iw) (v) (vi) satisfies the following: l~(iw)l< I~CO)I • lim ~(iw) "' 1 w :f: 0 0 w+±ClO The above conditions (v) and (vi) can be derived easily from definition ()C) ~(iw) • J exp(iwT) p(T) dT 0 and the Riemann-Lebesgue lemma. -12Chapter II FIRST ORDER STATISTICS 2.1 Introduction In this chapter we proceed to examine the first basic property of the scheme outlined in Section 1.2--the expected value of the estimator defined by equation 1.26. In particular, we will derive a general analytical expression for the spectral window Qm(w) and investigate how certain of its properties affect the resulting spectral estimates. The notion of aliasing for periodic sampling is reviewed and generalized to the non-periodic case. Testing criteria for alias-free estimation simpler than those of Beutler [5], and Shapiro and Silverman [9], are discussed and presented without proof. 2.2 Analytical Derivation of the Spectral Window To derive an expression for the spectral window we apply the expectation operator E E t E X A to S(f ) r in equation 1.26. Thus, E{S(f )} = 2 r (2.21) Interchange of expectation and summation is justified, since we are dealing with finite sums here. process By the stationarity of the sampling {t } , the summand in equation 2.21 does not depend on n k , -13so that Next, we take expectation with respect to X and introduce the Wiener- Khinchin relations to have: 00 E{S(fr)} a Il ; n• I S(f) Et{cos w(tk+n- tk)cos wr(tk+n- tk) ...co Now, (Re - "Real part of") Et {cos w(tk+n- tk) cos wr(tk+n- tk) ak('\+n+ aon 1 )} • ~ Et {"k("k+n+ aOnl) Re (exp[i(w+wr)(tk+n- tk)] + exp[i(w-wr)(tk+n- tk)] )} • ; [qn(w+wr) + qn(w-wr)] where The last equality follows from the assumption that the sampling intervals are independent and identically distributed. I Finally, we write 00 E{§(f )} r • ...co where S(f) Hm(w;w r ) df (2.22) -14b. 1 H (w;w ) ""- {'<m (w+wr) + ~ (w-wr)} m r 2 ~(Q) t, = m 2 Re{ a ~n I ¢n(Q) aH1 ¢ + a¢} n=l m = 2 Re {¢(1 - cp ) 1 - ¢ (2.22a) a ~n ¢ + a¢} aH2 (2. 23) Equation 2.23 gives an analytical expression for the so-called spectral window for non-uniform sampling in which the intervals are independent and identically distributed. For the uniform sampling scheme of Blackman and Tukey, the spectral window takes the form: Ou(n) ~ 2.3 ~ Re{ (1 + p)(l- pm) a t ~} 1 - ¢ ain n ~ (2.24) Desired Properties of the Spectral Window By equation 2.22 " E[S(f )] of the true spectral density H (w,w ). m r can be interpreted as a convolution r S(f ) r and the window function It is also customary to regard the mean value of " S(f ) r as a weighted average of the true spectrum over the bandwidth of the spectral window when the latter's main lobe is at f r For the purpose of useful-interpretation we can rewrite equation 2.22 in the following more general form: 00 f S(f) Hm(w;wr ) df E{S(f )} r - f -oo Hm(w;wr ) df --00 (2.31) -15Whence it is expected that (i) < ~(0) ~(w) = with equality iff w .. 0 00 (ii) I ~(w) dw = a O<a<oo , any m -oo For absolute convergence of the convolution integral (see, for example, Apostol [25] p. 490), and since S(f) is assumed bounded on (-eo,+oo), it is sufficient that (iii) J l~(w) I dw < oo Finally, it is desirable to have spectral windows ~(w) , indexed on m ' which form a defining sequence for the generalized function o(w) i.e. , (iv) lim Qm(w) = a o(w) This last requirement guarantees that biased. s"' (f r ) is asymptotically un- Throughout the preceding it has been taken for granted that . ~ (w) is continuous, real and even in 2.3.1 Bandwidth of Spectral Windows w • The bandwidth over which averaging of the true spectrum is done plays a major role in the stability of spectral estimates. For example, if the spectral window is "too wide", satisfactory resolution becomes rather difficult to achieve. In particular for spectra with bandwidths of the order of the window bandwidth, estimates at all frequencies will be very highly correlated. We are not concerned, at the moment, with -16how wide or how narrow the spectral window should be. Discussion on the trade off between resolution and stability will be taken up later. Rather, we are interested in establishing a working definition of bandwidth. Of the several definitions of bandwidth that abound in the litterature, the one that lends itself to easy calculation is that of Parzen [8]. It is simply that the bandwidth ~(w) is the base of the rectangle whose height is the peak of and whose area equals that of ~(w) Sm(Q) . of a spectral window ~(w) Symbolically we write: 00 f ~(w) dw . max -oo (2.32) ~(w) w The above definition assumes property (ii) of Section 2.3. Use will be made of equation 2.32 and property (i) of Section 2.3 to compute the bandwidths of periodic and Poisson windows later in this chapter. 2.4 Windows Resulting from Some Common Sampling Densities In Section 2.2 we derived general analytical formulae for ~(w) under the assumption that the sampling intervals are indepen- dent and identically distributed. istic function ~(w) By suitably defining the character- in equations 2.23 and 2.24, particular expressions can be obtained for spectral windows resulting from any sampling scheme for which our basic assumption is valid. In what follows, we shall derive particular expressions for the periodic, Poisson and rectangular sampling densit"ies, and exhibit some graphical plots as visual aids. -172.4 .1 Periodic Sampling: The Blackman-Tukey Window When sampling is done periodically without jitter as required for the Blackman-Tukey algorithm , we have ~t - sampling period, with sampling density O(T - ~t ) and characteristic function which, with puts equation 2.24 in the form: (1 + exp(in~t)) (1 - exp(ill~t) 1 - Multiplying the R.H.S. by exp (i~t)) 1 - exp(-in~t) 1 - exp(-in~t) before taking the real part leads to ~t n6t sin m n 6t cot -2- which is Blackman and Tukey's result 2.4.2 Poisson Sampling Process: (2.41) Q (w) (cf. [1], p.35). 0 The Poisson Window Here the sampling times are the occurrence times of events in a Poisson process (e.g., shot noise, radioactive decay). between successive events ak has the density The interval -18exp(-T/a) , p ~ (T) T < 0 and characteristic function = <P (in) -1 where a is the mean sampling rate and val. The joint density for n n-1 Pet a. · ·· ' a k+l k.+n' l:<.+n-1' the mean sampling inter- successive intervals is (T) .,. T exp(-T/a) n a (n-1) I from where we have the probability in an interval of length F(n,T) = (2.4la) (1- iS"2a)-l F(n,T) that there are n samples T , given as n T exp(- -r/a) an nl These last two expressions can be derived easily from the inverse Fourier transform of the joint characteristic function <P n (in) = (1 - H2a) -n Now , from equation 2.4la, we have and consequently, (from equation 2. 23) -19- (1- ~*)(1- ~*)-l Multiply R.H.S. by to get 1~1 :magnitude of the complex function where ~ • To simplify, note that we can write ~ • cos ~ exp(i~) with ~ • arc tan (na) so that QP (n) • 2a {cosm+ 2 ~ sin(m+l) ~/~in ~} m 2 . 4.3 Rectangular Sampling Process: As Rectangular Window a final example we consider the rectangular sampling process by which sampling intervals are uniformly distributed in mean (2.42) a , i • e. , l/2a p"' (T) • '""'k . { O otherwise The corresponding characteristic function is ~(n) • . (2illa) -1 [exp(i2na) - 1) and after carrying out the indicated differentiation, [0,2a) with -20a[l + i(~a- cot na)] whence equation 2.23 gives does not lend itself easily to The above expression for further analytical investigation. ~(n) - 2a ( s~;;,zna + However, we shall use the form 1 +sin nna(cot na - l_)]} na (2.43) for obtaining graphical plot of the rectangular window. Shown in Figure 1 are the periodic, rectangular and Poisson windows as given in equations 2.41, 2.43 and 2.42 respectively . They are plotted for m = 25 and a = 1 Based on the figure, one is . inclined to deduce that both Poisson and rectangular sampling processes are "alias free" for non-bandlimited spectra, while the periodic process is not. Analytical confirmation of this observation will proceed in subsequent sections for periodic aQd Poisson windows. As was pointed out earlier, the rectangular window is not tractable analytically. 2.5 Details on the Pe riodic a nd Poisson Windows In what follows we take a closer look at the periodic and Poisson windows (the rectangular window will not be investigated further due to lack of a suitable closed form expression for -21In particular, we will verify properties (i) - (iv) of Section 2. 3 and compute analy tically the bandwidth and sidelobes for each of the windows for comparative analysis. 2.5.1 Integrability and Absolute Integrability (a) To see tha t Periodic window. 2n/~t lobes at multiples of Qu (w) a ~t {1 + 2 m Qu(w) m has repeated major , we write m-1 L cos n~t + cos m~t} n=l cos 9 c cos (9 + 2nTI), n ~ 1,2,3,••• to have and apply (2.51) Equation 2.51 suggests r estricting the periodic window to the band lwl ~ n/~t otherwise integrable. On the other hand, Qu(w) m is neither integrable nor absolutely ·Qu(w) w m rect 2TI/~t is both integrable and absolutely integrable, It! < T rect 2t T "' { 01 t < T In particular, we have from tables TI/~t sin mw ~t J -TI/~t ~t dw "' 2TI , (cf Ref. [21], p. 366) tan - 2 Absolute integrability follows from the finiteness of the integration limits and ordinary integrability demonstrated above. Further, we note Fig. 1. 1.000 Spectra l HindoHs Arising Poisson Sampling Schemes 1.200 1.1100 1.600 SCII.f fFCTOR : 10 I .000 from Periodic Rectangular and FREQ.-RRD/SEC 0.800 POISSON WINDOW - M = 25 MEAN SAMPL. INT.<SECS) = RECT .-WINDOW---M = 25 MEAN SAHPL. INT.<SECS) = PERIODIC WINDOW-M I I N N -23that ~(w) (b) lwl From equation 2.2a, we have Poisson window. cos n+l \jJ cos (n+l) \jJ, . -1 \jJ • tan wet. whence it follows that ~(w) ~ ~(0) a is integrable, since Further, oo TI/2 -oo • 2 -TI/2 I ~(w)dw m I cosm\jl sin(m+l)\jl sin \jJ d\jl 2TI (cf. (21), p. 377) QP(w) m It follows too that interval. 2(m+l)et. · ~s ab so 1 ute1y i n t egrabl e over any fi n it e .Now,using the result . p lim w-..oo ~ (w) .. derived in the appendix, we can write for large L L • I IQ~(w)ldw 5 + 0(1/1 ) -L whence it follows that ~(w) is absolutely integrable. -242.5.2 Behavior of ~(w) , Q!(w) for Large m The main interest here is to see whether or not the window functions approach the delta function as m increases without bound. our purpose it will be sufficient if for F(w) For continuous, bounded and · absolutely integrable, 00 lim J ~(w) F(w)dw • a F(O) m-+oo -00 a where is the area under (a) Q (w) • m Consider the integral rr/6.t Iu • I Q~(w) F(w) dw and assume that -rr/6.t the sense described above. Then rr/6.t I u • J w6.t sin mw 6.t cot - 2- F(w) dw -rr/6.t mrr • I sin :t: F(y!m) t.an(y/2m) .£y_ m -mrr 00 i.e., lim I • m-+oo u 2F(O) I sin :t: dy • y 2rrF(O) -oo Consequently, lim m-+oo (b) u w q;<w) rect 2rr/6.t • In a similar manner 2rro(w) F(w) is "good" in -2500 J Q:(w) F(w) dw rr/2 m cos ~ sin(m+l)~ F(tan ~) d~ sin ~ --rr/2f (m+l)¥ I .. 2 cosm(~) sin y - - -- - - - F(tan ~) Tr sin~ -(m+l)2 !;1 00 lim RHS m+oo I sin y_ F(O) dy "' y - 2rr6(w) =2 27T F(O ) - 00 which implies that lim m+oo QP (w) m Additional evidence is provided by the observation that the major peaks of ~(w), ~(w) with m . given respectively by 2ma , 2(m+l)a tend to infinity Als o the bandwidths obtained from equation 2."32 as cps , l / 2(m+l}a cps respectively , tend to zero with increasing 2.5 . 3 l/2filf,t m • Computation of Side Lobes (a) Periodic. has zeros near ~t u Qm (w) ~ ~t sin mw ~t cot --2 wk = krr/ filf,t , k=l , 2, 3 so that its side lobes have peaks near wp = (2p+l) given approximately by Tr 2~t p ... 1 , 2 , · ·· C\1 . . 0 0 r- If) ...... 0 0 0 ..... . I 0 N In zo 0 0::: ~ ~1:: ~ 0 - N >--< w Cl 3: zc:) ON - O~.n I 3:o 0 '-" 3:0 ,-...1/) -0 O'o '-" 0 ,..... 0 If) I I \ -_ - - - /( I "\.c ./ -. -- 7 II I --- '"'"" II I .r 2- ::s..:A.__ =x. r* 6 :.,;;:>- lll""'o- -=:::..,~;-= ~ Period ic and Poi sson Sp e ctra l WindO\·JS fo r m "" 50, a: = 1 ~L.. -1"'\.~7.-- ~(w)! ~(0) Fig. 2. <I ~(w) I ~(0) t") I I N I (1\ . 0 . I 1/} 0 0 I 0 C\1 l/) zo 0 ~0 No w A 3: Zo ..._, D~ E)lf) I 3:o 0 '-" 0 ,.... Fig . 3. Periodic and Poisson Spec tral Windovrs f o r m I N I -...J -28(2.52) (b) Poisson. 2a cosm+ 2 ~ sin(m+l)~ sin ~ QP (w) = m has zeros near ~k = krr/m+l 1 or wk = a tan(kTI/m+l), k= 1,2,·. • or 1 (2p+l) 1T 1 2 wp "" a tan m+l 2 ' P.. ' '· • • and consequently, peaks near approximately equal to (2. 53) Using equations 2.52, 2.53 we have computed some coordinates for the first three side lobes (corresponding to m = 50,100,··· ,300. Results for those shown in Figures 2 and 3. m = 50,100 p=l,2,3) and agree quite well with It is observed from Table 1 that side lobes move closer to the origin but do not decrease with increasing as is expected. ~ather, the side lobes increase with m m but do not exceed a fixed fraction (namely, about 1/5 for the first side lobe) of the main lobe. Figures 2 and 3 lead to the conclusion that side lobe contributions are less pronounced with Poisson sampling than with periodic sampli~g. However, this discrepancy disappears if sufficiently large. Finally, we remark that the results of this sec- m is tion confirm the known result that Poisson sampling is alias free even -29(a) m 50 100 150 200 250 300 wn 1 0.0942 0.0471 0.0314 0.0235 0.0188 0.0156 (b) m 50 100 150 200 250 300 First Side-lobe p -0 . 21205 -0.21217 -0.21219 -0.21220 -0.21220 -0.21220 wl 0.0927 0.0467 0.0312 0.0234 0.0188 0.0156 Second side-lobe p = wu 2 0 . 1570 0 . 0784 0.0522 0 . 0392 0.0314 0.0261 (c) Q~(wl) p = 1 ~(w~) 0.12706 0.12726 0.12729 0.12731 0.12731 0.12732 Q! (wl) -0.16945 -0.18980 -0.19702 -0.20072 -0.20297 -0.20449 2 wr 2 0.1552 0.0779 0.0521 0.0391 0.0313 0.0261 ~(w~) 0 .06761 o. o·9329 0.10357 0.10908 0.11252 0.11486 Third Side-lobe p • 3 m u w3 Q~(w~) p w3 Q!(w~) 50 100 150 200 250 300 0. 2191 0.1099 0.0732 0.0549 0.0439 0.0366 -0.09058 -0.09085 -0.090.90 -0.09092 -0.09093 -0.09093 0.2190 0.1093 0.0729 0.0548 0.0438 0.0365 -0.02567 -0.04929 -0.06062 -0.06714 -0.07136 -0.07432 Table 1. _Some Coordinates for the first three side-lobes. -30for non-bandlimited [9]) 2.6 spectra (c.f., Beutler [5], Shapiro and Silverman 0 Aliasing of Spectral Estimates In the practica l estimation of spectra from uniformly sampled da ta, errors have been known to occur in the estimates due to the periodic nature of the sampling scheme. These errors come as unwanted addition of estimates of the true spectrum at certain integral multiples of the Nyquist frequency (f n = 1/~t , ~t = sampling interval). Furthermore, estimates at these multiples of the Nyquist frequency are indistinguishable from one another. For bandlimited spectra, this problem, often referred to as aliasing, is easily surmounted by taking samples at or above the Nyquist rate. For non-bandlimited spectra, random sampling schemes (namely, the Poisson sampli~g scheme) have been found to reduce or even eliminate aliasing effects. In this section we derive the basic concept of aliasing for periodic sampling directly from the corresponding window function and generalize to nonuniform sampling schemes. To this end, let us rewrite equation 2.22 ~as: 00 E {S (fr).} J S (f) Qm (f - fr) df = (2.61) 0 and assume for the moment that S(f) is non-bandlimited. form sampling interval we showed in Section 2.5.1 that periodic with period f p = 1/~t cps. ~(f) .. ~(f) rect 1 J~t@ With uni- ~(f) is Consequently, we may write 00 L c(f + vf p ) \)=-00 (2. 62) -31~ defines the convolution integral. where Now, using 2.62 in 2.61 and recalling (from Section 2.5.2) that J lim {~(f) rect 1 llt } "' 211'0 (f) m-+oo we get A lim E{S(f )}= 211' r m-+oo S(f + vf ) (2. 63) P r \)=..00 and oo+k A lim = 211' E {S(f + kf )} r p 2 s (f + \)f ) V=-oo+k (2. 64) P r Equations 2.63, 2.64 together define aliasing for uniform sampling scheme. Both equations show that on the average, estimates at A f , f + kf r r p k=l,2,··· are indistinguishable. is the true spectral density S(fr) plus magnitudes of fr + vfp S(f) at Each estimate S(f ) r at the frequency of interest v=1,2,···. f r The latter inter~ pretation becomes obvious when we write equation 2.63 (leaving out the 211' for convenience) as: . 00 lim "' E {S(f )} = S(fr>+ r L s (f r + vf P) ~ \)=..00 vi'O To recover S(f ) r from S(f ) above equations suggest that we make closer together) so that S(f) with periodic sampling, the r f p large enough (by sampling is zero outside (-f ,f ) • p p For band- limited spectra this is possible to within the capability of the sampling equipment. However, for non-bandlimited spectra nonuniform sampling patterns must be used,. since it is practically impossible -32to sample uniformly at an infinite rate. Definition 1: An estimate is said to be aliased in the ordin- ary sense iff for some integer lim m~~ k and sampling interval EjS(f)S(f r + k/~t)l r - ~t 0 In terms of definition 1 above, a sufficient condition for estimates to be aliased in the ordinary sense is that the spectral window a periodic function of ~(w) be w , with period w - 2~/~t . p 2. 7 Modification of Spectral Windows One peculiarity of the estimator under investigation is the term aonl filter which was included to insure that the resulting averaging ~(w) is integrable in (-~,~) • As we have shown earlier, integrability is only one of the essential features of a usable averaging window. The need for modification of this type did not arise for the periodic sampling, since the latter scheme is restricted to bandlimited spectra by which integrability of the spectral windows is virtually guaranteed. All efforts in the past have been directed primarily at de-emphasizing side lobes of windows and forcing Qm(w) to become positive definite--a necessary condition for realizability of filters and other transfer functions. Tailoring of the spectral window in this manner is usually accomplished by weighting the mean-lagged products with an even,positive definite function. Several weighting functions bearing the names of their innovators abound in the literature for the periodic sampling~ In this section we extend -33- these tailoring techniques to the nonuniform sampling scheme and examine in detail the effects on the Poisson window of Hann's and Bartlett's weighting functions. A general modification of the scheme under study can be written as: (2. 71) X where T m T m By appropriately defining D (T) , the resulting modified window can \) be derived analytically as was done in Section 2.2. - To illustrate, we wi l l derive the Poisson-Hanning and Poisson-Bartlett windows. 2.7.1 Application of the Bartlett Function to Nonperiodic Sampling Scheme The weighting function suggested by Bartlett is given as T m T m and the corresponding spectral estimator is SB(fn) ~ 2 m N-n 1 \ ---- \ (1 1.. N-n k!::.l n•l (2.72) -34- Proceeding as before, we get 00 m E {S(f )} r Cl 2 s (f) Et I n=1 J ( (1 - tk+n- tk ma ) cos w (t r k+n tk) -00 X cos w(tk+1- tk)~(ak+n+ acn1)} df 00 -J B S(f) Hm(w;wr ) df -oo with HB(w·w ) c m where • r ~(n) . 1 {QB(w+w ) + OB(w-w )} 2 m r 1n r is the unmodified window derived in Section 2.2 and with B tJ, qn(n) "'Et { ~(~+n+ acn 1 )(tk+n- tk) Re[exp H1(tk+n- tk))} . k+n k+n .. Et {Re[~c~+n+ ac 1 )( L a) exp(iQ L a n k+1 ~ ~=k+1 ~ Taking expectations as indicated, and noting that )J} -35- and it can be shown that ~ · ~} + aunl <Hn (2.73) Finally, we have the Bartlett window arising from nonuniform sampling given by m QB(n) = 2 Re L (2.74) n=l m The Poisson-Bartlett Window With Poisson sampling intervals, so that the Poisson-Bartlett window can be written as (substitute above into 2.73, 2.74): ~B(n) • m 2~ Re{ L ~n+l - a n=O m L (n+l) $n+2(iSG) (2. 74a) m n=O m = 2a Re{ L n=O Now the second summation on the R.H.S. of 2.74a can be carried out by rewriting it as· -36- s .. m m+l L n<j>n+l n=-1 and noting that m+2 <j>S m =- L (n-1) <l>n+l n=2 Subtracting the latter from the former, and manipulating accordingly, we obtain Next, substitute tjJ = tan -1na an d take real parts as before, to obtain 2 2 - (m+l) cos tjJ cos(m+2)tjl] - 2a cos tjl } m+2 1, -i { = (m sin2,~) 2a cos tjJ [m sin(m+l)tjl sin tjl- a(m+2)cos(m+l)tjl cos tjJ (2.75) Shown graphically in Figures 4 and 5 are the normalized Poisson window, Q!(w)/~(0) ~B(O) for respectively, m = 50 and and its Bartlett modification m = 100 -~ (0) • 2 (m+l)a and respectively. ~B(w) I The main peaks are --- 0 . r- lf) 0 0 ..... 0 0 N 0 I N . lf) :z:: 0 0 ~ ~ql <to N Cl LLJ 3: ....... zc) ON E)lf) I 3:o OS 3:0 '"'LI! 0 ,o CY 0 ........ ..... 0 lf) Fig. 4 . I I 4.000 FREQ.-RRO/SEC 3.200 1 1 I Poisson Spec tral Window and I t s Bar t lett Mocifica t ion \ I"-.... r/~J 1 1 m = 50 , a = 1 I I I --.! w '· 0 . . . . I 0 I 1./) 0 0 0 N z: 1./) o 0::: C> :E ex:: _JQ ...... No w 0 3: ....... Z o DC\! Dl.f) I :S: o C5F 3: 0 r-.~ 0 .......... 0 CY E C> - r-: -a 1./) 0 - 0 0 0 "! - 1./) Fig. 5 . LLOOO FREQ .-RRD/SEC 3 . 200 Poisson Spectral WindovT an d Its Bartl e tt Hodificat i on , m = 100, a= 1 ~(w) / ~ (0) 0:> I I w -39- 2ci ~~0) • ~ (0) - --;- m L (n+l) n•O 2 • 2 (m+l)a _ 2a m a [(m+l)(m+2)] 2 2 .2(m+l) a[l - m+ a] 2m While reducing the side lobes considerably, the Bartlett modifi'cation tends also to increase the bandwidth as the plots inFigures 4 and 5 indicate. 2.7.2 Again, results are improved with increased m • Application of Hanning, Hamming Functions to Nonperiodic Sampling Schemes A general expression for the weighting function of the type proposed by Julius Von Hann and R. W. Hamming can be written as: < T m D(T) • > T - m where a a 0 0 • a • 1 ... 0.5 for "Hanning" 0.54, a 1 • 0.46 -The corresponding estimator is for "Hamming" (w m • 7T/ma) Taking expectations as before and noting that -40- we have (2.73) where ~(w) has been defined in equations 2.22a and 2.23. The Poisson-Hanning Window Application of the Hanning function to the Poisson scheme leads to the Poisson-Hanning window defined as: where as before 2a cos~2~ sin(m+l) ~ sin iJJ p 0 - (w) • 1n ,,, , 't' • tan -1 wa Further, "v "' since for ~ "v ma .. 'IT/m for large and very large m 2 m ' cosm+ 'IT/m ~ 1 sin(m+l)'IT/m = -sin 'IT/m. HP HP P P The plots of Qm (w)/~ (0) , Qm(W)/Qm(O) and shown in Figures 6 and 7 for m = 50, 100 respectively, show a similarity between the Bartlett effect -41a nd the Hanning effect on the Poisson scheme. The Hanning modifica- tion leads to smaller side lob es but larger bandwidth than the Bartlett modification. 2.8 Aliasing : Extension to Nonuniform Sampling and Criteria for Alias-Free Estimation In refe rences [4] and [9] the authors have laid down some criteria for a sampling process to lead to alias-free estimates. Beutler [5] defines alias-free sampling in terms of the capability to recover the true spectrum of r(n) x(t) from the "correlation sequence" via A s (f r ) = r(n) - (2.81) where 00 .. --2'JT J S(w) ~:(iw) dw 1 (2.82) -co According to Beutler, alias-free recovery of spectra from correlation sequence requires that s(f r > be uniquely defined by equation 2.81. Equation 2.82 relates spectral recoverability from r(n) to the sampling process by means of the joint characteristic function ~n* (iw) of n successive intervals. From 2.82 it is inferred that sampling is alias-free if a one-to-one mapping between implied by equation 2.82. r(n) and S(w) is Some criteria for alias-free sampling have been laid down by the authors mentioned above. Whereas Shapiro and Silverman restrict their findings to spectra which are square 0 . . r- L/) 0 --< 0 g --< ........... 0 • . . 0 I L/) 0 0 I 0 N L/) zo 0:: 0 ::E c::cO _j 0 W ~0~ 3: ....... zc) ON E)L!) I 3:o o= 3:0 -......o ,.....~ cYo E 0 ........ . L/) N I 1 I Fig . 6. 0-ILOOO I FREQ.-RRO/SEC 3.200 1 I I I I Poisson Spectral Hi ndm·J and It s Hann ' s Mo dification , m = 5 0 , ~<:::: CL =1 I """" I I N ---. r-:0 Lf) 0 ...... ~ 0 0 ...... . I 0 0 ~ 0 I ~ 0 Lf) Zo 0 0::: ~ No ...... ....Jo <I: LLJ 0 3: ....... Zo 0~ 3:o DLn OS I 3:0 ""'~ ......._o c:ro E 0 ...._, . Lf) 0 N Fig. 7. 14.000 Poisson Spectral \\'indmv and I t s Hann ' s H.odifica t i on , m = 100, FREQ .-RRO/SEC 3.200 a = 1 I w I .t- -44integrable, Beutler consid ers a more general class of spectra. Their criteria provide sufficient conditions for a sampling process to yield alias-free recovery of spectra from the correlation sequence. words, they provide conditions on 2.82 a one-to-one mapping of ¢> r(n) In other (iw) sufficient to make equation into S(w). n this can be found in the cited references. Further details on In what follows we shall define alias-free estimation of spectra from the sampled data X(t ) n and establish some simple criteria in terms of the spectral window ~(f) • We refer to Section 2.6 and note that, whereas equations 2.63, A 2.64 together imply that on the average A S(f ) r and S(f identical, this does not necessarily hold in general. is conceivable to have a spectral window function r + kf p ) are For example, it ~(f) of the form: 00 ~(f) = Qm (f) rect f ~a) @ L p(v) o(f+vfp (a)) p where p (v) :f 1 and \1=-00 (2. 83) f (a) is fixed for each p a • Then we will have corresp'o nding to equation 2. 63, 00 lim m -+ oo E {s(f I p(v) s(f + vf' (a)) )} r \1=-oo r (2.84) p and to equation 2.64 oo+k lim m-+oo E {s<f r + kfp (a))}- L p(v+-k)S(f + vf (a)) v~-oo+k r P (2.85) Now suppose there exist constants then it follows that such that -45,.. ~k lim E {crkS[f + kf (a) J} • r m-+oo p L p(v) S[fr + kfp (a)] \P-~k (2.86) Thus we have Definition 2: A spectral estimate is said to be aliased in the generalized sense, if for some integer interval a , there exist O'k' fp(a) lim EjS(f) - crkS(f r + kf p (a))j r k and average sampling such that • 0 The above generalizes Definition 1 to include spectra which are identical to within a multiplicative constant. Thus a sufficient condition for aliasing in the generalized sense is that the spectral window ~(f) exhibits main lobes at fk= kfp(a); k .a ·,··,-2,-l,O,l,2,••• • Finally we have: An estimation scheme is said to be alias-free with respect to non-bandlimited spectra iff the characteristic window has one and only one main lobe at f • 0. -46Chapter III VARIABILITY AND COVARIABILITY 3.1 Introduction In the last chapter we investigated the first -order properties of our estimation scheme. One consequence of the subsequent analysis is that the question of whether an estimate is aliased or not can be . At answered directly from the characteristic spectral window ~(w) the first order level there is also the question of bias. Because of A the finiteness of data, S(f ) given in equation 1.26 is necessarily r However, it can be shown from property (iv) of Section 2.3 biased. A that S(f) is asymptotically unbiased, i.e., r lim m~oo E{S(f )} r S(f ) r Thus far, the quality of our estimation scheme has not been well related to the amount of data through m where N used. It has, however , been related m is chosen to be always less than N • In the last chapter we found that for certain classes of spectra (bandlimited for periodic, non-bandlimited for Poisson sampling) estimates , A S(f ) r on the average approach the true spectra without bound. S(f ) r as m increases We will see later in this chapter that unless N is increased accordingly, the stability of estimates will be adversely affected. Consequently the pertinent questions to be investigated in this chapter are: (a) How close to the true value is the estimate available from a finite data size, N ; and what is the covariability of the estimates -47" " S(f ),S(f ) 7 r q (b) What is the effect of increased data on the quality of our estimation scheme? In particular, we will derive analytical formulae for S(f )} and q 3.2 cov{S(f ), r var{S(f )} • r Analytical Derivation of cov[S(fr), S(fq)] , var {S(fr)} To derive a formula for the covariance of two spectral estimates " " S(f ),S(f ) r q obtained from the same record via equation 1.26, we will assume that an effective length Ne, of data has been selected. Details on how to choose Ne are given in [1] (see for example, [1], p.l02) but for our purposes here it suffices to point out that N- m < Ne < N . Also, we shall make use of (3.21) Equation 3.21 assumes that x1 ,x2 ,x 3 ,x4 are joint gaussian variates with zero means and is derived in Parzen [23], pp. 92-93. For analytical convenience we shall make us~ of the following approximately equivalent form of equation 1.26 wr - 27Tfr to write -48"' "' m Ne Tk Tkn 4 cov{S(fr) ,S( f ) } =- 22 L L E {cov[X(t - 2 n) X(tk+ - 2-), q Ne a j,n=l i,k=l t k (3.22) Now , assume that X(tk) 's are gaussian with zero means and apply equation 3.21 to get: R (T) = E (X(t+T) X(t)) where X process X X(t) • is the autocorrelation function of the Next substitute the Wiener-Khinchin relations to get 00 J J S(f1 ) S( f 2 ) {cos w1 [t i- tk- t<-rij- Tkn)] RHS a - 00 which, upon expanding the cosine functions, becomes 00 ~ J J S(f 1 ) S(f 2 ) {cos[(w1+ w2 )(ti- tk) - ~(w1- w2 )(Tij- Tkn)] -00 1 + cos[(w1 ~w2 )(ti- tk) - 2<w1+ w2 )(Tij- Tkn)] 1 + cos[(w1+ w2 )(ti- tk)- 2<w1 - w2 )(Tij+ Tkn)] -49- The last result suggests the following substitution w • w' + w , 1 w2 a w' - w and consequently df 1 df 2 = 2dfdf', whereby T. . Ti . X(t.; - ~) X(t + ~)] 2 ... i 2 00 • J I S(f+f') S(f-f') {cos[2w'(t 1 - tk)- w(Tij- Tkn)] -oo + cos[2w'(t . - tk) - w(T .. + Tk )] ~ n ~J (3. 23) Using equation 3.23 we can now write an expression for the covariability A of two spectral density estimates A S(f ), S(f) r q in the form: 00 cov{S(fn) S(fq)} • 4 II S(f'+ f) S(f'- f) {A 1 (w',w,wr ,w) q ..00 + A (w' ,w,w ,w )} df df' 2 r q where b. A (w' ,w,w ,w ) • 1 r q + cos[2w(t 1 - tk) -w'(Tij- Tkn)]) (cos wrTkncos wqTij) ai~(ai+j+ aojl)(ak+n+ aonl)} and (3.24) -50- (3.24a) Again, we expand the cosine functions using +cos (a-b)] 1 cos a cos b •2[cos(a+b) to rewrite defining equations 3.24, 3.24a compactly as: A (w' ,w,w ,w ) r ~ • q 1 {>.~ (w' ,w+wr , w+wq ) + A~ (w' ,w-wr , w+wq ) 4 + A (w' ,w+w ,w-w ) + A (w' ,w-w , w-w ) ~ r q r ~ q + A (w,w'+w ,w'+w ) + A (w,w'-w ,w'+w ) ~ r q + A~ (w,w'+wr ,w-w) q A1 (w', w+w, w+w) r . q A (w', w+w , w+w ) r q 2 Finally, we have ~ r q + A,(w, w'-w w'-w )}. ,.. r' q ~ •1,2 -51- (3.25) (3.25a) In taking expectation as indicated in equations 3.25, 3.25a, overlapping of the intervals into account. must be taken One way to get around this is to consider a given per- mutation of the times ti-tk, ti+j-ti' tk+n-tk ti' tk, ti+j' tk+n and break up the intervals into non-overlapping intervals so that we can use the assumed independence of disjoint intervals. Subject to -52the requirement that the sampling instants form an ordered sequence, it is shown in the appendix that there are only six permutations of the times (i) ti+j > ti ~ tk+n > tk (ii) tk+n > tk ~ ti+j > ti (iii) ti+j ~ tk+n > t.l.. > tk (iv) tk+n > ti+j > tk > ti (v) tk+n > ti+j > ti ~ tk (vi) ti+j ~ tk+n > tk ~ ti Now for each permutation given above, we can proceed to take expectation as stated earlier. for all six permutations. This is done in detail in the appendix The results are summarized in the following expression for the covariability of spectral estimates. 00 • II S(f'+f) S(f'-f) {A(w' ,w;w ,w ) +A(w,w' ,w ,w ) } df df' r q where A(w,w' ,w ,w ) r q • 1 {A(w' ,w+w ,w+w ) + A(w' ,w+w , w-w ) 4 r q r q + A(w', w-w , w+wq ) + A(w', w-wr , w-wq ) } . r r ~ q -53A(w', w+w, w+w) r q [~n (2w'+w+w) + ~ (2w'-w-w) + ao 1(~(2w'+w+w) + ~(2w'-w-w >)l r n r n r . r + I R2 (k , j , k , n) * * ~k_._.(2w')(~j(2w'+w+w )+ao. ~(2w'+w+w >) [~n(w+wr) 1 l. J q q J + ~ * (w+w) + ao 1(~(w+w) + ~ * (w+w )] n r n r r + \ ( . L . k R 3 l.,J, ,n ) ( * i-k-n+1· ) ~.+. k (w+w )+ao . ~ (w+w) [~. k(2w'+w+w) • l. J - -n q J1 q l.r k-i+1 ( ~k+n- 1 (w r -wq )+ao n 1 Q ) *' (w -w ) +~ k(2w'-w-w ) (~k+ 1 (2w+w +w ) 1r q r nr q .+ao + n1 Q*k-i+1< 2w+w +w ) ) 1 r q I 1 [C k(2w'+w+w )(~k+n. j(w+w )+ao Qk-i-j+ (w+w )) 1 R (i,j ,k,n) .l.r -1.r n r 4 X (~j(Wr-Wq)+aoj Q(Wr-Wq)) 1 + ~1.. k(2w'-w-wr ) (~k*+ . . (w+w ) +ao Q*k-i-j+1 (w+uJ )) (~j*(2w+w +w ) . n-l.-J r n r r q + * ( * i-k-n+1 ) ~-· (2w'+w+w ) ~i+ ' -k- (w+w ) +aojlQ (w+w ) l. q J n q q R6 (i , j , k , n) 1 \ L -54[1; (w -w) + n where r q ~;*(2w+w +w) + ao l(<ll (w -w )+<ll*(2w+w +w )) ) } n r q n r q r q I; (n) ~ <ll~-l(n) ~ ...b. <:Hn ' (0) D.= ~ (in) R (i,j ,k,n) "' (i,j,k,n) 1 b. > t i> tk+n > tk t.+. ~ J R (i,j ,k,n) b. (i,j ,k,n) 2 tk+n > tk ~ ti+j > t.~ R/i,j ,k,n) = (i,j,k,n) b. t ,+ , ~ tk+n > ti > tk R (i,j ,k,n) b. (i,j ,k,n) 4 > tk > tk tk+n > t.+. ~ J ~ D. D. R (i ,j, k,n) = (i,j,k,n) 6 i,k = 1,2, • • • ,Ne n,j = 1,2,···,m - Note: J tk+n > t ~.+.J > ti > tk R ( i , j , k, n) = (i,j ,k,n) 5 3.3 ~ (3.27) ti+j ~ tk+n > tk ~ ti } b. ::& 6 u lJml R ~ 0 Derivation of the Blackman-Tukey Result for Periodic Sampling The Blackman-Tukey result can be derived from the general result given in equation 3.27 by setting the terms containing m aan 1 equal to zero and replacing the ordinary sums r aojl' m I a a. with n=l n j=l J the "trapezoidal sums", m z a = n=O n a m-1 + 2 a +am 0 n=l n L in accordance wit;h the Blackman Tukey estimator. In addition, we will need the following properties of the periodic sampling process: -55- = e ~t (3.31) = i6t(r21+ $"22+ $"23) e (3.32) (3. 33) Substitute 3.31, 3.32, 3.33 into the expression so derived from 3.27 to get A.(w',w+w, w+w) r a q 6t~ Re { ZI ~i-k(2w') ~*j(w+w) [~n(w+wr) + ~*n(w+wr)] Ne R 1 q n* + Z L ~i-k (2w') ~*j(w+w )[ ~n(w+w ) + ~ (w+w ) ] r q r R2 + + ~*ncw+w )J z I ~i-k (2w ') ~*j(w+wq )[ ~n(w+w) r r R3 + z I ~i-k(2w') ~*j(w+wq ) [~n(w+wr ) + ~ R4 *n (w+w )] r + <l>*~w+w ) ] + Z L ~i-k(2w') ~*j(w+wq )[ ~n(w+w) r r R5 + Z L ~i-k(2w') $*j(w+w )[~n(w+w) + q r R6 Conseq uent1y: ~*(.,..,r) 1} -56A.(w' ,w;w ,w ) r q and J\(w,w' ;w ,w ) r q Now, with equation 3.33 it can be shown easily that K u (w' ,Ne) IJ. a 1 Ne . k -2 l: ~1- (2w'> .. Ne i ,k=l (sin New' IJ. t ) 2 Ne sin w'IJ.t is real, whence we write: 1\(w,w';wr ,w) • q 1 -2 Ku(w' ,Ne) Hm(w;w) Hm (w;w) q r where m ~(x) 2/J.t Re z ~n(x) n=O = x/J.t /J.t sin mx /J.t cot --2 Finally, we use .the evenness of S(f) (from 2.41) to write -5700 cov{S(f ),S(f )} r q -J H(w;w ) H(w;w ) f(w) df q r -00 with 00 r(w) !:, ... J s (f '+f) S(f-f') Ku (w,Ne) df' -oo Also, 00 var{S(fr)} J [H(w;wr)] 2 f(w) df = -oo Both of these are Blackman-Tukey's results (see, for example, p. 125 of [1)), Note further, that the periodic variance kernel periodic with period w' p .. 7T/llt fI or p a l/2flt and 7T /2 t.t J 7T /2 1 6t ... K (w' ,Ne) dw' u -7T/26t J (sin Ne x2 dx Ne sin x> -7T/2 .. 7T/Ne 6t 3.4 Results for the Poisson Sampling Scheme For the Poisson distributed sampling intervals 1 1 - Wa (1 - i!1a) ~ (!1) • ~ 2 a~~+ 1 (n) - K u (w' ,Ne) is -58so that t.(w' ·w+w , r, w+wq ) - a22 Re Ne r I l R1 (i,j ,k,n) ~i-k-n+l(2w') (~*j+l(w+w) +O.l~*(w+w )) q (~n+l(2w'+w+w) + ~n+l(2w'-w-w ) + o r r n q J + ~(2w'-w-w )) 1 (~(2w'+w+w) r r + o.l~i-k-n+l(w+w ))[~i-k+l(2w'+w+w) (~k+n-i+l(w -w) J + 0 nl . r q r q ~k-i+l(w -w )) + ~i-k+l( 2w'-w-w ) (~*k+n-i+l( 2w+w +w ) r r q r q + 0 ~*k-i+l (2w+w +w ) ) J nl r q \ [.~i-k+l ( 2 w'+w+w ) (~k+n-i-j+l (w+w ) + 0 ~k-i-j+l (w+w ) ) 'r r nl r R (i,j ,k,n) + 4 (~j+l(w -w ) + o r q jl ~(w -w )) + ~i-k+l(2w'•w-w ) (~*k+n-i-j+l(w+w ) r r q , r + o ~*k-i-j+l(w+w) (~*j+l(2w+w +w) + o. ~(w -w ))] nl .r r q Jl r q + I . k ,n ) R ( ~· ,J, 5 ~*k-i+l( 2w'+w+w) [ (~k+n-i-j+l(w+w) +o q r nl ~k-i-j+l(w+w )) r (~i+j-k+l(w -w ) + 0 ~i-k+l(w -w )) + (~*k+n-i-j+l(w+w ) r q jl r q r + 0 ~k-i-j+l (w+w )) (~*i+j-k+l( 2 w+w +w) +o ~*i-k+l( 2 w+w -w ))] nl . r r q jl r q -59- I + R (i,j,k,n) 6 ell*k-i+l( 2 w'+w+w) (ell*i+j-k-n+l(w+w) +o . ell*i-k-n+l(w+w >) q q Jl q ) + ell *n+l (2w+w +w ) + ao 1 ( ell(w -w ) +ell* (2w+w +w ) )]~ [ elln+l (wr-w q rq n rq rq R~ (i ,j, k ,n) The ranges of summation ~"' 1, • • ·, 6 (3.41) are those given fol- lowing equation 3.27. 3.4.1 An Exact Result for cov[S(fr),S(fq)], var{S(fr)} Equation 3.41 can be simplified a step further by noting that 0nl ell - By appropriately replacing m m I L with n=l n=O we have for the Poisson sampling scheme 00 cov{S(f ),S(f )} • r q If S(f+f') S(f-f') {A(w' ,w,w w ) r q + A(w,w';w ,w )} df' df r q where + A(w';w-w ,w+w) + A(w' ,w-w ,w-w )} r and q r q -60- A(W I ,w+w ,w+w ) r ... a~ Re( Ne q I l~(i,j,k,n) [<l>i-k-n+l( 2w') <I>*j+l(w+w ) (<I>n+l( 2w'+w+w ) r q + <I>n+l(2w'-w-w >)]+ I [<t>*k-i-j+l(2w') <I>*j+l(2w'+w+w) r R (i,j,k,n) q 2 I (<I>n+l(w+wr) + <I>*n+l(w+wr)) J + [ <I>*i+j-k_:n+l(w+w ) R (i,j,k,n) 3 (~i-k+l(2w'+w+w ) ~k+n-i+l(w -w ) + <I>i-k+l(2w'-w-w ) r r r q q <t>*k+n-i+l( 2w+w +w >)] + I [ <I>i-k+l( 2w'+w+w )<I>j+l(w -w ) r q R4 (i ,J, . k ,n ) r r q <I>k+n-i-j+l(w+w ) + <I>i-k+l( 2w'-w-w )<I>*j+l( 2w+w +w )<I>*k+n-i-j+l(w+w >] r r r q r L [<t>*k-i+l( 2w'+w+w) <I>*i+j -k-n+l(w+w) + R (i,j,k,n) q q 6 + w•n+l (2w+wr+wqll]} i,k .. 1,2, • • • ,Ne (3.42) n,j • 0,1,···,m To obtain an expression for the variance of spectral estimates, we simply note that var {S(f )} r • cov {S(f) S(f )} r'' r -61- 3.4.2 From the Blackman-Tukey result derived in Section 3.3, we can again relate the covariability of spectral density estimates s(f q > to the spectral window . ~(w) For example, if H (w,w ) m = q s(fr ) ' H (w,w ) m r and 0 • With Poisson sampling intervals, the results obtained so far do not make this so obvious, but intuitively we expect a similar correspondence between covariability of estimates and overlapping of spectral windows . Now consider A(w ' ,w+w ,w+w) r as given in equation 3.42 and q note that IA(w' ,w+wr ,w+wq >I I Re[ ] I+ •.. I[ ll + The last inequality follows, since - lz < IRe zl I Making use of (i) ~(w) ( ii ) ~(w+w - 'V r ) • 1/iwa w~oo ~(w+w +wr -wq ) q 0 0 o+ I ~6 I[ ll } -62and taking the leading term in each summation it can be shown that 1 .A(w' ,w+w ,w+w ) lim -+ oo r q lwr -wq I Ne 21 w -w 14a 2 r q and consequently, 0 ( cov{S(f ), S(f )} r q 1 4 (f -f ) r as If-f I r q 3.4.3 -+ ) q oo • An Approximation for var{S(fr)} In Section 3.3 we were able to derive a closed-form expression A for A cov{S(f ), S(f )} q r for the periodic sampling scheme. The proce- dure outlined in Section 3.3 exploited certain properties of the corresponding characteristic function (cf. equations 3.31, 3.32, 3.33) which do not hold for the general nonuniform sampling scheme. In particular, the characteristic function for the Poisson sampling intervals does not satisfy equations 3.31, 3.32 and 3.33 exactly. However, an approximate closed form expression for be derived from equation 3.41 for values of near 4l .. ± w and 1/Ne w' = 0 and n~ar ~(w), respectively, where m, Ne K(w' ,Ne) w and are of the order of are typically very large. Specifically, we will derive an approximation for lw' I ~ 1/Nea can r Such an approximation is the most natural one to seek r since the effective bandwidths of 1/m w' vadS(f ) } var{S(f )} r valid in lw ± wr I ..5. l/2ma • We begin .by assuming the following form for .A(w',w+wr,w+wr) -63A. (w' ,w+w ,w+w ) a r r + ~*n+l(w+w ))] + R(w' ,w+w ) r (3.43) r and proceed to derive an approximation for the remainder function, R(w',w+w) r valid in lw' I~ 1/Nea and lw + w I <l/2ma r From equations 3.42, 3.43 we can write the remainder f unction exactly as: R(w' ,w+w ) r + L [(~i-k+l(Zw'+w+wr)~k+n-i-j+l(w+wr) · R4 _ ~i+l(Zw')~*k+.l(Zw')~*j+l(wtw )~n+l(w+w )) r r + -64- + L [ (cl>*k-i+l( 2w'+w+w ) cl>k+n-i-j+l(w+w ) _ cl>i+l( 2w')cl>*k+l( 2w') R r r . 5 cl>*j+l(w+w )$n+l(w+w )) + ($*k-i+l( 2w'+w+w )$*i+j-k+l( 2w+ 2w ) r r r r (3.44) where the ranges R (i,j,k,n) ~ =1,2,··· ,6 ~ Now consider the sum over we rewrite as and note that are as previously specified. ' R1 (i,j,k,n) which, for convenience, -65- - - 1 (1- 2iw'~)(l- i(w+w )~) 1 1 - i(2w'+w+w )~ r r 2 -2w' (w+wr)~ ------------~------------------------------- [1- 2w' (w+w )~2- i(2w'+w+w )~] [1 - i(2w'+w+w )~] r r r By maximizing the numerator and minimizing the denominator of the above I w' I ~ 1/Ne~ , expression in Iw'+w I ~ 1/m~ , obtain that. r I~C2w'+w+w)- ~C2w')~Cw+w )I r r 2w' (w+w )~ 2 • ----------------------~r ________~------------------~ w'(w+w )~2)2+ ( w'+w+w )2~2 112 [l+ ( w'+w+w )2~2]1/2 [(1- 2 < 2 r 2 (1 2 )-1 - Nem Nem • r 1 2 r _2_(1 + _2_ + 4 + •• ·) Nem Nem N 2 2 em ·whence we have (3.45a) Similarly, __.~i-k+2( :;. .___ 2w') ___::o=~--- ~i+1(2w')~*k+l(2w') ~ • 1 , near w' • 0 Next, we write, using 3.45a and 3.45b, the sum over (3.45b) R (i,j,k,n) 1 approximately as L ~*j+l(w+w )~i+l(2w')~*k+l(2w') [ (Hw+w ) ~ r . 2 ) n+l r - Nem H2w') _ ~n+l(w+w ) + ' (~*(w+w ) _ 2 )n+l _ ~*n+l(w+w r r Nem ~(2w') r >] -66- ~ ___ 2 2 [~*j+l(w+wr )~i+l( 2 w')~*k( 2 w') (~n+l(w+wr ) Ne R . 1 + ~*n+l(w+wr)) + 0(1/Ne)J (3.46) The latter is obtained from the former by taking the first term in the binomial expansion of 2 (~(w+wr) - Nem ~(2w') and noting that m+l n+l/m ~--m ~ )n+l 1. Approximations similar to equation 3.46 can be derived for each of the composite summations in equation 3.44 by first establishing approximations similar to equations 3.45a, 3.45b and then substituting into appropriate terms. Summarized below are some useful approximations whose derivations follow closely that of equation 3.45a: ~(2w'+w+w )~*(w+w ) r r ~*(2w+2w ) r ~(2w'+w+w ) ~*(w+w ) r r ,v • + 0(1/Nem) (3.45c) "'• ~(2w') + 0(1/Nem) (3.45d) ~(2w') (3.45e) *(2w+-2wr ) ~ *(w+w ) r ~ (3.45f) They are valid in quite easily. lw' I ~ 1/Nea , lw ± wr I ~ 1/ma and can be derived To illustrate further, we derive an approximation for -67the sum over I R (i,j ,k,n) 3 R (i,j,k,n) 3 which, for convenience, we write as ~*j+l(w+w ) ~i+l(2w') ~*k+l(2w') r Introducing equations 3.45, we have I ~* j +l (w+w ) ~ R3(i,j ,k,n) i+l k+l ( 2w' ) ~ * ( 2w' ) r [( (~(2w') +0(1/Nem))i-k+l -.,..., . -=------:-,.-:----- ~~+l(2w')~*k+l(2w') (~(w+wr ) + O(l/m2)) n+l _ ~n+l(w+wr >) (~(2w') + 0(1/Nem))i-k+l + ~i+l(2w') ~*k+l(2w') which, on taking the first terms in the binomial expansions a·s before, reduces approximately to l 1 m R 3 ~*j+l(w+w )~i(2w')~*k(2w')[~n+l(w+W) r (i,j ,k,n) r + ~*n+l(w+wr) + 0(1/Ne)] Except for the multiplicative constants and 3.46 are identical. (3. 4) 1/m, -2/Ne, equations 3.47 In a similar manner, it can be shown that approximations to the other sums are identical to equation 3.46 to within multiplic?tive constants and that these constants are of the order of 1/m, 1/Ne • m • O(Ne)) Since m is typically a fraction of we can write finally: Ne (i.e., -682 = 0(1/Ne) { a 2 Re R(w' ,w+wr) Ne 3 L R (i,j~k,n) ~*k+l(2w') (~n+l(w+w) + ~*n+l(w+w )) r r ~*j+l(w+w )~i+l(2w') r + 0(1/Ne)} Consequently: A. (w' ,w+w ,w+w ) a r r which, by the method of Section 3.3, leads to 2 A (w' ·w w ) • l K (w' Ne)H (w·w ) (1 + 0(1/Ne)) a ''r 2p' m'r where, as before: l/Ne 2 Kp (w' ' Ne) . • H (w·w ) • m ' r Ne L ~i+l(2w') ~*k+l(2w') i,k=l qP(w+w ) + oP(w-w ) m r r 'In m ~(w) • 2a Re L ~n+l(w) n=O Thus we have 00 A var{S(f ) r } 1\J • II S(f+f') S(f-f')K (f',Ne)H 2 (f;f )[1+0(1/Ne)]df df' p m r (3.48) Substitute ~(2w') • cos cp e icp , cp = tan -1 2w'a. to obtain the following closed-form expression for the Poisson variance kernel K (w',Ne): p -69K (w , Ne ) • cos 4 cp (l - 2 cos Ne cp cos Ne cp + Cos 2Ne cp ) p , 2 2 Ne sin cp Obtai n further, co J K(w' ,Ne)dw' TT/2 Ne J i , k=l - TT/2 -oo I cosi+k cp[cos(i-k)cp + i sin(i-k)cp] dp and i ntegrate term by t erm using (cf. tables, [21]) rr/2 J 0 cosp+q- 2x cos(p-q)x dx • ----~--~TT~-----2p+q-l(p+q-l)B(p,q) and the fact that K(w' ,Ne) is real to get Ne co I where K (f',N )df' • p e 1 I 2 N a i,k=l e TT 2i+k+l(i+k+l) B(i+l,k+l) (3.49) B(p,q) is the well known beta function, viz., B(p,q) • f(p) f(q ) f (p+q) Shown graphically i.n Figures 8a,b,c are the variance kernels K (w',Ne) , K (w' ,Ne ) u p for Ne • 200, 500, 1000 respectively. Observe that for sufficiently large data (namely, Ne > 1000) both kernels coincide. . 0 . CD 0 0 0 ~ 0 0 ..... 0 0 0 ..... ~ . I 0 I C'! 0 0 0 0 0 0 0 C'! 0 0 3:o '-" zo .. :;1' Wo ,... 0 0 C'! Fig. Sa. Periodic and Poisson Variance Kernels, Ne 0.500 FREQ.-RRD/SEC 0.400 200, a = 1 I -....J 0 I . 0 0 0 0 ~ 0 0 0 ~ 0 0 _, ~ '-" . 0 . I 0 0 N I 0 ~ 0 0 0 0 0 N 3:o · =': Z o W o ,.... . _, 0 0 N Fig. 8b. 0.500 Periodic and Pois son Variance Kernels, Ne FREQ.-RAD/SEC 0.1!00 500 , (l 1 I -...J I I-' ,.... ~ ......, I 0 I ~ 0 0 0 0 0 0 . ~ 0 0 0 ~0 •::J' Wo Z o 0 ~ 0 0 0 ~ 0 0 ~ - 0 0 ~ - 0 0 Fig . Sc. 0 .500 Pe riodic and Pois son Variance Ke rnels , Ne FREQ. - RRO/ SEC 0.1100 1000, CL = 1 I -...,J I N -73- 3 .5 Special Cas es of Practical Inter est From the results of the preceeding sections, a number of useful results bearing on assumptions about the spectrum S(f), can be extracted, In (1] pp, 104-106, we find results based on four such assumptions. Here we shall examine the two most useful of these, namely: (i) Slow-varying spectra, and (ii) Spectral spikes. Usually in practice, we are faced with estimating spectra which are fairly smooth except for one or two jumps; in which case (i) and (ii) can be applied piecewise. Analytically we shall be deriving approximations to the integrals: Cov {~(f ) ~(f )} = ~H (f,f )H (f,f )r(f) df r q _oo m r m q J S(f+f' )S(f-f' )K(f',Ne) df' 00 f(f) = _oo 3. 5.1 Estimation of Spectral Spikes of Suppose S (f) consists of very sharp peaks (of widths K(f',Ne) ) at f = ±f and S(f+f')S(f-f') 0 with area= A, << width so that we can write: -74- = A2{6 ( f+f '+f0 )6 (f-f '+f0 ) + 6 (f+f' -f 0 ) 6 (f-f'+f 0 ) + 6(f+f'+f0 ) o(f+f'-f0 ) + 6(f+f'-f0 )6(f-f'-f 0 )} Now, let J I(f,f 0 ) ~ _ 6(f+f'+f0 )6(f-f'+f0 )K(f',Ne) df' 00 00 = 6(2f+2f0 )K(-f-f',Ne) 0 = i 6(f+f0 )K(f+f',Ne) 0 The last equality follows from the evennes:: of K(f,Ne) and its boundedness and absolute integrability (cf. Lighthill), so that, 2 2 f(f) ~ ~ {6(f+f0 )K(f+f',Ne) + 6(f)K(f-f0 ,Ne) 0 + 6(f).K(f+f0 ,Ne) + 6(f-f0 )K(f-f0 1 Ne)} A2{6(f+f )K(0 Ne) + 6(f)K(-f Ne) ~ 2 0 1 01 + 6(f)K(f0 1 Ne) + 6(f-f0 )K(-O,Ne)Jl Where the last equality follows from ~ (x) 6 (x) and the observation that hill p. 42) Next~ = ¢ (0) 6 (x) K(f,Ne) substitute is a fairly good function (cf Light- K(01 Ne) = l to have -75- Whence cov{s<fr),s(fq)} = J~H (f;fr )Hm(f;fq )f(f) df ~ m A2{ = -2 Hm(-fo·;frm )H (-foq ;f ) + Hmo· (f ;frmo· )H (f ;f q ) + 2Hm(O,· f r )Hm(O,·fq )K(fo ,.Ne)} By the evenness of Hm(f ;f ) 1 2 cov{s (fr) we can write: ,s(fq)} = A2{Hm (f o ;fr)H m(fo;f q) + For Hm(O;frm )H (O;fq )K(fo ;Ne)} so that cov{S(fr),S(fq)} ~ A2{Hm(fo ;fr )Hm(f o;fq )} and In particular, -76Thus, we see that even with large data size reliable estimates of spectral spikes are difficult to obta~ regardless of the sampling scheme being used. 3.5.2 Estimation of Slow-varying Spectra Here we are considering S (f) which vary slowly enough so that within the bandwidth of K(f',Ne) the quadratic terms in its Taylor series expansion may be neglected. Then for f'~ 1 / Ne~ ~ NeCX ' we write: f'2 s (f+f') ::! S(f) + f'Sf' (f) + 2: sf'f'(f) + ......... f' '2 s (f-f') ::! s (f) - f'Sf' (f) + 2T sf'f' (f) + ......... and consequently 2 2 s(f+f')S(f-f')::! s (f) + f' {s(f)Sf'f'(f)- s~,(f)} + ..... ::! s2(f) + o (l/Ne2a2) The last approximate equality follows since and continuous and bandlimited to lfl~ B , f(f) and in particular, where f' ~ 1 INca , we have Ne~ ~ S (f) is assumed bounded Consequently for S(f) -77o(Ne) ~ (B+f) CX) J K(f',Ne)df ~ JK(f',Ne)df - (B+f) 1 for B>> 1 /Nea. _oo Note that, so that = + where the equivalent width, We 1 O(m '2) ' Ne is defined as: Now, for both the periodic and Poisson sampling schemes, we saw in Chapter I I that E{~(fr)} is approximately the same when the spectrum is restricted to the Nyquist band, and m is very large. Also, the equivalent width of estimate can be shown to approximate the corresponding window bandwidths. Since the bandwidths of the windows are for all practical purposes equal, we look to the quantity o(N ) e of the variability of estimates which w.e define as for a measure -78cr(Ne) + We where for both schemes, We~ 1 ; m:x For the periodic scheme we have from section 3.3 and for the Poisson (cf eqn. 3.49) Ne cr (Ne) where 1 L: --:--:--"::"'1_ _ _ _ _ __ 2i+k+l( i+k+l) B(i+l, k+l) P = rf!ea .i, k = 1 B(p,q) is the beta function. Figure 9 exhibits graphical plots of crp (Ne) and cru (Ne) in the computation of cr (Ne) for large values of range 10 $ Ne $ 100 p Ne consumes excessive machine time. However, to see that crp (Ne) monotone decreasing, note from figures 8 a,b,c, that positive definite and narrows with increasing Ne equal to 1 for all· Ne K(f',Ne) u is while its peak remains Consequently we conclude that crp(Ne) asymptotically bounded from above by cr (Ne) is is Tl}us we can ( cf section 4.2 below) say that it is possible to achieve a smaller v~iability with Poisson sampling than with the method of periodic samp~ing. . . . a: > . u cr (Ne) crp(Ne) . I 0 8N Fig. 9. Scaled Variance versus Data Size NE--DATA SIZE 9p~.-1-oo----oL:l-s-o----o~:-2a-o----o-:L3-70 _____oL:ij-S-o----o~:-ss-o----~·~-----L•------~·------~·~----_J• 0 0 0 0 0 0 N -~ CCo w u Zo a:o 0 (D 0 0 0 ~ 0 0 ..... 0 0 0 ..... 0 0 N I .....a I 1.0 -80CHAPI'ER IV SUMMARY AND CONCLUSIONS 4 .1 S umma.ry and Conclu.s ions For readers unfamiliar with the work of Blackman and Tukey [1], the research reported in this thesis is by no means complete. Details on 'planning for measurement' - 'prewhitening', 'post greening', etc. have been purposely left out of this account to avoid undue repetition of the main reference cited above. Primarily this thesis has been concerned with: (i) Outlining an estimation scheme admitting of data with random sampling intervals, and, (ii) Analyzing the mean and variance of such a scheme. The assumption of independent, identically distributed random variables which governs most of our analysis, is a useful one since it includes most of the practical sampling schemes. Although the analytical results contained in Chapters II and III seem to take advantage of this assumption, similar expressions can be derived for the most general sampling scheme. For this case, closed-form expressions for windows and kernels will present great analytical difficulties. On the question of aliasing, we found in this thesis a more practical way of testing for aliasing - obtain a plot of the spectral window, and check for maxima. ~(rn) In particular, we verified via an estimating algorithm that the Poisson sampling process is alias-free even for nonbandlimited spectra. One shortcoming of our algorithm is its non- general nature, since, in fact, it was tailored for the Poisson sampling process. As was pointed out in Chapter I, it is quite possible using -81techniques of numerical calculus to design algorithms to suit some particular sampling schemes, Our algorithm which invokes the rectan- gular approximation (instead of Blackman and Tukey's trapezoidal approximation) to the fourier cosine integral was found unsuitable for a detailed investigation of the rectangular sampling process. The gaussian assumption to which most of the analysis in Chapter III has been subjected, is not merely for analytical convenience since we are dealing with very large data size and the central limit theorem validates this assumption even when the process is not gaussian. On variability and covariability of estimates we were able to obtain closed form analytical expressions of the Blackman-Tukey type for the case when sampling intervals are Poisson distributed. These results, though approximate, are just as useful as any others derived in this thesis since they make use of the governing assumption of very large data size. Making use of these results we found further, that the Poisson sampling process achieves a smaller variability than the periodic sampling process for spectra which are very smooth, On the other hand, for rapidly varying spectra the Poisson is just as unreliable as the periodic scheme. However, this is not to say that the algorithm of Blackman and Tukey should be discarded as we will see presently. 4.2 Poisson Vs Periodic Sampling In the last section we inferred that for the same number of samples of the process the .Poisson sampling scheme achieves a smaller variance than the periodic scheme. This presumes that data are available for as long as -82- we want. In some cases of practical interest, the duration of data is limited somewhat, so that one is constrained to sample more closely to increase the number of samples or estimate from fewer record samples than needed. Already, we know that increasing the data size achieves smaller variance, so that it is quite possible that in a given time T we can obtain more samples by sampling periodically than with Poisson sampling. Even when signal duration is not limited, the question of how much longer one waits for say, N samples by sampling one way instead of the other, is worth looking into. Answers to these questions can be sought only in a probabilistic sense such as outlined below. Let TN = ai + 1 + ai + 2 + • • • + ai + N be the length of time it takes to obtain N data samples in a general non-uniform sampling scheme with mean sampling interval equal to a. uniformly every a sees., TN = N:X sees. For large totically normal with mean = E{Ti} cf- ~N(.n.) i.n. = N:X I .t\.=0 a2 N = aw ~ (.n.) .t\.=0 = N(N+l) a 2 and variance, = E(Tn2 ) - ~ 2 (Tn ) = N a 2 From statistical tables we have If sampling is done N1 TN is asymp- -83- i.e. Na - 3CX /N < T < Na + 3CX /N n with probability 0.99. Alternately, we write 3 Tn 3 1--< - < 1 + - /N IN N:l with probability 0.99. Now, for the case when the waiting time T, is fixed, the prob- ability of obtaining between N-K and N+K samples by Poisson sampling wi th mean rate 1/a is given by Pr{ 80 s NT .s 120} = 0.954319 Pr{ 70 4.3 S NT S 130} = 0.997057 Practical Implications of the Random Sampling Scheme In certain fields of application, prior statistical information on the sample times is not readily available, in which case, it may be necessary to record these sample times simultaneously. In sane cases this will require modification of the existing recording hardware and possibly more computation time. On the other hand, the scheme allows -84a greater flexibility in data acquisition and errors due to small jitter are less significant in this scheme than they are in the Blackman and Tukey scheme. In what follows we outline the relevance of random sampling in some areas of application. 4.3.1 Communications Technology In some coded communication systems the time of channel availability are random. By sampling a time sequence in a random manner corresponding to the random availability of the channel rather than sampling periodically, the need for buffer storage can be eliminated. 4. 3.2 Seismic Data Processing In future design of field experiments to investigate random seismic noise by means of correlation or spectral techniques, geophone arrays may now be set up randomly spaced in a manner that will optimize seismic data acquisition. It may also be possible to filter out, using appropriate geophone distributions, the propagating modes or coherent noise as is sometimes called. This will eliminate the need for delay lines as is the usual practice (cf [19]). 4.3.3 Oceanography . Whether it be swells from distant storms or ocean waves, the times of arrival are a random phenomenon and it will no doubt be more expedient and space-saving to take readings whenever they are available rather than periodically as is. usually done. In this application arrival times have been known to follow certain well known distributions • -85- 4.3.4 Structural Design In the design of tall structures and aircraft, destructive power of atmospheric turbulence is of very great concern. Engineers must design to accomodate oscillations due to wind gusts whose magnitude and direction are random phenomena with random occurence times. Estimation of spectra of wind velocities and of ground motion (velocity and acceleration) due to earthquakes normally preceeds actual structural design. -86Appendix A A-1 Behavior of ~(ill) for large = arctan ill a and write Recall that lim ~ ~ (ill) = 2 a m lim ~ ... lt/2 ill , n+l ~ cos ~ cos(n+l) ~ n=o 2 a ~l'~lt/ { cos 2~(1 + cos 2~) + cos 3~cos3~ 2 + cos ~ cos4~ + ........ ,+ cosm+ ~cos(m+l)~} 4 lim RHS ~ 2 a 1 1 im { 2 cos 4~ + cos 3 ~cos3~ + cos 4~cos4~' L ~ .... lt/2 J Expand the cosine functions noting that sinlt/ 2 = 1, coslt/ 2 and get lim ~(ill) ill ... co = 2 a .l im { 2 sin4€ + sin4 € cos4€ - sin3 €sin3€} € ... 0 = - 16 a € 6 Now, ~- € = tan-1 ill a ~ € ~ 1 I o:a so that = 0 etc., -87- ~ (w) = 0( -€ 6 ) Behavior of ~(w) A-2 PB ~ (w) = aa for large { m w • n+l ~ (1-i) cos ~ cos(n+l) ~ + cos 2~ n=l Note from A-1 that we can write Proceeding as in A-1 we write PB lim ~ (w) (l) ..... CXI = lim ~P (w) - -aa • (l) ..... CXI m lim ..... ~ { 2 - /2 !Cos ~cos2~ 1{ ' + 2'Cos 3cp-cos3~ + 3cos 4~.cos~ + •••••} Now, lim {cos ~cos2~ + 2cos3~cos3~ + 3Cos ~cos~} ~ ..... 7!/ 2 = 2 lim {cos (1f/ -e)cos(7!-2e) + 2 cos 3 2 € ..... 0 2 4 (¥- e)cos(~l! 3e) 4 + 3cos (1!/ -e)cos(27! - 4e)} 2 = 4 lim {-sin2e cos2e - 2 sin3esin3e + 3 sin e cos4e} € ..... 0 -88- 2 2 4 4 --- 2 - € (l-2€ + ••• ) - 6€ + 3€ (l-8€ + ••• ) = ~ 2 - € Thus, whence we say PB ~·(oo) = 0 ( l 2 2) ~m>ex A-3 ~(oo) ~(oo) &p max .1\. = 2CXR:l { ncx cos n+l r: (1-m) cp cos ( n+l )cp m n=o = 2 ex Re{ r: (l - :XHn+l)[ -cosncp cos (n+l) cp sincp = 0 when cp = 0 , ~ . 2 1 2 2~ 1 3~ } - (mtl) (mt2 ). ex cos~ cos (mt2) cp , -89- It can be deduced from above that < 0 q> = o, 2kn Further, we note that which implies max J\. = ~(0) -90Appendix B B-1 Permutations of ti, tk, tk+n' ti+j subject to: First, assume also that and so that only the following four permutations are possible (subject to (1) and (2) above): (3) ti ~ .tk < ti+j < tk+n (v) (4) ti ~ tk < tk+n ~ ti+j (vi) (5) tk < ti ~ ti+j < tk+n (iv) (6) tk < .ti < tk+n ~ ti+j (iii) Now relax restrictions in· (2) to have from (3) and (6) respectively, (7) ti < ti+j ~ tk < tk+n (ii) (8) tk < tk+n ~ ti < ti+j (i) An interval tree for the above permutations is illustrated in Figure Bl. Note that the permutations are unique only as far as absolute inequalities are concerned. - I . t ... ·- · t; < tk+n t Fig. B1 tk < ti-tj_ J t k !> t i JLk+n t > t. . > t i 2: t kj (v) I+J 1 lt L k+n > t k ->t·I+ J· >t-j I ( i i) - {ti+j ~ tk+n > tk ~ t; (vi) tk > ti J t k+n > ti +j > tk > td- (iv) tk~ti+j t; ~ tk 1I t·I+J· > tk+n > t; > tkj (iii) A representative interval tree for the permutat ions der ived in appendix B-1. t · .j < tk+n . L - I+ fi+j ~ t k+n . ti > tk -f t i+j > ti ~ tk+n > tkj ( i ) I I ~ \0 -92B-2 Derivation of . A(W' ,w+wr,w+wq) subject to inequalities (i) through (vi) of Section B-1 so that equation 3.25 becomes: A(i)(w' w+wr, w+w) 1 , q - 21 2 L Ne a R1 (i,j,k,n) Etraiak(~+n+ aonl)(ai+j+ aojl) l cos[2w' (ti-tk+n) + (2w'+w+wr) (tk+n-tk)- (w+wq) (ti+j-ti) 1 i 2w e (<~+n+ ~nl) e i(2w'+w+w) r Define: and ~*(n) and take expectation as indicated to get i E. k+n+l a ) l.1 1} -93A(i)(w' w+w w+w) 1 r' ' q ___,...;;1~ L Re N! a R (i,j,k,n) 1 • ~i k - -n (2w') ~ (2w'+w+w ) + ae 1 ~(2w'+w+w >) n r n r X Similarly, from equation 3.25a, obtain A(i)(w' w+w 2 ' r' w+w) q To obtain equation 3.27, follow the above procedure and derive A~t), A~i) i • (ii),•••,(vi), then note that equation 3.27 is given by 6 A(w',w+wr ,w+w) • q . L L i•l R ( i,j ,k,n ) 1 {A1 <w',w+wr ,w+wq )+A12 (w',w+wr ,w+wq)} 1 -94References 1. Blackman, R. B., and Tukey, J. W., The Measurement of Power Spectra From the Point of View of Communications Engineering, New York 1953. 2. Beutler, F. J., and Leneman, O.A.Z., "Random Sampling of Random Processes: Stationary Point Processes", Information and Control, 2 (1966), 325-346. 3. Beutler, F. J ., and Lenema.n, O.A.Z., "The Theory of Stationary Point Process", Acta Math., ll6 (1966). 4. Beutler, F. J., "Alias-Free Randomly Timed Sampling of Stochastic Processes", Information and Control Engineering Program, The University of Michigan- (to be published). 5. Beutler, F. J., and Lenema.n, O.A.Z., "Spectral Analysis of Impulse Process", Information and Control Engineering Program, The University of Michigan- (to be published). 6. Jenkins, G. M., "General Considerations in the Analysis of Spectra", Technometrics, 3, No. 2, May 1961. 7. Rice, S. o., "Mathematical Analysis of Random Noise", from Selected Papers on Noise and Stochastic Processes, Nelson Wax, Editor; Dover Publications, Inc., New York, N. Y. 1954. 8. Parzen, E., "Mathematical Considerations in the Estimation of Spectra", Technometrics 2_, No. 2, May 1961. 9. Shapiro, H. S ., and Silverman, R. A., "Alias-Free Sampling of Random Noise", J. Soc. Indus, Appl. Math., 8, (1960), 225-248. 10. Lorena, C. s., "Recovery of Randomly Sampled Time Sequences", IRE Trans. · on Comm. Sys., June 1962 . ll. Pa.nofsky, H. A., "Meteorological Applications of Power Spectrum", -95- Bulletin American Meteorological Society, 36, No.4, April 1955. 12. Levin, M. J.: "Power Spectrum Parameter Estimation", PGIT, Jan. '65. 13. Bello, Phillip: "Joint Estimation of Delay, Doppler, and Doppler Rate", PGIT, June 1960. 14. Mathematical Methods of Statistics , Princeton Cramer, Harold: Mathematical Series, No. 9. 15. Lighthill, M. J.: Fourier Analysis and Generalized Functions, Cambridge University Press, 1964. 16. Feller, William: Probability Theory and its Applications , Vol.I, John Wiley and Sons, Inc., New York, N.Y. 1950. 17. Davenport, W. B., and Root, W. L., Random Signal and Noise, McGraw-Hill Book Co., Inc., New York, 1958. 18. Whittle, P. "Estimation and Information in Stationary Time Series", Arkiv for Mathematic Ban 2 nr 23, 1952 •. 19. Adegbola, M. 0.: "A Theoretical Model for Seismic Noise", Basic Geophysics Division, Esso Production Research Co., unpublished Report (Summer, 1967). 20. Zadeh, L. A.: "Correlation Functions Power Spectra in Variable Networks", Proceedings of the I.R.E., Nov. 1950. 21. Gradshteyn, I. s., and Ryzhik, I.M., Tables of Integrals, Series and Products , Academic Press, New York and London, 1965. 22. Todd, John: A Survey of Numerical Anal.ysis , McGraw-Hill Book Co. Inc., New York, 1962 . 23. Parzen, E., dtochastic Processes , Holden-Day, Inc., San Francisco 1962. 24. Molina, E. c., Poisson's Exponential Binomial Limit , D. Van Nostrand Co., Inc., New York, 1942. -9625. Apostol, T.M., }fathematical Analysis, Addison-Wesley Publishing Co., Inc. Palo Alt, 1964.