Estimation of Ambiguity Functions With Limited Spread

DEP AR TMENT OF ST A TISTICAL SCIENCE RESEARCH REPOR T 293 1 Estimation of Ambiguity Functions W ith Limited Spread Heidi Hindberg & Sofia C. Olhede Abstract This paper proposes a new estimation procedure for the ambiguity function of a non-stationary time series. The stochastic properties of the empirical ambiguity function calculated from a single sample in time are deriv ed. Different thresholding procedures are introduced for the estimation of the ambiguity function. Such estimation methods are suitable if the ambiguity function is only non-negligible in a limited region of the ambiguity plane. The thresholds of the procedures are formally derived for each point in the plane, and methods for the estimation of nuisance parameters that the thresholds depend on are proposed. The estimation method is tested on several signals, and reductions in mean square error when estimating the ambiguity function by factors of ov er a hundred are obtained. An estimator of the spread of the ambiguity function is proposed. Index T erms Ambiguity function, estimation, harmonizable process, non-stationary process, thresholding, underspread process. I . I N T RO D U C T I O N W e propose a new estimation procedure for the Ambiguity Function (AF) of a non-stationary Gaussian process. The suitable definition of a ‘locally stationary process’, and the development of inference methods for a given sample from a non-stationary process is still an open question, see [1]. The AF forms an essential component to this task, as it provides a characterisation of dependency between a given time series and its translates in time and frequency . The AF of a non-stationary process can often be assumed to be mainly limited in support to a small region of the ambiguity plane, and this corresponds to correlation in the process of interest being limited in support. If the support is additionally centred at time and frequency lags of zero, then the process is underspread, see Matz and Hlawatsch [2]. An important class of non-stationary processes, namely semi-stationary processes [3], exhibit limited essential spread in the ambiguity plane. Non-stationary processes can also be constructed from the viewpoint of time-variant linear filtering. Pfander and W alnut [4] hav e shown that if the spreading of a filtering operator is sufficiently limited then the operator may be identified from its measurements. It is not necessary in this case to assume that the spread is centered at the origin [4]. The key to the estimation or characterization of the generating mechanism of non-stationary processes, in all of these cases, is the compression of the AF of the process. The AF has also been popularly used in radar and sonar applications [5], where the echo of a known signal is recorded. The delay and frequency shift of the echo can be determined from the AF of the signal [6]. In such applications the emitted signal is deterministic, ev en if the echo has been contaminated by noise. By using the compression of the emitted signal in the AF domain, the properties of the process of reflection can be determined, see Ma and Goh [7]. The estimation of the AF of a non-zero mean signal immersed in noise is of interest. Unlike Ma and Goh, we shall not assume the compression of the AF is known, but propose an estimator of the AF , suitable if the AF satisfies some (unknown) compression constraints. Giv en the importance of the AF it is surprising to find that existing methods for its estimation are quite naiv e. Methods are usually based on calculating the EMpirical AF (EMAF), or some smoothed version of the EMAF . It is clearly unpalatable to implement a smoothing procedure on the EMAF , as the compression of the EMAF should be preserved, or ev en preferably increased, by any proposed estimation procedure. In this spirit Jachan et al. [8] have used shrinkage methods to estimate the EMAF with an assumed knowledge of its support, using a multiplier that attenuates larger local time and frequency lags. Often the support of the AF is not known, and arbitrary shrinkage at larger lags will not be an admissible estimation procedure. Giv en the compression or sparsity of the AF , shrinkage or threshold estimators are a natural choice. In this paper we propose a set of threshold estimators for the AF . Thresholding has been implemented before in the global time-frequency plane for wa velet spectra, or ev olutionary spectra, see Fryzle wicz and Nason [9], or von Sachs and Schneider [10]. The estimation H. Hindberg thanks the Research Council of Norway for financial support under project 162831/V00. H. Hindberg is with the Department of Physics and T echnology , Univ ersity of T romsø, NO-9037 T romsø, Norway (heidih@phys.uit.no) T el:+4777645118, Fax+4777645580. S. C. Olhede is with the Department of Statistical Science, Univ ersity College London, Gower Street, London, WC1 E6BT , UK (s.olhede@ucl.ac.uk). T el: +44 (0)20 7679 8321, Fax: +44 (0)20 7383 4703. DEP AR TMENT OF ST A TISTICAL SCIENCE RESEARCH REPOR T 293 2 problem in this case is quite different to estimating the AF , as we cannot assume the raw wavelet/e volutionary spectra exhibit compression. Hedges and Suter [11] calculated the spread of the AF , based on exceedances of the magnitude of the EMAF ov er a (fixed and user tuned) threshold in a giv en direction in the local time-frequency plane. The stochastic properties of the EMAF were not treated by Hedges and Suter , and instead their in vestigation focused on some effects from different boundary treatments. T o be able to select a suitable threshold procedure at each local time-frequency point we must establish the distribution of the EMAF , and this requires modelling the observed process. W e discuss the classes of non-stationary processes that we intend to treat in Section II-A, and initial method of moments estimators of second order structure in Section II-B. In Section III we determine the first and second order properties of the EMAF . For strictly underspread processes satisfying constraints on their degree of uniformity we show that the EMAF is asymptotically Gaussian proper , if calculated for local frequency and time lags outside the support of the AF . W e introduce the proposed estimation procedure in Section IV for deterministic signals immersed in noise and stochastic processes. In the first estimation procedure we threshold the EMAF based on a variance estimated from the entire plane, this yielding the Thresholded EMAF (TEAF). F or stochastic processes, we also introduce a Local TEAF (L TEAF), where we estimate the variance locally in regions of the ambiguity plane. Such procedures are suitable as smooth time-frequency spectra will correspond to decaying ambiguity functions and are a natural extension to treating wa velet coefficients differently depending on the level of the wa velet transform, so called lev el-dependent thresholding, see e.g. Johnstone and Silverman [12]. The EMAF of deterministic signals immersed in analytic white noise will be biased. W e propose a method to remove this bias prior to local thresholding, thus obtaining the Bias-corrected L TEAF (BL TEAF). W e define an estimator of the total spread of a signal in the ambiguity plane in Section IV -A, to numerically measure the support of the signal in the ambiguity plane. The proposed TEAFs will produce sparse approximations to processes that are not strictly underspread. Whenev er the EMAF is small in comparison to the v ariance of the EMAF , the AF will be estimated as zero. The AF can provide essential information as to the natural bandwidth of a process via the spread of the TEAF . In Section V we provide simulation studies of the proposed estimators. The Mean Square Errors (MSE) of our procedures reduce in many cases to around a hundredth of the MSE of the EMAF . Plots of single realisations show the accuracy of the threshold estimators, based on a single sample. The proposed methodology thus introduces a new method of characterising the features of structured non-stationary processes with high accuracy . I I . S E C O N D O R D E R S T R U C T U R E D E S C R I P T I O N S A. Modelling the Observations W e assume that the process under consideration, X [ t ] , is an analytic, zero-mean, discrete-time, Gaussian harmonizable random process. Note that [ · ] is used to indicate a discrete argument, and ( · ) is used for a continuous argument. As X [ t ] is harmonizable it admits the spectral representation [13] of X [ t ] = Z 1 2 0 e j 2 π tf d e X ( f ) , t ∈ Z , (1) where n d e X ( f ) o is a Gaussian increment process, and X [ t ] ∈ C . W e note that n d e X ( f ) o has a correlation structure specified by the dual-frequency spectrum S X X ∗ ( ν, f ) dν d f = E n d e X ( ν + f ) d e X ∗ ( f ) o . The AF of X [ t ] is defined as a Fourier T ransform (FT) pair with the dual-frequency spectrum in v ariable f . The AF forms an FT pair with the dual-time second moment of the process, or M X X ∗ [ t, τ ] = E { X [ t ] X ∗ [ t − τ ] } , for t and τ ∈ Z , now in variable t . W e refer to [14], [15] for a more detailed exposition of the AF and forms of AFs of common processes. Estimation of the AF of a real-valued process is known to display artifacts from aliasing and interference from negati ve frequencies [16]. For this reason we consider analytic processes solely . Note that since we are working with the analytic signal of a non-stationary process, it may (in general) be improper [17] and thus may exhibit complementary correlation [18]. W e leav e the estimation of the complimentary AF outside the scope of this paper , noting that our proposed methods can in a straight forward fashion be extended to include the estimation of such objects. Furthermore, we will in this paper assume that we are working with proper processes, and any terms that will only be non-zero for improper processes will be omitted. B. Initial Estimation W e observe x [ t ] , a finite sample of a realization of X [ t ] . The Sample AF (SAF) of a sampled deterministic signal, { g [ t ] } N − 1 t =0 , is giv en by A g g ∗ ( ν, τ ] = N − 1+min(0 ,τ ) X t =max(0 ,τ ) g [ t ] g ∗ [ t − τ ] e − j 2 π ν t , ν ∈  − 1 2 , 1 2  , τ = [ − ( N − 1) , . . . , ( N − 1)] . (2) DEP AR TMENT OF ST A TISTICAL SCIENCE RESEARCH REPOR T 293 3 The properties of A g g ∗ ( ν, τ ] follow directly from properties of AFs of deterministic sequences, see [19], [20]. The ef fects of discretization on deterministic structure, i.e. trying to infer properties of g ( t ) from g [ t ] ≡ g ( t ∆ t ) for some appropriate sampling period ∆ t > 0 , will also be left outside the scope of this article [21], [22]. The Empirical Second Moment (ESM) of a stochastic process based on the realization x [ t ] , t = 0 , . . . , N − 1 is giv en by c M X X ∗ [ t, τ ] = x [ t ] x ∗ [ t − τ ] , which has an expectation E n c M X X ∗ [ t, τ ] o = M X X ∗ [ t, τ ] and a v ariance v ar n c M X X ∗ [ t, τ ] o = M X X ∗ [ t, 0] M ∗ X X ∗ [ t − τ , 0] . Note that the ESM, its expected value and variance is defined as above for 0 ≤ t ≤ N − 1 and t − ( N − 1) ≤ τ ≤ t , and they are zero otherwise. Clearly the variance of c M X X ∗ [ t, τ ] is very large and the estimator should not be put to direct use. W e form the EMAF of X [ t ] by b A X X ∗ ( ν, τ ] = N − 1+min(0 ,τ ) X t =max(0 ,τ ) c M X X ∗ [ t, τ ] e − j 2 π ν t = N − 1+min(0 ,τ ) X t =max(0 ,τ ) x [ t ] x ∗ [ t − τ ] e − j 2 π ν t . (3) The EMAF corresponds to an estimator of the AF of X [ t ] . Howe ver , it is not a suitable estimator of A X X ∗ ( ν, τ ] , as its variance will be extremely large. When implementing (3) we have chosen zero-padding as our choice of e xtending the estimated cov ariance matrix. W e found that periodic (circular) extension created unpalatable mixing effects of a fundamentally non-stationary object. Edge-effects are discussed more carefully in [11]. The best choice of edge treatment will depend on the statistical and deterministic properties of the chosen EMAF . I I I . P RO P E RT I E S O F T H E E M P I R I C A L A M B I G U I T Y F U N C T I O N A. Repr esentation and Moments of the EMAF W e shall construct estimators of the AF , based on the EMAF . The EMAF of an analytic zero-mean harmonizable process X [ t ] admits the representation of b A X X ∗ ( ν, τ ] = Z 1 / 2 0 Z 1 / 2 0 e j π ( ν 0 − ν )( N + τ − 1) e j 2 π f τ D N −| τ | ( ν 0 − ν ) d e X ( ν 0 + f ) d e X ∗ ( f ) , where D N ( f ) = sin( πf N ) / sin( π f ) . W e respectiv ely denote µ A ( ν, τ ] , σ 2 A ( ν, τ ] and r A ( ν, τ ] as the mean, variance and relation of the EMAF . W e denote an analytic process corresponding to a real-valued stationary white process as an analytic white process, but keep in mind that this process is not actually white. Thus, an analytic white process is a stationary process that has a power spectral density that is constant for nonnegati ve frequencies, and zero for negati ve frequencies. Pr oposition 3.1: Moments of the EMAF for noisy Deterministic Signals For a process that is the sum of a deterministic analytic signal and analytic white noise with v ariance σ 2 W , X [ t ] = g [ t ] + W [ t ] , the expected value, variance and relation of the EMAF take the form µ A ( ν, τ ] = b A g g ∗ ( ν, τ ] + σ 2 W 2 e − j π ν ( N + τ − 1) D N −| τ | ( ν ) e j π τ / 2 sinc( τ / 2) (4) σ 2 A ( ν, τ ] = σ 2 W h ( ν, τ ] + σ 2 W h ( − ν, − τ ] + σ 4 W ( N − | τ | )(1 / 2 − | ν | ) + O (1) (5) r A ( ν, τ ] = − σ 4 W | ν | ( N − | τ | ) L ( N − τ , ν ) e − 2 j π ν ( N − 1) sinc(2 | ν | τ ) + 2 σ 2 W h 0 ( ν, τ ] + O (1) L ( N − | τ | , ν ) = Z ∞ −∞ sinc( f )sinc( f + 2( N − | τ | ) ν ) d f , h ( ν, τ ] = Z 1 / 2 − max(0 ,ν ) max( − ν, 0) | G ( f , τ ] | 2 d f h 0 ( ν, τ ] = Z 1 / 2+min(0 , − ν ) max(0 , − ν ) G ∗ ( f , τ ] G ( f + 2 ν, − τ ] e j 2 π ( f − sign ( τ ) ν ) τ d f where sinc( x ) = sin( π x ) / ( π x ) and G ( f , τ ] = N − 1 −| τ | P t =0 g [ t − | τ | I ( τ < 0)] e − j 2 π f t . The indicator function is defined as I ( τ < 0) = 1 if τ < 0 and I ( τ < 0) = 0 otherwise. Pr oof: See appendix A. Pr oposition 3.2: Moments of the EMAF for a Stationary Process For a zero-mean, analytic stationary stochastic process X [ t ] the expected value, variance and relation of the EMAF take the form µ A ( ν, τ ] = D N −| τ | ( ν ) e − j π ν ( N + τ − 1) f M X X ∗ [ τ ] (6) σ 2 A ( ν, τ ] = ( N − | τ | )(1 / 2 − | ν | ) A X X ∗ ( − ν, 0] + O (1) (7) r A ( ν, τ ] = e − 2 j π ν ( N + τ − 1) ( N − | τ | )(1 / 2 − | ν | ) L ( N − | τ | , ν ) A X X ∗ ( ν, τ ] + O (1) (8) DEP AR TMENT OF ST A TISTICAL SCIENCE RESEARCH REPOR T 293 4 where f M X X ∗ [ τ ] is the autocorrelation sequence of the process, e S X X ∗ ( f ) is the spectral density of the process, and A X X ∗ ( ν, τ ] = Z 1 / 2+min( ν, 0) max(0 , − ν ) e S X X ∗ ( f − ν ) e S X X ∗ ( f ) e j 4 π f τ d f / (1 / 2 − | ν | ) . (9) Pr oof: See appendix B. The AF of a stationary process will be nonzero only on the stationary manifold ν = 0 , and takes the form A X X ∗ ( ν, τ ] = f M X X ∗ [ τ ] δ ( ν ) . For a fixed value of N , D N −| τ | ( ν ) does not correspond to a delta function in ν . Ideally howe ver the estimated EMAF should not exhibit spreading in ν . Pr oposition 3.3: Moments of the EMAF for a Uniformly Modulated White Noise Process For a zero-mean, analytic process X [ t ] corresponding to a real-valued uniformly modulated white noise process, the expected value, variance and relation of the EMAF take the form µ A ( ν, τ ] = (1 / 2 − | ν | )Σ X X ∗ ( ν ) e j π (1 / 2 −| ν | ) τ sinc((1 / 2 − | ν | ) τ ) + O (1) (10) σ 2 A ( ν, τ ] = (1 / 2 − | ν | )( N − | τ | ) Z 1 / 2 −| ν | − 1 / 2+ | ν | | Σ X X ∗ ( f ) | 2 N − | τ | e j 2 π f τ d f + O (1) (11) r A ( ν, τ ] = e − 4 j π ν τ Z 1 / 2+min(0 ,ν ) max(0 ,ν ) Z 1 / 2+min(0 ,ν ) max(0 ,ν ) Σ X X ∗ ( f − α + ν )Σ ∗ X X ∗ ( f − ν − α ) e j 2 π ( f + α ) τ d f dα + O (1) . where Σ X X ∗ ( ν ) is the DFT of the time-varying v ariance of X [ t ] , σ 2 X [ t ] . Pr oof: See appendix C. W e note that the result of the integral in the expression for the variance decreases with | ν | and also with | τ | . B. Distribution of the EMAF of an Underspr ead Pr ocess T o determine suitable estimation procedures for the AF we need to deriv e the distribution of the EMAF . Theor em 3.1: Distribution of the EMAF of an underspread process Assume X [ t ] is the analytic process corresponding to a harmonizable real-v alued zero-mean Gaussian process whose EMAF is strictly underspread. A strictly underspread process has an ambiguity function that is only non-zero for ( ν, τ ] ∈ D , where there exists some finite non-negati ve T and Ω such that D ⊂ [ − Ω , Ω] × [ − T , T ] . Then the EMAF of X [ t ] ev aluated at ( ν , τ ] has a mean, variance and relation of µ A ( ν, τ ] = 1 / 2 Z 0 1 / 2 − f Z − f S X X ∗ ( ν 0 , f ) e j 2 π f τ e j π ( N + τ − 1)( ν 0 − ν ) D N −| τ | ( ν 0 − ν ) dν 0 d f (12) σ 2 A ( ν, τ ] = N − 1+min(0 ,τ ) X t =max(0 ,τ ) T − 1 X τ 0 = − ( T − 1) e − j 2 π ν τ 0 M X X ∗ [ t, τ 0 ] M ∗ X X ∗ [ t − τ , τ 0 ] + O  log  N T  (13) r A,N ( ν, τ ] = ( O  log  N T  if | τ | > T 0 if | ν | > Ω , (14) while if | ν | < Ω and | τ | < T then the relation in noted in (A-8). Fix ( ν, τ ) / ∈ D . Let τ ≥ 0 and take Q N −| τ | [ t ] = e − j 2 π ν ( t + τ ) ( X [ t + τ ] X ∗ [ t ] − M X X ∗ [ t + τ , τ ]) , t = 0 , . . . , N − 1 − τ . (15) Assume that the triangular array  Q N −| τ | [ t ] , t = 0 , . . . , N − 1 − τ  is strongly mixing and satisfies for any  > 0 1 σ 2 A ( ν, τ ] N − τ − 1 X t =0 v ar { Q N −| τ | [ t ] } I  | Q N −| τ | [ t ] | > σ A ( ν, τ ]  − → 0 , (16) as N − | τ | → ∞ , as well as the condition: sup N 1 σ 2 A ( ν, τ ] N − τ − 1 X t =0 M X X ∗ [ t + τ , 0] M X X ∗ [ t, 0] < ∞ . (17) Then, with L − → denoting con ver gence in law [23], b A X X ∗ ( ν, τ ] − µ A ( ν, τ ] σ A ( ν, τ ] L − → N C (0 , 1 , 0) , ( ν , τ ) / ∈ D , t ≥ 0 , (18) DEP AR TMENT OF ST A TISTICAL SCIENCE RESEARCH REPOR T 293 5 where N C (0 , 1 , 0) denotes a complex Gaussian distribution with mean 0, variance 1 and relation 0. The result follows mutatis mutandis for τ < 0 . Pr oof: See appendix D. It is tempting to attempt to deduce the distributional results directly from the first and second order structure established in the first part of the theorem. In fact, if X [ t ] corresponds to a stationary process, where T and τ are both ev en, the result follows directly . For more general classes of processes the result is slightly more in volved, and a condition on the variance needs to be combined with a suitable mixing assumption, as abov e. Note that the CL T theorem that we use does not provide rates of conv ergence. Also, in some degenerate cases the joint distribution of n b A X X ∗ ( ν, τ ] o ( ν,τ ) may not be asymptotically multiv ariate Gaussian. For ease of exposition, such cases are not treated here. Finally , whilst (asymptotically) retrieving the Gaussian distribution for the EMAF by Theorem 3.1, we fail to obtain simple interpretable forms for σ 2 A ( ν, τ ] and r A ( ν, τ ] . For this reason, it is con venient to derive the first and second order structure directly from the modelling assumptions of some commonly used processes, as done in Section III. I V . E S T I M A T I O N P R O C E D U R E T o determine the approximate distribution of the EMAF of a deterministic signal immersed in analytic white noise, we need to note the distributions of the four components of (A-1). Here, A g g ∗ ( ν, τ ] is constant, whilst A W W ∗ ( ν, τ ] is asymptotically Gaussian, and its form is determined by Theorem 3.1. W e note that since W [ t ] is Gaussian analytic white noise, and so from (A-3) and (A-6) A W g ∗ ( ν, τ ] + A g W ∗ ( ν, τ ] d = N C  0 , σ 2 W ( h ( ν, τ ] + h ( − ν , − τ ]) , 2 σ 2 W h 0 ( ν, τ ]  . From (A-2) and (A-5) we can note that A W g ∗ ( ν, τ ] + A g W ∗ ( ν, τ ] is independent of A W W ∗ ( ν, τ ] . Combining (A-1) and Theorem 3.1, as well as using proposition 3.1, it follows for ( ν, τ ) 6 = (0 , 0) b A X X ∗ ( ν, τ ] − µ A ( ν, τ ] σ A ( ν, τ ] L − → N C  0 , 1 , r A ( ν, τ ] /σ 2 A ( ν, τ ]  (19) (20) Hence we may note | µ A ( ν, τ ] | 2  σ 2 A ( ν, τ ] in most of the ambiguity plane. Furthermore it follows that    b A X X ∗ ( ν, τ ] − µ A ( ν, τ ]    2 σ 2 A ( ν, τ ] d = 1 2  1 + | r A ( ν, τ ] | σ 2 A ( ν, τ ]  U 2 1 + 1 2  1 − | r A ( ν, τ ] | σ 2 A ( ν, τ ]  U 2 2 + o (1) , (21) where U 1 and U 2 are independent standard Gaussian random variables. Thus for such points in the ambiguity plane it follows that    b A X X ∗ ( ν, τ ] − µ A ( ν, τ ]    2 is a weighted sum of χ 2 1 ’ s. These are equally weighted if r A ( ν, τ ] = 0 . For most points in the ambiguity plane this statement will be roughly valid as apparent from proposition 3.1. Then if ( ν, τ ] / ∈ D , where D denotes the support region of the ambiguity function of the real part of X [ t ] , we may note that    b A X X ∗ ( ν, τ ] − µ A ( ν, τ ]    2 σ 2 A ( ν, τ ] d = " 1 2 χ 2 2 + O | log( N ) | 2 N − | τ | !# + o (1) . (22) A suitable estimation procedure for such random variables with σ 2 A ( ν, τ ] known a priori, is then the following thresholding procedure, for some giv en threshold λ 2 , b A (ht) X X ∗ ( ν, τ ] =      b A X X ∗ ( ν, τ ] if    b A X X ∗ ( ν, τ ]    2 > λ 2 σ 2 A ( ν, τ ] 0 if    b A X X ∗ ( ν, τ ]    2 ≤ λ 2 σ 2 A ( ν, τ ] . (23) If we normalize     b A X X ∗ ( ν, τ ]    2  via di viding the sequence, element by element, by σ 2 A ( ν, τ ] , then we retriev e a set of corr elated positi ve random variables. For any such collection, we may note from [24], that using a threshold λ 2 N X ( C ) = 2 log( N X [log( N X )] C ) for C ≥ 1 is suitable. In our example N X is chosen to be twice the number of observ ations, as we threshold the ambiguity function frequency by frequency , across the total collection of all time lags. Olhede [24][p. 1529] has calculated the risk of this non-linear estimator for sums of unequally weighted χ 2 1 ’ s. Of course σ 2 A ( ν, τ ] is not known. W e standardize the EMAF by b A ( S ) X X ∗ ( ν, τ ] = b A X X ∗ ( ν, τ ] p ( N − | τ | )(1 / 2 − | ν | ) , (24) DEP AR TMENT OF ST A TISTICAL SCIENCE RESEARCH REPOR T 293 6 and note that v ar n b A ( S ) X X ∗ ( ν, τ ] o = σ 4 W + σ 2 W h ( ν, τ ] + h ( − ν , − τ ] ( N − | τ | )(1 / 2 − | ν | ) + O  1 N − | τ |  . (25) W e ha ve assumed that g [ t ] ∈ ` 2 , and so | h ( ν, τ ] | = O (1) . Thus σ 2 W h ( ν,τ ]+ h ( − ν, − τ ] ( N −| τ | )(1 / 2 −| ν | ) = o (1) , where the form of h ( ν, τ ] at high values of | ν | and | τ | ensure that the statement is still true for such values. W e shall use the Median Absolute Deviation (MAD) estimator of the v ariance. MAD has been used for estimating the scale of correlated data before, see [12]. The MAD estimator will need an adjustment factor that is different for 1 2 χ 2 2 random variables compared to χ 2 1 random variables. W e note the median of a 1 2 χ 2 2 is ln[2] . W e define an estimator of the analytic white noise variance for any region { ( ν, τ ] } = M ⊂ [ − 1 / 2 , 1 / 2] × {− ( N − 1) , . . . , ( N − 1) } by b σ 4 W ( M ) = median     b A ( S ) X X ∗ ( ν, τ ]    2  ( ν,τ ) ∈M ln[2] . (26) The imprecision of this procedure will depend on the lack of compression of the representation of { X [ t ] } in the ambiguity domain. W e note that MAD has a breakdown point of 50 %, and so with quite sev ere contamination the estimator will still be useful, if somewhat inefficient. W e then take b σ 2 A ( ν, τ |M ] = b σ 4 W ( M )( N − | τ | )(1 / 2 − | ν | ) . (27) A suitable threshold procedure (now with an unknown σ 2 A ( ν, τ ] ) simply corresponds to using (23) with the variance replaced by its estimated value from (27). If the variance is estimated from the entire plane, i.e., M = {− N / (2 N ) , . . . , N / (2 N ) } × {− ( N − 1) , . . . , ( N − 1) } , we denote b A (ht) X X ∗ ( ν, τ ] the Thresholded EMAF (TEAF). The EMAF of a deterministic signal immersed in analytic white noise is biased, see proposition 3.1. Giv en an estimator of σ 2 W , we can remove this bias prior to thresholding. W e define the bias-corrected EMAF by b A ( B ) X X ∗ ( ν, τ ] = b A X X ∗ ( ν, τ ] − 1 2 b σ 2 W ( M 1 ) e − j π ν ( N + τ − 1) D N −| τ | ( ν ) e j π τ / 2 sinc ( τ / 2) . (28) The noise variance is estimated from some gi ven region M 1 , which is chosen to only include the rims of the ambiguity plane. W e normalize b A ( B ) X X ∗ ( ν, τ ] as in (24). For smaller sample sizes, it may be unreasonable to assume that the contributions of h ( ν, τ ] are negligible compared to p N − | τ | . W e define σ 4 X ( ν, τ ] = σ 2 A ( ν, τ ] / [( N − | τ | )(1 / 2 − | ν | )] . This is explicitly subscripted by X , as it is only non-constant from contributions due to h ( ν, τ ] . W e propose to segment the plane into regions, and implement the proposed procedure after taking a separate estimate of σ 4 X ( ν, τ ] in each region. The optimal choice of region is given by (25). Since h ( ν, τ ] is unknown, we propose to use a centre square around (0 , 0] , and square annuli, see Fig. 1(a), to estimate the local variance. This is based on assuming σ 2 A ( ν, τ ] smooth and decaying from (0 , 0) . In general we separate the ambiguity plane into K regions {M k } , so that σ 4 X ( ν, τ ] ≈ σ 4 X ( M k ] , a constant for all ( ν , τ ] ∈ M k . W e estimate the variance in each region as b σ 4 X ( M k ) = median     b A ( S ) X X ∗ ( ν, τ ]    2  ( ν,τ ) ∈M k ln[2] . (29) W e threshold b A ( B ) X X ∗ ( ν, τ ] according to (23), but now using b σ 4 X ( M k ) in each defined region instead of b σ 4 W ( M ) in (27), this yielding b A (lbht) X X ∗ ( ν, τ ] , or the Local Bias-corrected TEAF (LBTEAF). Next, we deriv e results for estimating the EMAF of an underspread zero-mean stochastic process. This is strictly speaking a correct treatment for processes satisfying the constraints of Theorem 3.1, but we expect that the distributional result is valid under less constrained scenarios. The conditions are suf ficient, but by no means necessary , for the distributional result to hold. W e note from Theorem 3.1 that for most points in the ambiguity plane    b A X X ∗ ( ν, τ ] − µ A ( ν, τ ]    2 σ 2 A ( ν, τ ] d = 1 2 χ 2 2 + O | log( N ) | 2 N − | τ | ! + o (1) . (30) The question then arises of yet again estimating σ 2 A ( ν, τ ] appropriately . W e note from Theorem 3.1 that σ 2 A ( ν, τ ) takes different forms depending on the spreading of the process. The more X [ t ] is like an analytic white noise process, the less variation will σ 2 A ( ν, τ ] exhibit from the form of the variance of the EMAF of an analytic white noise process. W e again define σ 4 X ( ν, τ ] = σ 2 A ( ν, τ ] / [( N − | τ | )(1 / 2 − | ν | )] . Then    b A X X ∗ ( ν, τ ] − µ A ( ν, τ ]    2 ( N − | τ | )(1 / 2 − | ν | ) d = σ 4 X ( ν, τ ] " 1 2 χ 2 2 + O | log( N ) | 2 N − | τ | ! + o (1) # . (31) DEP AR TMENT OF ST A TISTICAL SCIENCE RESEARCH REPOR T 293 7 ! " ! 0.4 ! 0.2 0 0.2 0.4 ! 200 ! 100 0 100 200 ! 0.# ! 0.2 0 0.2 0.# ! 200 ! 100 0 100 200 ! " Fig. 1. (a) Partitioning of the ambiguity plane. (b) Radial lines over which the integrand in (32) exhibits stationary phase. Using in verse FTs we note from Theorem 3.1 that the variance of the EMAF of a strictly underspread process can in fact be rewritten for τ > 0 σ 2 A ( ν, τ ] = T − 1 X τ 0 = − ( T − 1) Z Ω − Ω Z Ω − Ω e − j 2 π ( ν τ 0 − ν 0 τ ) A X X ∗ ( ν 0 , τ 0 ] A ∗ X X ∗ ( ν 00 , τ 0 ] e j π ( N − τ − 1)( ν 0 − ν 00 ) D N − τ ( ν 0 − ν 00 ) dν 0 dν 00 + O  log  N T  = T − 1 X τ 0 = − ( T − 1) Z Ω − Ω e − j 2 π ( ν τ 0 − ν 0 τ ) | A X X ∗ ( ν 0 , τ 0 ] | 2 dν 0 + O  log  N T  + O (1) . (32) Deriving this requires that A X X ∗ ( ν, τ ] is sufficiently smooth in ν . If this is not the case, use (13). The v ariance of the EMAF corresponds to aggregating the total magnitude squared of the AF over the plane and implementing phase shifts. If the function is smooth in ν with a stationary phase approximation to the integral we mainly pick up contributions on the line ν 0 = ν τ 0 /τ , which we later sum over τ 0 . This is exemplified for some choices of ν and τ in Figure 1(b), where we plot the lines summed ov er . Despite this geometrical intuition, the form in (32) is not extremely informativ e. Equation (32) does indicate that σ 2 A ( ν, τ ] should be a smooth function of ν and τ , as we only integrate over D . W e rely on the forms of propositions 3.2 and 3.3 to giv e more precise understanding of the variance of the EMAF in these special cases. W e note from these set of results that the variance is smooth in ν and τ , and often exhibits a decay in ( ν, τ ) that resembles that of the variance of analytic white noise. If σ 4 X ( ν, τ ] is exactly constant across the ambiguity plane we can propose an estimator of σ 4 X ( ν, τ ] given by equating b σ 4 X ( ν, τ ] to the estimator from (26). Propositions 3.2 and 3.3 argue that in regions of the ambiguity plane σ 4 X ( ν, τ ] is smooth, and can be approximated as constant over a given region, again using (26). Noting that as the process is underspread | µ A ( ν, τ ] |  σ A ( ν, τ ] for most ( ν, τ ] , and so thresholding is an admissible estimation procedure. W e therefore propose a Local TEAF (L TEAF), by estimating the variance in giv en regions, just like the LBTEAF , but without removing the bias. W e will divide the ambiguity plane into regions according to Fig. 1(a) to obtain both the L TEAF and LBTEAF . Using separate regions in the local time- frequency domain to estimate the variance of the EMAF is a natural procedure also for stochastic processes. It is reasonable to assume that the time-varying spectrum is smooth in global frequency and global time. If the time-varying spectrum is a member of a Sobolev space of order k , then the modulus of the Fourier transform must decay faster than ν − ( k +1) in local frequency , and τ − ( k +1) in local time, see for example [25]. It would therefore be natural to use regions defined in terms of distance from ( ν, τ ) = (0 , 0) as M k , but to more naturally meld with digital sampling we propose to use the set of square annuli, similar to sampling discussed in [26]. Johnstone and Silverman in a similar spirit [12] proposed level dependent thresholding using a different variance for each wav elet le vel, in 1-D. As a final note the proposed estimators may not correspond to valid AFs. Only AFs that are the FTs of valid autocovariance sequences are valid. W e do not think this corresponds to a major flaw of the procedure, as we are mainly concerned with functions whose support will inform us of the resolution of the autocov ariance. A. Estimated Spread Using thresholding methods, we hav e produced an estimate of the AF in the entire ambiguity plane. W e also propose an estimator of the spread of the estimated AF , see [2], [11]. Note that for processes that are not strictly underspread, our thresholding procedure will identify regions where the mean of the EMAF is sufficiently distinct in magnitude from the variance DEP AR TMENT OF ST A TISTICAL SCIENCE RESEARCH REPOR T 293 8 of the EMAF . This in essence corresponds to comparing the magnitude of the AF at a point, with the magnitude of the AF at other points, see (32). W e additionally note from (32) that if the AF is mainly centered at other ( ν 0 , τ 0 ] then the variance of the EMAF will be large compared to the mean of the EMAF . In this instance the contribution at ( ν, τ ] is not significantly contributing to the structure of X [ t ] . W e define the ambiguity indicator cell to be ζ X ( ν, τ ] = I ( | A X X ∗ ( ν, τ ] | 6 = 0) , ν ∈  − 1 2 , 1 2  , τ ∈ Z . (33) The total spread of the AF of X [ t ] over a time-frequency region S is gi ven by 0 ≤ ξ X ( S ) = R P S ζ Z ( ν, τ ] dν R P S dν ≤ 1 . (34) A process X [ t ] is strictly AF Compressible if ξ X ( S ) << 1 . W e define b ζ X ( ν, τ ] = I     b A (ht) X X ∗ ( ν, τ ]    6 = 0  , ν ∈  − 1 2 , 1 2  , τ ∈ [ − ( N − 1) , ( N − 1)] , (35) and thus b ξ X ( S ) = P P S b ζ X ( ν, τ ] P P S 1 (36) is an estimator of the total spread of the process. Matz and Hlawatsch define extended underspread processes, as those whose spread is essentially limited. Such processes will be approximated by the TEAF to the region of their essential support, i.e., where their magnitude is non-negligible in comparison to the rest of the ambiguity plane. This permits the quantification of the degree of variability of the time series at each lag, and the band of variability supported at a giv en lag τ 0 , can be determined by b ξ X ( { ( ν, τ ] : τ = τ 0 } ) . V . E X A M P L E S W e estimate the mean squared error (MSE) by comparing the estimated AFs to a theoretically based quantity that is the expected v alue of the EMAF , for a finite v alue of N , limited to the support ρ of the AF of the process of interest, µ A ( ν, τ ] I (( ν , τ ] ∈ ρ ) , and refer to this as the N -AF . If the process is the analytic process constructed from a real-valued process, then ρ is the support of the AF of the real-valued process. Note that A X X ∗ ( ν, τ ] is based on the infinite sequence { M X X ∗ [ t, τ ] } t,τ ∈ Z and does not exhibit any finite sample issues with spreading in local frequency and time lag due to finite sample effects. W e can nev er observe such v alues in a finite digital sample, because the maximum concentration of energy will behave like N (our number of sample points), and thus will be finite for finite N . Therefore, we ideally insist on the concentration of the AF to where A X X ∗ ( ν, τ ] is supported, b ut only letting the N -AF take finite sample length realizable values. A. Linear chirp immersed in white noise. W e first consider the case of a deterministic linear chirp, g [ t ] = exp[ j π (2 αt + β t 2 )] , immersed in analytic white noise. Chirps are commonly characterised in the ambiguity plane [7]. The chirp has a starting frequency α = 0 . 1 , with chirp rate β = 9 . 0196 · 10 − 4 , and the noise has variance σ 2 W = 0 . 6 . W e generate data sets of length N = 256 from this process. W e find the N -AF of the chirp b A ( N ) X X ∗ ( ν, τ ] = exp  j π  2 ατ − β τ 2 + ( β τ − ν )( N + | τ | − 1)  D N −| τ | ( β τ − ν ) (37) for β τ = ν and zero otherwise. W e have not inserted β τ = ν in this equation, since we are dealing with discrete values and β τ will not be identically equal to ν anywhere. Instead, we will for each τ i find ν i = min | β τ i − ν | . Thus, the N -AF will be zero except for at the points ( ν i , τ i ) , where it is defined by (37). Fig. 2 (a)–(c) shows the EMAF , TEAF and LBTEAF for one realisation of the process. The EMAF is substantially corrupted by noise, but both thresholding methods hav e removed most of the noise. It should be noted that the line has exhibited some increasing thickness especially tow ards high | τ | . This will explain why the MSE results are not quite as good as might be anticipated from this figure. T o quantify the performance of the estimators, we will estimate the MSE by Monte Carlo simulation, where we compare with the N -AF . W e generate K = 5 , 000 realisations of the process. The estimated MSE is shown in Fig. 2(d)–(f). Thresholding has reduced the MSE a great deal. T o observe differences between the two thresholding procedures, note that the TEAF suf fers from estimating a single noise variance, in that the regions where h ( ν , τ ] inflates the variance, noise remains in the estimation (see for example the band around the chirp in Fig. 2(e)). The LBTEAF looks like a considerable improvement, as it removes more noise, but does suffer from difficulties in estimating the noise around the rim of the domain. This leaves a “badly cleaned window” ef fect. T o quantify our visual impression of these procedures, we look at the total MSE of each estimator , which is obtained by summing the MSE over the ν - τ -plane, averaged ov er the fiv e thousand realisations. The resulting total MSEs and their DEP AR TMENT OF ST A TISTICAL SCIENCE RESEARCH REPOR T 293 9 ! " ! 0.4 ! 0.2 0 0.2 0.4 ! 200 ! 100 0 100 200 ! 20 ! 10 0 10 20 (a) ! " ! 0.4 ! 0.2 0 0.2 0.4 ! 200 ! 100 0 100 200 ! 20 ! 10 0 10 20 (b) ! " ! 0.4 ! 0.2 0 0.2 0.4 ! 200 ! 100 0 100 200 ! 20 ! 10 0 10 20 (c) ! " ! 0.4 ! 0.2 0 0.2 0.4 ! 200 ! 100 0 100 200 ! 40 ! 20 0 20 40 (d) ! " ! 0.4 ! 0.2 0 0.2 0.4 ! 200 ! 100 0 100 200 ! 40 ! 20 0 20 40 (e) ! " ! 0.4 ! 0.2 0 0.2 0.4 ! 200 ! 100 0 100 200 ! 40 ! 20 0 20 40 (f) Fig. 2. Chirp immersed in white noise. (a) EMAF (b) TEAF (c) LBTEAF for one realisation. Estimation of the MSE for (d) EMAF (e) TEAF (f) LBTEAF . All are on a dB scale. standard deviations are shown in T able I. W e see that the local bias corrected thresholding yields the lowest MSE, but both thresholding methods hav e quite significantly reduced the MSE of the EMAF . The reduction, is as mentioned, not quite as large as might have been anticipated, as the chirp has been broadened near its true support. But the reduction in MSE is respectable corresponding to a factor of almost four . W e estimate the total spread for the chirp immersed in analytic white noise by Monte Carlo simulation. The estimated total spread and the standard deviation of the total spread is gi ven in T able II. The EMAF will not be zero anywhere, so we do not need to estimate the total spread for this estimate, but equate this to 1 . Theoretically , the AF of the chirp should be nonzero for only 511 of 512 × 511 cells, which gi ves us a spread of 0 . 002 . Our estimated spread is larger than this number , but reflects the sparsity of the TEAF . B. Stationary pr ocess W e consider the analytic process X [ t ] corresponding to the real-valued stationary Moving A verage (MA) process R [ t ] = P L i =0 w i ξ [ t − i ] , where ξ ∼ N (0 , σ 2 ξ ) and E { ξ [ t ] ξ [ t − τ ] } = σ 2 ξ δ [ τ ] . The autocorrelation of this process is f M RR ∗ [ τ ] = E { R [ t ] R [ t − τ ] } = σ 2 ξ L −| τ | X i =0 w i w i + | τ | for | τ | ≤ L. (38) W e generate realisations of the process with w =  1 0 . 33 0 . 266 0 . 2 0 . 133 0 . 066  T , length N = 256 , L = 5 , and σ 2 ξ = 1 . Fig. 3(a)–(c) show the EMAF of one realisation of the process and the result of thresholding the EMAF . Both DEP AR TMENT OF ST A TISTICAL SCIENCE RESEARCH REPOR T 293 10 Chirp MA UM TVMA T otal MSE of EMAF 7 . 5382 · 10 7 2 . 0337 · 10 8 3 . 3225 · 10 7 5 . 0842 · 10 7 T otal MSE of TEAF 1 . 9622 · 10 7 9 . 7961 · 10 5 1 . 7409 · 10 5 3 . 1209 · 10 5 T otal MSE of L TEAF - 9 . 5598 · 10 5 1 . 7076 · 10 5 2 . 951 · 10 5 T otal MSE of LBTEAF 1 . 8398 · 10 7 - - - STD of total MSE of EMAF 9 . 9790 · 10 6 4 . 6983 · 10 7 7 . 3221 · 10 6 1 . 3809 · 10 7 STD of total MSE of TEAF 4 . 0799 · 10 6 4 . 9168 · 10 5 5 . 1692 · 10 4 1 . 7156 · 10 5 STD of total MSE of L TEAF - 4 . 1841 · 10 5 4 . 6430 · 10 4 1 . 1522 · 10 5 STD of total MSE of LBTEAF 4 . 0407 · 10 6 - - - T ABLE I A V E R AG E T OTA L M S E A N D T H E S TA N DA R D D E V I A T I O N ( S T D ) O F T OTA L M S E . Chirp MA UM TVMA T otal spread of TEAF 0 . 013 0 . 0013 6 . 5513 · 10 − 4 0 . 0012 T otal spread of L TEAF - 0 . 001 6 . 3735 · 10 − 4 0 . 001 T otal spread of LBTEAF 0 . 0156 - - - STD of total spread of TEAF 0 . 0028 0 . 0016 8 . 9984 · 10 − 4 0 . 0018 STD of total spread of L TEAF - 0 . 0012 8 . 1374 · 10 − 4 0 . 0014 STD of total spread of LBTEAF 0 . 0039 - - - T ABLE II A V E R AG E T OTA L S P R E A D A N D S T D O F T O TAL S P R E A D . thresholding procedures have zeroed out most of the points for ν 6 = 0 , and correspond to substantial improv ements to the EMAF . Fig. 4(a) shows the EMAF and the thresholded EMAF for one realisation with the N -AF for ν = 0 and τ ∈ [ − 10 , 10] . W e see that the thresholded EMAF and the N -AF are very similar , and only when the v ariance of the EMAF become comparable to the modulus square of the AF is the process thresholded. This exactly corresponds to the support of the real-valued process’ autocorrelation. The results of the MSE estimation for these estimators are shown in Fig. 3(d)–(f). Again, in most regions of the ambiguity plane we have a satisfying small value of the MSE, but at the ends, due to the lo w value of the threshold, suffer from the previously noted effects. Thresholding has reduced the MSE, and we see from T able I that the MSE is actually reduced by more than a factor of 100, which is a substantial reduction. Note also that the L TEAF has lower MSE than the TEAF . W e estimate the total spread of this process, see T able II. The theoretical spread is in this case equal to 11 cells out of 512 × 511 which is 4 . 2044 × 10 − 5 . Our estimated spread is larger , due to the analytic signal spreading energy in τ and the points at the ends that hav e not been successfully thresholded, but still reflects the geometry of the support. Using the discrete analytic signal and problems with edge effects makes the inflated spread inevitable. Whilst propositions 3.1 – 3.3 give the decay of the variance of the EMAF , at the rim of the ambiguity plane such expressions are not accurate as the O(1) error term will be of comparable magnitude. C. Uniformly modulated pr ocess Next, we consider the analytic process X [ t ] corresponding to the real-valued, non-stationary , Uniformly Modulated (UM) process, R [ t ] = σ X [ t ] ξ [ t ] , where ξ ∼ N (0 , 1) and E { ξ [ t ] ξ [ t − τ ] } = δ [ τ ] . The time-varying variance is defined as σ 2 X [ t ] = sin 2 (2 π f 0 t ) , which has an FT of Σ X X ∗ ( ν ) = 1 4 (2 δ ( ν ) − δ ( ν − 2 f 0 ) − δ ( ν + 2 f 0 )) . The N -AF is A ( N ) X X ∗ ( ν, τ ] = N for ν = 0 and τ = 0 , A ( N ) X X ∗ ( ν, τ ] = − N (1 / 2 − | 2 f 0 | ) for ν = ± 2 f 0 , τ = 0 , and zero otherwise. W e generate a realisation of the process of length N = 256 with f 0 = 0 . 09 . The EMAF and the results of the two thresholding procedures based on this realisation are shown in 5(a) – (c). The EMAF is very noisy , but both the TEAF and the L TEAF are substantially cleaner . In Fig. 4(b) we see the EMAF and TEAF for τ = 0 , and we observe that the thresholding has kept only the three points specified by the N -AF . Estimates of the MSE of these methods are shown in Fig. 5(d)–(f). Again we see that the thresholding has greatly reduced the MSE. From T able I we see that also for this process the MSE is reduced with a factor over one hundred from the EMAF to the TEAF and L TEAF , with the L TEAF doing a bit better than the TEAF . The theoretical spread of this process is 1 . 1466 × 10 − 5 , and our estimated spread is giv en in T able II. D. T ime-varying MA Finally , we will combine the MA with the uniformly modulated process, thus defining a real-valued T ime-V arying MA (TVMA) as R [ t ] = σ X [ t ] P L i =0 w i ξ [ t − i ] , where the w ’ s and ξ [ t ] is as defined in Section V -B. W e use σ X [ t ] = sin(2 π f 0 t ) DEP AR TMENT OF ST A TISTICAL SCIENCE RESEARCH REPOR T 293 11 ! " ! 0.4 ! 0.2 0 0.2 0.4 ! 200 ! 100 0 100 200 ! 10 ! 5 0 5 10 15 20 25 (a) ! " ! 0.4 ! 0.2 0 0.2 0.4 ! 200 ! 100 0 100 200 ! 10 ! 5 0 5 10 15 20 25 (b) ! " ! 0.4 ! 0.2 0 0.2 0.4 ! 200 ! 100 0 100 200 ! 10 ! 5 0 5 10 15 20 25 (c) ! " ! 0"4 ! 0"2 0 0"2 0"4 ! 200 ! 100 0 100 200 ! 20 ! 10 0 10 20 30 40 50 (d) ! " ! 0.4 ! 0.2 0 0.2 0.4 ! 200 ! 100 0 100 200 ! 20 ! 10 0 10 20 30 40 50 (e) ! " ! 0.4 ! 0.2 0 0.2 0.4 ! 200 ! 100 0 100 200 ! 20 ! 10 0 10 20 30 40 50 (f) Fig. 3. MA process. (a) EMAF (b) TEAF (c) L TEAF for one realisation. Estimation of the MSE for (d) EMAF (e) TEAF (f) L TEAF . All are on a dB scale. with f 0 = 0 . 042 , and we generate samples of length N = 256 . As always, we work with the analytic process corresponding to the real-valued process. W e find the N -AF of this process as b A ( N ) X X ∗ ( ν, τ ] =        ( N − | τ | ) R 1 / 2 0 h e S ( f − f 0 ) + e S ( f + f 0 ) i e j 2 π f τ d f for ν = 0 , | τ | ≤ L − ( N − | τ | ) R 1 / 2 − 2 f 0 0 e S ( f + f 0 ) e j 2 π f τ d f for ν = 2 f 0 , | τ | ≤ L − ( N − | τ | ) R 1 / 2 2 f 0 e S ( f − f 0 ) e j 2 π f τ d f for ν = − 2 f 0 , | τ | ≤ L (39) where e S ( f ) is the FT of the autocorrelation in (38). W e show the EMAF , TEAF and L TEAF of one realisation in Fig. 6(a)–(c), and we see that the thresholding has worked well, retaining only a few points. Fig. 4(c) shows the N -AF , the EMAF and the L TEAF for ν = 0 . Likewise, Fig. 4(d) shows the EMAF and L TEAF at τ = 0 . These two figures demonstrate that the geometry of the non-stationary process is retained by the thresholding estimator , even if the spread of the EMAF is cut short when the variance of the EMAF becomes too large. The estimated MSEs are shown in 6(d) – (f), and the total MSE in T able I. The MSE for the L TEAF shows a distinct improv ement, especially near the region (0 , 0) . The thresholding methods have again reduced the MSE with a factor over one hundred, and the total estimated spread demonstrates that the AF of this process is extremely sparse. V I . C O N C L U S I O N S In this paper we have introduced a new class of estimators for the AF of a non-stationary process that exhibits sparsity in the ambiguity plane. The AF is a fundamental characterisation of a non-stationary process, and many important properties of a DEP AR TMENT OF ST A TISTICAL SCIENCE RESEARCH REPOR T 293 12 ! 10 ! 5 0 5 10 0 100 200 300 400 500 600 700 ! ! 0.5 0 0.5 0 50 100 150 200 250 300 ! ! 10 ! 5 0 5 10 0 50 100 150 200 250 300 350 ! ! 0.4 ! 0.2 0 0.2 0.4 0 50 100 150 200 250 t x[t] Fig. 4. (a) TEAF (dashed), EMAF (dotted) and N -AF (solid) for ν = 0 for one realisation of the MA process. (b) TEAF (solid) and EMAF (dotted) for one realisation of the UM process for τ = 0 . (c) L TEAF (dashed), EMAF (dotted) and N -AF (solid) for ν = 0 for one realisation of the TVMA process. (d) L TEAF (solid) and EMAF (dotted) for one realisation of the TVMA process for τ = 0 . process can be deduced from its AF . The inherent resolution of the process in global time and frequency is an example of such. The AF has been used in estimating the generating mechanism of a non-stationary process [8]. It is also a popular tool in radar and sonar [5]. Despite this fact little efforts hav e been focused on the estimation of the AF , especially to produce estimators amenable to the determination of spread (size of the support). Characterisation of limited support is vital in designing the best analysis tools for generic second order non-stationary processes, a problem that remains open [1]. Based on the assumption of compression of the AF (small support), we proposed different threshold procedures for estimating the AF . The advantage of using the threshold estimator is that the size of the magnitude of the AF is compared with the magnitude of the AF at other time-frequency cells. Only if the local magnitude is suf ficiently large is the value not thresholded. This should enforce a strict support from for example a process whose AF is only extended underspread, see [2][p. 1072]. Unlike Matz and Hlawatsch [2] we do not calculate the spread of the process by fitting the support of the AF into a box centered at (0 , 0) , but simply count the number of non-zero AF coefficients. Conceptually this can be thought of as determining a finite number of cells that would represent most of the structure of the observ ed process. If a finite number of cells represent a full process then estimation of its generating mechanism is possible. Hence whilst it clearly is not necessary to constrain the process to be underspread, i.e. be concentrated in support around the origin in the AF domain [4], it is necessary to constrain the possible degree of dependency in the process to ensure that the generating mechanism of the process could be estimated. Piv otal to thresholding procedures is determining a suitable threshold. W e implemented the thresholding with a global variance estimate as well as a local variance estimate, both adjusted for each point in the ambiguity plane. For deterministic signals in analytic white noise we also proposed a procedure which remov es bias in the estimator caused by the noise. W e demonstrated the superior properties of the threshold estimator in finite sample sizes for a variety of common types of non- stationary and stationary processes. The reductions in MSE compared to a pre viously used, if naiv e, estimator for the AF were considerable, up to factors of over a hundred. Formally our proposed threshold procedure for zero-mean processes is only valid if the process is underspread in which case we determine the asymptotic distribution of the EMAF . W e conjecture that asymptotic normality can be shown for a larger class of processes, and that in the case of any process with compressed AF , thresholding is an appropriate estimation procedure. The estimated spread does not fully reflect the degree of sparsity of the AF: an inflation is accrued due to the usage of the discrete analytic signal and edge effects. Resolving such effects fully remains an outstanding issue. The definition of spread in this paper has the clear advantage of interpretation, corresponding to the fraction of points where the mean of the EMAF dominates the variance of the EMAF . The AF remains the least studied of time-frequency representations of non-stationary signals, perhaps because its arguments lack direct global interpretability . A large class of processes that can be estimated exhibit compression in this domain. W e anticipate that further study of the AF of non-stationary processes will yield understanding into what generating mechanisms can be determined from a giv en sample with fixed sampling rates and sample size. DEP AR TMENT OF ST A TISTICAL SCIENCE RESEARCH REPOR T 293 13 ! " ! 0.4 ! 0.2 0 0.2 0.4 ! 200 ! 100 0 100 200 ! 10 ! 5 0 5 10 15 20 (a) ! " ! 0.4 ! 0.2 0 0.2 0.4 ! 200 ! 100 0 100 200 ! 10 ! 5 0 5 10 15 20 (b) ! " ! 0.4 ! 0.2 0 0.2 0.4 ! 200 ! 100 0 100 200 ! 10 ! 5 0 5 10 15 20 (c) ! " ! 0.4 ! 0.2 0 0.2 0.4 ! 200 ! 100 0 100 200 ! 30 ! 20 ! 10 0 10 20 30 40 (d) ! " ! 0.4 ! 0.2 0 0.2 0.4 ! 200 ! 100 0 100 200 ! 30 ! 20 ! 10 0 10 20 30 40 (e) ! " ! 0.4 ! 0.2 0 0.2 0.4 ! 200 ! 100 0 100 200 ! 30 ! 20 ! 10 0 10 20 30 40 (f) Fig. 5. UM process. (a) EMAF (b) TEAF (c) L TEAF for one realisation. Estimation of the MSE for (d) EMAF (e) TEAF (f) L TEAF . All are on a dB scale. A P P E N D I X A P RO O F O F P RO P O S I T I O N 3 . 1 The Empirical Cross Ambiguity Function for two finite equal length samples of signals { g 1 [ t ] } t and { g 2 [ t ] } t of length N is giv en by b A g 1 g ∗ 2 ( ν, τ ] = P N − 1+min( τ , 0) t =max(0 ,τ ) g 1 [ t ] g ∗ 2 [ t − τ ] e − j 2 π ν t . Then b A X X ∗ ( ν, τ ] = N − 1+min(0 ,τ ) X t =max(0 ,τ ) c M X X ∗ [ t, τ ] e − j 2 π ν t = N − 1+min(0 ,τ ) X t =max(0 ,τ ) ( g [ t ] + W [ t ]) ( g ∗ [ t − τ ] + W ∗ [ t − τ ]) + e − j 2 π ν t = b A g g ∗ ( ν, τ ] + b A W g ∗ ( ν, τ ] + b A g W ∗ ( ν, τ ] + b A W W ∗ ( ν, τ ] . (A-1) DEP AR TMENT OF ST A TISTICAL SCIENCE RESEARCH REPOR T 293 14 ! " ! 0.4 ! 0.2 0 0.2 0.4 ! 200 ! 100 0 100 200 ! 20 ! 1& ! 10 ! & 0 & 10 1& 20 (a) ! " ! 0.4 ! 0.2 0 0.2 0.4 ! 200 ! 100 0 100 200 ! 20 ! 15 ! 10 ! 5 0 5 10 15 20 (b) ! " ! 0"4 ! 0"2 0 0"2 0"4 ! 200 ! 100 0 100 200 ! 20 ! 15 ! 10 ! 5 0 5 10 15 20 (c) ! " ! 0.4 ! 0.2 0 0.2 0.4 ! 200 ! 100 0 100 200 ! 40 ! 20 0 20 40 (d) ! " ! 0.4 ! 0.2 0 0.2 0.4 ! 200 ! 100 0 100 200 ! 40 ! 20 0 20 40 (e) ! " ! 0.4 ! 0.2 0 0.2 0.4 ! 200 ! 100 0 100 200 ! 40 ! 20 0 20 40 (f) Fig. 6. TVMA process. (a) EMAF (b) TEAF (c) L TEAF for one realisation. Estimation of the MSE for (d) EMAF (e) TEAF (f) L TEAF . All are on a dB scale. As g [ t ] is deterministic, and W [ t ] is zero-mean, it follo ws E n b A W g ∗ ( ν, τ ] o = 0 , and E n b A g W ∗ ( ν, τ ] o = 0 . Using (4) for τ ≥ 0 and the fact that S X X ∗ ( ν, f ) = σ 2 W δ ( ν ) for f ≥ 0 and zero otherwise for the noise process, we find for τ > 0 E n b A W W ∗ ( ν, τ ] o = Z 1 / 2 0 Z 1 / 2 0 e j π ( ν 0 )( N − 1 − τ ) E n d e X ( ν 0 + f ) d e X ∗ ( f ) o e j 2 π ( ν 0 + f ) τ D N − τ ( ν 0 − ν ) e − j π ν ( N + τ − 1) = Z 1 / 2 0 Z 1 / 2 0 e j π ν 0 ( N − 1 − τ ) σ 2 W δ ( ν 0 ) e j 2 π ( ν 0 + f ) τ D N − τ ( ν 0 − ν ) e − j π ν ( N + τ − 1) d f 1 d f 2 = e − j π ν ( N + τ − 1) D N − τ ( ν ) Z 1 / 2 0 σ 2 W e j 2 π f τ d f = σ 2 W 2 e − j π ν ( N + τ − 1) D N − τ ( ν ) e j π τ / 2 sinc( π τ / 2) , where (as usual) sinc( x ) = sin( π x ) /π x . Thus taking expectations of (A-1) and using the linearity of E {·} , the result follows. Mutatis mutandis the calculations are implemented for τ < 0 . W e start from (A-1) and note that v ar n b A X X ∗ ( ν, τ ] o can be written in terms of covariances of A g g ∗ ( ν, τ ] , A g W ∗ ( ν, τ ] etc. DEP AR TMENT OF ST A TISTICAL SCIENCE RESEARCH REPOR T 293 15 Here, g [ t ] is deterministic and W [ t ] is Gaussian proper and so v ar n b A g g ∗ ( ν, τ ] o = 0 , cov n b A g g ∗ ( ν, τ ] , A W g ∗ ( ν, τ ] o = 0 co v n b A g g ∗ ( ν, τ ] , b A g W ∗ ( ν, τ ] o = 0 , cov n b A g g ( ν, τ ] , b A W W ∗ ( ν, τ ] o = 0 co v n b A W g ∗ ( ν, τ ] , A g W ∗ ( ν, τ ] o = 0 , cov n b A W g ∗ ( ν, τ ] , A W W ∗ ( ν, τ ] o = 0 co v n b A g W ∗ ( ν, τ ] , A W W ∗ ( ν, τ ] o = 0 . (A-2) W e next define G ( f , τ ] = N − 1 −| τ | P t =0 g [ t − | τ | I ( τ < 0)] e − j 2 π f t , where I ( · ) is the indicator function, and note that G ( f , τ ] is supported on negati ve frequencies as well ev en if g [ t ] is analytic. If N − | τ | → ∞ , G ( f , τ ] will tend toward only being supported on positiv e frequencies. Then, for τ > 0 , v ar n b A W g ∗ ( ν, τ ] o = v ar ( N − 1 − τ X t =0 W [ t + τ ] g ∗ [ t ] e − j 2 π ν ( t + τ ) ) = v ar ( Z 1 / 2 0 d e X ( f 1 ) N − 1 − τ X t =0 g [ t ] e − j 2 π ( f 1 − ν ) t ! ∗ e j 2 π τ f 1 d f 2 ) = v ar ( Z 1 / 2 0 d e X ( f 1 ) G ∗ ( f 1 − ν , τ ] e j 2 π τ f 1 d f 2 ) = Z 1 / 2 0 Z 1 / 2 0 E n d e X ( f 1 ) d e X ∗ ( f 1 ) o G ∗ ( f 1 − ν , τ ] G ( f 2 − ν , τ ] e j 2 π τ ( f 1 − f 2 ) = σ 2 W Z 1 / 2 0 | G ( f − ν, τ ] | 2 d f (A-3) After a change of variables f 0 = f − ν , we are giv en an outer integral of R 1 / 2 − ν − ν ov er f 0 , rather than R 1 / 2 0 . For ν > 0 we can rewrite this as R 1 / 2 − ν 0 plus a contribution that will hav e the inner integrals integrating to negligible contributions. For ν < 0 we can rewrite this as R 1 / 2 − ν plus a contribution that will have the inner integrals integrating to negligible contributions. W e denote h ( ν, τ ] = R 1 / 2 − max(0 ,ν ) max( − ν, 0) | G ( f , τ ] | 2 d f . For a generic harmonizable process X [ t ] , τ > 0 giv es v ar n b A X X ∗ ( ν, τ ] o = Z 1 / 2 0 Z 1 / 2 0 Z 1 / 2 0 Z 1 / 2 0 e j 2 π ( f 1 − α 1 ) τ e j π ( f 1 − f 2 − α 1 + α 2 )( N − 1 − τ ) D N − τ ( f 1 − f 2 − ν ) D N − τ ( α 1 − α 2 − ν )E n d e X ( f 1 ) d e X ∗ ( f 2 ) d e X ∗ ( α 1 ) d e X ( α 2 ) o −    E n b A X X ∗ ( ν, τ ] o    2 . This follows by direct calculation starting from (4).Using Isserlis’ theorem [27] we write the variance as v ar n b A X X ∗ ( ν, τ ] o = Z 1 / 2 0 Z 1 / 2 0 Z 1 / 2 0 Z 1 / 2 0 S X X ∗ ( f 1 − α 1 , α 1 ) S X X ∗ ( α 2 − f 2 , f 2 ) e j π ( f 1 − f 2 − α 1 + α 2 )( N − 1 − τ ) e j 2 π ( f 1 − α 1 ) τ D N − τ ( f 1 − f 2 − ν ) D N − τ ( α 1 − α 2 − ν ) d f 1 d f 2 dα 1 dα 2 , (A-4) which for analytic white noise reduces to v ar n b A X X ∗ ( ν, τ ] o = σ 4 W Z 1 / 2 0 Z 1 / 2 0 D 2 N − τ ( f 1 − f 2 − ν ) d f 1 d f 2 = σ 4 W ( N − τ ) Z 1 / 2+min(0 , − ν ) max(0 , − ν ) Z g 2 ( N − τ ) − (1 / 2 − g 2 )( N − τ ) sinc 2 ( ξ ) dξ dg 2 = σ 4 W ( N − τ )(1 / 2 − | ν | ) + O (1) . Implementing the same calculations for τ < 0 deriv es analogous expressions with τ replaced by | τ | . Putting (A-2), (A-3) and v ar { A g W ∗ ( ν, τ ] } = v ar { A W g ∗ ( − ν, − τ ] } together , and mutatis mutandis implementing the calculations for τ < 0 , thus yields v ar n b A X X ∗ ( ν, τ ] o = 0 + σ 2 W h ( ν, τ ] + σ 2 W h ( − ν, − τ ] + σ 4 W ( N − | τ | )(1 / 2 − | ν | ) + O (1) . DEP AR TMENT OF ST A TISTICAL SCIENCE RESEARCH REPOR T 293 16 W e start from (A-1) and note that rel n b A X X ∗ ( ν, τ ] o can be written as an aggregation of relations. As g [ t ] is deterministic and W [ t ] is zero-mean Gaussian, rel n b A g g ∗ ( ν, τ ] o = 0 , rel n b A g g ∗ ( ν, τ ] , A W g ∗ ( ν, τ ] o = 0 (A-5) rel n b A g g ∗ ( ν, τ ] , b A g W ∗ ( ν, τ ] o = 0 , rel n b A g g ∗ ( ν, τ ] , b A W W ∗ ( ν, τ ] o = 0 rel n b A W g ∗ ( ν, τ ] , b A W W ∗ ( ν, τ ] o = 0 , rel n b A g W ∗ ( ν, τ ] , b A W W ∗ ( ν, τ ] o = 0 . For τ ≥ 0 rel n b A W g ∗ ( ν, τ ] , b A g W ∗ ( ν, τ ] o = E n b A W g ∗ ( ν, τ ] b A g W ∗ ( ν, τ ] o = E ( Z 1 / 2 0 Z 1 / 2 0 d e X ( f 1 ) d e X ∗ ( f 2 ) e j 2 π ( f 1 − ν ) τ G ∗ ( f 1 − ν , τ ] G ( f 2 + ν , − τ ] e j 2 π ν τ ) = σ 2 W Z 1 / 2 0 G ∗ ( f − ν, τ ] G ( f + ν, − τ ] e j 2 π ( f − 2 ν ) τ d f . (A-6) The factor e − j 4 π ν τ is not present for τ < 0 . W e define h 0 ( ν, τ ] = R 1 / 2+min(0 , − ν ) max(0 , − ν ) G ∗ ( f , τ ] G ( f + 2 ν, − τ ] e j 2 π ( f − sign ( τ ) ν ) τ d f , where sign ( τ ) = 1 for τ ≥ 0 and sign ( τ ) = − 1 for τ < 0 . If G ( f , τ ] is mainly supported only ov er a range of frequencies we would assume with increasing ν that h 0 ( ν, τ ] → 0 . Furthermore we note that rel n b A W g ∗ ( ν, τ ] o = 0 , and rel n b A g W ∗ ( ν, τ ] o = 0 . W e can note that for a generic harmonizable process X [ t ] , rel n b A X X ∗ ( ν, τ ] o = Z 1 / 2 0 Z 1 / 2 0 Z 1 / 2 0 Z 1 / 2 0 e j π ( f 1 − f 2 + α 1 − α 2 )( N − 1 − τ ) e j 2 π ( f 1 + α 1 ) τ e − 2 j π ν ( N + τ − 1) D N − τ ( f 1 − f 2 − ν ) D N − τ ( α 1 − α 2 − ν )E n d e X ( f 1 ) d e X ∗ ( f 2 ) d e X ( α 1 ) d e X ∗ ( α 2 ) o − E n b A X X ∗ ( ν, τ ] o 2 . This follows by direct calculation starting from (4). Using Isserlis’ theorem and duplicating the calculations for the variance, we write the relation as rel n b A X X ∗ ( ν, τ ] o = Z 1 / 2 0 Z 1 / 2 0 Z 1 / 2 0 Z 1 / 2 0 S X X ∗ ( f 1 − α 2 , α 2 ) S X X ∗ ( α 1 − f 2 , f 2 ) e j π ( f 1 − f 2 + α 1 − α 2 )( N − 1 − τ ) e j 2 π ( f 1 + α 1 ) τ D N − τ ( f 1 − f 2 − ν ) D N − τ ( α 1 − α 2 − ν ) e − 2 j π ν ( N + τ − 1) d f 1 d f 2 dα 1 dα 2 . (A-7) For an analytic white-noise process this corresponds to rel n b A X X ∗ ( ν, τ ] o = Z 1 / 2 0 Z 1 / 2 0 Z 1 / 2 0 Z 1 / 2 0 σ 2 W δ ( f 1 − α 2 ) σ 2 W δ ( α 1 − f 2 ) e j π ( f 1 − f 2 + α 1 − α 2 )( N − 1 − τ ) e j 2 π ( f 1 + α 1 ) τ D N − τ ( f 1 − f 2 − ν ) D N − τ ( α 1 − α 2 − ν ) e − 2 j π ν ( N + τ − 1) d f 1 d f 2 dα 1 dα 2 = Z 1 / 2 0 Z 1 / 2 0 σ 4 W e j π ( f 1 − f 2 + f 2 − f 1 )( N − 1 − τ ) e j 2 π ( f 1 + f 2 ) τ D N − τ ( f 1 − f 2 − ν ) D N − τ ( f 2 − f 1 − ν ) e − 2 j π ν ( N + τ − 1) d f 1 d f 2 . If ν = O (1 / ( N − τ )) then this corresponds to (up to order ( N − τ ) / N ) rel n b A X X ∗ ( ν, τ ] o = Z 1 / 2 0 Z f ( N − τ ) ( f − 1 / 2)( N − τ ) σ 4 W e j 4 πf τ D N − τ f 0 ( N − τ ) − ν ! D N − τ − f 0 ( N − τ ) − ν ! e − 2 j πν ( N + τ − 1) N − τ d f d f 0 = Z 1 / 2 0 Z ( f − ν )( N − τ ) ( f − 1 / 2 − ν )( N − τ ) σ 4 W e j 4 πf τ D N − τ „ f 0 ( N − τ ) « D N − τ „ − f 0 ( N − τ ) − 2 ν « e − 2 j πν ( N + τ − 1) d f d f 0 N − τ . If ν is not O (1 / ( N − τ )) , then rel n b A X X ∗ ( ν, τ ] o is negligible. If ν = O (1 / ( N − τ )) then with ξ ( ν ) = ( N − τ ) ν = O (1) we obtain that rel n b A X X ∗ ( ν, τ ] o = ( N − τ ) σ 4 W Z ∞ −∞ sinc( f 0 )sinc( f 0 + 2 ξ ( ν )) e − 2 j π ν ( N + τ − 1) d f 0 Z 1 / 2+min(0 ,ν ) max(0 ,ν ) e j 4 π f τ d f ! = ( N − τ ) σ 4 W L ( N − τ , ν ) e − 2 j π ν ( N + τ − 1) e j 4 π (1 / 2+min(0 ,ν )) τ − e j 4 π max(0 ,ν ) τ j 4 π τ + O (1) = − ( N − τ ) σ 4 W L ( N − τ , ν ) e − 2 j π ν ( N − 1)  sin(2 π | ν | τ ) 2 π τ I ( τ 6 = 0) − 1 2 I ( τ = 0)  + O (1) . DEP AR TMENT OF ST A TISTICAL SCIENCE RESEARCH REPOR T 293 17 This defines the quantity L ( N − τ , ν ) , which decays like O  1 ξ ( ν ) 2  . Thus if τ 6 = 0 rel n b A X X ∗ ( ν, τ ] o = − ( N − τ ) | ν | σ 4 W L ( N − τ , ν ) e − 2 j π ν ( N − 1) sinc(2 | ν | τ ) + 2 σ 2 W h 0 ( ν, τ ] + O (1) . Mutatis mutandis calculations can be replicated for τ < 0 . A P P E N D I X B P RO O F O F P RO P O S I T I O N 3 . 2 For τ ≥ 0 we note that E n b A X X ∗ ( ν, τ ] o = Z 1 / 2 0 Z 1 / 2 0 e j π ( f 1 − f 2 )( N − 1 − τ ) E n d e X ( f 1 ) d e X ∗ ( f 2 ) o e j 2 π f 1 τ D N − τ ( f 1 − f 2 − ν ) e − j π ν ( N + τ − 1) = D N − τ ( ν ) e − j π ν ( N + τ − 1) f M X X ∗ [ τ ] . Note that D N − τ ( ν ) is an even function, and the expression is valid for ν < 0 . Calculations can mutatis mutandis be replicated for negati ve τ . W e note that with τ ≥ 0 v ar n b A X X ∗ ( ν, τ ] o = Z 1 / 2 0 Z 1 / 2 0 e S X X ∗ ( f 1 ) e S X X ∗ ( f 2 ) D 2 N − τ ( f 1 − f 2 − ν ) d f 1 d f 2 = ( N − τ )(1 / 2 − | ν | ) A X X ∗ ( − ν, 0] + O (1) , with A X X ∗ ( ν, τ ] = R 1 / 2+min( ν, 0) max(0 , − ν ) e S X X ∗ ( f − ν ) e S X X ∗ ( f ) e j 4 π f τ d f / (1 / 2 − | ν | ) . Mutatis mutandis the calculations can be replicated for τ < 0 . Finally we start from (A-7) and note that for stationary processes rel n b A X X ∗ ( ν, τ ] o = Z 1 / 2 0 Z 1 / 2 0 e S X X ∗ ( f 1 ) e S X X ∗ ( f 2 ) e j 2 π ( f 1 + f 2 ) τ D N − τ ( f 1 − f 2 − ν ) e − 2 j π ν ( N + τ − 1) D N − τ ( f 1 − f 2 − ν ) d f 1 d f 2 = e − 2 j π ν ( N + τ − 1) Z 1 / 2 0 Z ( N − τ )( g 2 − ν ) ( N − τ )( g 2 − 1 / 2 − ν ) e S X X ∗ ( g 2 ) e S X X ∗  g 2 − g 1 N − τ − ν  e j 2 π (2 g 2 − g 1 N − τ − ν ) τ D N − τ  g 1 N − τ  D N − τ  g 1 + 2( N − | τ | ) ν ) N − τ  dg 1 dg 2 / ( N − τ ) = e − 2 j π ν ( N + τ − 1) ( N − τ )(1 / 2 − | ν | ) L ( N − τ , ν ) A X X ∗ ( ν, τ ] + O (1) . mutatis mutandis results may be deriv ed for τ < 0 . A P P E N D I X C P RO O F O F P RO P O S I T I O N 3 . 3 Using (4) and by direct calculation we note that E n b A X X ∗ ( ν, τ ] o = Z 1 / 2 0 Z 1 / 2 0 e j π ( f 1 − f 2 )( N − 1 − τ ) E n d e X ( f 1 ) d e X ∗ ( f 2 ) o e j 2 π f 1 τ D N − τ ( f 1 − f 2 − ν ) e − j π ν ( N + τ − 1) = Z 1 / 2 0 Z f f − 1 / 2 e j π ν 0 ( N − 1 − τ ) X τ 0 A X X ∗ ( ν 0 , τ 0 ] e − j 2 π ( f − ν 0 ) τ 0 e j 2 π f τ D N − τ ( ν 0 − ν ) e − j π ν ( N + τ − 1) dν 0 d f = Z 1 / 2 − max( ν, 0) max(0 , − ν ) e j 2 π f τ Z f ( N − τ ) ( f − 1 / 2)( N − τ ) e j π ν 0 ( N − 1 − τ ) / ( N − τ ) Σ X X ∗  ν 0 N − τ + ν  D N − τ  ν 0 N − τ  dν 0 / ( N − τ ) d f = Z 1 / 2 − max( ν, 0) max(0 , − ν ) e j 2 π f τ Σ X X ∗ ( ν ) Z ∞ −∞ e j π ν 0 ( N − 1 − τ ) / ( N − τ ) sinc( ν 0 ) dν 0 d f + O (1) = Σ X X ∗ ( ν ) e j 2 π (1 / 2 − max( ν, 0)) τ − e j 2 π max( − ν, 0)) τ j 2 π τ + O (1) = (1 / 2 − | ν | )Σ X X ∗ ( ν ) e j π (1 / 2 − ν ) τ sinc((1 / 2 − | ν | ) τ ) + O (1) . DEP AR TMENT OF ST A TISTICAL SCIENCE RESEARCH REPOR T 293 18 Mutatis mutandis we may deriv e the results for τ < 0 . W e start from (A-4), and again by direct calculation v ar n b A X X ∗ ( ν, τ ] o = Z 1 / 2 0 Z 1 / 2 0 Z 1 / 2 0 Z 1 / 2 0 S X X ∗ ( f 1 − α 1 , α 1 ) S X X ∗ ( f 2 − α 2 , α 2 ) e j 2 π ( f 1 − α 1 ) τ e j π ( f 1 − f 2 − α 1 + α 2 )( N − 1 − τ ) D N − τ ( f 1 − f 2 − ν ) D N − τ ( α 1 − α 2 − ν ) d f 1 d f 2 dα 1 dα 2 = Z 1 / 2 0 Z 1 / 2 0 Z f f − 1 / 2 Z α α − 1 / 2 Σ X X ∗ ( f − α )Σ ∗ X X ∗ ( f − ν 0 − α + α 0 ) e j π ( ν 0 − α 0 )( N − 1 − τ ) e j 2 π ( f − α ) τ D N − τ ( ν 0 − ν ) D N − τ ( α 0 − ν ) dν 0 dα 0 d f dα = Z 1 / 2 − ν − ν Z 1 / 2 − ν − ν Z f ( N − τ ) ( f − 1 / 2)( N − τ ) Z α ( N − τ ) ( α − 1 / 2)( N − τ ) Σ X X ∗ ( f − α )Σ ∗ X X ∗  f − ν 0 − α 0 N − τ − α  e j π ( ν 0 − α 0 )( N − 1 − τ ) / ( N − τ ) e j 2 π ( f − α ) τ D N − τ ( ν 0 / ( N − τ )) D N − τ ( α 0 / ( N − τ )) dν 0 dα 0 d f dα/ ( N − τ ) 2 = Z 1 / 2 − max( ν, 0) max(0 , − ν ) Z 1 / 2 − max( ν, 0) max(0 , − ν ) Σ X X ∗ ( f − α )Σ ∗ X X ∗ ( f − α ) e j 2 π ( f − α ) τ dαd f + O (1) = Z 1 / 2 −| ν | − 1 / 2+ | ν |  1 2 − | ν |  | Σ X X ∗ ( f ) | 2 e j 2 π y τ d f + O (1) = (1 / 2 − | ν | )( N − τ ) Z 1 / 2 −| ν | − 1 / 2+ | ν | | Σ X X ∗ ( f ) | 2 ( N − τ ) e j 2 π f τ d f + O (1) . W e start from (A-7) rel n b A X X ∗ ( ν, τ ] o = Z 1 / 2 0 Z 1 / 2 0 Z 1 / 2 0 Z 1 / 2 0 S X X ∗ ( f 1 − α 2 , α 2 ) S X X ∗ ( α 1 − f 2 , f 2 ) e j π ( f 1 − f 2 + α 1 − α 2 )( N − 1 − τ ) × e j 2 π ( f 1 + α 1 ) τ D N − τ ( f 1 − f 2 − ν ) D N − τ ( α 1 − α 2 − ν ) e − 2 j π ν ( N + τ − 1) d f 1 d f 2 dα 1 dα 2 = Z 1 / 2 0 Z 1 / 2 0 Z 1 / 2 0 Z 1 / 2 0 Σ X X ∗ ( f 1 − α 2 )Σ ∗ X X ∗ ( f 2 − α 1 ) e j π ( f 1 − f 2 + α 1 − α 2 )( N − 1 − τ ) e j 2 π ( f 1 + α 1 ) τ D N − τ ( f 1 − f 2 − ν ) D N − τ ( α 1 − α 2 − ν ) e − 2 j π ν ( N + τ − 1) d f 1 d f 2 dα 1 dα 2 W e implement the change of variables ν 00 = ( ν 0 − ν )( N − τ ) and α 00 = ( α 0 − ν )( N − τ ) . rel n b A X X ∗ ( ν, τ ] o = Z 1 / 2 0 Z 1 / 2 0 Z ( f − ν )( N − τ ) ( f − 1 / 2 − ν )( N − τ ) Z ( α − ν )( N − τ ) ( α − 1 / 2 − ν )( N − τ ) Σ X X ∗  f − α + α 00 N − τ + ν  × Σ ∗ X X ∗  f − ν 00 N − τ − ν − α  e j π ( ν 00 + α 00 )( N − 1 − τ ) / ( N − τ )+ j π ( N − 1 − τ )2 ν e j 2 π ( f + α ) τ D N − τ  ν 00 N − τ  D N − τ  α 00 N − τ  e − 2 j π ν ( N + τ − 1) 1 ( N − τ ) 2 dα 00 dν 00 d f dα = e − 4 j π ν τ Z 1 / 2+min(0 ,ν ) max(0 ,ν ) Z 1 / 2+min(0 ,ν ) max(0 ,ν ) Σ X X ∗ ( f − α + ν )Σ ∗ X X ∗ ( f − ν − α ) e j 2 π ( f + α ) τ d f dα + O (1) . Mutatis mutandis the expressions may be deriv ed for τ < 0 . A P P E N D I X D P RO O F O F T H E O R E M 3 . 1 W e outline the proof for τ ≥ 0 , the same arguments hold mutatis mutandis for τ < 0 . W e let { X [ t ] } be the analytic extension of a real-v alued process { Y [ t ] } or equiv alently X [ t ] = Y [ t ] + j U [ t ] . As { Y [ t ] } is a real-valued underspread process there exists a constant T = O (1) such that ∀ | τ | ≥ T , cov { Y [ t ] , Y [ t + τ ] } = 0 . This implies that ∀ | τ | ≥ T co v { Y [ t ] , U [ t + τ ] } = O (1 / ( τ − T )) , where { U [ t ] } is the discrete Hilbert transform of { Y [ t ] } . For conv enience assume N − τ = T n . This is of course not a necessity but simplifies the proofs. Otherwise for T m such that T m < N − τ and T ( m + 1) > N − τ we can split the sum in the EAF into two parts, and show the latter part has negligible contributions to the total sum. The first part, summing from zero ov er t up to T m , can be treated as if N − τ = T m with m = n . W e now take a full length observation Y = [ Y 0 , . . . , Y N − 1 ] and construct n subv ectors by the following mechanism Y 1 = [ Y [0] , Y [ T ] , . . . , Y [( n − 1) T ]] , Y 2 = [ Y [1] , Y [ T + 1] , . . . , Y [( n − 1) T + 1]] , . . . . DEP AR TMENT OF ST A TISTICAL SCIENCE RESEARCH REPOR T 293 19 The vectors are constructed so that the elements of Y k are uncorrelated and if we define U k and X k , as the obvious extensions to Y k , then cov ( Y k [ t ] , U l [ s ]) = O (1 / | k − l + T ( t − s ) | ) , cov ( X k [ t ] , X l [ s ]) = O (1 / | k − l + T ( t − s ) | ) if k 6 = l . Before constructing the CL T type of argument let us study how we shall represent b A X X ∗ ( ν, τ ] . W e write b A X X ∗ ( ν, τ ] = T − 1 X u =0 n − 1 X v =0 x [ v T + u + τ ] x ∗ [ v T + u ] e − j 2 π ν ( vT + u ) e − j 2 π ν τ . The expectation of b A X X ∗ ( ν, τ ] follows directly from this expression, and we additionally note that E n c M X X ∗ [ t, τ ] o = M X X ∗ [ t, τ ] . Using Isserlis’ formula [27] and the assumption that the process is proper , we find the variance of b A X X ∗ ( ν, τ ] as v ar n b A X X ∗ ( ν, τ ] o = n − 1 X v 1 =0 n − 1 X v 2 =0 T − 1 X u 1 =0 T − 1 X u 2 =0 c N ( ν, τ ; u 1 , u 2 , v 1 , v 2 ] , where c N ( ν, τ ; u 1 , u 2 , v 1 , v 2 ] = e − j 2 π ν ( T ( v 1 − v 2 )+ u 1 − u 2 ) M X X ∗ [ v 1 T + u 1 + τ , T ( v 1 − v 2 ) + u 1 − u 2 ] M ∗ X X ∗ [ v 1 T + u 1 , T ( v 1 − v 2 ) + u 1 − u 2 ] . As { Y [ t ] } is exactly underspread if v 1 6 = v 2 and v 1 6 = v 2 ± 1 we hav e that with t l = v l T + u l , c N ( ν, τ ; u 1 , u 2 , v 1 , v 2 ] = O  1 / ( t 1 − t 2 ) 2  . Furthermore, P | v 1 − v 2 | > 1 c N ( ν, τ ; u 1 , u 2 , v 1 , v 2 ] = O (log ( n )) . Otherwise c N ( ν, τ ; u 1 , u 2 , v 1 , v 2 ] = M X X ∗ [ v 1 T + u 1 + τ , u 1 − u 2 + ( v 1 − v 2 ) T ] M ∗ X X ∗ [ v 1 T + u 1 , u 1 − u 2 + ( v 1 − v 2 ) T ] . Thus we find that (up to O (log ( n )) ) v ar n b A X X ∗ ( ν, τ ] o = T − 1 X u 1 =0 u 1 +( T − 1) X u 2 = u 1 − ( T − 1) n − 1 X v =0 e − j 2 π ν [ u 1 − u 2 ] M X X ∗ [ v T + u 1 + τ , u 1 − u 2 ] M ∗ X X ∗ [ v T + u 1 , u 1 − u 2 ] . Let x = v T + u 1 and τ 0 = u 1 − u 2 . Then this expression is re written for τ ≥ 0 as v ar n b A X X ∗ ( ν, τ ] o = N − τ − 1 X x =0 T − 1 X τ 0 = − ( T − 1) e − j 2 π ν τ 0 M X X ∗ [ x + τ , τ 0 ] M ∗ X X ∗ [ x, τ 0 ] + O (log ( n )) . Thus v ar n b A X X ∗ ( ν, τ ] o = O ( N ) . W e note that rel n b A X X ∗ ( ν, τ ] o = n − 1 X v 1 =0 n − 1 X v 2 =0 T − 1 X u 1 =0 T − 1 X u 2 =0 r N ( ν, τ ; u 1 , u 2 , v 1 , v 2 ]( ν, τ ] , r N ( ν, τ ; u 1 , u 2 , v 1 , v 2 ] = e − j 2 π ν ( T ( v 1 + v 2 )+ u 1 + u 2 +2 τ ) M X X ∗ [ v 1 T + u 1 + τ , T ( v 1 − v 2 ) + u 1 − u 2 + τ ] × M ∗ X X ∗ [ v 1 T + u 1 , T ( v 1 − v 2 ) + u 1 − u 2 − τ ] If | t 1 − t 2 + τ | > T and | t 1 − t 2 − τ | > T r N ( ν, τ ; u 1 , u 2 , v 1 , v 2 ] = O  1 ( t 1 − t 2 + τ )( t 1 − t 2 − τ )  , whilst if | t 1 − t 2 + τ | > T and | t 1 − t 2 − τ | < T r N ( ν, τ ; u 1 , u 2 , v 1 , v 2 ] = O  1 t 1 − t 2 + τ  . Thus if | τ | > T , P v 1 ,v 2 r N ( ν, τ ; u 1 , u 2 , v 1 , v 2 ] = O (log ( n )) . If | τ | < T P | v 1 − v 2 | > 1 r N ( ν, τ ; u 1 , u 2 , v 1 , v 2 ] = O (log ( n )) . If τ < T and | v 1 − v 2 | ≤ 1 we have that r N ( ν, τ ; u 1 , u 2 , v 1 , v 2 ] = e − j 2 π ν ( T ( v 1 + v 2 )+ u 1 + u 2 +2 τ ) M X X ∗ [ v 1 T + u 1 + τ , T ( v 1 − v 2 ) + u 1 − u 2 + τ ] M ∗ X X ∗ [ v 1 T + u 1 , T ( v 1 − v 2 ) + u 1 − u 2 − τ ] . W e deduce that with x = v T + u 1 and τ 0 = u 1 − u 2 that the relation of b A X X ∗ ( ν, τ ] is =  O (log ( n )) if | τ | > T P N − τ − 1 x =0 P T − 1 τ 0 = − ( T − 1) e − j 2 π ν [2 x +2 τ − τ 0 ] M X X ∗ [ x + τ , τ 0 + τ ] M ∗ X X ∗ [ x, τ 0 − τ ] if | τ | < T . (A-8) DEP AR TMENT OF ST A TISTICAL SCIENCE RESEARCH REPOR T 293 20 This completes the proof of the order structure of the EAF of an underspread process. T o prove a CL T for the EAF we add the extra assumptions stated in the theorem. W e need to use results for random variables that are both non-stationary and correlated to determine large same results. W e note that for τ > 0 b A X X ∗ ( ν, τ ] = N − τ − 1 X t =0 Q N −| τ | [ t ] + N − τ − 1 X t =0 e − j 2 π ν ( t + τ ) c M X X ∗ [ t + τ , τ ] . W e note that  Q N −| τ | [ t ] , t = 0 , . . . , N − τ − 1  is a triangular array of centred random variables, that these by assumption are strongly mixing, and have finite second moments. The constraint Peligrad [28] makes on the correlation is satisfied because X [ t ] is strictly underspread, and the Hilbert transform only induces suitably decaying correlation from such a process. W e hav e assumed the Lindeberg condition holds, and note (17) from our assumptions. Hence from Theorem 2.1 of Peligrad [28] the theorem follows. The condition in (17) does not have to be stated for the real and imaginary part separately as the relation divided by the variance goes to zero for ( ν, τ ) / ∈ D . W e can deduce the exact form of the mean, variance and correlation from the previous part of the theorem. R E F E R E N C E S [1] D. L. Donoho, S. Mallat, and R. von Sachs, “Estimating cov ariances of locally stationary processes: rates of con vergence of best basis methods, ” Statistics, Stanford Univ ersity , Standford, California, USA, T ech. Rep., 1998. [2] G. Matz and F . Hlawatsch, “Nonstationary spectral analysis based on time-frequenc y operator symbols and underspread approximations, ” IEEE T rans. Inf. Theory , vol. 52, pp. 1067–1086, March 2006. [3] M. B. Priestley , “Evolutionary spectra and non-stationary processes, ” J. Roy . Stat. Soc. B , vol. 47, no. 2, pp. 204–237, 1965. [4] G. E. Pfander and D. F . W alnut, “Measurement of time-variant linear channels, ” IEEE T rans. Inf. Theory , vol. 52, pp. 4808–4820, Nov ember 2006. [5] R. E. Blahut, W . Miller, and C. H. Wilcox, Eds., Radar and Sonar , P art I . New Y ork, USA: Springer V erlag, 1991. [6] P . Flandrin, T ime-F requency/T ime-Scale Analysis . San Diego, CA: Academic Press, 1999. [7] N. Ma and J. T . Goh, “ Ambiguity-function-based techniques to estimate doa of broadband chirp signals, ” IEEE T rans. Signal Process. , vol. 54, pp. 1826–1839, May 2006. [8] M. Jachan, G. Matz, and F . Hlawatsch, “T ime-frequency arma models and parameter estimators for underspread nonstationary random processes, ” IEEE T rans. Signal Pr oc. , vol. 55, pp. 4366–4381, 2007. [9] P . Fryzle wicz and G. P . Nason, “Haar-Fisz estimation of ev olutionary wav elet spectra, ” J. Roy . Stat. Soc. B , vol. 68, pp. 611–634, 2006. [10] R. von Sachs and K. Schneider, “W avelet smoothing of ev olutionary spectra by nonlinear thresholding, ” ACHA , vol. 3, pp. 268–282, 1996. [11] R. A. Hedges and B. W . Suter, “Numerical spread: Quantifying local stationarity , ” Digital Signal Pr ocessing , vol. 12, pp. 628–643, 2002. [12] I. M. Johnstone and B. W . Silverman, “W avelet threshold estimators for data with correlated noise, ” J. Roy . Stat. Soc. B , vol. 59, pp. 319–351, 1997. [13] H. Cram ´ er , “On the theory of stationary random processes, ” Ann. Math. , vol. 41, pp. 215–230, January 1940. [14] W . Kozek, “Matched W eyl-Heisenber g expansions of nonstationary en vironments, ” Ph.D. dissertation, Universit ¨ at W ien, V ienna, Austria, 1996. [15] G. Matz, “ A time-frequency calculus for time-varying systems and non-stationary processes with applications, ” Ph.D. dissertation, Univ ersit ¨ at W ien, V ienna, Austria, 2000. [16] J. Jeong and W . J. W illiams, “ Alias-free generalized discrete-time time-frequency distributions, ” IEEE T rans. Signal Pr ocess. , vol. 40, pp. 2757–2765, November 1992. [17] F . D. Neeser and J. L. Massey , “Proper complex random processes with applications to information theory , ” IEEE T rans. Inf. Theory , vol. 39, pp. 1293–1302, July 1993. [18] P . J. Schreier and L. L. Scharf, “Stochastic time-frequency analysis using the analytic signal: Why the complementary distribution matters, ” IEEE T rans. Signal Pr ocess. , vol. 51, pp. 3071–3079, December 2003. [19] L. Auslander and R. T olimieri, “Characterizing the radar ambiguity functions, ” IEEE T rans. Inf. Theory , vol. 30, pp. 832–836, November 1984. [20] ——, “Radar ambiguity functions and group-theory , ” SIAM Journal on Mathematical Analysis , vol. 16, pp. 577–601, 1985. [21] E. Bekir , “Unaliased discrete-time ambiguity function, ” JOSA , vol. 94, pp. 817–826, 1993. [22] R. T olimieri and R. S. Orr , “Poisson summation, the ambiguity function and the theory of weyl-heisenberg frames, ” The Journal of F ourier Analysis and Applications , vol. 1, pp. 233–247, 1995. [23] T . S. Ferguson, A Course in Lar ge Sample Theory . London, UK: Chapman & Hall, 1996. [24] S. Olhede, “Hyperanalytic denoising, ” IEEE T rans. Image Pr ocess. , vol. 16, pp. 1522–1537, June 2007. [25] L. W asserman, All of Nonparametric Statistics . New Y ork, USA: Springer V erlag, 2007. [26] A. A verb uch, R. R. Coifman, D. L. Donoho, M. Elad, and M. Israeli, “Fast and accurate polar fourier transform, ” Applied and Computational Harmonic Analysis , vol. 21, pp. 145–167, 2006. [27] L. Isserlis, “On a formula for the product-moment coef ficient of any order of a normal frequenc y distrib ution in any number of variables, ” Biometrika , vol. 12, pp. 134–139, 1918. [28] M. Peligrad, “On the asymptotic normality of sequences of weak dependent random variables, ” J. Theo. Prob . , vol. 9, pp. 703–715, 1996.