Minimum bias multiple taper spectral estimation

Minim um bias m ultiple tap er sp ectral estimation ∗ Kurt S. Riedel and Alexander Sidorenk o Couran t Institute of Mathematical Sciences, New Y ork Univ ersity New Y ork, New Y ork 10012-1185 EDICS: SP 3.1.1 Abstract Tw o families of orthonormal tap ers are prop osed for m ultitap er sp ectral analysis: minimum bias tap ers, and sinusoidal tap ers { v ( k ) } , where v ( k ) n = q 2 N +1 sin π kn N +1 , and N is the n umber of p oints. The re- sulting sin usoidal multitaper sp ectral estimate is ˆ S ( f ) = 1 2 K ( N +1) P K j =1 | y ( f + j 2 N +2 ) − y ( f − j 2 N +2 ) | 2 , where y ( f ) is the F ourier transform of the sta- tionary time series, S ( f ) is the sp ectral densit y , and K is the num b er of tap ers. F or fixed j , the sinusoidal tap ers con verge to the minimum bias tap ers lik e 1 / N . Since the sin usoidal tap ers hav e analytic ex- pressions, no n umerical eigenv alue decomp osition is necessary . Both the minimum bias and sinusoidal tap ers hav e no additional parame- ter for the sp ectral bandwidth. The bandwidth of the j th tap er is simply 1 N cen tered ab out the frequencies ± j 2 N +2 . Th us the bandwidth of the multitaper sp ectral estimate can b e adjusted lo cally b y simply adding or deleting tap ers. The band limited sp ectral concentration, R w − w | V ( f ) | 2 d f , of b oth the minimum bias and sin usoidal tap ers is very close to the optimal concentration achiev ed by the Slepian tap ers. In con trast, the Slepian tapers can ha ve the lo cal bias, R 1 / 2 − 1 / 2 f 2 | V ( f ) | 2 d f , m uch larger than of the minim um bias tap ers and the sin usoidal ta- p ers. ∗ The authors thank D. J. Thomson and the referees for useful commen ts. Researc h funded b y the U.S. Departmen t of Energy . 1 1 In tro duction W e consider a stationary time series, { x n , n = 1 . . . N } with a sp ectral den- sit y , S ( f ). A common estimator of the sp ectral density is to smo oth the square of the discrete F ourier transform (DFT) lo cally: ˆ S ( f ) = 1 (2 L + 1) N L X j = − L | y ( f + j N ) | 2 , (1) where y ( f ) is the F ourier transform (FT) of the stationary time series: y ( f ) ≡ P N n =1 x n e − i 2 π nf . Since (1) is quadratic in the FT, y ( f ), it is natural to consider a more general class of quadratic sp ectral estimators. W e examine quadratic estimators where the underlying self-adjoin t matrix has rank K , where K is prescrib ed. Using the eigen vector represen tation, the resulting quadratic sp ectral estimator can b e recast as a weigh ted sum of K orthonormal rank one sp ectral estimators. This class of sp ectral estimators w as originally prop osed by Thomson [19] under the name of multiple tap er sp ectral analysis (MTSA). W e refer the reader to [3, 8, 11, 12, 16, 18, 19] for excellen t exp ositions and generalizations of Thomson’s theory . In MTSA, a rank K quadratic sp ectral estimate is constructed b y c ho os- ing an orthonormal family of tap ers/sp ectral windo ws and then a v eraging the K estimates of the sp ectral density . In practice, only the Slepian tapers (also known as discrete prolate spheroidal sequences [17]) are routinely used for MTSA. In the present pap er, we prop ose and analyze tw o new orthonormal fam- ilies of tap ers: minimum bias (MB) tap ers and sin usoidal tap ers. The MB tap ers minimize the lo cal frequency bias, R f 2 | V ( f ) | 2 d f , sub ject to orthonor- malit y constrain ts, where V ( f ) is the DFT of the tap er. F or con tin uous time, the MB tap ers ha ve simple analytic expressions. The first tap er in the fam- ily is Papoulis’ optimal taper [9]. F or discrete time, the MB tap ers satisfy a selfadjoin t eigenv alue problem and may b e computed n umerically . In the case of discrete time, we define the k th sin usoidal taper, v ( k ) , as v ( k ) n = q 2 N +1 sin π kn N +1 , where N is the sequence length. The sinusoidal tap ers are an orthonormal family that conv erge to the MB tap ers with rate 1 / N as N → ∞ . These results are given in Section 3. Section 4 compares the lo cal bias, R 1 / 2 − 1 / 2 f 2 | V ( f ) | 2 d f , and the sp ectral concen tration, R w − w | V ( f ) | 2 d f , of the MB tapers, the sin usoidal tap ers and the Slepian tap ers. 2 In Section 5, w e sho w that the quadratic sp ectral estimator which mini- mizes the exp ected square lo cal error is w eighted m ultitap er estimate using the MB tap ers. A lo cal error analysis is given and the optimal num b er of tap ers is determined. A t frequencies where the sp ectral densit y is c hang- ing rapidly , few er tap ers should b e used. In Section 6, w e show that k ernel smo other sp ectral estimates [4, 10] are m ultitap er estimates and w e show that smo othing the logarithm of the multitaper estimate significan tly re- duces the v ariance in comparison with smoothing athe logarithm of a single tap er estimate. W e also describ e our data adaptiv e metho d for estimating the sp ectrum. In Section 7, we apply our sp ectral estimation tec hniques to real data and show that our tap ers outp erform the Slepian tapers whenever a v ariable bandwidth is needed. In the App endix, we sho w that the leading principal comp onen ts of kernel smoother sp ectral estimates resemble the MB tap ers. 2 Quadratic Estimators of the P o w er Sp ec- trum Let N discrete measuremen ts, x 1 , x 2 , . . . , x N , b e given as a realization of a stationary sto c hastic pro cess. W e normalize the time interv al b etw een measuremen ts to unit y . The Cramer representation of a discrete stationary sto c hastic process [4, 12] is x ( t ) = Z 1 / 2 − 1 / 2 e 2 π inf d Z ( f ) , where d Z has indep enden t sp ectral increments: E [ d Z ( f ) d Z ( g )] = S ( f ) δ ( f − g ) d f dg . W e assume that the sp ectral densit y , S ( f ), is t wice con tinuously differen tiable. The sp ectral in verse problem is to estimate the sp ectral densit y , S ( f ), giv en { x n } . As shown in [3, 8], every quadratic, mo dulation-inv ariant p o wer sp ectrum estimator has the form: b S ( f ) = N X n,m =1 q nm e 2 π i ( m − n ) f x n x m , (2) where Q = [ q nm ] is a symmetric matrix of order N and do es not dep end on frequency . Consider the eigenv ector decomp osition: Q = P K k =1 µ k v ( k )  v ( k )  T , 3 where K is the rank of Q , and v (1) , v (2) , . . . , v ( K ) is an orthogonal system of eigenv ectors. The multitaper represen tation of the quadratic sp ectral es- timator is b S ( f ) = K X k =1 µ k      N X n =1 v ( k ) n x n e − 2 π inf      2 . (3) In the case K = 1, estimator (3) turns out a tap er e d p erio do gr am estimator : b S v ( f ) =      N X n =1 v n x n e − 2 π inf      2 , (4) with a tap er v = ( v 1 , v 2 , . . . , v N ) T . If the tap ering is uniform (i.e. v 1 = v 2 = . . . = v N = 1 √ N ), w e name (4) the p erio do gr am estimator . The estimator (3) is a linear com bination of K orthogonal tap ered p erio dogram estimators. In MTSA, K is normally chosen to b e muc h less than N . The multiple tap er sp ectral estimate can b e thought of as a low rank, “principal comp onen ts” appro ximation of a general quadratic estimator. Multiple tap er analysis has also been applied to nonstationary sp ectral analysis [13, 1]. In practice, one does not b egin the analysis with a given quadratic es- timator, Q . Instead, one usually sp e cifies a family of orthonormal tap ers { v (1) , . . . , v ( K ) } with desir able pr op erties. Previously , only the family of Slepian tap ers were used in practice. The goal of this article is to intro- duce other families of tap ers. W e define the k th sp ectral windo w, V ( k ) , to b e the FT of the k th tap er: V ( k ) ( f ) = N X n =1 v ( k ) n e − i 2 π nf . (5) The tap ers are normally c hosen to hav e their sp ectral densit y localized near zero frequency . W e define t wo common measures of frequency lo calization. The lo c al bias of a sp ectral windo w V is R 1 / 2 − 1 / 2 f 2 | V ( f ) | 2 d f . The term “lo cal bias” is used b ecause it is prop ortional to the leading order term in the bias error of a taper estimate as N → ∞ . The sp e ctr al c onc entr ation in band [ − w , w ] is defined as R w − w | V ( f ) | 2 d f . The bandwidth, w , is a free parameter. The Slepian tapers are the unique se- quences which maximize the sp ectral concentration sub ject to the constrain t that they form an orthonormal family . Detailed analysis of the Slepian se- quences is given in [17]. W e stress that the Slepian tap ers depend on the 4 bandwidth parameter, w , and that the first 2 N w spectral windows are con- cen trated in the band [ − w , w ] while the remaining windows are concen trated outside. 3 Minim um Bias T ap ers 3.1 Con tin uous Time Case W e consider time-limited signals; the time in terv al is normalized to [0 , 1]. In the time domain, the tap er ν ( t ) is a function in L 2 [0 , 1] whic h we normalize to R 1 0 ν 2 ( t ) dt = 1. The functions { sin( π k t ) , k = 1 , 2 , . . . } form a complete orthogonal basis on [0 , 1]. (Completeness can b e prov en b y extending ν ( t ) to b e an o dd function on [ − 1 , 1] and using the completeness of the complex exp onen tials on [ − 1 , 1]. See [5].) Setting a k = 2 R 1 0 ν ( t ) sin( π k t ) dt , then P K k =1 a k sin( π k t ) dt con verges to ν ( t ) in L 2 [0 , 1] as K → ∞ . The tap er normalization is equiv alen t to 1 2 P ∞ k =1 a 2 k = 1. The F ourier transform of the tap er is the complex-v alued, sp ectral windo w function: V ( f ) = Z 1 0 ν ( t ) e − i 2 π f t dt . (6) V ( f ) is defined on the frequency domain [ −∞ , ∞ ], b elongs to L 2 [ −∞ , ∞ ], and satisfies R ∞ −∞ | V ( f ) | 2 d f = 2 π R 1 0 ν 2 ( t ) dt = 2 π . The local bias of a tap er spectral estimate [4, 9, 10] is E [ b S ( f )] − S ( f ) = Z ∞ −∞ | V ( g − f ) | 2 ( S ( g ) − S ( f )) dg ≈ S 00 ( f ) 2 Z ∞ −∞ | V ( h ) | 2 h 2 dh . (7) W e consider tap ers which minimize the leading order term: Z ∞ −∞ | V ( f ) | 2 f 2 d f = Z 1 0 1 2 π d dt ∞ X k =1 a k sin( π k t ) ! 2 dt = 1 4 Z 1 0 ∞ X k =1 a k k sin( π k t ) ! 2 dt = 1 8 ∞ X k =1 a 2 k k 2 . (8) The last expression attains the global minimum when a 2 1 = 2 , a 2 = a 3 = . . . = 0. Hence, the leading order term in the bias expression (7) is minimal for the tap er √ 2 sin( π t ) (This result was obtained b y P ap oulis [9]). In [9], P ap oulis extends ν ( t ) to b e zero outside of (0 , 1), and therefore has a F ourier 5 in tegral represen tation of ν ( t ). W e ha v e extended ν ( t ) to b e p erio dic and v anish at each in teger v alue. Since ν ( t ) is optimized for t ∈ [0 , 1], both represen tations are v alid. Equation (8) implies the more general result: Theorem 3.1 v k ( t ) = √ 2 sin( π k t ) ( k = 1 , 2 , . . . ) is the only system of func- tions in L 2 [0 , 1] which satisfy the r e quir ements: (i) R 1 0 v 2 k ( t ) dt = 1 , and (ii) v k minimizes R ∞ −∞ | V ( k ) ( f ) | 2 f 2 d f in the subsp ac e of func- tions ortho gonal to v 1 , . . . , v k − 1 . The kth minimum value is R ∞ −∞ | V ( k ) ( f ) | 2 f 2 d f = k 2 4 . W e name v k ( t ) = √ 2 sin( π k t ) ( k = 1 , 2 , . . . ) the c ontinuous time minimum bias tap ers . The F ourier transform of v k ( t ) is V ( k ) ( f ) = e − iπ ( f − k 2 ) i √ 2    sin h π  f − k 2 i π  f − k 2  − ( − 1) k sin h π  f + k 2 i π  f + k 2     = e − iπ ( f − k − 1 2 ) · k 2 sin  π f − π k 2  4 √ 2 π  f 2 −  k 2  2  . Th us, | V ( k ) ( f ) | deca ys as f − 2 for large frequencies. 3.2 Discrete Time Case W e no w consider the discrete time domain { 1 , 2 , . . . , N } with the corre- sp onding normalized frequency domain [ − 1 2 , 1 2 ]. A taper is a vector, ν = ( ν 1 , . . . , ν N ), normalized b y P N n =1 ( ν n ) 2 = 1. By the same argument as in the previous section, the leading order term of the bias of a tap er sp ectral esti- mate is prop ortional to the local bias, R 1 2 − 1 2 | V ( f ) | 2 f 2 d f , where the frequency windo w, V ( f ), is defined in Eq. (5). Lemma 3.2 F or the fr e quency window of a discr ete time tap er, Z 1 2 − 1 2 | V ( f ) | 2 f 2 d f = ν A ν ∗ , 6 wher e A = [ a nm ] with a nm = Z 1 2 − 1 2 e i 2 π ( n − m ) f f 2 d f = ( 1 12 if n = m ; ( − 1) n − m 2 π 2 ( n − m ) 2 if n 6 = m . Corollary 3.3 The tap ers ν (1) , ν (2) , . . . , ν ( N ) , define d by the r e quir ements (i) P N n =1 ( ν ( k ) n ) 2 = 1 , (ii) ν ( k ) minimizes R 1 2 − 1 2 | V ( k ) ( f ) | 2 f 2 d f in the subsp ac e of ve ctors ortho gonal to ν (1) , . . . , ν ( k − 1) , ar e the eigenve ctors of the matrix A sorte d in the incr e asing or der of the eigenvalues. The inte gr al R 1 2 − 1 2 | V ( k ) ( f ) | 2 f 2 d f is e qual to the kth eigenvalue. W e name ν (1) , ν (2) , . . . , ν ( N ) the discr ete minimum bias tap ers . They can b e approximated b y the sinusoidal tap ers , v (1) , v (2) , . . . , v ( N ) , which are discrete analogs of the contin uous time minimum bias tapers. Namely , w e define v ( k ) = ( v ( k ) 1 , . . . , v ( k ) N ) T with v ( k ) n = q 2 N +1 sin π kn N +1 , k = 1 , 2 , . . . , N . Lemma 3.4 The sinusoidal tap ers, v (1) , v (2) , . . . , v ( N ) , form an orthonormal b asis in R N and the lo c al bias of v ( k ) is k 2 4 N 2  1 + O ( 1 N )  . Corollary 3.5 The multitap er estimate (3) using K sinusoidal tap ers has lo c al bias e qual to P K k =1 µ k k 2 4 N 2 + O ( K 2 N 3 ) and has the fol lowing r epr esentation: ˆ S ( f ) = K X j =1 µ j 2( N + 1) | y ( f + j 2 N + 2 ) − y ( f − j 2 N + 2 ) | 2 . (9) The uniformly weighte d estimate, µ k = 1 K , has lo c al bias e qual to K 2 12 N 2 + O ( K 2 N 3 ) . F rom (9), the reason for the low bias of the sinusoidal tap ers is apparent: the fr e quency sidelob e fr om y ( f + j 2 N +2 ) c anc els the sidelob e of y ( f − j 2 N +2 ) . As a result, the sidelob e of y ( f + j 2 N +2 ) minus y ( f − j 2 N +2 ) is muc h smaller than that of the p erio dogram. Our preferred w eighting is the parab olic weigh ting: µ j = C (1 − j 2 /K 2 ) b ecause the parab olic weigh ting minimizes the exp ected square error in kernel smo others as K and N tend to infinity . Since the w eights decrease smo othly to zero, the resulting estimate is smo oth in frequency . 7 Corollary 3.6 The uniformly weighte d multitap er estimate using K sinu- soidal tap ers c an b e c ompute d in O ( N ln N ) + O ( K N ) op er ations while the generic multitap er estimate r e quir es O ( K N ln N ) op er ations plus the c ost of c omputing the K tap ers. The follo wing result demonstrates that the k th sp ectral window is con- cen trated on h k − 1 2( N +1) , k +1 2( N +1) i ∪ h − k +1 2( N +1) , − k − 1 2( N +1) i . Corollary 3.7 The F ourier tr ansform of v ( k ) e quals V ( k ) ( f ) = e − iπ ( ( N +1) f − k 2 ) i q 2( N + 1)    sin h N π  f − k 2( N +1) i sin h π  f − k 2( N +1) i − ( − 1) k sin h N π  f + k 2( N +1) i sin h π  f + k 2( N +1) i    = e − iπ ( ( N +1) f − k − 1 2 ) · sin π k N +1 q 2( N + 1) · sin h ( N + 1) π f − π k 2 i sin 2 ( π f ) − sin 2 π k 2( N +1) . Th us | V ( k ) ( f ) | = q N +1 2 for | f | = k 2( N +1) , and | V ( k ) ( f ) | ∼ 1 N 3 / 2 for f = O (1). In particular, | V ( k ) ( f ) | ≈ π k √ 2 N 3 / 2 when f → 1 2 . In the in termediate frequencies, k 2( N +1) < f  1 2 , | V ( k ) ( f ) | deca ys as 1 f 2 . Numerical ev aluation (see T able 1) shows that k v ( k ) − ν ( k ) k L 2 < k 4( N +2) for all k . The same rate of conv ergence is observ ed in L ∞ norm:     v ( k ) k v ( k ) k L ∞ − ν ( k ) k ν ( k ) k L ∞     L ∞ < k 2( N +2) . ( k · k L ∞ is the suprem um norm in the time domain.) Figure 1 plots the en velopes of | V ( k ) ( f ) | 2 for b oth the minimum bias and sinusoidal tapers with N = 200 , k = 1. The sp ectral energy of b oth tap ers is nearly iden tical for | f | < . 25. Near the Nyquist frequency , the sp ectral energy of the sin usoidal tap er is roughly three times larger than that of the minim um bias taper. In the time domain, the sin usoidal tap ers are virtually indistinguishable from the MB tap ers. 4 Comparison of Sp ectral Lo calizations W e now compare the lo cal bias and the sp ectral concen tration of three fam- ilies of orthonormal tap ers: minim um bias (MB) tap ers, sin usoidal tap ers and Slepian tap ers. Both the lo cal bias and the sp ectral concentration of the 8 Slepian tap ers dep end on the bandwidth parameter, w . F or prop erties of the Slepian tapers, we refer the reader to [12, 16, 17, 18, 19]. Since the MB tap ers minimize the lo cal bias, clearly the sinusoidal tap ers and the Slepian tap ers hav e larger lo cal bias. The only question is whether the difference is large or small. T able 2 giv es the lo cal bias, P K k =1 R 1 / 2 − 1 / 2 f 2 | V ( k ) ( f ) | 2 d f , of the three families of tap ers for N = 50. The sinusoidal tap ers c ome within 0.2% of achieving the optimal lo c al bias. In con trast, the lo cal bias of the Slepian tapers can b e man y times larger. W e compute the lo cal bias for three different v alues of the bandwidth, w . The general pattern is that the k th Slepian tap er has roughly the same lo cal bias as the MB tap er do es when N w < k < 2 N w . The ratio of the lo cal bias of the Slepian tapers to that of the MB tap ers is smallest at k ≈ 1 . 2 N w . As | k − 1 . 2 N w | increases, the lo cal bias rapidly departs from the optimal v alue. T able 3 compares the spectral concentration of the tapers for N = 50 and w = . 08 . Both the MB tap ers and the sin usoidal tap ers are within 1.7% of the optimal v alue, except for k = 8 , 9. Notice that 2 N w = 8. Although the ratio of the sp ectral concentration, R w − w | V ( f ) | 2 d f , for the MB and sinusoidal tap ers to that of the Slepian tap ers is usually v ery close to one, the ratio of the sp ectral energy outside of the frequency band | f | < w can b e quite large. Th us our conclusions dep end on using R w − w | V ( f ) | 2 d f and not 1 − R w − w | V ( f ) | 2 d f as the measure of frequency concen tration. Figures 1-3 compare | V ( k ) ( f ) | 2 of the Slepian and MB tap ers for N = 200. F or Figures 1 and 2, we select the Slepian parameter, w , equal to .01 so that K = 2 N w equals four. Figure 2 plots | V 1 ( f ) | 2 for the frequencies up to f = 0 . 14. The central p eak of the MB tap er is more concentrated around f = 0 than the Slepian tap er is. The first sidelob e of the MB tap er is visible while the first Slepian sidelob e is muc h smaller. Figure 1 plots the logarithm of | V 1 ( f ) | 2 o ver the en tire frequency range. The MB tap er has smaller range bias in the frequency range | f | < 0 . 3 w and in the frequency range | f | > 0 . 13. In the middle frequency range, the Slepian tap er is clearly better. The Slepian penalty function maximizes the energy inside the frequency band, [ − w , w ], and thus it is natural that the Slepian tap ers do b etter for f ∼ w . By using a discontin uous p enalt y function, the Slepian sp ectral windows exp erience Gibbs phenomenon and decay only as 1 f , ( | V ( f ) | 2 ∼ 1 f 2 ). The MB sp ectral windows deca y as 1 f 2 , and th us, it is natural that the MB tap ers ha ve low er bias for f ∼ O (1). Figure 3 plots P 3 k =1 | V ( k ) ( f ) | 2 for | f | < 0 . 14 and on this scale, the MB 9 tap ers are clearly preferable to the Slepian tapers. F or larger frequencies, the energy of the multitaper estimate, P K k =1 | V ( k ) ( f ) | 2 , is v ery similar to Fig. 2 on the logarithmic scale pro vided that K << N . In summary , the sinusoidal tap ers p erform nearly as well as the MB tap ers while the Slepian tap ers ha ve several times larger lo cal bias (except when k ≈ 1 . 2 N w ). F or k  2 N w , the Slepian tap ers ha ve b etter broad-band bias protection than the minim um bias tap ers do. F or k ∼ 2 N w , the minim um bias tap ers pro vide b oth smaller lo cal bias and b etter broad-band protection due to the Gibbs phenomena whic h the Slepian tap ers exp erience. 5 Lo cal Error Analysis and Optimal Multita- p ering W e no w give a lo cal error analysis of MTSA and determine the optimal n umber of tap ers. Our results are the m ultitap er analog of the lo cal error analysis of the smo othed p erio dogram [4, 10]. W e assume the time series is a Gaussian processes and do not consider frequencies near f = 0 and f = 1 / 2. In this case, the v ariance of the multitaper estimate is appro ximately V ariance[ ˆ S ( f )] ≈ S ( f ) 2 P K k =1 µ 2 k due to the orthonormality of the tap ers Asymptotically , the lo cal bias of the multitaper estimate of Eq. (3) is Bias[ ˆ S ] = S ( f ) K X k =1 µ k − 1 ! + 1 2 S 00 ( f ) K X k =1 λ k µ k , where λ k = R 1 / 2 − 1 / 2 f 2 | V ( k ) ( f ) | 2 d f . The second term is the MT generalization of (7). When P K k =1 µ k 6 = 1, the MT estimate has bias ev en in white noise. When w e require P K k =1 µ k = 1, the lo cal exp ected loss simplifies: Theorem 5.1 F or a Gaussian pr o c ess, away fr om f = 0 and f = 1 / 2 , the exp e cte d squar e err or of the multitap er sp e ctr al estimate (3) with P K k =1 µ k = 1 is asymptotic al ly (to le ading or der in K / N ) Bias 2 + V ariance ≈ " 1 2 S 00 ( f ) K X k =1 λ k µ k # 2 + S ( f ) 2 K X k =1 µ 2 k . (10) Theorem 5.2 The multitap er estimate which minimizes the lo c al loss (10) (with µ k ≥ 0 ) is c onstructe d with the minimum bias tap ers. 10 Pro of: W e order the µ k suc h that µ 1 ≥ µ 2 ≥ . . . ≥ µ K and define µ K +1 = 0. Since the weigh ts, µ k are fixed, w e need to minimize P K k =1 µ k u k A u ∗ k o ver all sets of K orthonormal tap ers, u 1 , . . . , u K . W e split the series in K subseries and minimize each subseries separately: min u 1 ,..., u K K X k =1 µ k u k A u ∗ k = min u 1 ,..., u K K X k =1 ( µ k − µ k − 1 )   k X j =1 u j A u ∗ j   ≥ K X k =1 ( µ k − µ k − 1 )   min u ( k ) 1 ,..., u ( k ) k k X j =1 u ( k ) j A u ( k ) ∗ j   = K X k =1 ( µ k − µ k − 1 )   k X j =1 λ A,j   = K X k =1 µ k λ A,j , (11) where the λ A,j are the eigen v alues of A , given in increasing order. The u ( k ) j are sub ject to orthonormalit y constrain ts that u ( k ) j · u ( k ) j 0 = δ j,j 0 , but are otherwise indep endent and minimized separately . In the last line of (11), w e use F an’s Theorem [6]: min u ( k ) 1 ,..., u ( k ) k P k j =1 u ( k ) j A u ( k ) ∗ j = P k j =1 λ A,j , where the u ( k ) j are again sub ject to orthonormality constraints. The theorem is no w pro ved b ecause the MB tapers are precisely the eigenv ectors of A . Theorem 5.3 The uniformly weighte d multitap er estimate using K sinu- soidal tap ers has an asymptotic lo c al loss of Bias 2 + V ariance ' " S 00 ( f ) K 2 24 N 2 # 2 + S ( f ) 2 K . (12) Corollary 5.4 The asymptotic lo c al loss of (12) is minimize d when the num- b er of tap ers is chosen as K opt ∼ " 12 S ( f ) N 2 S 00 ( f ) # 2 / 5 . (13) Th us, the optimal num b er of tap ers is prop ortional to N 4 / 5 and v aries with the ratio of S ( f ) to S 00 ( f ). Intuitiv ely (13) sho ws that few er tap ers should b e used when the sp ectrum v aries more rapidly . A key adv an tage of the MB and sin usoidal tap ers is that the tap ers need not b e recomputed as K is c hanged. In contrast, the Slepian tap ers are most efficient when the bandwidth parameter, w , is chosen such that K ∼ 2 N w . Thus, when the n umber of tap ers is c hanged, as in (13), the Slepian tap ers should be recomputed. 11 6 Smo othed Multitap er Estimates In our o wn comparison of kernel smo othing and multitaper estimation [16], w e found that a smo othed m ultiple tap er estimate work ed b est. W e no w ev aluate the exp ected error of the k ernel smo othed multitaper estimator and sho w that smo othing the logarithm of the m ultitap er estimate is useful for estimating the l og ar ithm of the sp ectrum. W e b egin by ev aluating that the quadratic estimator (2) whic h is equiv- alen t to a kernel smoother estimates of the spectrum [4, 10]. Let ˆ S ( f ) be the quadratic sp ectral estimator (2), and smo oth it with a kernel κ ( · ) of halfwidth w : b b S ( f ) = Z w − w κ ( g w ) b S ( f + g ) dg , (14) where w is the bandwidth parameter and κ ( · ) is a k ernel smo other with domain [ − 1 , 1]. This can b e rewritten as b b S ( f ) = N X n,m =1 ˜ q nm e i 2 π ( m − n ) f x n x m , (15) where ˜ q nm = q nm ˆ κ m − n with ˆ κ m = R w − w κ ( g /w ) e 2 π img dg . Th us smo oth- ing replaces the original quadratic estimator with matrix [ q nm ] by another quadratic estimator with matrix ˜ Q = [ ˜ q nm ]. By Theorem 5.2, this h ybrid metho d cannot outp erform the pure m ultitap er metho d with minimum bias tap ers. W e now sho w that com bining k ernel smo othing with multitapering do es impro ve the estimation of the l og ar ithm of the sp ectral density , θ ( f ) = log[ S ( f )]. One standard approach is to k ernel smo oth the logarithm of the tap ered p erio dogram. This approach has the disadv an tage that | y ( f ) | 2 has a χ 2 2 distribution and log[ χ 2 2 ] has a long low er tail of its distribution. As a re- sult, log [ | y ( f ) | 2 ] has an appreciable bias and its v ariance is inflated by π 2 / 6. A common alternativ e is to estimate the spectrum either by kernel smo othing or b y multitapering and then to tak e logarithms. This approach has the dis- adv an tage that the smo othed sp ectral estimate tends to b e more sensitive to nonlo cal bias effects than the corresp onding smo othed log-sp ectral estimate. T o robustify the log-spectral estimate while reducing the v ariance infla- tion from the long tail, w e propose the follo wing h ybrid estimate: 1) compute the m ultitap er estimate using the sin usoidal tap ers with µ k = 1 K and then 2) smo oth ˆ θ M T ( f ) ≡ ln[ ˆ S M T ( f )] − B K /K , where B K is the bias of ln[ χ 2 2 K ]. 12 ( B K ≡ ψ ( K ) − ln K where ψ ( K ) is the digamma function). F or white noise, the v ariance of ˆ θ M T ( f ) = ψ 0 ( K ) ∼ 1 K + 1 2 K 2 , so the v ariance enhancemen t from the logarithm tends rapidly to zero. In [15], we show that the asymptotic error for this sc heme is θ 00 ( f ) 2 " b k w 2 + K 2 24 N 2 # 2 + C κ N w  1 + 1 2 K  2 , (16) where b κ and C κ are constan ts whic h dep end on the k ernel shap e. In (16), w e assume uniformly weigh ted sinusoidal tap ers are used and 1  K  N w . In (16), one factor of (1 + 1 2 K ) is the v ariance enhancemen t from the logarithmic transformation and one factor of (1 + 1 2 K ) arises in the v ariance calculation of (15) with sin usoidal tap ers. Optimizing (16) with resp ect to b oth w and K yield w ∼ N − 1 / 5 and K ∼ N 8 / 15 , thus the smo othing halfwidth w is m uch larger than K / N . The exp ected error (16) dep ends w eakly on K provided that 1  K  N w . F or simplicity , we set K = N 8 / 15 and optimize (16) with resp ect to the halfwidth w . The resulting halfwidth dep ends on θ 00 ( f ): w opt ( θ 00 ( f )) with w opt ∼ | θ 00 ( f ) | − 2 / 5 N − 1 / 5 . Th us when the log-sp ectrum v aries rapidly , the halfwidth should b e reduced as | θ 00 ( f ) | − 2 / 5 . Since θ 00 ( f ) is unknown, w e consider t wo stage estimators which b egin b y making a preliminary estimate of θ 00 ( f ) prior to estimating θ ( f ). W e then insert the estimate d θ 00 ( f ) in to the expression for w opt : w ( f ) = w opt ( c θ 00 ( f )) and use a v ariable halfwidth k ernel smo other with halfwidth ˆ w ( f ) to estimate θ ( f ). Multiple stage k ernel estimators are described in [2, 7, 13, 14, 15]. These m ultiple stage schemes hav e a conv ergence rate of N − 4 / 5 and ha v e a relativ e conv ergence rate of at least N − 2 / 9 . A more detailed description is giv en in [2, 13, 15]. 7 Application W e no w compare spectral estimates on an actual data series. W e use the micro wa ve scattering data set which is describ ed in [16]. The data measures turbulen t plasma fluctuations in the T ok amak F usion T est Reactor at Prince- ton. The sp ectrum is dominated by a 1 MHz p eak which is quasicoherent. The spectral density v aries b y o ver five orders of magnitude. 13 The bias v ersus v ariance trade-off of Sec. 5 sho ws that fewer tap ers should b e used near the p eak. T o make the sp ectral estimate smo oth, a parab olic w eighting of the tap ers is used as describ ed in Sec. 3. T o determine how man y tap ers to use lo cally , we use the multiple stage “plug-in” method as describ ed in the previous section; i.e. w e determine the n umber of tapers using a pre-estimate on the same data. T o reduce the fluctuations from the estimate of the optimal n umber, we use a longer data segment to determine the n umber of tap ers at each frequency . W e find the optimal n umber of sin usoidal tap ers is roughly 24 for frequencies in the 200 to 800 kHz range. Near the 1 MHz p eak, as few as 12 tap ers are used to minimize the lo cal bias error. Betw een 1300 and 2400 kHz, the sp ectrum is flatter and w e use up to 40 tapers. The dotted line is the sinusoidal m ultitap er estimate, and the solid, more wiggly , curv e is the corresp onding Slepian estimate using 24 tap ers with w = 60 kHz. The 1 MHz p eak is p o orly resolved in the Slepian estimate, and the regions of high curv ature are artificially flattened. F or f ≥ 1 . 5 MHz, the Slepian estimate is artificially bumpy due to statistical noise. The v ariable tap er n umber estimate suppresses these bumps by av eraging o ver a larger frequency halfwidth. W e ha ve also used a v ariable tap er num b er estimate with the Slepian tapers. Since the Slepian parameter, w w as fixed at 100 kHz to allo w for fort y tap ers, the artificial broadening w as ev en more exreme. Comparing with a con verged estimate of the sp ectrum based on N = 45 , 000 sho ws that the sinuosidal tap er estimate is more accurate. Another significan t difference is that the Slepian multitaper estimate requires muc h more CPU time than the sin usoidal m ultitap er estimate. 8 Conclusion W e ha ve proposed and analyzed the minim um bias and the sin usoidal tap ers, v ( k ) n = q 2 N +1 sin π kn N +1 , for multitaper sp ectral estimation. The resulting sin u- soidal multitaper sp ectral estimate is ˆ S ( f ) = 1 2 K ( N +1) P K j =1 | y ( f + j 2 N +2 ) − y ( f − j 2 N +2 ) | 2 . The sin usoidal tap ers hav e lo w bias b ecause the fr e quency sidelob e fr om y ( f + j 2 N +2 ) c anc els the sidelob e of y ( f − j 2 N +2 ) . The minim um bias tapers minimize the lo cal bias, R 1 / 2 − 1 / 2 f 2 | V ( k ) b ( f ) | 2 d f , and ha v e goo d broad-band bias protection as w ell. Asymptotically , the quadratic sp ectral estimate whic h minimizes the exp ected lo cal square er- 14 ror is a multiple tap er estimate using the minimal bias tapers. The sinusoidal tap ers hav e a simple analytic form and appro ximate the minim um bias tap ers to O  1 N  . The k th sinusoidal tap er has its spectral energy concen trated in the frequency bands k − 1 2( N +1) ≤ | f | ≤ k +1 2( N +1) . The minim um bias and sin usoidal tap ers ha ve no auxiliary bandwidth pa- rameter, and the bandwidth of the sp ectral estimate is determined solely b y the n um b er of used tap ers. By adaptively adding and deleting tap ers, a m ul- titap er estimate with the optimal con vergen ce prop erties of k ernel smo others can b e constructed. In contrast, the Slepian tap ers need to b e recomputed with a differen t bandwidth. Th us the Slepian tap ers are only practical for fixed bandwidth estimation and this is inherently inefficient. 9 App endix: Multitap er decomp osition of k er- nel estimates In Sec. 6, w e show ed that k ernel smo other estimators (14) hav e an equiv alen t m ultitap er representation (3) W e now sho w that the equiv alen t m ultitap ers of some p opular k ernel smoother estimates of the sp ectrum strongly resem ble the MB/sinusoidal tap ers. In one special case, this corresondence is exact; i.e. the smo othed p erio dogram can b e exactly decomp osed into MB tap ers. Theorem 9.1 L et w = 1 2 and κ ( f ) b e the p ar ab olic kernel, κ ( f ) = 3 2 − 6 f 2 . The eigenve ctors of the kernel smo othe d p erio do gr am ar e exactly the discr ete minimum bias tap ers. Pro of: The ˜ Q matrix in (15) can b e calculated explicitly for this case. W e find ˜ Q = [ b nm ] = 1 N ( 3 2 I − 6 A ) where A is the matrix from Lemma 3.2. Th us ˜ Q and A hav e the same eigenv ectors. T o illustrate that this result is typical ev en when w e apply a taper and smo oth ov er a small band, we consider a smo othed tap ered p erio dogram with N = 200. W e use T ukey’s split-cosine tap er [16] and then smo oth the estimate with a square b o x k ernel with a halfwidth of . 01. W e then ev aluate the corresp onding ˜ Q matrix and compute its eigenv ectors. Figure 5 displa ys the first 4 eigen v ectors. They are v ery close to the sinusoids q 2 N +1 sin π kn N +1 . T able 4 shows that this sp ectral estimate is virtually a K = 4 multiple tap er 15 sp ectral estimate. After k > 4, the eigenv alues decrease sharply , and these higher eigen vectors con tribute v ery little to the ov erall estimate. References [1] M. Amin, “Optimal estimation of evolutionary spectra,” I.E.E.E. T r ans. on Signal Pr o c essing v ol. 42, 2?, 1994. [2] T. Bro c kman, Th. Gasser and E. Hermann, “Lo cally adaptiv e bandwidth choice for kernel regression estimators,” J. A mer. Stat. Asso c. 88 , 1302-1309 (1994). [3] T. P . Bronez, Nonp ar ametric Sp e ctr al Estimation of Irr e gularly Sample d Multidimensional R andom Pr o c esses , PhD Thesis, Ari- zona State Universit y , 1985. [4] U. Grenander and M. Rosen blatt, Statistic al Analysis of Station- ary Time Series , New Y ork: Wiley , 1957. [5] A. N. Kolmogoro v and S. V. F omin, R e ele F unktionen und F unk- tionalanalysis , Section 7.3.2, Berlin: VEB Deutsc her V erlag der Wissensc haften, 1975. [6] A. W. Marshall and I. Olkin, Ine qualities: The ory of Majorization and its Applic ations , p. 511, New Y ork: Academic Press 1979. [7] H.-G. M ¨ uller and U. Stadtm ¨ uller, “V ariable bandwidth kernel estimators of regression curv es,” A nnals of Statistics , vol. 15, pp. 182-201, 1987. [8] C. T. Mullis and L. L. Scharf, in A dvanc es in Sp e ctrum A nalysis , S. Ha ykin, Ed., New Y ork: Pren tice-Hall, 1990, Chapter 1, pp. 1-57. [9] A. P ap oulis, “Minimum bias windo ws for high resolution sp ectral estimates,” IEEE T r ans. Information The ory , vol. 19, pp. 9-12, 1973. [10] E. P arzens, “On asymptotically efficien t consisten t estimates of the sp ectral densit y of a stationary time series,” J. R oyal Stat. So c. , v ol. 19, pp. 303-322, 1958. 16 [11] J. Park, C. R. Lindberg, and F. L. V ernon, “Multitap er Sp ectral analysis of high frequency seismograms,” J. Ge ophys. R es. , v ol. 92B, pp. 12765-12684, 1987. [12] D. Perciv al and A. W alden, Sp e ctr al A nalysis for Physic al Ap- plic ations: Multitap er and Conventional Univariate T e chniques . Cam bridge: Cam bridge Universit y Press, 1993. [13] K. S. Riedel, “Optimal k ernel estimation of ev olutionary sp ectra,” I.E.E.E. T r ans. on Signal Pr o c essing v ol. 41, 2439-2447, 1993. [14] K.S. Riedel and A. Sidorenko, “F unction estimation using data adaptiv e kernel smo others- How muc h smo othing?” Computers in Physics vol. 8, 402-409, 1994. [15] K. S. Riedel and A. Sidorenk o, “Smo othed log-m ultitap er sp ec- tral estimation and data adaptive implementation,” submitted for publication. [16] K. S. Riedel, A. Sidorenk o, and D. J. Thomson, “Sp ectral densit y estimation for plasma fluctuations I: Comparison of metho ds,” Physics of Plasmas v ol. 1, pp. 485-500. [17] D. Slepian, “Prolate spheroidal w av e functions, F ourier analysis, and uncertain ty - V: the discrete case,” Bel l System T e ch. J. , vol. 5, pp. 1371-1429, 1978. [18] D. J. Thomson, “Spectrum estimation and harmonic analysis,” Pr o c. IEEE , v ol. 70, pp. 1055-1096, 1982. [19] D. J. Thomson, “Quadratic inv erse sp ectrum estimates: appli- cations to paleo climatology ,” Phil. T r ans. R. So c. L ond. A , v ol. 332, p. 539-597, 1990. 17 T able Captions: T able 1: Con v ergence of the sin usoidal tap ers to the minimum bias tap ers. T able 2: Normalized bias term, 4( N + 1) 2 P K k =1 R 1 / 2 − 1 / 2 f 2 P K k =1 f 2 | V ( k ) ( f ) | 2 d f , for N = 50. T able 3: Sp ectral concen tration, R w − w P K k =1 | V ( k ) ( f ) | 2 d f , for N = 50. T able 4: Eigen v ectors of the smooth tap ered p erio dogram estimator. Figure Captions: Figure 1: Sp ectral energy of the minimum bias, sin usoidal and Slepian tap ers, N = 200 , k = 1 , w = . 01. Figure 2: Sp ectral energy of the minim um bias and sinusoidal tap ers, N = 200 , k = 1. Figure 3 Spectral energy , P 3 k =1 R 1 / 2 − 1 / 2 f 2 | V ( k ) ( f ) | 2 d f , of the minim um bias and Slepian tapers, N = 200 , K = 3 , w = . 01. Figure 4: Estimated sp ectral densit y of the plasma fluctuations. Dashed line is sinusoidal multitaper estimate and solid line is estimate using Slepian tap ers with w = 60 kHz. Because the Slepian tap ers hav e a fixed bandwidth, the corresp onding estimate sp ectral densit y at 1 MHz is artificially broadened while being undersmo othed for f ≥ 1 . 5 MHz. Figure 6: First eigen vectors of the smo oth tap ered p erio dogram estima- tor. 18 19 20 21 22 23 24 T able 1. Con v ergence of the sin usoidal tap ers to the minimum bias tap ers N max k n N +2 k k v ( k ) − ν ( k ) k L 2 o max k ( N +2 k     v ( k ) k v ( k ) k L ∞ − ν ( k ) k ν ( k ) k L ∞     L ∞ ) 20 0.24750 0.4602 50 0.24844 0.4760 200 0.24852 0.4829 800 0.24844 0.4844 T able 2. Normalized bias term, 4( N + 1) 2 P K k =1 R 1 / 2 − 1 / 2 f 2 | V ( k ) ( f ) | 2 d f , for N = 50 K Minim um Sin usoidal Slepian tapers bias tapers tap ers w =0.04 w =0.08 w =0.16 1 1.0095 1.0116 1.3439 2.6316 5.1039 2 5.0475 5.0580 5.7724 10.5484 20.4670 3 14.1328 14.1622 18.0651 23.8086 46.1953 4 30.2846 30.3475 58.9520 42.5181 82.4018 5 55.5217 55.6366 154.4818 66.9996 129.2087 6 91.8634 92.0528 305.4382 99.1800 186.7507 7 141.3284 141.6185 496.7959 150.0103 255.1797 8 205.9362 206.3570 721.1743 251.8833 334.6717 9 287.7056 288.2899 976.5088 437.5993 425.4379 10 388.6562 389.4409 1262.4251 702.1523 527.7433 25 T able 3. Sp ectral concen tration, R w − w | V ( f ) | 2 d f , for N = 50 , w = 0 . 08 k Minim um bias tap ers Sin usoidal tap ers Slepian tapers 1 .9997 .9997 1. 2 .9988 .9988 .9999999 3 .9972 .9972 .9999989 4 .9940 .9937 .99997 5 .9888 .9887 .9995 6 .9760 .9753 .9928 7 .9381 .9417 .9380 8 .6084 .6247 .7002 9 .1688 .1780 .2981 10 .0637 .0624 .0628 T able 4. Eigen v ectors of the smooth tap ered p erio dogram estimator k W eigh t of the eigen vector Normalized lo cal bias Lo cal bias in comparison with λ k ( B ) / tr( B ) 4( N + 1) 2 R 1 / 2 − 1 / 2 f 2 | V ( f ) | 2 d f the minimum bias tap er (ratio) 1 .2856 1.5138 1.509 2 .2828 4.7371 1.181 3 .2519 9.6254 1.067 4 .1416 19.2095 1.198 5 .0340 33.7118 1.345 6 .0037 51.3616 1.423 7 .0002 72.9747 1.486 26