arXiv:1803.04078v2 [stat.ME] 30 Mar 2018

Minimum bias multiple taper spectral
estimation ∗
Kurt S. Riedel and Alexander Sidorenko
Courant Institute of Mathematical Sciences, New York University
New York, New York 10012-1185
EDICS: SP 3.1.1

Abstract
Two families of orthonormal tapers are proposed for multitaper
(k)
spectral analysis:
q minimum bias tapers, and sinusoidal tapers {v },
(k)
where vn = N 2+1 sin Nπkn
+1 , and N is the number of points. The resulting sinusoidal multitaper spectral estimate is Ŝ(f ) =
j
2N +2 ) − y(f

1
2K(N +1)

PK

j=1 |y(f +

j
2
2N +2 )| ,

−
where y(f ) is the Fourier transform of the stationary time series, S(f ) is the spectral density, and K is the number
of tapers. For fixed j, the sinusoidal tapers converge to the minimum
bias tapers like 1/N . Since the sinusoidal tapers have analytic expressions, no numerical eigenvalue decomposition is necessary. Both
the minimum bias and sinusoidal tapers have no additional parameter for the spectral bandwidth. The bandwidth of the jth taper is
simply N1 centered about the frequencies 2N±j+2 . Thus the bandwidth
of the multitaper spectral estimate can be adjusted locally by simply
adding
or deleting tapers. The band limited spectral concentration,
Rw
|V
(f
)|2 df , of both the minimum bias and sinusoidal tapers is very
−w
close to the optimal concentration achieved by the Slepian tapers. In
R 1/2
contrast, the Slepian tapers can have the local bias, −1/2 f 2 |V (f )|2 df ,
much larger than of the minimum bias tapers and the sinusoidal tapers.
∗

The authors thank D. J. Thomson and the referees for useful comments. Research
funded by the U.S. Department of Energy.

1

1

Introduction

We consider a stationary time series, {xn , n = 1 . . . N } with a spectral density, S(f ). A common estimator of the spectral density is to smooth the
square of the discrete Fourier transform (DFT) locally:
L
X
1
j
Ŝ(f ) =
|y(f + )|2 ,
(2L + 1)N j=−L
N

(1)

where y(f ) is the Fourier transform (FT) of the stationary time series:
PN
−i2πnf
y(f ) ≡
. Since (1) is quadratic in the FT, y(f ), it is
n=1 xn e
natural to consider a more general class of quadratic spectral estimators. We
examine quadratic estimators where the underlying self-adjoint matrix has
rank K, where K is prescribed. Using the eigenvector representation, the
resulting quadratic spectral estimator can be recast as a weighted sum of K
orthonormal rank one spectral estimators. This class of spectral estimators
was originally proposed by Thomson [19] under the name of multiple taper
spectral analysis (MTSA). We refer the reader to [3, 8, 11, 12, 16, 18, 19] for
excellent expositions and generalizations of Thomson’s theory.
In MTSA, a rank K quadratic spectral estimate is constructed by choosing an orthonormal family of tapers/spectral windows and then averaging
the K estimates of the spectral density. In practice, only the Slepian tapers
(also known as discrete prolate spheroidal sequences [17]) are routinely used
for MTSA.
In the present paper, we propose and analyze two new orthonormal families of tapers: minimum bias (MB) tapers
and sinusoidal tapers. The MB
R 2
tapers minimize the local frequency bias, f |V (f )|2 df , subject to orthonormality constraints, where V (f ) is the DFT of the taper. For continuous time,
the MB tapers have simple analytic expressions. The first taper in the family is Papoulis’ optimal taper [9]. For discrete time, the MB tapers satisfy a
selfadjoint eigenvalue problem and may be computed numerically.
(k)
In the
q case of discrete time, we define the kth sinusoidal taper, v , as
vn(k) = N2+1 sin Nπkn
, where N is the sequence length. The sinusoidal tapers
+1
are an orthonormal family that converge to the MB tapers with rate 1/N as
N → R∞. These results are given in Section 3. Section 4R compares the local
1/2
w
bias, −1/2 f 2 |V (f )|2 df , and the spectral concentration, −w
|V (f )|2 df , of the
MB tapers, the sinusoidal tapers and the Slepian tapers.

2

In Section 5, we show that the quadratic spectral estimator which minimizes the expected square local error is weighted multitaper estimate using
the MB tapers. A local error analysis is given and the optimal number of
tapers is determined. At frequencies where the spectral density is changing rapidly, fewer tapers should be used. In Section 6, we show that kernel
smoother spectral estimates [4, 10] are multitaper estimates and we show
that smoothing the logarithm of the multitaper estimate significantly reduces the variance in comparison with smoothing athe logarithm of a single
taper estimate. We also describe our data adaptive method for estimating
the spectrum. In Section 7, we apply our spectral estimation techniques to
real data and show that our tapers outperform the Slepian tapers whenever
a variable bandwidth is needed. In the Appendix, we show that the leading
principal components of kernel smoother spectral estimates resemble the MB
tapers.

2

Quadratic Estimators of the Power Spectrum

Let N discrete measurements, x1 , x2 , . . . , xN , be given as a realization of
a stationary stochastic process. We normalize the time interval between
measurements to unity. The Cramer representation of a discrete stationary
stochastic process [4, 12] is
x(t) =

Z 1/2

e2πinf dZ(f ) ,

−1/2

where dZ has independent spectral increments: E[dZ(f )dZ(g)] = S(f )δ(f −
g)df dg. We assume that the spectral density, S(f ), is twice continuously
differentiable.
The spectral inverse problem is to estimate the spectral density, S(f ),
given {xn }. As shown in [3, 8], every quadratic, modulation-invariant power
spectrum estimator has the form:
b )
S(f

=

N
X

qnm e2πi(m−n)f xn xm ,

(2)

n,m=1

where Q = [qnm ] is a symmetric matrix of order N and does not depend on
frequency. Consider the eigenvector decomposition: Q =
3

PK

k=1



µk v(k) v(k)

T

,

where K is the rank of Q , and v(1) , v(2) , . . . , v(K) is an orthogonal system
of eigenvectors. The multitaper representation of the quadratic spectral estimator is
b )
S(f

=

K
X

µk

N
X

2

vn(k) xn e−2πinf

.

(3)

n=1

k=1

In the case K = 1, estimator (3) turns out a tapered periodogram estimator:
Sbv (f ) =

N
X

2

vn xn e

−2πinf

,

(4)

n=1

with a taper v = (v1 , v2 , . . . , vN )T . If the tapering is uniform (i.e. v1 = v2 =
. . . = vN = √1N ), we name (4) the periodogram estimator. The estimator (3)
is a linear combination of K orthogonal tapered periodogram estimators. In
MTSA, K is normally chosen to be much less than N . The multiple taper
spectral estimate can be thought of as a low rank, “principal components”
approximation of a general quadratic estimator. Multiple taper analysis has
also been applied to nonstationary spectral analysis [13, 1].
In practice, one does not begin the analysis with a given quadratic estimator, Q . Instead, one usually specifies a family of orthonormal tapers
{v(1) , . . . , v(K) } with desirable properties. Previously, only the family of
Slepian tapers were used in practice. The goal of this article is to introduce other families of tapers.
We define the kth spectral window, V (k) , to be the FT of the kth taper:
V (k) (f ) =

N
X

vn(k) e−i2πnf .

(5)

n=1

The tapers are normally chosen to have their spectral density localized near
zero frequency. We define two common measures
of frequency localization.
R 1/2 2
The local bias of a spectral window V is −1/2 f |V (f )|2 df . The term
“local bias” is used because it is proportional to the leading order term in
the bias error of a taper estimate as N → ∞.
Rw
The spectral concentration in band [−w, w] is defined as −w
|V (f )|2 df .
The bandwidth, w, is a free parameter. The Slepian tapers are the unique sequences which maximize the spectral concentration subject to the constraint
that they form an orthonormal family. Detailed analysis of the Slepian sequences is given in [17]. We stress that the Slepian tapers depend on the
4

bandwidth parameter, w, and that the first 2N w spectral windows are concentrated in the band [−w, w] while the remaining windows are concentrated
outside.

3
3.1

Minimum Bias Tapers
Continuous Time Case

We consider time-limited signals; the time interval is normalized to [0, 1]. In
theR time domain, the taper ν(t) is a function in L2 [0, 1] which we normalize
to 01 ν 2 (t)dt = 1. The functions {sin(πkt), k = 1, 2, . . .} form a complete
orthogonal basis on [0, 1]. (Completeness can be proven by extending ν(t)
to be an odd function on [−1, 1] and using the completeness of the complex
exponentials on [−1,R 1]. See [5].)
P
to
Setting ak = 2 01 ν(t) sin(πkt)dt, then K
k=1 ak sin(πkt)dt converges
P
2
a
ν(t) in L2 [0, 1] as K → ∞. The taper normalization is equivalent to 21 ∞
k=1 k =
1. The Fourier transform of the taper is the complex-valued, spectral window
function:
Z 1
ν(t)e−i2πf t dt .
(6)
V (f ) =
0

V (f ) is defined
on the frequencyR domain [−∞, ∞], belongs to L2 [−∞, ∞],
R∞
and satisfies −∞ |V (f )|2 df = 2π 01 ν 2 (t)dt = 2π.
The local bias of a taper spectral estimate [4, 9, 10] is
S 00 (f ) Z ∞
|V (h)|2 h2 dh .
2
−∞
−∞
(7)
We consider tapers which minimize the leading order term:

b )]−S(f ) =
E[S(f

Z ∞
−∞

Z ∞

2 2

|V (g−f )|2 (S(g)−S(f ))dg ≈

|V (f )| f df =

Z 1
0

!2

∞
1 d X
ak sin(πkt)
2π dt k=1

∞
1Z 1 X
=
ak k sin(πkt)
4 0 k=1

dt

!2

dt =

∞
1X
a2 k 2 .
8 k=1 k

(8)

The last expression attains the global minimum when a21 = 2, a2 = a3 =
. . . = 0. Hence,
√ the leading order term in the bias expression (7) is minimal
for the taper 2 sin(πt) (This result was obtained by Papoulis [9]). In [9],
Papoulis extends ν(t) to be zero outside of (0, 1), and therefore has a Fourier
5

integral representation of ν(t). We have extended ν(t) to be periodic and
vanish at each integer value. Since ν(t) is optimized for t ∈ [0, 1], both
representations are valid.
Equation (8) implies the more general result:
√
Theorem 3.1 vk (t) = 2 sin(πkt) (k = 1, 2, . . .) is the only system of functions in L2 [0, 1] which satisfy the requirements:
(i)

R1 2
v (t)dt
0

k

= 1, and

∞
(ii)
vk minimizes −∞
|V (k) (f )|2 f 2 df in the subspace of functions orthogonal to v1 , . . . , vk−1 .

R

The kth minimum value is

R∞

−∞

|V (k) (f )|2 f 2 df =

k2
4

.

√
We name vk (t) = 2 sin(πkt) (k = 1, 2, . . .) the continuous time minimum
bias tapers. The Fourier transform of vk (t) is


h 

k
k
e−iπ(f − 2 )  sin π f − 2
(k)


√
V (f ) =
i 2  π f − k2



f − k−1
2

= e−iπ(

i

h 

− (−1)

sin π f +
k

k
2





π f+

k
2

i 





k 2 sin πf − πk
2
)·
 2  .
√ 
k
4 2π f 2 − 2

Thus, |V (k) (f )| decays as f −2 for large frequencies.

3.2

Discrete Time Case

We now consider the discrete time domain {1, 2, . . . , N } with the corresponding normalized frequency domain [− 12 , 12 ]. A taper is a vector, ν =
P
2
(ν1 , . . . , νN ), normalized by N
n=1 (νn ) = 1. By the same argument as in the
previous section, the leading order term of the bias of a taper spectral estiR1
mate is proportional to the local bias, −2 1 |V (f )|2 f 2 df , where the frequency
2
window, V (f ), is defined in Eq. (5).
Lemma 3.2 For the frequency window of a discrete time taper,
Z

1
2

− 12

|V (f )|2 f 2 df = ν A ν
6

∗

,

where A = [anm ] with
anm =

Z

1
2

− 12

( 1
i2π(n−m)f

e

2

f df =

if n = m ;
if n 6= m .

12
(−1)n−m
2π 2 (n−m)2

Corollary 3.3 The tapers ν (1) , ν (2) , . . . , ν (N ) , defined by the requirements
(i)

PN

(ii)

ν (k) minimizes

(k) 2
n=1 (νn )

orthogonal to

= 1,
R

1
2

− 12

|V (k) (f )|2 f 2 df in the subspace of vectors

ν (1) , . . . , ν (k−1) ,

are the eigenvectors of the matrix A sorted in the increasing order of the
R1
eigenvalues. The integral −2 1 |V (k) (f )|2 f 2 df is equal to the kth eigenvalue.
2

We name ν (1) , ν (2) , . . . , ν (N ) the discrete minimum bias tapers. They
can be approximated by the sinusoidal tapers, v(1) , v(2) , . . . , v(N ) , which are
discrete analogs of the continuous time minimum
bias tapers. Namely, we
q
(k)
(k) T
2
(k)
(k)
define v = (v1 , . . . , vN ) with vn = N +1 sin Nπkn
, k = 1, 2, . . . , N .
+1
Lemma 3.4 The sinusoidal tapers, v(1) , v(2), . . . , v(N ) , form an orthonormal
k2
1 + O( N1 ) .
basis in RN and the local bias of v(k) is 4N
2
Corollary 3.5 The multitaper estimate (3) using K sinusoidal tapers has
P
k2
K2
local bias equal to K
k=1 µk 4N 2 + O( N 3 ) and has the following representation:
Ŝ(f ) =

K
X

µj
j
j
|y(f +
) − y(f −
)|2 .
2(N
+
1)
2N
+
2
2N
+
2
j=1

The uniformly weighted estimate, µk =
2
O( K
).
N3

1
,
K

has local bias equal to

(9)
K2
12N 2

+

From (9), the reason for the low bias of the sinusoidal tapers is apparent:
the frequency sidelobe from y(f + 2Nj+2 ) cancels the sidelobe of y(f − 2Nj+2 ).
As a result, the sidelobe of y(f + 2Nj+2 ) minus y(f − 2Nj+2 ) is much smaller
than that of the periodogram.
Our preferred weighting is the parabolic weighting: µj = C(1 − j 2 /K 2 )
because the parabolic weighting minimizes the expected square error in kernel
smoothers as K and N tend to infinity. Since the weights decrease smoothly
to zero, the resulting estimate is smooth in frequency.
7

Corollary 3.6 The uniformly weighted multitaper estimate using K sinusoidal tapers can be computed in O(N ln N ) + O(KN ) operations while the
generic multitaper estimate requires O(KN ln N ) operations plus the cost of
computing the K tapers.
The following
result demonstrates
that the kth
i h
i spectral window is conh
k+1
k−1
k+1
k−1
,
,
−
∪
−
.
centrated on 2(N
+1) 2(N +1)
2(N +1)
2(N +1)
Corollary 3.7 The Fourier transform of v(k) equals
h



V (k) (f ) =



k
e−iπ((N +1)f − 2 )  sin N π f −

k
2(N +1)

i

h



i 

sin N π f + 2(Nk+1) 
k
h 
i − (−1)
h 
i
q

sin π f + 2(Nk+1)
i 2(N + 1)  sin π f − 2(Nk+1)
(N +1)f − k−1
2

= e−iπ(

sin Nπk
+1

)·q

h

·
2
2(N + 1) sin (πf ) −

i

πk
2
2
πk
sin 2(N +1)

sin (N + 1)πf −

.

q

Thus |V (k) (f )| = N2+1 for |f | = 2(Nk+1) , and |V (k) (f )| ∼ N 13/2 for f =
πk
1
O(1). In particular, |V (k) (f )| ≈ √2N
3/2 when f → 2 . In the intermediate
frequencies, 2(Nk+1) < f  12 , |V (k) (f )| decays as f12 .
Numerical evaluation (see Table 1) shows that kv(k) − ν (k) kL2 < 4(Nk+2)
(k)
− k νν(k) k
L∞
k
. (k·kL∞ is the supremum norm in the time domain.) Figure 1 plots the
2(N +2)
envelopes of |V (k) (f )|2 for both the minimum bias and sinusoidal tapers with
N = 200, k = 1. The spectral energy of both tapers is nearly identical for
|f | < .25. Near the Nyquist frequency, the spectral energy of the sinusoidal
taper is roughly three times larger than that of the minimum bias taper. In
the time domain, the sinusoidal tapers are virtually indistinguishable from
the MB tapers.

for all k. The same rate of convergence is observed in L∞ norm:

4

v(k)
kv(k) kL∞

Comparison of Spectral Localizations

We now compare the local bias and the spectral concentration of three families of orthonormal tapers: minimum bias (MB) tapers, sinusoidal tapers
and Slepian tapers. Both the local bias and the spectral concentration of the
8

<
L∞

Slepian tapers depend on the bandwidth parameter, w. For properties of the
Slepian tapers, we refer the reader to [12, 16, 17, 18, 19].
Since the MB tapers minimize the local bias, clearly the sinusoidal tapers
and the Slepian tapers have larger local bias. The only question
is whether
R 1/2 2 (k)
P
(f )|2 df ,
the difference is large or small. Table 2 gives the local bias, K
k=1 −1/2 f |V
of the three families of tapers for N = 50. The sinusoidal tapers come within
0.2% of achieving the optimal local bias. In contrast, the local bias of the
Slepian tapers can be many times larger.
We compute the local bias for three different values of the bandwidth, w.
The general pattern is that the kth Slepian taper has roughly the same local
bias as the MB taper does when N w < k < 2N w. The ratio of the local bias
of the Slepian tapers to that of the MB tapers is smallest at k ≈ 1.2N w. As
|k − 1.2N w| increases, the local bias rapidly departs from the optimal value.
Table 3 compares the spectral concentration of the tapers for N = 50 and
w = .08 . Both the MB tapers and the sinusoidal tapers are within 1.7% of
the optimal value, except for k = 8,R 9. Notice that 2N w = 8. Although the
w
ratio of the spectral concentration, −w
|V (f )|2 df , for the MB and sinusoidal
tapers to that of the Slepian tapers is usually very close to one, the ratio of
the spectral energy outside of the frequency
band |f | < w can be
quite large.
Rw
Rw
2
Thus our conclusions depend on using −w |V (f )| df and not 1− −w |V (f )|2 df
as the measure of frequency concentration.
Figures 1-3 compare |V (k) (f )|2 of the Slepian and MB tapers for N = 200.
For Figures 1 and 2, we select the Slepian parameter, w, equal to .01 so that
K = 2N w equals four. Figure 2 plots |V1 (f )|2 for the frequencies up to
f = 0.14. The central peak of the MB taper is more concentrated around
f = 0 than the Slepian taper is. The first sidelobe of the MB taper is visible
while the first Slepian sidelobe is much smaller.
Figure 1 plots the logarithm of |V1 (f )|2 over the entire frequency range.
The MB taper has smaller range bias in the frequency range |f | < 0.3w and
in the frequency range |f | > 0.13. In the middle frequency range, the Slepian
taper is clearly better. The Slepian penalty function maximizes the energy
inside the frequency band, [−w, w], and thus it is natural that the Slepian
tapers do better for f ∼ w. By using a discontinuous penalty function, the
Slepian spectral windows experience Gibbs phenomenon and decay only as
1
, (|V (f )|2 ∼ f12 ). The MB spectral windows decay as f12 , and thus, it is
f
natural that the MB tapers have lower bias for f ∼ O(1).
P
Figure 3 plots 3k=1 |V (k) (f )|2 for |f | < 0.14 and on this scale, the MB

9

tapers are clearly preferable to the Slepian tapers. For larger frequencies, the
P
(k)
(f )|2 , is very similar to Fig. 2
energy of the multitaper estimate, K
k=1 |V
on the logarithmic scale provided that K << N .
In summary, the sinusoidal tapers perform nearly as well as the MB tapers
while the Slepian tapers have several times larger local bias (except when
k ≈ 1.2N w). For k  2N w, the Slepian tapers have better broad-band bias
protection than the minimum bias tapers do. For k ∼ 2N w, the minimum
bias tapers provide both smaller local bias and better broad-band protection
due to the Gibbs phenomena which the Slepian tapers experience.

5

Local Error Analysis and Optimal Multitapering

We now give a local error analysis of MTSA and determine the optimal
number of tapers. Our results are the multitaper analog of the local error
analysis of the smoothed periodogram [4, 10]. We assume the time series is
a Gaussian processes and do not consider frequencies near f = 0 and f =
1/2. In this case, the variance of the multitaper estimate is approximately
P
2
Variance[Ŝ(f )] ≈ S(f )2 K
k=1 µk due to the orthonormality of the tapers
Asymptotically, the local bias of the multitaper estimate of Eq. (3) is
Bias[Ŝ] = S(f )

K
X
1
λk µk ,
µk − 1 + S 00 (f )
2
k=1
k=1
K
X

!

R 1/2

where λk = −1/2 f 2 |V (k) (f )|2 df . The second term is the MT generalization
P
of (7). When K
µ 6= 1, the MT estimate has bias even in white noise.
k=1
P k
When we require K
k=1 µk = 1, the local expected loss simplifies:
Theorem 5.1 For a Gaussian process, away from f = 0 and f = 1/2, the
P
expected square error of the multitaper spectral estimate (3) with K
k=1 µk = 1
is asymptotically (to leading order in K/N )
"
2

Bias + Variance ≈

K
X
1 00
S (f )
λk µk
2
k=1

#2

+ S(f )2

K
X

µ2k .

(10)

k=1

Theorem 5.2 The multitaper estimate which minimizes the local loss (10)
(with µk ≥ 0) is constructed with the minimum bias tapers.
10

Proof: We order the µk such that µ1 ≥ µ2 ≥ . . . ≥ µK and define µK+1 =
P
∗
0. Since the weights, µk are fixed, we need to minimize K
k=1 µk uk A uk
over all sets of K orthonormal tapers, u1 , . . . , uK . We split the series in K
subseries and minimize each subseries separately:
min

u1 ,...,uK

K
X

µk uk A u∗k =

min

u1 ,...,uK

k=1

K
X

≥

K
X



(µk − µk−1 ) 

K
X



uj A u∗j 

k
X
(k)

(µk − µk−1 )  min

(k)
(k)
u1 ,...,uk



(µk − µk−1 ) 

k
X


(k)∗ 

uj A uj

j=1



K
X

λA,j  =

j=1

k=1



j=1

k=1

k=1

=

k
X

µk λA,j , (11)

k=1

where the λA,j are the eigenvalues of A , given in increasing order. The
(k)
(k)
(k)
uj are subject to orthonormality constraints that uj · uj 0 = δj,j 0 , but are
otherwise independent and minimized separately. In the last line of (11), we
P
P
(k)
(k)∗
use Fan’s Theorem [6]: minu(k) ,...,u(k) kj=1 uj A uj = kj=1 λA,j , where the
1

k

(k)
uj

are again subject to orthonormality constraints. The theorem is now
proved because the MB tapers are precisely the eigenvectors of A .
Theorem 5.3 The uniformly weighted multitaper estimate using K sinusoidal tapers has an asymptotic local loss of
"
2

Bias + Variance '

S 00 (f )K 2
24N 2

#2

+

S(f )2
.
K

(12)

Corollary 5.4 The asymptotic local loss of (12) is minimized when the number of tapers is chosen as
12S(f )N 2
∼
S 00 (f )
"

Kopt

#2/5

.

(13)

Thus, the optimal number of tapers is proportional to N 4/5 and varies
with the ratio of S(f ) to S 00 (f ). Intuitively (13) shows that fewer tapers
should be used when the spectrum varies more rapidly. A key advantage
of the MB and sinusoidal tapers is that the tapers need not be recomputed
as K is changed. In contrast, the Slepian tapers are most efficient when
the bandwidth parameter, w, is chosen such that K ∼ 2N w. Thus, when
the number of tapers is changed, as in (13), the Slepian tapers should be
recomputed.
11

6

Smoothed Multitaper Estimates

In our own comparison of kernel smoothing and multitaper estimation [16],
we found that a smoothed multiple taper estimate worked best. We now
evaluate the expected error of the kernel smoothed multitaper estimator and
show that smoothing the logarithm of the multitaper estimate is useful for
estimating the logarithm of the spectrum.
We begin by evaluating that the quadratic estimator (2) which is equivalent to a kernel smoother estimates of the spectrum [4, 10]. Let Ŝ(f ) be
the quadratic spectral estimator (2), and smooth it with a kernel κ(·) of
halfwidth w :
Z w
g b
bb
+ g)dg ,
(14)
S(f ) =
κ( )S(f
w
−w
where w is the bandwidth parameter and κ(·) is a kernel smoother with
domain [−1, 1]. This can be rewritten as
b ) =
S(f
b

N
X

q̃ nm ei2π(m−n)f xn xm ,

(15)

n,m=1
w
where q̃ nm = qnm κ̂m−n with κ̂m = −w
κ(g/w)e2πimg dg. Thus smoothing replaces the original quadratic estimator with matrix [qnm ] by another
quadratic estimator with matrix Q̃ = [q̃ nm ]. By Theorem 5.2, this hybrid
method cannot outperform the pure multitaper method with minimum bias
tapers.
We now show that combining kernel smoothing with multitapering does
improve the estimation of the logarithm of the spectral density, θ(f ) =
log[S(f )]. One standard approach is to kernel smooth the logarithm of the
tapered periodogram. This approach has the disadvantage that |y(f )|2 has a
χ22 distribution and log[χ22 ] has a long lower tail of its distribution. As a result, log[|y(f )|2 ] has an appreciable bias and its variance is inflated by π 2 /6.
A common alternative is to estimate the spectrum either by kernel smoothing
or by multitapering and then to take logarithms. This approach has the disadvantage that the smoothed spectral estimate tends to be more sensitive to
nonlocal bias effects than the corresponding smoothed log-spectral estimate.
To robustify the log-spectral estimate while reducing the variance inflation from the long tail, we propose the following hybrid estimate: 1) compute
the multitaper estimate using the sinusoidal tapers with µk = K1 and then
2) smooth θ̂M T (f ) ≡ ln[ŜM T (f )] − BK /K, where BK is the bias of ln[χ22K ].

R

12

(BK ≡ ψ(K) − ln K where ψ(K) is the digamma function). For white noise,
the variance of θ̂M T (f ) = ψ 0 (K)∼ K1 + 2K1 2 , so the variance enhancement
from the logarithm tends rapidly to zero.
In [15], we show that the asymptotic error for this scheme is
K2
bk w 2 +
24N 2

"
00

θ (f )

2

#2

1
Cκ
1+
+
Nw
2K


2

,

(16)

where bκ and Cκ are constants which depend on the kernel shape. In (16), we
assume uniformly weighted sinusoidal tapers are used and 1  K  N w. In
1
) is the variance enhancement from the logarithmic
(16), one factor of (1 + 2K
1
transformation and one factor of (1 + 2K
) arises in the variance calculation
of (15) with sinusoidal tapers. Optimizing (16) with respect to both w and
K yield w ∼ N −1/5 and K ∼ N 8/15 , thus the smoothing halfwidth w is much
larger than K/N . The expected error (16) depends weakly on K provided
that 1  K  N w.
For simplicity, we set K = N 8/15 and optimize (16) with respect to the
halfwidth w. The resulting halfwidth depends on θ00 (f ): wopt (θ00 (f )) with
wopt ∼ |θ00 (f )|−2/5 N −1/5 . Thus when the log-spectrum varies rapidly, the
halfwidth should be reduced as |θ00 (f )|−2/5 .
Since θ00 (f ) is unknown, we consider two stage estimators which begin by
making a preliminary estimate of θ00 (f ) prior to estimating θ(f ). We then
insert the estimate θ00d
(f ) into the expression for wopt : w(f ) = wopt (θc00 (f )) and
use a variable halfwidth kernel smoother with halfwidth ŵ(f ) to estimate
θ(f ). Multiple stage kernel estimators are described in [2, 7, 13, 14, 15].
These multiple stage schemes have a convergence rate of N −4/5 and have a
relative convergence rate of at least N −2/9 . A more detailed description is
given in [2, 13, 15].

7

Application

We now compare spectral estimates on an actual data series. We use the
microwave scattering data set which is described in [16]. The data measures
turbulent plasma fluctuations in the Tokamak Fusion Test Reactor at Princeton. The spectrum is dominated by a 1 MHz peak which is quasicoherent.
The spectral density varies by over five orders of magnitude.

13

The bias versus variance trade-off of Sec. 5 shows that fewer tapers should
be used near the peak. To make the spectral estimate smooth, a parabolic
weighting of the tapers is used as described in Sec. 3. To determine how
many tapers to use locally, we use the multiple stage “plug-in” method as
described in the previous section; i.e. we determine the number of tapers
using a pre-estimate on the same data. To reduce the fluctuations from the
estimate of the optimal number, we use a longer data segment to determine
the number of tapers at each frequency. We find the optimal number of
sinusoidal tapers is roughly 24 for frequencies in the 200 to 800 kHz range.
Near the 1 MHz peak, as few as 12 tapers are used to minimize the local bias
error. Between 1300 and 2400 kHz, the spectrum is flatter and we use up to
40 tapers.
The dotted line is the sinusoidal multitaper estimate, and the solid, more
wiggly, curve is the corresponding Slepian estimate using 24 tapers with
w = 60 kHz. The 1 MHz peak is poorly resolved in the Slepian estimate, and
the regions of high curvature are artificially flattened. For f ≥ 1.5 MHz, the
Slepian estimate is artificially bumpy due to statistical noise. The variable
taper number estimate suppresses these bumps by averaging over a larger
frequency halfwidth. We have also used a variable taper number estimate
with the Slepian tapers. Since the Slepian parameter, w was fixed at 100 kHz
to allow for forty tapers, the artificial broadening was even more exreme.
Comparing with a converged estimate of the spectrum based on N =
45, 000 shows that the sinuosidal taper estimate is more accurate. Another
significant difference is that the Slepian multitaper estimate requires much
more CPU time than the sinusoidal multitaper estimate.

8

Conclusion

We have
q proposed and analyzed the minimum bias and the sinusoidal tapers,
(k)
, for multitaper spectral estimation. The resulting sinuvn = N2+1 sin Nπkn
+1
j
soidal multitaper spectral estimate is Ŝ(f ) = 2K(N1 +1) K
j=1 |y(f + 2N +2 ) −
y(f − 2Nj+2 )|2 . The sinusoidal tapers have low bias because the frequency
sidelobe from y(f + 2Nj+2 ) cancels the sidelobe of y(f − 2Nj+2 ).

P

R 1/2

The minimum bias tapers minimize the local bias, −1/2 f 2 |V (k) b(f )|2 df ,
and have good broad-band bias protection as well. Asymptotically, the
quadratic spectral estimate which minimizes the expected local square er14

ror is a multiple taper estimate using the minimal bias tapers.
The sinusoidal tapers have
 a simple analytic form and approximate the
minimum bias tapers to O N1 . The kth sinusoidal taper has its spectral
k+1
k−1
≤ |f | ≤ 2(N
.
energy concentrated in the frequency bands 2(N
+1)
+1)
The minimum bias and sinusoidal tapers have no auxiliary bandwidth parameter, and the bandwidth of the spectral estimate is determined solely by
the number of used tapers. By adaptively adding and deleting tapers, a multitaper estimate with the optimal convergence properties of kernel smoothers
can be constructed. In contrast, the Slepian tapers need to be recomputed
with a different bandwidth. Thus the Slepian tapers are only practical for
fixed bandwidth estimation and this is inherently inefficient.

9

Appendix: Multitaper decomposition of kernel estimates

In Sec. 6, we showed that kernel smoother estimators (14) have an equivalent
multitaper representation (3) We now show that the equivalent multitapers
of some popular kernel smoother estimates of the spectrum strongly resemble
the MB/sinusoidal tapers. In one special case, this corresondence is exact;
i.e. the smoothed periodogram can be exactly decomposed into MB tapers.
Theorem 9.1 Let w = 21 and κ(f ) be the parabolic kernel, κ(f ) = 23 − 6f 2 .
The eigenvectors of the kernel smoothed periodogram are exactly the discrete
minimum bias tapers.
Proof: The Q̃ matrix in (15) can be calculated explicitly for this case.
We find Q̃ = [bnm ] = N1 ( 32 I − 6A ) where A is the matrix from Lemma 3.2.
Thus Q̃ and A have the same eigenvectors.
To illustrate that this result is typical even when we apply a taper and
smooth over a small band, we consider a smoothed tapered periodogram
with N = 200. We use Tukey’s split-cosine taper [16] and then smooth the
estimate with a square box kernel with a halfwidth of .01. We then evaluate
the corresponding Q̃ matrix and compute its eigenvectors. Figure
q 5 displays
the first 4 eigenvectors. They are very close to the sinusoids N2+1 sin Nπkn
.
+1
Table 4 shows that this spectral estimate is virtually a K = 4 multiple taper
15

spectral estimate. After k > 4, the eigenvalues decrease sharply, and these
higher eigenvectors contribute very little to the overall estimate.

References
[1] M. Amin, “Optimal estimation of evolutionary spectra,”
I.E.E.E. Trans. on Signal Processing vol. 42, 2?, 1994.
[2] T. Brockman, Th. Gasser and E. Hermann, “Locally adaptive
bandwidth choice for kernel regression estimators,” J. Amer. Stat.
Assoc. 88, 1302-1309 (1994).
[3] T. P. Bronez, Nonparametric Spectral Estimation of Irregularly
Sampled Multidimensional Random Processes, PhD Thesis, Arizona State University, 1985.
[4] U. Grenander and M. Rosenblatt, Statistical Analysis of Stationary Time Series, New York: Wiley, 1957.
[5] A. N. Kolmogorov and S. V. Fomin, Reele Funktionen und Funktionalanalysis, Section 7.3.2, Berlin: VEB Deutscher Verlag der
Wissenschaften, 1975.
[6] A. W. Marshall and I. Olkin, Inequalities: Theory of Majorization
and its Applications, p. 511, New York: Academic Press 1979.
[7] H.-G. Müller and U. Stadtmüller, “Variable bandwidth kernel
estimators of regression curves,” Annals of Statistics, vol. 15, pp.
182-201, 1987.
[8] C. T. Mullis and L. L. Scharf, in Advances in Spectrum Analysis,
S. Haykin, Ed., New York: Prentice-Hall, 1990, Chapter 1, pp.
1-57.
[9] A. Papoulis, “Minimum bias windows for high resolution spectral
estimates,” IEEE Trans. Information Theory, vol. 19, pp. 9-12,
1973.
[10] E. Parzens, “On asymptotically efficient consistent estimates
of the spectral density of a stationary time series,” J. Royal
Stat. Soc., vol. 19, pp. 303-322, 1958.
16

[11] J. Park, C. R. Lindberg, and F. L. Vernon, “Multitaper Spectral
analysis of high frequency seismograms,” J. Geophys. Res., vol.
92B, pp. 12765-12684, 1987.
[12] D. Percival and A. Walden, Spectral Analysis for Physical Applications: Multitaper and Conventional Univariate Techniques.
Cambridge: Cambridge University Press, 1993.
[13] K. S. Riedel, “Optimal kernel estimation of evolutionary spectra,”
I.E.E.E. Trans. on Signal Processing vol. 41, 2439-2447, 1993.
[14] K.S. Riedel and A. Sidorenko, “Function estimation using data
adaptive kernel smoothers- How much smoothing?” Computers
in Physics vol. 8, 402-409, 1994.
[15] K. S. Riedel and A. Sidorenko, “Smoothed log-multitaper spectral estimation and data adaptive implementation,” submitted
for publication.
[16] K. S. Riedel, A. Sidorenko, and D. J. Thomson, “Spectral density
estimation for plasma fluctuations I: Comparison of methods,”
Physics of Plasmas vol. 1, pp. 485-500.
[17] D. Slepian, “Prolate spheroidal wave functions, Fourier analysis,
and uncertainty - V: the discrete case,” Bell System Tech. J., vol.
5, pp. 1371-1429, 1978.
[18] D. J. Thomson, “Spectrum estimation and harmonic analysis,”
Proc. IEEE, vol. 70, pp. 1055-1096, 1982.
[19] D. J. Thomson, “Quadratic inverse spectrum estimates: applications to paleoclimatology,” Phil. Trans. R. Soc. Lond. A, vol.
332, p. 539-597, 1990.

17

Table Captions:
Table 1: Convergence of the sinusoidal tapers toR the minimum bias tapers.
P
1/2
2 PK
2
(k)
(f )|2 df ,
Table 2: Normalized bias term, 4(N + 1)2 K
k=1 −1/2 f
k=1 f |V
for N = 50.
R w PK
(k)
(f )|2 df , for N = 50.
Table 3: Spectral concentration, −w
k=1 |V
Table 4: Eigenvectors of the smooth tapered periodogram estimator.

Figure Captions:
Figure 1: Spectral energy of the minimum bias, sinusoidal and Slepian tapers,
N = 200, k = 1, w = .01.
Figure 2: Spectral energy of theRminimum bias and sinusoidal tapers, N = 200, k = 1.
P
1/2
Figure 3 Spectral energy, 3k=1 −1/2 f 2 |V (k) (f )|2 df , of the minimum bias and
Slepian tapers, N = 200, K = 3, w = .01.
Figure 4: Estimated spectral density of the plasma fluctuations. Dashed
line is sinusoidal multitaper estimate and solid line is estimate using Slepian
tapers with w = 60 kHz. Because the Slepian tapers have a fixed bandwidth,
the corresponding estimate spectral density at 1 MHz is artificially broadened
while being undersmoothed for f ≥ 1.5 MHz.
Figure 6: First eigenvectors of the smooth tapered periodogram estimator.

18

19

20

21

22

23

24

Table 1. Convergence of the sinusoidal tapers to the minimum bias tapers
N

maxk

n

20
50
200
800

N +2
kv(k)
k

− ν (k) kL2

o

1
2
3
4
5
6
7
8
9
10

maxk

0.24750
0.24844
0.24852
0.24844

Table 2. Normalized bias term, 4(N + 1)2
N = 50
K

(

Minimum
bias tapers
1.0095
5.0475
14.1328
30.2846
55.5217
91.8634
141.3284
205.9362
287.7056
388.6562

Sinusoidal
tapers
1.0116
5.0580
14.1622
30.3475
55.6366
92.0528
141.6185
206.3570
288.2899
389.4409

w=0.04
1.3439
5.7724
18.0651
58.9520
154.4818
305.4382
496.7959
721.1743
976.5088
1262.4251

25

R 1/2
k=1 −1/2

PK

N +2
k

(k)
− k νν(k) k
L∞
0.4602
0.4760
0.4829
0.4844

v(k)
kv(k) kL∞

f 2 |V (k) (f )|2 df , for

Slepian tapers
w=0.08
w=0.16
2.6316
5.1039
10.5484
20.4670
23.8086
46.1953
42.5181
82.4018
66.9996
129.2087
99.1800
186.7507
150.0103
255.1797
251.8833
334.6717
437.5993
425.4379
702.1523
527.7433

)
L∞

Table 3. Spectral concentration,
k Minimum bias tapers
1
.9997
2
.9988
3
.9972
4
.9940
5
.9888
6
.9760
7
.9381
8
.6084
9
.1688
10
.0637

Rw

−w

|V (f )|2 df , for N = 50, w = 0.08

Sinusoidal tapers
.9997
.9988
.9972
.9937
.9887
.9753
.9417
.6247
.1780
.0624

Slepian tapers
1.
.9999999
.9999989
.99997
.9995
.9928
.9380
.7002
.2981
.0628

Table 4. Eigenvectors of the smooth tapered periodogram estimator
k

1
2
3
4
5
6
7

Weight of the eigenvector
Normalized local bias
R 1/2
λk (B)/tr(B)
4(N + 1)2 −1/2 f 2 |V (f )|2 df
.2856
.2828
.2519
.1416
.0340
.0037
.0002

1.5138
4.7371
9.6254
19.2095
33.7118
51.3616
72.9747

26

Local bias in comparison with
the minimum bias taper (ratio)
1.509
1.181
1.067
1.198
1.345
1.423
1.486

