On the Underspread/Overspread Classification of Random Processes

QUADRA TIC TIME–V AR YING SPECTRAL ESTIMA T ION F OR UNDERSPREAD PR OCESSES ∗ W erner Kozek NUHA G, Dept. of Mathematics Univ ersit y of Vienna Strudlhofg. 4, A–1090 Vienna, Austria (k ozek@t yc he.mat.univie.ac.at) Kurt Riedel Couran t Institute of Mathematical Sciences New Y ork Univ ersit y New Y o rk, New Y ork 10012 (riedel@cims.n yu.edu) Abstract: Time–v a rying sp ectral estima tion is studied for nonstationary pro cesses with restricted time–frequency (TF) correlation. Explicit bias and v ariance expressions are giv en for quadr atic TF–in v arian t (QTFI) estimators of an exp ected real–v a lued QTFI representa tion based on a sin gle noisy observ ation. Un biased theoretical estimators with globally minimal v ariance are derived and appro ximately realized b y a matc hed m ulti-windo w metho d. 1 INTR ODUCTION Time-frequency distrib ution are w id ely used to searc h for h idden structur e in the signal. When the signal consists of a small n umb er of slo wly v arying sinusoids, the W entz el-Kramer-Brillioun represent ation reduces the signal to curv es in the time-frequency plane [16]. W e consider the case of n onstationary s to chastic p ro cesses with u nderlying time-frequ en cy structure in the correlation op erator. The ev olutionary sp ectrum is one common represen tation of nons tationary pro cesses. In [15], we prop ose estimating the ev olutionary sp ectrum b y smo othing the log-sp ectrogram using a data-adaptiv e kernel s mo other in the time-frequency p lane. The evo lutionary sp ectrum has t wo adv an tage s: it is alw a ys p ositiv e and it con v erges to the sp ectrum as the r atio of th e characte ristic time s cale to the sampling rate b ecomes large. Its disadv an tages are its lac k of uniqueness and its relativ ely p o or time f requency resolutions. W e consider a d ifferen t class of represen tations of n onstationary pro cesses: qu adratic Cohen’s class s p ectra. These represen tations c orresp ond to the exp ected v alue of Cohen’s class time- frequency represen tations. An imp ortan t mem b er of this class is th e W igner-Ville sp ectrum. This class of sp ectral represent ations p ossess useful op erator prop erties and a repro d ucing k ernel Hilb ert space stru cture. In this article, we consider th e estimation problem: ho w to estimate the C ohen’s class sp ectra. This same problem has b een considered by Sa y eed and Jones as w ell. In [19], a complete knowle dge of the correlation op erator is assumed. This assumption is appr opriate f or the signal classificiation problem of recognizing one or more sp ecific signals. Our app roac h assumes m uc h we ake r a priori kno wledge. W e assu me only that the signal is und er s pread whic h corresp on d s to b eing double band limited in th e am biguity plane. In Section 2, w e review time-frequency representa tions of deterministic signal s. Section 3 present s the analogous theory for time v arying sp ectra. S ection 4 d efi nes and motiv ate s u nder- spread pr o cesses. Section 5 analyzes the bias and v ariance of a sp ecial class of quadratic estimators of Cohen’s class sp ectra. Secti on 6 d etermines the minimum v ariance unbiase d estimato r of an undersp read ∗ F unding b y g rant 4913 of the Jubil¨ aumsfonds der ¨ osterr eichisc hen Nationalb a nk and t he US D ept. of Ener gy. pro cess. Section 7 describ es a related estimation us ing m ultiple win do ws. Section 8 pr esen ts a biomedical example. 2 QUADRA TIC TIME-FREQUENCY DISTRIBU- TIONS Ev ery real–v alued quadratic time–frequency (TF) sh ift–in v arian t (QTFI) representa tion of a signal x ( t ) can b e r ep resen ted as a quad r atic form [1 ]. Compr eh ensiv e reviews of Coh en ’s class are giv en in [2, 4]. W e no w cast Cohen ’s class of time-frequency representa tions in an op erator th eoretic framew ork. Let P b e a self–adjoin t Hilb ert-Sc hmid t (H-S) “pr otot yp e” op erator. W e d efine the quadratic time-frequency shift in v arian t (QTFI) distribution T x ( t, f ) = h P ( t,f ) x, x i , (1) where P ( t,f ) is a TF shifted v ersion of the self–adjoin t protot yp e op erator P . The choic e of the kernel, p ( s, t ), determines a particular repr esentati on in Cohen ’s class. The TF–shiftin g of op erators is defined as P ( t,f ) def = S ( t,f ) PS ( t,f )+ , w here S ( τ ,ν ) is a u nitary TF shift op erator, acting as  S ( τ ,ν ) x  ( t ) = x ( t − τ ) e − j 2 πν t . The standard H-S inner pro du ct, < R , P > , is < R , P > = Z Z r ( t, s ) p ( s, t ) dtds , where r ( t, s ) and p ( s, t ) are the resp ectiv e kernels of th e H.-S. op erators R and P . Throu gh ou t this article, w e assume an infinite time domain and suppress replace R ∞ −∞ dt with R dt . W e n o w r eview the basic u n itary TF r epresen tations of HS op erators [6]. The generalized W eyl symb ol is defined as L ( α ) H ( t, f ) def = Z t h  t +  1 2 − α  τ , t −  1 2 + α  τ  e − j 2 πf τ dτ , where | α | ≤ 1 / 2. The W eyl corresp ondence is giv en by α = 0, and the Kohn-Nirenberg corr e- sp onden ce (time–v a ryin g transfer fun ction) by α = 1 / 2 [20]. (When we su ppress the su p erscript, this means v alidit y for any α .) Th e TF shifting of op erators corresp onds to a shift of the s ym- b ol, L P ( τ ,ν ) ( t, f ) = L P ( t − τ , f − ν ) , which shows that wh en ev er P is a TF lo calization op erator that selects signals cen tered in the origin of the TF plane then P ( t,f ) lo calizes signal comp onents cen tered aroun d ( t, f ). The gener alize d spr e ading fu nction (GSF) of a linear op erator [6 ] is S ( α ) H ( τ , ν ) def = Z t h  t +  1 2 − α  τ , t −  1 2 + α  τ  e − j 2 πν t dt. The GSF is the symplectic F ourier transform of the generalized W eyl symb ol L ( α ) H ( t, f ): S ( α ) H ( τ , ν ) = Z t Z f L ( α ) H ( t, f ) e − j 2 π ( ν t − τ f ) dtd f , (2) L ( α ) H ( t, f ) = Z τ Z ν S ( α ) H ( τ , ν ) e − j 2 π ( − ν t + τ f ) dτ dν . (3) When the W eyl sy mb ol is sm o othly v arying in time and frequ ency , then the generalizing spr eading function deca ys in τ and ν . 2 3 TIME V AR YING SPECTR UM F or a n onstationary pro cess, a time–v ary in g sp ectrum ma y b e d efined as the exp ectation of (1) P x ( t, f ) def = E n h P ( t,f ) x, x i o , (1) A pr ominen t example for P x ( t, f ) is th e Wigner–Vil le sp ectrum [13]. Priestley’s ev olutionary sp ectrum [14, 15] is a differen t, p opular d efinition of a sto c hastic time–v arying s p ectrum that cannot b e b rough t into the form of (1). W e consider cir cu lar complex, zero–mean Gaussian pro cesses with trace–cla ss correlation ker- nel ( R x ) ( t, t ′ ) = r x ( t, t ′ ) = E  x ( t ) x ∗ ( t ′ )  , tr R x < ∞ . The tr ace–c lass conv entio n imp lies a HS inn er p ro duct representa tion of P x ( t, f ), alternativ ely written as the trace of the pr o duct op erator: P x ( t, f ) = h R x , P ( t,f ) i = tr n R x P ( t,f ) o . (2) The exp e cte d ambiguity function is defi ned as the GSF of the correlation op erator [8] E A ( α ) x ( τ , ν ) def = S ( α ) R x ( τ , ν ) (3) With the generalized Wigner–Vill e sp ectrum, defined as E W ( α ) x ( t, f ) def = L ( α ) R x ( t, f ) , Eq. (2) carries o v er to a “nonstationary Wiener–Kh in tc hine relation”: E W ( α ) x ( t, f ) = Z τ Z ν E A ( α ) x ( τ , ν ) e j 2 π ( ν t − τ f ) dτ dν , (4) E A ( α ) x ( τ , ν ) = Z t Z f E W ( α ) x ( t, f ) = e j 2 π ( τ f − ν t ) dt d f . (5) These relation s h ips are summ arized in T able 1. L ( α ) R ( t, f ) = E W ( α ) ( t, f ) W eyl symb ol of correlation Generalized W-V sp ectrum ⇑ t ↔ ν ⇑ t ↔ ν ⇓ f ↔ τ ⇓ f ↔ τ S ( α ) R ( t, ν ) = E A ( τ , ν ) GSF of correlation Exp ected ambiguit y function T able 1 As an example, the real–v alued generalize d Wigner–Ville sp ectrum can b e written as Re n E W ( α ) x ( t, f ) o = E n h P ( t,f ) ( α ) x, x i o , where the α –dep en d en t p rotot yp e op erator is giv en b y: S (0) P ( α ) ( τ , ν ) = cos(2 π τ ν α ) . (6) 3 Since b oth the W eyl symb ol and the spreading f unction are unitary representa tions of HS op erators we can rewrite the general time–v arying sp ectrum, h R x , P ( t,f ) i = h E W x , L P ( t,f ) i = h E A x , S P ( t,f ) i . Note furtherm ore th at th e GS F o f the TF shifted pr otot yp e op erator is just a mo d ulated v ersion of the GSF of th e original v ersion: S P ( t,f ) ( τ , ν ) = S P ( τ , ν ) e j 2 π ( ν t − τ f ) , th us in particular | S P ( t,f ) ( τ , ν ) | = | S P ( τ , ν ) | . 4 UNDERSPREAD PR OCESSES The bias-v ariance analysis of Sec. 5 is v alid for an y circular Gaussian pr o cess with a trace class co v ariance. W e no w restict our consideration to th e case wh er e the pro cess’ exp ected ambiguit y function E A ( α ) x ( τ , ν ) is zero ou tsid e a rectangle in the am biguit y p lane. Ou r r equiremen t that the exp ected ambiguit y is double b and-limited implies th at the W eyl sym b ol is smo oth in time and frequency . W e d enote the maxim um temp oral correlation width τ max and the maximum sp ectral corre- lation width ν max ; i.e., we assume that the exp ected am biguit y function satisfies E A ( α ) x ( τ , ν ) = E A ( α ) x ( τ , ν ) χ x ( τ , ν ) , (1) where χ x ( τ , ν ) is the 0 / 1–v alued indicator fun ction of a cen tered rectangle with area s x = 4 τ max ν max . Ac cording to the recen tly introd uced t erminology w e c all a pr o cess with s x < 1 underspr e ad and in the conv erse case overspr e ad [7]. F or asymptotics we assum e th at the un- derspread parameter is very small: s x ≪ 1. The undersp read parameter, s x , corresp onds to the expansion parameter 1 / ( τ λ f ), which is u sed in the analysis of ev olutionary sp ectra [15]. As to the relev ance and realiza bility of the un derspread pro cesses we n ote that practically imp ortant linear time–v ary in g (L TV) systems, as e .g. the mobile radio c hann el or underw a- ter a coustic c hannel [21], are c haracterized b y an (at least in goo d approximat ion) restricted spreading function (this is the field where the un derspread/o verspread terminology w as originally in tro du ced). No w, we apply sta tionary white n oise n ( t ) with E { n ( t ) n ∗ ( t ′ ) } = δ ( t − t ′ ) to an undersp read L TV system H c haracterized by S ( α ) H ( τ , ν ) = S ( α ) H ( τ , ν ) χ H ( τ , ν ) where χ H ( τ , ν ) co vers a cent ered rectangle with h alfwidths τ max,H and ν max,H . Then the output pro cess x ( t ) = ( H n )( t ) is n on s tationary with correlation R x = HH + . Applying the triangle inequalit y to the spreading fun ction of the pro duct op erator [11] give s | E A x ( τ , ν ) | ≤ | S H ( τ , ν ) | ∗ ∗| S H ( τ , ν ) | , where the ∗∗ denotes double con vo lution. The output pr o cess is th us u n derspread with τ max,x = 2 τ max,H and ν max,x = 2 ν max,H . Hence, w e ha v e shown that und erspread pro cesses are realizable and relev ant. In view of the “nonstationary Wiener–Khin tc hine relation” (4), the ov ersp read/undersp read classification ma y b e in terpreted as a smo othness condition for th e time–v ary in g sp ectrum of the pro cess. Applying th e sampling theorem on the symbol lev el leads to a discr ete Weyl–Heisenb er g exp ansion of the correlatio n op er ator [11]: R x = X l X m E W ( α ) x ( lT , mF ) P ( lT ,mF ) ( α ) 4 v alid for a sampling grid w ith T ≤ 1 2 ν max and F ≤ 1 2 τ max and the protot yp e op erator defined b y S ( α ) P ( α ) ( τ , ν ) = χ x ( τ , ν ) . (2) The critical spr ead s x = 1 corresp onds to the Nyquist sampling density T F = 1. Hence, consider- ing bandlimited pro cesses, for s x = 1 the rate of inno v ation in the pro cess second order statistics is equ al to the s amp ling rate of th e realization [8]. Ho w ev er, a robu st estimation pro cedu re maps a time series with N samples on a mo del with less than N co efficien ts suc h that th e cr itical spread is a treshold for robu st estimatio n of the generalized Wigner–Ville sp ectrum. It is furtherm ore remark able that the ev olutionary sp ectrum of an undersp read pro cess is 2D band limited in exactly the same manner as the generalized Wigner–Ville sp ectrum [11]. It should b e noted that one can view th e stat ionarit y assumption underlying any time– in v ariant sp ectrum estimation as a limit case of (1 ) since the exp ected am biguity function of a wide–sense stationary pro cess is c haracterized by ideal concen tr ation on th e τ –axis: E A x ( τ , ν ) = r x ( τ ) δ ( ν ) , where r x ( τ ) is the auto correlation function. The Wiener–Kh in tc hine relation requires strict band-limiting th e am biguity plane. The r e- mainder of our analysis requires only a concen tration in the am biguity plane with characte ristic spread, s x ≪ 1, but not complete band limitation. 5 REPR ODUCI NG KERNEL HILBER T SP A CE W e n o w sh o w that time-frequency distrib utions are a repro ducing k ernel Hilb ert sp aces (RKHS) [5] using the Wigner-Ville k ernel. A RK HS is Hilb ert space H of complex v alued fun ctions, defined on a set S , that has a r epro ducing k ernel K ( s, t ) defined on S × S with tw o prop erties: (i) for eac h t, th e fun ction K(s,t) lies in H and (ii) for eac h x ∈ H and eac h t ∈ S one has the repro d u cing prop erty: x ( t ) = h x, K ( ., t ) i = Z t ′ K ∗ ( t ′ , t ) x ( t ′ ) dt ′ . In our case, the Hilb ert space, H , is th e set of W eyl sym b ols of un derspread op erators H satisfying a giv en spreadin g constraint (1). The repro du cing k ernel is giv en by the W eyl symbol of the protot yp e op erator: K ( t ′ , f ′ , t, f ) = L P ( t,f ) ( t ′ , f ′ ) . This is in f act a repro ducing k ernel as (i) for eac h ( t, f ) P ( t,f ) remains und erspread sin ce | S P ( t,f ) ( τ , ν ) | = | S P ( τ , ν ) | , and (ii) one has the repr o ducing form ula as follo ws: L H ( t, f ) = h L H , L P ( t,f ) i = Z t ′ Z f ′ L H ( t ′ , f ′ ) L P ( t ′ − t, f ′ − f ) dt ′ d f ′ (1) WERNER: DO YOU MEAN 1 for th e Wigner Ville kernel or for the k ernel in2 or WHA T ? 5 6 QTFI ESTIMA TION W e no w consider QTFI estimators of the time v arying sp ectrum of the signal pr o cess x ( t ) when it is con taminated with noise. W e are giv en a single noisy observ ation, y ( t ) of the signal pro cess x ( t ): y ( t ) = x ( t ) + n ( t ) with E  n ( t ) n ∗ ( t ′ )  = σ 2 n δ ( t − t ′ ) , where n ( t ) is statistically in dep endent, zero–mean, cir cular complex Gaussian w hite noise. T o estimate P x ( t, f ) we u se a generally different QTFI representat ion of the observ ation: b P x ( t, f ) = h b P ( t,f ) y , y i . W e define the “bias op erator” e P as e P def = b P − P . The QTFI estimator is consistent with classical, “non–parametric” time–in v ariant sp ectrum estimation wh ere th e p redominant class of estimators [22, 18] can b e basically written as a frequency parametrized quadratic form: ˆ S x ( f ) = h b P (0 ,f ) y , y i . (1) 7 BIAS AND V ARIANCE ANAL YSIS With the statistic al indep endence of signal and noise and usin g (2) we h a v e the follo wing exp ec- tation of the estimate: E n b P x ( t, f ) o = E n h b P ( t,f ) x, x i o + σ 2 n tr b P , (1) suc h that the bias is give n by B ( t, f ) def = E n b P x ( t, f ) o − P x ( t, f ) = tr n e P ( t,f ) R x o + σ 2 n tr b P . Using the S ch w arz inequalit y for op erator inner pro ducts and tr iangle inequalit y , we imm ediately get a tigh t b oun d for the maxim um b ias: | B ( t, f ) | ≤ k e P kk R x k + σ 2 n | tr b P | , (2) where the op erator norm is the HS norm. W e assume kno wledge of the noise leve l σ 2 n suc h that w e can trivially correct the T F–indep endent bias term: b P ′ x ( t, f ) = b P x ( t, f ) − σ 2 n tr b P , where b P ′ x ( t, f ) denotes the corrected estimate. The v ariance, V ( t, f ) def = E n b P 2 x ( t, f ) o −  E n b P x ( t, f ) o 2 , is ev aluated us ing of Isser lis’ fourth order momen t form ula (for the sp ecial case of circular complex v ariables), E { x ( t 1 ) x ∗ ( t 2 ) x ( t 3 ) x ∗ ( t 4 ) } = r x ( t 1 , t 2 ) r x ( t 3 , t 4 ) + r x ( t 1 , t 4 ) r x ( t 3 , t 2 ) , one h as: V ( t,f ) = tr   b P ( t,f ) R x  2  + 2 σ 2 n tr   b P ( t,f )  2 R x  + σ 4 n k b P k 2 . The Sch warz inequalit y f or the op erator inner pro d uct leads to a b ound on the maxim um v ariance, V max ≤ k b P k 2  k R x k + σ 2 n  2 , (3) prop ortional to the HS norm of the protot yp e op erator b P . 6 Global Mean Square E rror. The bias an d v ariance results are complicated TF–dep endent ex- pressions. Due to our restriction to QTFI estimators we need TF–in v ariant, th us global ind icators for the estimator p erformance. After correcting for the TF ind ep endent b ias term, B 0 def = σ 2 n tr b P , w e charact erize the global square bias as follo ws: B 2 tot def = Z t Z f  B ( t, f ) − σ 2 n tr b P  2 dt d f = h    S e P    2 , | E A x | 2 i . (4) Just as for the bias we giv e a global c haracterization of the v ariance. Th e TF indep end en t term is giv en b y: V 0 = σ 4 n k b P k 2 . W e d efine a total v ariance as the in tegral o ver the TF dep endent v ariance terms, one has: V tot def = Z t Z f ( V ( t, f ) − V 0 ) dt d f = k b P k 2  tr R 2 x + 2 σ 2 n tr R x  . (5) Equations (4) and (5) are derived in the app en dix. Observe that any of the glob al v arianc e c onstants; i.e., the maximum varianc e V max , the TF– indep endent varianc e term V 0 , and the total varianc e V tot ar e pr op ortional to the HS norm of the pr ototyp e op er ator: V 0 , V max , V tot ∝ k b P k 2 . (6) 8 ESTIMA TOR OPTIMIZA TION Classical sp ectrum estimation pro du ces smo oth sp ectra since — du e to the absence of a m o del — smo othing is the actual tool for v ariance reduction. The prop osed estimators usually are the result of mean–squared error considerations. In the present w ork, w e d eviate from this p oint of view in a pragmatic wa y: we restrict ours elves to u n derspr ead pro cesses wh ose true sp ectra are itself smo oth (in the sense of 2D band limitation) suc h that there exist a whole class of u n biased estimators. While suc h a mo delling ingredient ma y b e questionable for time–in v ariant sp ectrum analysis we feel that it is necessary for time–v arying sp ectral estimation. T he reason lies in the often o ve rlo oked p oin t that fr equency parametrization is matc hed to an y stationary p ro cess (the F ourier trans f orm d iagonaliz es th e correlation op erator) while TF parametrization is not matc hed to a ge ne r al nonstationary pro cess. F rom the p oint of view of op erator d iagonalization it is the class of u nderspr ead pro cesses wh ere TF–parametrizatio n is appropriate [11]. Un biased estimation w ithout fu rther assump tion on the signal pro cess x ( t ) requires a v anishing “bias op erator”, i.e., P = b P . In the case of the generalized Wigner–Ville sp ectrum, the protot yp e op erator (cf. (6)) is not HS since k P k 2 = Z τ Z ν | S P ( τ , ν ) | 2 dτ dν , (1) so that one c an exclude finite–varianc e unbi ase d e stimation of the gener alize d Wigner–Vil le sp e c- trum without a priori know le dge on the pr o c ess. This is w ell–kno wn [13]. Based up on the kno wn supp ort of E A x ( τ , ν ) one has a large cl ass of non trivial unbiased estimators (with n on v anish ing “bias op erator”, ) S ( α ) b P U B ( τ , ν ) = ( S ( α ) P ( τ , ν ) , where E A x ( τ , ν ) 6 = 0 arbitrary , where E A x ( τ , ν ) = 0 . W e interpret minimum v ariance in the sen s e of th e com bin ed consideration of the global v ariance constan ts V 0 , V max , V tot . Due to (6) one h as to select the u n biased estimator with minimum HS 7 norm p rotot yp e op erator. Using (1) this turns out to b e trivial: the m inim um–v ariance u n biased (MVUB) QTFI estimator is obtained by setting the spreading fu nction of the protot yp e op er ator zero whereve r p ossible: S ( α ) b P M V U B ( τ , ν ) = ( S ( α ) P ( τ , ν ) , where E A x ( τ , ν ) 6 = 0 0 , where E A x ( τ , ν ) = 0 . When χ x ( τ , ν ) is th e smallest in dicator function conta ining th e sup p ort of E A x ( τ , ν ), then the MVUB QTFI estimator can b e written as: S ( α ) b P M V U B ( τ , ν ) = S ( α ) P ( τ , ν ) χ x ( τ , ν ) . In particular, for the α –parametrized real–v alued generalized Wigner–Ville sp ectrum one h as: S (0) b P M V U B ( α ) ( τ , ν ) = cos(2 π τ ν α ) χ x ( τ , ν ) . This estimator is optimal among all QTFI estimators thus in the sense of global v ariance minimiza- tion. The estimate is lo cally stable since it minimizes a b ound on the maxim um v ariance ( V max ) and it is unbiase d for arb itrary time and fr equ ency , but it deviates from the lo cal TF–dep end en t MVUB estimate. Mean–Squared Error. The theoretical MVUB estimator serve s wel l as a s tarting p oin t for obtaining practical estimators w ith go o d mean–squared err or p erformance. The mean squared error is giv en b y E ( t, f ) = V ( t, f ) + B 2 ( t, f ). F or any pro cess that satisfies th e spr eading constrain t (1) one can f ormally redefin e the estimation target via th e protot yp e op erator of an y unbiased estimator: P x ( t, f ) = tr n R x P ( t,f ) o = tr n R x b P ( t,f ) U B o , so that one can obtain a u seful b ound on the integrat ed mean–squared error E tot < k e P k 2 tr 2 R x + k b P k 2  k R x k 2 + 2 σ 2 n tr R x  , (2) with e P = b P − b P M V U B . This b oun d is based on (4), (5) and h    S e P    2 , | E A x | 2 i < k e P k 2 tr 2 R x . 9 MA TCHED MUL TI-WINDO W ESTIMA TOR The eigenfunction decomp osition of the protot yp e op erator P shows that P ( t,f ) is a we ight ed su m of rank one pro jections. Equiv alen tly , an y QTFI represen tation can b e written as a w eight ed sum of sp ectrograms with orthonormal windo ws [20 ]. F or practicalit y , we require our estimator to b e based on a fi nite–rank protot yp e op erator with finite–length eigenfun ctions. The MVUB estimator of S ec. 6 do es not s atisfy these requirements. Thus, we choose the finite–rank, time– limited estimator whic h minimizes the up p er b oun d on th e in tegrated mean–squared error as giv en b y (2). When we imp ose the additional requir emen t that the protot yp e op erator b e pro jection t yp e with normalized trace, b P has the represen tation: b P N = 1 N N X k =1 γ k ⊗ γ k (1) where γ k ⊗ γ k denotes the rank–one pro jection on the orthonormal w indo w fu nctions γ k ( t ) and N is th e rank. I n th is case, k b P N k 2 = 1 / N , and the optimization of (2) reduces to minimizing 8 k b P M V U B − b P N k 2 sub ject to orthonormalit y constr aints on the γ k . W e define the matc hed m ulti- windo w estimator as the q u adratic f orm b ased on a pr otot yp e op erator the min imizes E tot sub ject to (1). The optimizatio n is p erf orm ed in a t wo step p ro cedure: w e optimize the windo ws sub ject to a fi xed r an k and then w e optimize the r ank. F or practicalit y , we imp ose that the γ k ha v e supp ort on [ − T / 2 , T / 2]. T o imp ose this time lo calization on the optimization of b P N , w e define T as the pro jection onto the cen tered interv al and require b P N = T b P N T . Minimizing (2) yields the optimal win do ws equation: T b P M V U B T γ k ,opt = λ k γ k ,opt (2) The optim um window set is in dep end ent of N . F or the sp ecific case wh ere b P M V U B is an ideal bandpass (whic h ma y b e consider ed as a th eoretically optimal estimator f or s tationary pro cesses) (2) yields the time–limited and optimally band limited pr olate spherio d al wa ve fu nctions consistent with [22, 18]. A more realistic and simpler family of tap ers are the discrete sinusoidal tap er s , { v ( k ) } , where v ( k ) n = q 2 N +1 sin π k n N +1 , and N is the num b er of p oint s [17]. The resulting sinusoidal multi -tap er sp ectral estimate is ˆ S ( t, f ) = 1 2 K ( N + 1) P K j =1 | y ( t, f + j 2 N +2 ) − y ( t, f − j 2 N +2 ) | 2 , where y ( t, f ) is the lo cal F ourier transf orm cen tered at time t w ith length N . S ( t, f ) is the instanta nteo us sp ectral densit y , and K is the n umb er of tap er s . The sin usoidal tap ers are asymptotically optimal wh en the bias error is lo cal. 10 FREE P ARAMETER OPTIMIZA TION F or a strongly un derspread pro cess s x ≪ 1, S P ( τ , ν ) is approximate ly constan t in the supp ort of E A x ( τ , ν ). Using the optimal windo w fun ctions of (2), w e approximat e b P M V U B with 1 /s x suc h rank one pro jections. In this case, k b P M V U B − b P N k 2 reduces to ( s x − 1 / N ). Optimizing (2) with resp ect to N f or mo derate noise lev el yields N opt ≈ 1 s x , for σ 2 n tr R x < 1 − s x 2 . WERNER: YOUR ES T IMA TE of s x and T ... 11 CONCLUSIONS W e h a v e studied time–v aryin g sp ectral estimation via quadratic TF–in v ariant estimators. F or circular complex Gauss ian signal and noise pro cesses w e ha v e presente d explicit (lo cal an d global) bias and v ariance results. F or th e sp ecific case of an und er s pread pro cess the design of matc h ed m ulti-window estimators has b een b ased on approximat ing a theoretical MVUB estimator. The theoretical MVUB est imator a s deriv ed in Section 6 is a sp ecific case of the rece ntly prop osed optim um ke rn el d esign for Wigner–Ville sp ectrum estimation [19]. W e emph asize that [19] requir es a complete kn owledge of a second order statistic what mak es th is appr oac h p urely theoretical wh ile our p rop osed estimator uses a more realistic, incomplete a pr iori kno wledge of the pr o cess statistics. F or Cohen’s class time v arying sp ectra. Using the repro du cing ke rn el Hilb ert space formalism, w e d eriv e expressions for the leading order bias and v ariance. Unders p read pro cesses are band limited in the am biguit y plane and smo oth in the time fr equency domain. F or u nderspr ead pro cesses, we give unbiase d minim um v ariance estimators. APPENDIX: PR OOFS 9 W e n o w derive (4) wh ich equates the int egral square bias with the inner pr o duct of the squared GSF of the “bias op erator” and the pro cess’ exp ected am biguit y function: Z t Z f tr 2 n e P ( t,f ) R x o dt d f = Z t Z f    h S e P ( t,f ) , E A x i    2 dt d f = Z t Z f Z τ 1 Z ν 1 Z τ 2 Z ν 2 S e P ( τ 1 , ν 1 ) E A x ( τ 1 , ν 1 ) S ∗ e P ( τ 2 , ν 2 ) E A ∗ x ( τ 2 , ν 2 ) · e − j 2 π [ ( ν 1 − ν 2 ) t − ( τ 1 − τ 2 ) f ] dt d f dτ 1 dν 1 dτ 2 dν 2 = = Z τ Z ν    S e P ( τ , ν )    2 | E A x ( τ , ν ) | 2 dτ dν = h    S e P    2 , | E A x | 2 i . Deriv ation of (5): P ( t,f ) P ( t,f ) = S ( t,f ) P 2 S ( t,f )+ =  P 2  ( t,f ) , together with Z t Z f P ( t,f ) dt d f = tr { P } I (whic h follo ws directly f rom the trace f ormula of the W eyl corresp ondence [6]); as w ell as  P ( t,f ) R  ( s, s ′ ) = Z s ′′ p ( s − t, s ′′ − t ) e j 2 πf ( s − s ′′ ) r ( s ′′ , s ′ ) ds ′′ , whence Z t Z f tr   P ( t,f ) R  2  dt d f = Z t Z f Z s Z s ′ Z s 1 Z s 2 p ( s − t, s 1 − t ) r ( s 1 , s ′ ) · p ∗ ( s − t, s 2 − t ) r ∗ ( s 2 , s ′ ) e j 2 πf ( s 1 − s 2 ) dt d f ds ds ′ ds 1 ds 2 = Z t Z s Z s ′ Z s 1 | p ( s − t, s 1 − t ) | 2   r ( s 1 , s ′ )   2 dt ds ds ′ ds 1 = k P k 2 k R k 2 . References [1] L. Cohen, “Generalized ph ase–space distrib ution functions,” J. Math. Phys. , V ol. 7, 1966, pp. 781-786. [2] Cohen, L. (1989). “Time-frequency d istributions - a review.” Pr o c . I.E.E. E. V ol. 77, 941-9 81. [3] F. Hla watsc h and W. K ozek, “Time–frequency w eigh ting and displacemen t effects in linear time–v arying systems,” Pr o c. IEEE ISCAS–92 , Mai 1992 , p p. 1455 –1458. [4] Hla wa tsc h, F. and Boudreaux-Bartels G.F. (1992). “Linear and quadr atic time-frequency represent ations.” I.E.E.E. Signal Pr o c essing M ag. 9 , 21-67. [5] T. K ailath, “RKHS Approac h to Detectio n and Estimation Problems—Pa rt V: Pa rameter Estimation,” IEEE T r ans. Info. The ory , S ept. 1971, pp . 530–548. [6] W. Kozek, “On the generalized W eyl corresp ondence and its application to time–frequency analysis of linear time–v arying systems,” Pr o c. IE EE Int. Symp. on Time–F r e quency and Time–Sc ale Analysis , Victoria, C anada, Oct. 1992, p p. 167– 170. 10 [7] W. Kozek, “Optimally Karh un en–Lo ev e–lik e STFT expansion of nonstationary p ro cesses,” Pr o c. IEEE ICASSP–93 , Minn eap olis, MN, April 1993, V ol. 4, pp . 428–4 31. [8] W. Kozek, F. Hla watsc h, H. Kirchauer and U. T r aut w ein, “Correlativ e time–frequency anal- ysis and classification of non s tationary r andom pr o cesses,” Pr o c. IEEE Workshop on Time– F r e quency/Time–Sc ale Analysis , Ph iladelphia, P A, (199 4), 417–420 . [9] W. Kozek, H. G. F eic h tinger a nd T. Strohmer, “Time–frequency synthesis of statistic ally matc hed W eyl–Heisen b er g p rotot yp e s ignals,” this issue. [10] W. Kozek, “Matc hed generalized Gab or expansion of nonstationary p r o cesses,” Pr o c. IEE E Int. Conf. Signals, Systems, and Computers , P acific Gro ve , CA, Nov. 1993, pp. 499–503. [11] W. K ozek, “Matc hed W eyl–Heisen b erg Expansions of Nonstationary Environmen ts,” Ph.D. Th esis, Un iv ersit y of T ec hnology Vienna, und er review. [12] W. K ozek and K.S. Riedel, Quadratic time–v aryin g sp ectral estimation for underspread pro cesses; Pr o c. IEE E Worksho p on Time–F r e quency/Time–Sc ale A nalysis , Philadelphia, P A, (1994), 460–463. [13] W. Martin and P . Flandrin,“Wigner–Ville sp ectral analysis of n onstationary p ro cesses,” IEEE T r ans. A c oust., Sp e e ch, Signal Pr o c essing , V ol. 33, No. 6, Dec. 1985, pp. 1461–1 470. [14] M. Priestley , “Ev olutionary sp ectra and non–sta tionary pro cesses,” J. Ro y . Statist. S o c. S er. B, V ol. 27, 1965, pp . 204–2 37. [15] K. S. Riedel, “Optimal data based k ernel estimation of evol utionary sp ectra,” IE EE T r ans. Signal Pr o c essing , V ol. 41, No. 7, July 1993, pp. 2439– 2447. [16] K. S. Riedel, “Optimal kernel estimation of th e instanta neous frequency .” I.E.E .E. T r ans. Sig- nal Pr o c essing, V ol. 42 p . 2644-2 649, (Oct. 1994 ). [17] K.S. Riedel, A. Sidorenk o, “Minim um bias m ultiple tap er sp ectral estimati on.” I.E.E.E. T r ans. Signal Pr o c essing V ol. 43, p. 188-195 , (Jan. 1995). [18] K.S. Riedel, A. S idorenk o, D.J. Thomson, “Sp ectral estimation of plasma fluctuations I: Comparison of metho d s.” Physics of Plasmas , V ol. 1, Marc h 1994, pp. 485-500. [19] A. M. Say eed and D. L. Jones, “Optimal k ernels for Wigner–Ville s p ectral estimation,” Pr o c. IEEE ICASSP–94 , Sy d ney , Australia, Apr il 1994. [20] R. G. Sh eno y and T.W. P arks, “The W eyl corresp ondence and time–frequency analysis,” IEEE T r ans. Signal Pr o c essing , V ol. 42, No. 2, F eb. 1994 , pp . 318–331. [21] K. A. S ostrand, “Mathematics o f t he t ime–v aryin g c hann el,” Pr o c. of NA TO A dv. Study Inst. on Sign. Pr o c. V ol. 2 , Enschede (The Netherlands), 1968. [22] D. J. T hompson,“Sp ectrum estimation and harmonic analysis,” Pr o c. IE EE , V ol. 70, No. 9, Sept. 1982, p p.1055–1 096. 11