Thomson’s Multitaper Method Revisited

Santhosh Karnik, Justin Romberg, Mark A. Davenport

Abstract

Thomson’s multitaper method estimates the power spectrum of a signal from $N$ equally spaced samples by averaging $K$ tapered periodograms. Discrete prolate spheroidal sequences (DPSS) are used as tapers since they provide excellent protection against spectral leakage. Thomson’s multitaper method is widely used in applications, but most of the existing theory is qualitative or asymptotic. Furthermore, many practitioners use a DPSS bandwidth $W$ and number of tapers that are smaller than what the theory suggests is optimal because the computational requirements increase with the number of tapers. We revisit Thomson’s multitaper method from a linear algebra perspective involving subspace projections. This provides additional insight and helps us establish nonasymptotic bounds on some statistical properties of the multitaper spectral estimate, which are similar to existing asymptotic results. We show using $K=2NW-O(\log(NW))$ tapers instead of the traditional $2NW-O(1)$ tapers better protects against spectral leakage, especially when the power spectrum has a high dynamic range. Our perspective also allows us to derive an $\epsilon$ -approximation to the multitaper spectral estimate which can be evaluated on a grid of frequencies using $O(\log(NW)\log\tfrac{1}{\epsilon})$ FFTs instead of $K=O(NW)$ FFTs. This is useful in problems where many samples are taken, and thus, using many tapers is desirable.

^†^†S. Karnik, J. Romberg, and M. A. Davenport are with the School of Electrical and Computer Engineering, Georgia Institute of Technology, Atlanta, GA, 30332 USA (e-mail: [email protected], [email protected], [email protected]). This work was partially supported by NSF grant CCF-1409406, a grant from Lockheed Martin, and a gift from the Alfred P. Sloan Foundation. A preliminary version of this paper highlighting some of the key results also appeared in [1].

1 Introduction

Perhaps one of the most fundamental problems in digital signal processing is spectral estimation, i.e., estimating the power spectrum of a signal from a window of $N$ evenly spaced samples. The simplest solution is the periodogram, which simply takes the squared-magnitude of the discrete time Fourier transform (DTFT) of the samples. Obtaining only a finite number of samples is equivalent to multiplying the signal by a rectangular function before sampling. As a result, the DTFT of the samples is the DTFT of the signal convolved with the DTFT of the rectangular function, which is a slowly-decaying sinc function. Hence, narrow frequency components in the true signal appear more spread out in the periodogram. This phenomenon is known as “spectral leakage”.

The most common approach to mitigating the spectral leakage phenomenon is to multiply the samples by a taper before computing the periodogram. Since multiplying the signal by the taper is equivalant to convolving the DTFT of the signal with the DTFT of the taper, using a taper whose DTFT is highly concentrated around $f=0$ will help mitigate the spectral leakage phenomenon. Numerous kinds of tapers have been proposed [2] which all have DTFTs which are highly concentrated around $f=0$ . The Slepian basis vectors, also known as the discrete prolate spheroidal sequences (DPSSs), are designed such that their DTFTs have a maximal concentration of energy in the frequency band $[-W,W]$ subject to being orthonormal [3]. The first $\approx 2NW$ of these Slepian basis vectors have DTFTs which are highly concentrated in the frequency band $[-W,W]$ . Thus, any of the first $\approx 2NW$ Slepian basis vectors provides a good choice to use as a taper.

In 1982, David Thomson [4] proposed a multitaper method which computes a tapered periodogram for each of the first $K\approx 2NW$ Slepian tapers, and then averages these periodograms. Due to the spectral concentration properties of the Slepian tapers, Thomson’s multitaper method also does an excellent job mitigating spectral leakage. Furthermore, by averaging $K$ tapered periodograms, Thomson’s multitaper method is more robust than a single tapered periodogram. As such, Thomson’s multitaper method has been used in a wide variety of applications, such as cognitive radio [5, 6, 7, 8, 9], digital audio coding [10, 11], as well as to analyze EEG [12, 13] and other neurological signals [14, 15, 16, 17, 18], climate data [19, 20, 21, 22, 23, 24, 25], breeding calls of Adélie penguins [26] and songs of other birds [27, 28, 29, 30], topography of terrestrial planets [31], solar waves[32], and gravitational waves [33].

The existing theoretical results regarding Thomson’s multitaper method are either qualitative or asymptotic. Here, we provide a brief summary of these results. See [4, 34, 35, 8, 36, 37] for additional details. Suppose that the power spectral density $S(f)$ of the signal is “slowly varying”. Let $\widehat{S}^{\text{mt}}_{K}(f)$ denote the multitaper spectral estimate. Then, the following results are known.

•

The estimate is approximately unbiased, i.e., $\mathbb{E}\widehat{S}^{\text{mt}}_{K}(f)\approx S(f)$ .
•

The variance is roughly $\operatorname{Var}[\widehat{S}^{\text{mt}}_{K}(f)]\approx\tfrac{1}{K}S(f)^{2}$ .
•

For any two frequencies $f_{1},f_{2}$ that are at least $2W$ apart, $\widehat{S}^{\text{mt}}_{K}(f_{1})$ and $\widehat{S}^{\text{mt}}_{K}(f_{2})$ are approximately uncorrelated.
•

The multitaper spectral estimate $\widehat{S}^{\text{mt}}_{K}(f)$ has a concentration behavior about its mean which is similar to a scaled chi-squared random variable with $2K$ degrees of freedom.
•

If the power spectral density $S(f)$ is twice differentiable, then choosing a bandwidth of $W=O(N^{-1/5})$ and using $K\approx 2NW=O(N^{4/5})$ tapers minimizes the mean squared error of the multitaper spectral estimate $\mathbb{E}|\widehat{S}^{\text{mt}}_{K}(f)-S(f)|^{2}$ .

These asymptotic results demonstrate that using more than a small constant number of tapers improves the quality of the multitaper spectral estimate. However, using more tapers increases the computational requirements. As a result, many practitioners often use significantly fewer tapers than what is optimal.

The main contributions of this work are as follows:

•

We revisit Thomson’s multitaper method from a linear algebra based perspective. Specifically, for each frequency $f$ , the multitaper spectral estimate computes the $2$ -norm of the projection of the vector of samples onto a $K$ -dimensional subspace. The subspace chosen can be viewed as the result of performing principle component analysis on a continuum of sinusoids whose frequency is between $f-W$ and $f+W$ .
•

Using this perspective, we establish non-asymptotic bounds on the bias, variance, covariance, and probability tails of the multitaper spectral estimate. These non-asymptotic bounds are comparable to the known asymptotic results which assume that the spectrum is slowly varying. Also, these bounds show that using $K=2NW-O(\log(NW))$ tapers (instead of the traditional choice of $K=\left\lfloor 2NW\right\rfloor-1$ or $\left\lfloor 2NW\right\rfloor-2$ tapers) provides better protection against spectral leakage, especially in scenarios where the spectrum has a large dynamic range.
•

We also use this perspective to demonstrate a fast algorithm for evaluating an $\epsilon$ -approximation of the multitaper spectral estimate on a grid of $L$ equally spaced frequencies. The complexity of this algorithm is $O(L\log L\log(NW)\log\tfrac{1}{\epsilon})$ while the complexity of evaluating the exact multitaper spectral estimate on a grid of $L$ equally spaced frequencies is $O(KL\log L)$ . Computing the $\epsilon$ -approximation is faster than the exact multitaper provided $K\gtrsim\log(NW)\log\tfrac{1}{\epsilon}$ tapers are used.

The rest of this work is organized as follows. In Section 2 we formally introduce Thompson’s multitaper spectral estimate, both from the traditional view involving tapered periodograms as well as from a linear algebra view involving projections onto subspaces. In Section 3, we state non-asymptotic bounds regarding the bias, variance, and probability concentration of the multitaper spectral estimate. These bounds are proved in Appendix A. In Section 4, we state our fast algorithm for evaluating the multitaper spectral estimate at a grid of frequencies. The proofs of the approximation error and computational requirements are in Appendix B. In Section 5, we show numerical experiments which demonstrate that using $K=2NW-O(\log(NW))$ tapers minimizes the impact of spectral leakage on the multitaper spectral estimate, and that our $\epsilon$ -approximation to the multitaper spectral estimate can be evaluated at a grid of evenly spaced frequencies significantly faster than the exact multitaper spectral estimate. We finally conclude the paper in Section 6.

2 Thomson’s Multitaper Method

2.1 Traditional view

Let $x(n)$ , $n\in\mathbb{Z}$ be a stationary, ergodic, zero-mean, Gaussian process. The autocorrelation and power spectral density of $x(n)$ are defined by

R_{n}=\mathbb{E}\left[x(m)\overline{x(m+n)}\right]\quad\text{for}\quad m,n\in\mathbb{Z},

and

S(f)=\sum_{n=-\infty}^{\infty}R_{n}e^{-j2\pi fn}\quad\text{for}\quad f\in\mathbb{R}

respectively. The goal of spectral estimation is to estimate $S(f)$ from the vector $\boldsymbol{x}\in\mathbb{C}^{N}$ of equispaced samples $\boldsymbol{x}[n]=x(n)$ for $n\in[N]$ . ¹¹1We use the notation $[N]$ to denote the set $\{0,\ldots,N-1\}$ . We will also use $[K]$ and $[L]$ in a similar manner.

One of the earliest, and perhaps the simplest estimator of $S(f)$ is the periodogram [38, 39]

\widehat{S}(f)=\dfrac{1}{N}\left|\sum_{n=0}^{N-1}\boldsymbol{x}[n]e^{-j2\pi fn}\right|^{2}.

This estimator can be efficiently evaluated at a grid of evenly spaced frequencies via the FFT. However, the periodogram has high variance and suffers from the spectral leakage phenomenon [40].

A modification to the periodogram is to pick a data taper $\boldsymbol{w}\in\mathbb{R}^{N}$ with $\|\boldsymbol{w}\|_{2}=1$ , and then weight the samples by the taper as follows

\widehat{S}_{\boldsymbol{w}}(f)=\left|\sum_{n=0}^{N-1}\boldsymbol{w}[n]\boldsymbol{x}[n]e^{-j2\pi fn}\right|^{2}.

If $\boldsymbol{w}[n]$ is small near $n=0$ and $n=N-1$ , then this “smoothes” the “edges” of the sample window. Note that the expectation of the tapered periodogram is given by a convolution of the true spectrum and the spectral window of the taper,

\mathbb{E}[\widehat{S}_{\boldsymbol{w}}(f)]=S(f)\circledast\left|\widetilde{\boldsymbol{w}}(f)\right|^{2}

where

\widetilde{\boldsymbol{w}}(f)=\sum_{n=0}^{N-1}\boldsymbol{w}[n]e^{-j2\pi fn}.

Hence, a good taper will have its spectral window $\left|\widetilde{\boldsymbol{w}}(f)\right|^{2}$ concentrated around $f=0$ so that $\mathbb{E}\left[\widehat{S}_{\boldsymbol{w}}(f)\right]=S(f)\circledast\left|\widetilde{\boldsymbol{w}}(f)\right|^{2}\approx S(f)$ , i.e., the tapered periodogram will be approximately unbiased.

For a given bandwidth parameter $W>0$ , we can ask “Which taper maximizes the concentration of its spectral window, $\left|\widetilde{\boldsymbol{w}}(f)\right|^{2}$ , in the interval $[-W,W]$ ?” Note that we can write

\int_{-W}^{W}\left|\widetilde{\boldsymbol{w}}(f)\right|^{2}\,df=\boldsymbol{w}^{*}\boldsymbol{B}\boldsymbol{w}

where $\boldsymbol{B}$ is the $N\times N$ prolate matrix [41, 42] whose entries are given by

\boldsymbol{B}[m,n]=\begin{cases}\dfrac{\sin[2\pi W(m-n)]}{\pi(m-n)}&\text{if}\ m\neq n,\\ 2W&\text{if}\ m=n.\end{cases}

Hence, the taper whose spectral window, $\left|\widetilde{\boldsymbol{w}}(f)\right|^{2}$ , is maximally concentrated in $[-W,W]$ is the eigenvector of $\boldsymbol{B}$ corresponding to the largest eigenvalue.

The prolate matrix $\boldsymbol{B}$ was first studied extensively by David Slepian [3]. As such, we will refer to the orthonormal eigenvectors $\boldsymbol{s}_{0},\boldsymbol{s}_{1},\ldots,\boldsymbol{s}_{N-1}\in\mathbb{R}^{N}$ of $\boldsymbol{B}$ as the Slepian basis vectors, where the corresponding eigenvalues $\lambda_{0}\geq\lambda_{1}\geq\cdots\geq\lambda_{N-1}$ are sorted in descending order. Slepian showed that all the eigenvalues of $\boldsymbol{B}$ are distinct and strictly between $0$ and $1$ . Furthermore, the eigenvalues of $\boldsymbol{B}$ exhibit a particular clustering behavior. Specifically, the first slightly less than $2NW$ eigenvalues are very close to $1$ , and the last slightly less than $N-2NW$ eigenvalues are very close to $0$ .

While $\boldsymbol{s}_{0}$ is the taper whose spectral window is maximally concentrated in $[-W,W]$ , any of the Slepian basis vectors $\boldsymbol{s}_{k}$ for which $\boldsymbol{s}_{k}^{*}\boldsymbol{B}\boldsymbol{s}_{k}=\lambda_{k}\approx 1$ will also have a high energy concentration in $[-W,W]$ , and thus, make good tapers. Thomson [4] proposed a multitaper spectral estimate by using each of the first $K\approx 2NW$ Slepian basis vectors as tapers, and taking an average of the resulting tapered periodograms, i.e.,

\widehat{S}^{\text{mt}}_{K}(f)=\dfrac{1}{K}\sum_{k=0}^{K-1}\widehat{S}_{k}(f)\quad\text{where}\quad\widehat{S}_{k}(f)=\left|\sum_{n=0}^{N-1}\boldsymbol{s}_{k}[n]\boldsymbol{x}[n]e^{-j2\pi fn}\right|^{2}.

The expectation of the multitaper spectral estimate satisfies

\mathbb{E}\left[\widehat{S}^{\text{mt}}_{K}(f)\right]=S(f)\circledast\psi(f)

where

\psi(f)=\dfrac{1}{K}\sum_{k=0}^{K-1}\left|\sum_{n=0}^{N-1}\boldsymbol{s}_{k}[n]e^{-j2\pi fn}\right|^{2}

is known as the spectral window of the multitaper spectral estimate. It can be shown that when $K\approx 2NW$ , the spectral window $\psi(f)$ approximates $\tfrac{1}{2W}\mathbbm{1}_{[-W,W]}(f)$ on $f\in[-\tfrac{1}{2},\tfrac{1}{2}]$ . Thus, the multitaper spectral estimate behaves in expectation like a smoothed version of the true spectrum $S(f)$ .

It can be shown that if the spectrum $S(f)$ is slowly varying around a frequency $f$ , then the tapered spectral estimates $\widehat{S}_{k}(f)$ are approximately uncorrelated, and $\operatorname{Var}[\widehat{S}_{k}(f)]\approx S(f)^{2}$ . Hence, $\operatorname{Var}[\widehat{S}^{\text{mt}}_{K}(f)]\approx\tfrac{1}{K}S(f)^{2}$ . Thus, Thomson’s multitaper method produces a spectral estimate whose variance is a factor of $K\approx 2NW$ smaller than the variance of a single tapered periodogram.

As we increase $W$ , the width of the spectral window $\psi(f)$ increases, which causes the expectation of the multitaper spectral estimate to be further smoothed. However, increasing $W$ also allows us to increase the number of tapers $K\approx 2NW$ , which reduces the variance of the multitaper spectral estimate. Intuitively, Thomson’s multitaper method introduces a tradeoff between resolution and robustness.

2.2 Linear algebra view

Here we provide an alternate perspective on Thomson’s multitaper method which is based on linear algebra and subspace projections. Suppose that for each frequency $f\in\mathbb{R}$ , we choose a low-dimensional subspace $\mathcal{S}_{f}\subset\mathbb{C}^{N}$ , and form a spectral estimate by computing $\|\text{proj}_{\mathcal{S}_{f}}(\boldsymbol{x})\|_{2}^{2}$ , i.e., the energy in the projection of $\boldsymbol{x}$ onto the subspace $\mathcal{S}_{f}$ . One simple choice is the one-dimensional subspace $\mathcal{S}_{f}=\text{span}\{\boldsymbol{e}_{f}\}$ where

\boldsymbol{e}_{f}=\begin{bmatrix}1&e^{j2\pi f\cdot 1}&e^{j2\pi f\cdot 2}&\cdots&e^{j2\pi f\cdot(N-1)}\end{bmatrix}^{T}

is a vector of equispaced samples from a complex sinusoid with frequency $f$ . For this choice of $\mathcal{S}_{f}$ , we have

\left\|\text{proj}_{\mathcal{S}_{f}}(\boldsymbol{x})\right\|_{2}^{2}=\dfrac{\left|\left<\boldsymbol{e}_{f},\boldsymbol{x}\right>\right|^{2}}{\|\boldsymbol{e}_{f}\|_{2}^{2}}=\dfrac{1}{N}\left|\sum_{n=0}^{N-1}\boldsymbol{x}[n]e^{-j2\pi fn}\right|^{2},

which is exactly the classic periodogram.

We can also choose a low-dimensional subspace $\mathcal{S}_{f}$ which minimizes the average representation error of sinusoids $\boldsymbol{e}_{f^{\prime}}$ with frequency $f^{\prime}\in[f-W,f+W]$ for some small $W>0$ , i.e.,

\operatorname*{\text{minimize}}_{\begin{subarray}{c}\mathcal{S}_{f}\subset\mathbb{C}^{N}\\ \dim(\mathcal{S}_{f})=K\end{subarray}}\int_{f-W}^{f+W}\left\|\boldsymbol{e}_{f^{\prime}}-\text{proj}_{\mathcal{S}_{f}}(\boldsymbol{e}_{f^{\prime}})\right\|_{2}^{2}\,df^{\prime},

where the dimension of the subspace $K$ is fixed. Using ideas from the Karhunen-Loeve (KL) transform [43], it can be shown that the optimal $K$ -dimensional subspace is the span of the top $K$ eigenvectors of the covariance matrix

\boldsymbol{C}_{f}:=\dfrac{1}{2W}\int_{f-W}^{f+W}\boldsymbol{e}_{f^{\prime}}\boldsymbol{e}_{f^{\prime}}^{*}\,df^{\prime}.

The entries of this covariance matrix are

	$\displaystyle\boldsymbol{C}_{f}[m,n]$	$\displaystyle=\dfrac{1}{2W}\int_{f-W}^{f+W}\boldsymbol{e}_{f^{\prime}}[m]\overline{\boldsymbol{e}_{f^{\prime}}[n]}\,df^{\prime}$
		$\displaystyle=\dfrac{1}{2W}\int_{f-W}^{f+W}e^{j2\pi f^{\prime}(m-n)}\,df^{\prime}$
		$\displaystyle=\dfrac{\sin[2\pi W(m-n)]}{2\pi W(m-n)}e^{j2\pi f(m-n)}$
		$\displaystyle=\dfrac{1}{2W}\boldsymbol{e}_{f}[m]\boldsymbol{B}[m,n]\overline{\boldsymbol{e}_{f}[n]},$

where again $\boldsymbol{B}$ is the $N\times N$ prolate matrix. Hence, we can write

\boldsymbol{C}_{f}=\dfrac{1}{2W}\boldsymbol{E}_{f}\boldsymbol{B}\boldsymbol{E}_{f}^{*},

where $\boldsymbol{E}_{f}=\operatorname{diag}(\boldsymbol{e}_{f})\in\mathbb{C}^{N\times N}$ is a unitary matrix which modulates vectors by pointwise multiplying them by the sinusoid $\boldsymbol{e}_{f}$ . Therefore, the eigenvectors of $\boldsymbol{C}_{f}$ are the modulated Slepian basis vectors $\boldsymbol{E}_{f}\boldsymbol{s}_{k}$ for $k\in[N]$ , and the corresponding eigenvalues are $\tfrac{\lambda_{k}}{2W}$ . Hence, we can choose $\mathcal{S}_{f}=\text{span}\{\boldsymbol{E}_{f}\boldsymbol{s}_{0},\ldots,\boldsymbol{E}_{f}\boldsymbol{s}_{K-1}\}$ , i.e., the span of the first $K$ Slepian basis vectors modulated to the frequency $f$ . Since $\boldsymbol{s}_{0},\ldots,\boldsymbol{s}_{K-1}$ are orthonormal vectors, and $\boldsymbol{E}_{f}$ is a unitary matrix, $\boldsymbol{E}_{f}\boldsymbol{s}_{0},\ldots,\boldsymbol{E}_{f}\boldsymbol{s}_{K-1}$ are orthonormal vectors. Hence, $\text{proj}_{\mathcal{S}_{f}}(\boldsymbol{x})=\boldsymbol{E}_{f}\boldsymbol{S}_{K}\boldsymbol{S}_{K}^{*}\boldsymbol{E}_{f}^{*}\boldsymbol{x}$ where $\boldsymbol{S}_{K}=\begin{bmatrix}\boldsymbol{s}_{0}&\ldots&\boldsymbol{s}_{K-1}\end{bmatrix}$ , and thus,

	$\displaystyle\left\\|\text{proj}_{\mathcal{S}_{f}}(\boldsymbol{x})\right\\|_{2}^{2}$	$\displaystyle=\left\\|\boldsymbol{E}_{f}\boldsymbol{S}_{K}\boldsymbol{S}_{K}^{}\boldsymbol{E}_{f}^{}\boldsymbol{x}\right\\|_{2}^{2}$
		$\displaystyle=\left\\|\boldsymbol{S}_{K}^{}\boldsymbol{E}_{f}^{}\boldsymbol{x}\right\\|_{2}^{2}$
		$\displaystyle=\sum_{k=0}^{K-1}\left\|(\boldsymbol{S}_{K}^{}\boldsymbol{E}_{f}^{}\boldsymbol{x})[k]\right\|^{2}$
		$\displaystyle=\sum_{k=0}^{K-1}\left\|\boldsymbol{s}_{k}^{*}\boldsymbol{E}_{f}\boldsymbol{x}\right\|^{2}$
		$\displaystyle=\sum_{k=0}^{K-1}\left\|\sum_{n=0}^{N-1}\boldsymbol{s}_{k}[n]\boldsymbol{x}[n]e^{-j2\pi fn}\right\|^{2}.$

Up to a constant scale factor, this is precisely the multitaper spectral estimate. Hence, we can view the the multitaper spectral estimate $\widehat{S}^{\text{mt}}_{K}(f)=\tfrac{1}{K}\|\boldsymbol{S}_{K}^{*}\boldsymbol{E}_{f}^{*}\boldsymbol{x}\|_{2}^{2}$ as the energy in $\boldsymbol{x}$ after it is projected onto the $K$ -dimensional subspace which best represents the collection of sinusoids $\{\boldsymbol{e}_{f^{\prime}}:f^{\prime}\in[f-W,f+W]\}$ .

2.3 Slepian basis eigenvalues

Before we proceed to our main results, we elaborate on the clustering behavior of the Slepian basis eigenvalues, as they are critical to our analysis in the rest of this paper. For any $\epsilon\in(0,\tfrac{1}{2})$ , slightly fewer than $2NW$ eigenvalues lie in $[1-\epsilon,1)$ , slightly fewer than $N-2NW$ eigenvalues lie in $(0,\epsilon]$ , and very few eigenvalues lie in the so-called “transition region” $(\epsilon,1-\epsilon)$ . In Figure 1, we demonstrate this phenomenon by plotting the first $1000$ Slepian basis eigenvalues for $N=10000$ and $W=\tfrac{1}{100}$ (so $2NW=200$ ). The first $194$ eigenvalues lie in $[0.999,1)$ and the last $9794$ eigenvalues lie in $(0,0.001]$ . Only $12$ eigenvalues lie in $(0.001,0.999)$ .

Refer to caption — Figure 1: A plot of the first $1000$ Slepian basis eigenvalues for $N=10000$ and $W=\tfrac{1}{100}$ . These eigenvalues satisfy $\lambda_{193}\approx 0.9997$ and $\lambda_{206}\approx 0.0003$ . Only $12$ of the $10000$ Slepian basis eigenvalues lie in $(0.001,0.999)$ .

Slepian [3] showed that for any fixed $W\in(0,\tfrac{1}{2})$ and $b\in\mathbb{R}$ ,

\lambda_{\left\lfloor 2NW+(b/\pi)\log N\right\rfloor}\to\dfrac{1}{1+e^{b\pi}}\quad\text{as}\quad N\to\infty.

From this result, it is easy to show that for any fixed $W\in(0,\tfrac{1}{2})$ and $\epsilon\in(0,\tfrac{1}{2})$ ,

\#\{k:\epsilon<\lambda_{k}<1-\epsilon\}\sim\dfrac{2}{\pi^{2}}\log N\log\left(\dfrac{1}{\epsilon}-1\right)\quad\text{as}\quad N\to\infty.

In [44], the authors of this paper derive a non-asymptotic bound on the number of Slepian basis eigenvalues in the transition region. For any $N\in\mathbb{N}$ , $W\in(0,\tfrac{1}{2})$ , and $\epsilon\in(0,\tfrac{1}{2})$ ,

\#\{k:\epsilon<\lambda_{k}<1-\epsilon\}\leq\dfrac{2}{\pi^{2}}\log(100NW+25)\log\left(\dfrac{5}{\epsilon(1-\epsilon)}\right)+7.

This shows that the width of the transition region grows logarithmically with respect to both the time-bandwidth product $2NW$ and the tolerance parameter $\epsilon$ . Note that when $\epsilon$ and $W$ are small, this result is a significant improvement over the prior non-asymptotic bounds in [45], [46], and [47], which scale like $O(\tfrac{\log N}{\epsilon})$ , $O(\tfrac{\log(NW)}{\epsilon})$ , and $O(\log N\log\tfrac{1}{\epsilon})$ respectively. In Section 4, we will exploit the fact that the number of Slepian basis eigenvalues in the transition region grows like $O(\log(NW)\log\tfrac{1}{\epsilon})$ to derive a fast algorithm for evaluating $\widehat{S}^{\text{mt}}_{K}(f)$ at a grid of evenly spaced frequencies.

Furthermore, in [44], the authors of this paper combine the above bound on the width of the transition region with the fact from [45] that $\lambda_{\left\lfloor 2NW\right\rfloor-1}\geq\tfrac{1}{2}\geq\lambda_{\left\lceil 2NW\right\rceil}$ to obtain the following lower bound on the first $\approx 2NW$ eigenvalues

\lambda_{k}\geq 1-10\exp\left(-\dfrac{\left\lfloor 2NW\right\rfloor-k-7}{\tfrac{2}{\pi^{2}}\log(100NW+25)}\right)\quad\text{for}\quad 0\leq k\leq\left\lfloor 2NW\right\rfloor-1.

This shows that as $k$ decreases from $\approx 2NW$ , the quantity $1-\lambda_{k}$ decays exponentially to $0$ . In particular, this means that the first $2NW-O(\log(NW))$ Slepian basis tapers $\boldsymbol{s}_{k}$ have spectral windows $\widetilde{\boldsymbol{s}}_{k}(f)$ which have very little energy outside $[-W,W]$ . However, the spectral windows of the Slepian basis tapers for $k$ near $2NW$ will have a significant amount of energy outside $[-W,W]$ . In Section 3, we show that the statistical properties of the multitaper spectral estimate are significantly improved when using $K=2NW-O(\log(NW))$ tapers instead of the traditional $K=\left\lfloor 2NW\right\rfloor-1$ tapers.

3 Statistical Properties and Spectral Leakage

For a given vector of signal samples $\boldsymbol{x}\in\mathbb{C}^{N}$ , using Thomson’s multitaper method for spectral estimation requires selecting two parameters: the half-bandwidth $W$ of the Slepian basis tapers and the number of tapers $K$ which are used in the multitaper spectral estimate. The selection of these parameters can greatly impact the accuracy of the multitaper spectral estimate.

In some applications, a relatively small number of samples $N$ are taken, and the desired frequency resolution for the spectral estimate is $O(\tfrac{1}{N})$ , i.e., a small multiple of the fundamental Rayleigh resolution limit. In such cases, many practitioners [14, 17, 19, 20, 22, 25, 28, 29] choose the half-bandwidth parameter $W$ such that $2NW$ is between $3$ and $10$ , and then choose the number of Slepian basis tapers to be between $K=\left\lfloor 2NW\right\rfloor$ and $K=\left\lfloor 2NW\right\rfloor-2$ . However, in applications where a large number of samples $N$ are taken, and some loss of resolution is acceptable, choosing a larger half-bandwidth parameter $W$ can result in a more accurate spectral estimate. Furthermore, if the power spectral density $S(f)$ has a high dynamic range (that is $\max_{f}S(f)\gg\min_{f}S(f)$ ), we aim to show that choosing $K=2NW-O(\log(NW)\log\tfrac{1}{\delta})$ tapers for some small $\delta>0$ (instead of $K=2NW-O(1)$ tapers) can provide significantly better protection against spectral leakage.

For all the theorems in this section, we assume that $\boldsymbol{x}\in\mathbb{C}^{N}$ is a vector of samples from a complex Gaussian process whose power spectral density $S(f)$ is bounded and integrable. Note that the analogous results for a real Gaussian process would be similar, but slightly more complicated to state. To state our results, define

M=\max_{f\in\mathbb{R}}S(f),

i.e., the global maximum of the power spectral density, and for each frequency $f\in\mathbb{R}$ we define:

m_{f}=\min_{f^{\prime}\in[f-W,f+W]}S(f^{\prime}),

M_{f}=\max_{f^{\prime}\in[f-W,f+W]}S(f^{\prime}),

A_{f}=\dfrac{1}{2W}\int_{f-W}^{f+W}S(f^{\prime})\,df^{\prime},

R_{f}=\sqrt{\dfrac{1}{2W}\int_{f-W}^{f+W}S(f^{\prime})^{2}\,df^{\prime}},

i.e. the minimum, maximum, average, and root-mean-squared values of the power spectral density over the interval $[f-W,f+W]$ . We also define the quantities

\Sigma_{K}^{(1)}=\dfrac{1}{K}\sum_{k=0}^{K-1}(1-\lambda_{k})

\Sigma_{K}^{(2)}=\sqrt{\dfrac{1}{K}\sum_{k=0}^{K-1}(1-\lambda_{k})^{2}}.

Before we proceed to our results, we make note of the fact that $m_{f}$ , $M_{f}$ , $A_{f}$ , and $R_{f}$ are all “local” properties of the power spectral density, i.e., they depend only on values of $S(f^{\prime})$ for $f^{\prime}\in[f-W,f+W]$ , whereas $M$ is a “global” property. Note that if the power spectral density is “slowly varying” over the interval $[f-W,f+W]$ , then $m_{f}\approx M_{f}\approx A_{f}\approx R_{f}\approx S(f)$ . However, $M$ could be several orders of magnitude larger than $m_{f}$ , $M_{f}$ , $A_{f}$ , and $R_{f}$ if the power spectral density has a high dynamic range.

By using the bound on the Slepian basis eigenvalues from Section 2.3, we can obtain $\lambda_{K-1}\geq 1-\delta$ for some suitably small $\delta>0$ by choosing the number of tapers to be $K=2NW-O(\log(NW)\log\tfrac{1}{\delta})$ . This choice of $K$ guarantees that $0\leq\Sigma^{(1)}_{K}\leq\Sigma^{(2)}_{K}\leq 1-\lambda_{K-1}\leq\delta$ , i.e., $\Sigma^{(1)}_{K}$ , $\Sigma^{(2)}_{K}$ , and $1-\lambda_{K-1}$ are all small, and thus, the global property $M=\max_{f}S(f)$ of the power spectral density will have a minimal impact on the non-asymptotic results below. In other words, using $K=2NW-O(\log(NW)\log\tfrac{1}{\delta})$ tapers mitigates the ability for values of the power spectral density $S(f^{\prime})$ at frequencies $f^{\prime}\not\in[f-W,f+W]$ to impact the estimate $\widehat{S}^{\text{mt}}_{K}(f)$ . However, if $K=2NW-O(1)$ tapers are used, then the quantities $\Sigma^{(1)}_{K}$ , $\Sigma^{(2)}_{K}$ , and $1-\lambda_{K-1}$ could be large enough for the global property $M=\max_{f}S(f)$ of the power spectral density to significantly weaken the non-asymptotic results below. In other words, energy in the power spectral density $S(f^{\prime})$ at frequencies $f^{\prime}\not\in[f-W,f+W]$ can “leak” into the estimate $\widehat{S}^{\text{mt}}_{K}(f)$ .

We begin with a bound on the bias of the multitaper spectral estimate under the additional assumption that the power spectral density is twice differentiable. Note this assumption is only used in Theorem 1.

Theorem 1.

For any frequency $f\in\mathbb{R}$ , if $S(f^{\prime})$ is twice continuously differentiable in $[f-W,f+W]$ , then the bias of the multitaper spectral estimate is bounded by

\operatorname{Bias}\left[\widehat{S}^{\text{mt}}_{K}(f)\right]=\left|\mathbb{E}\widehat{S}^{\text{mt}}_{K}(f)-S(f)\right|\leq\dfrac{M^{\prime\prime}_{f}NW^{3}}{3K}+(M+M_{f})\Sigma_{K}^{(1)},

where

M^{\prime\prime}_{f}=\max_{f^{\prime}\in[f-W,f+W]}|S^{\prime\prime}(f^{\prime})|.

If $K=2NW-O(\log(NW)\log\tfrac{1}{\delta})$ tapers are used for some small $\delta>0$ , then this upper bound is slightly larger than $\tfrac{1}{6}M^{\prime\prime}_{f}W^{2}$ , which is similar to the asymptotic results in [4, 35, 48, 36, 37] which state that the bias is roughly $\tfrac{1}{6}S^{\prime\prime}(f)W^{2}$ . However, if $K=2NW-O(1)$ tapers are used, the term $(M+M_{f})\Sigma_{K}^{(1)}$ could dominate this bound, and the bias could be much larger than the asymptotic result.

If the power spectral density is not twice-differentiable, we can still obtain the following bound on the bias of the multitaper spectral estimate.

Theorem 2.

For any frequency $f\in\mathbb{R}$ , the bias of the multitaper spectral estimate is bounded by

\operatorname{Bias}\left[\widehat{S}^{\text{mt}}_{K}(f)\right]=\left|\mathbb{E}\widehat{S}^{\text{mt}}_{K}(f)-S(f)\right|\leq(M_{f}-m_{f})(1-\Sigma_{K}^{(1)})+M\Sigma_{K}^{(1)}.

If $K=2NW-O(\log(NW)\log\tfrac{1}{\delta})$ tapers are used for some small $\delta>0$ , then this upper bound is slightly larger than $M_{f}-m_{f}$ . This guarantees the bias is small when the power spectral density is “slowly varying” over $[f-W,f+W]$ . However, if $K=2NW-O(1)$ tapers are used, the term $M\Sigma_{K}^{(1)}$ could dominate this bound, and the bias could be much larger than the asymptotic result.

Next, we state our bound on the variance of the multitaper spectral estimate.

Theorem 3.

For any frequency $f\in\mathbb{R}$ , the variance of the multitaper spectral estimate is bounded by

\operatorname{Var}\left[\widehat{S}^{\text{mt}}_{K}(f)\right]\leq\dfrac{1}{K}\left(R_{f}\sqrt{\dfrac{2NW}{K}}+M\Sigma^{(2)}_{K}\right)^{2}.

If $K=2NW-O(\log(NW)\log\tfrac{1}{\delta})$ tapers are used for some small $\delta>0$ , then this upper bound is slightly larger than $\tfrac{1}{K}R_{f}^{2}$ , which is similar to the asymptotic results in [4, 35, 48, 36, 37] which state that the variance is roughly $\tfrac{1}{K}S(f)^{2}$ . However, if $K=2NW-O(1)$ tapers are used, the term $M\Sigma_{K}^{(2)}$ could dominate this bound, and the variance could be much larger than the asymptotic result.

We also note that if the frequencies $f_{1},f_{2}$ are more than $2W$ apart, then the multitaper spectral estimates at those frequencies have a very low covariance.

Theorem 4.

For any frequencies $f_{1},f_{2}\in\mathbb{R}$ such that $2W<|f_{1}-f_{2}|<1-2W$ , the covariance of the multitaper spectral estimates at those frequencies is bounded by

0\leq\operatorname{Cov}\left[\widehat{S}^{\text{mt}}_{K}(f_{1}),\widehat{S}^{\text{mt}}_{K}(f_{2})\right]\leq\left((R_{f_{1}}+R_{f_{2}})\sqrt{\dfrac{2NW}{K}\Sigma^{(1)}_{K}}+M\Sigma^{(1)}_{K}\right)^{2}.

If $K=2NW-O(\log(NW)\log\tfrac{1}{\delta})$ tapers are used for some small $\delta>0$ , then the covariance is guaranteed to be small. However, if $K=2NW-O(1)$ tapers are used, the upper bound is no longer guaranteed to be small, and the covariance could be large.

Finally, we also provide a concentration result for the multitaper spectral estimate.

Theorem 5.

For any frequency $f\in\mathbb{R}$ , the multitaper spectral estimate satisfies the concentration inequalities

\mathbb{P}\left\{\widehat{S}^{\text{mt}}_{K}(f)\geq\beta\mathbb{E}\widehat{S}^{\text{mt}}_{K}(f)\right\}\leq\beta^{-1}\mathrm{e}^{-\kappa_{f}(\beta-1-\ln\beta)}\quad\text{for}\quad\beta>1,

and

\mathbb{P}\left\{\widehat{S}^{\text{mt}}_{K}(f)\leq\beta\mathbb{E}\widehat{S}^{\text{mt}}_{K}(f)\right\}\leq\mathrm{e}^{-\kappa_{f}(\beta-1-\ln\beta)}\quad\text{for}\quad 0<\beta<1,

where the frequency dependent constant $\kappa_{f}$ satisfies

\kappa_{f}\geq\dfrac{K\left(1-\Sigma^{(1)}_{K}\right)M_{f}-2NW(M_{f}-A_{f})}{M_{f}+(M-M_{f})(1-\lambda_{K-1})}.

We note that these are identical to the concentration bounds for a chi-squared random variable with $2\kappa_{f}$ degrees of freedom. If $K=2NW-O(\log(NW)\log\tfrac{1}{\delta})$ tapers are used for some small $\delta>0$ and the power spectral density is “slowly varying” over $[f-W,f+W]$ , then this lower bound on $\kappa_{f}$ is slightly less than $K$ . Hence, $\widehat{S}^{\text{mt}}_{K}(f)$ has a concentration behavior that is similar to a chi-squared random variable with $2K$ degrees of freedom, as the asymptotic results in [4, 34] suggest. However, if $K=2NW-O(1)$ tapers are used, then $\kappa_{f}$ could be much smaller, and thus, the multitaper spectral estimate would have significantly worse concentration about its mean.

The proofs of Theorems 1-5 are given in Appendix A. In Section 5, we perform simulations demonstrating that using $K=2NW-O(1)$ tapers results in a multitaper spectral estimate that is vulnerable to spectral leakage, whereas using $K=2NW-O(\log(NW)\log\tfrac{1}{\delta})$ tapers for a suitably small $\delta>0$ significantly reduces the impact of spectral leakage on the multitaper spectral estimate.

4 Fast Algorithms

Given a vector of $N$ samples $\boldsymbol{x}\in\mathbb{C}^{N}$ , evaluating the multitaper spectral estimate $\widehat{S}^{\text{mt}}_{K}(f)$ at a grid of $L$ evenly spaced frequencies $f\in[L]/L$ (where we assume $L\geq N$ ) can be done in $O(KL\log L)$ operations and using $O(KL)$ memory via $K$ length- $L$ fast Fourier transforms (FFTs). In applications where the number of samples $N$ is small, the number of tapers $K$ used is usually a small constant, and so, the computational requirements are a small constant factor more than that of an FFT. However in many applications, using a large number of tapers is desirable, but the computational requirements make this impractical. As mentioned in Section 1, if the power spectrum $S(f)$ is twice-differentiable, then the MSE of the multitaper spectral estimate is minimized when the bandwidth parameter is $W=O(N^{-1/5})$ and $K=O(N^{4/5})$ tapers are used [36]. For medium to large scale problems, precomputing and storing $O(N^{4/5})$ tapers and/or performing $O(N^{4/5})$ FFTs may be impractical.

In this section, we present an $\epsilon$ -approximation $\widetilde{S}^{\text{mt}}_{K}(f)$ to the multitaper spectral estimate $\widehat{S}^{\text{mt}}_{K}(f)$ which requires $O(L\log L\log(NW)\log\tfrac{1}{\epsilon})$ operations and $O(L\log(NW)\log\tfrac{1}{\epsilon})$ memory. This is faster than the exact multitaper spectral estimation provided the number of tapers satisfies $K\gtrsim\log(NW)\log\tfrac{1}{\epsilon}$ .

To construct this approximation, we first fix a tolerance parameter $\epsilon\in(0,\tfrac{1}{2})$ , and suppose that the number of tapers, $K$ , is chosen such that $\lambda_{K-1}\geq\tfrac{1}{2}$ and $\lambda_{K}\leq 1-\epsilon$ . Note that this is a very mild assumption as it only forces $K$ to be slightly less than $2NW$ . Next, we partition the indices $[N]$ into four sets:

	$\displaystyle\mathcal{I}_{1}$	$\displaystyle=\{k\in[K]:\lambda_{k}\geq 1-\epsilon\}$
	$\displaystyle\mathcal{I}_{2}$	$\displaystyle=\{k\in[K]:\epsilon<\lambda_{k}<1-\epsilon\}$
	$\displaystyle\mathcal{I}_{3}$	$\displaystyle=\{k\in[N]\setminus[K]:\epsilon<\lambda_{k}<1-\epsilon\}$
	$\displaystyle\mathcal{I}_{4}$	$\displaystyle=\{k\in[N]\setminus[K]:\lambda_{k}\leq\epsilon\}$

and define the approximate estimator

\widetilde{S}^{\text{mt}}_{K}(f):=\dfrac{1}{K}\Psi(f)+\dfrac{1}{K}\sum_{k\in\mathcal{I}_{2}}(1-\lambda_{k})\widehat{S}_{k}(f)-\dfrac{1}{K}\sum_{k\in\mathcal{I}_{3}}\lambda_{k}\widehat{S}_{k}(f),

where

\Psi(f):=\sum_{k=0}^{N-1}\lambda_{k}\widehat{S}_{k}(f).

Both $\widehat{S}^{\text{mt}}_{K}(f)$ and $\widetilde{S}^{\text{mt}}_{K}(f)$ are weighted sums of the single taper estimates $\widehat{S}_{k}(f)$ for $k\in[N]$ . Additionally, it can be shown that the weights are similar, i.e., the first $K$ weights are exactly or approximately $\tfrac{1}{K}$ , and the last $N-K$ weights are exactly or approximately $0$ . Hence, it is reasonable to expect that $\widetilde{S}^{\text{mt}}_{K}(f)\approx\widehat{S}^{\text{mt}}_{K}(f)$ . The following theorem shows that is indeed the case.

Theorem 6.

The approximate multitaper spectral estimate $\widetilde{S}^{\text{mt}}_{K}(f)$ defined above satisfies

\left|\widetilde{S}^{\text{mt}}_{K}(f)-\widehat{S}^{\text{mt}}_{K}(f)\right|\leq\dfrac{\epsilon}{K}\|\boldsymbol{x}\|_{2}^{2}\quad\text{for all}\quad f\in\mathbb{R}.

Furthermore, we show in Lemma 8 that $\Psi(f)=\boldsymbol{x}^{*}\boldsymbol{E}_{f}\boldsymbol{B}\boldsymbol{E}_{f}^{*}\boldsymbol{x}$ . This formula doesn’t involve any of the Slepian tapers. By exploiting the fact that the prolate matrix $\boldsymbol{B}$ is Toeplitz, we also show in Lemma 8 that if $L\geq 2N$ , then

\begin{bmatrix}\Psi(\tfrac{0}{L})&\Psi(\tfrac{1}{L})&\cdots&\Psi(\tfrac{L-2}{L})&\Psi(\tfrac{L-1}{L})\end{bmatrix}^{T}=\boldsymbol{F}^{-1}\left(\boldsymbol{b}\circ\boldsymbol{F}\left|\boldsymbol{F}\boldsymbol{Z}\boldsymbol{x}\right|^{2}\right),

where $\boldsymbol{Z}\in\mathbb{R}^{L\times N}$ is a matrix which zero-pads length- $N$ vectors to length- $L$ , $\boldsymbol{F}\in\mathbb{C}^{L\times L}$ is a length- $L$ FFT matrix, $\boldsymbol{b}\in\mathbb{R}^{L}$ is the first column of the matrix formed by extending the prolate matrix $\boldsymbol{B}$ to an $L\times L$ circulant matrix, $|\cdot|^{2}$ denotes the pointwise magnitude-squared, and $\circ$ denotes a pointwise multiplication. Hence, $\Psi(f)$ can be evaluated at a grid of $L$ evenly spaced frequencies $f\in[L]/L$ in $O(L\log L)$ operations via three length- $L$ FFTs/inverse FFTs. Evaluating the other $\#(\mathcal{I}_{2}\cup\mathcal{I}_{3})=O(\log(NW)\log\tfrac{1}{\epsilon})$ terms in the expression for $\widetilde{S}^{\text{mt}}_{K}(f)$ at the $L$ grid frequencies can be done in $O(L\log L\log(NW)\log\tfrac{1}{\epsilon})$ operations via $\#(\mathcal{I}_{2}\cup\mathcal{I}_{3})$ length- $L$ FFTs. Using these results, we establish the following theorem which states how quickly $\widetilde{S}^{\text{mt}}_{K}(f)$ can be evaluated at the grid frequencies.

Theorem 7.

For any vector of samples $\boldsymbol{x}\in\mathbb{C}^{N}$ and any number of grid frequencies $L\geq N$ , the approximate multitaper spectral estimate $\widetilde{S}^{\text{mt}}_{K}(f)$ can be evaluated at the $L$ grid frequencies $f\in[L]/L$ in $O(L\log L\log(NW)\log\tfrac{1}{\epsilon})$ operations and using $O(L\log(NW)\log\tfrac{1}{\epsilon})$ memory.

Note if $N\leq L<2N$ , we can apply the method briefly described above to evaluate $\Psi(f)$ at $f\in[2L]/2L$ , and then downsample the result. The proofs of Theorems 6 and 7 are given in Appendix B. In Section 5, we perform simulations comparing the time needed to evaluate $\widehat{S}^{\text{mt}}_{K}(f)$ and $\widetilde{S}^{\text{mt}}_{K}(f)$ at a grid of frequencies.

5 Simulations

In this section, we show simulations to demonstrate three observations. (1) Using $K=2NW-O(\log(NW))$ tapers instead of the traditional choice of $K=\left\lfloor 2NW\right\rfloor-1$ tapers significantly reduce the effects of spectral leakage. (2) Using a larger bandwidth $W$ , and thus, more tapers can produce a more robust spectral estimate. (3) As the number of samples $N$ and the number of tapers $K$ grows, our approximation $\widetilde{S}^{\text{mt}}_{K}(f)$ becomes significantly faster to use than the exact multitaper spectral estimate $\widehat{S}^{\text{mt}}_{K}(f)$ .

5.1 Spectral leakage

First, we demonstrate that choosing $K=2NW-O(\log(NW))$ tapers instead of the traditional choice of $K=\left\lfloor 2NW\right\rfloor-1$ tapers significantly reduces the effects of spectral leakage. We fix a signal length of $N=2000$ , a bandwidth parameter of $W=\tfrac{1}{100}$ (so $2NW=40$ ) and consider four choices for the number of tapers: $K=39$ , $36$ , $32$ , and $29$ . Note that $K=39=\left\lfloor 2NW\right\rfloor-1$ is the traditional choice as to how many tapers to use, while $36$ , $32$ , and $29$ are the largest values of $K$ such that $\lambda_{K-1}$ is at least $1-10^{-3}$ , $1-10^{-6}$ , and $1-10^{-9}$ respectively.

In Figure 2, we show three plots of the spectral window

\psi(f)=\dfrac{1}{K}\sum_{k=0}^{K-1}\left|\sum_{n=0}^{N-1}\boldsymbol{s}_{k}[n]e^{-j2\pi fn}\right|^{2}

of the multitaper spectral estimate for each of those values of $K$ . At the top of Figure 2, we plot $\psi(f)$ over the entire range $[-\tfrac{1}{2},\tfrac{1}{2}]$ using a logarithmic scale. The lines outside appear thick due to the highly oscillatory behavior of $\psi(f)$ . This can be better seen in the middle of Figure 2, where we plot $\psi(f)$ over $[-2W,2W]$ using a logarithmic scale. The behavior of $\psi(f)$ inside $[-W,W]$ can be better seen at the bottom of Figure 2, where we plot $\psi(f)$ over $[-2W,2W]$ using a linear scale.

All four spectral windows have similar behavior in that $\psi(f)$ is small outside $[-W,W]$ and large near $0$ . However, outside of $[-W,W]$ the spectral windows using $K=\left\lfloor 2NW\right\rfloor-O(\log(NW))$ tapers are multiple orders of magnitude smaller than the spectral window using $K=\left\lfloor 2NW\right\rfloor-1$ tapers. Hence, the amount of spectral leakage can be reduced by multiple orders of magnitude by trimming the number of tapers used from $K=\left\lfloor 2NW\right\rfloor-1$ to $K=2NW-O(\log(NW))$ .

We further demonstrate the importance of using $K=2NW-O(\log(NW))$ tapers to reduce spectral leakage by showing a signal detection example. We generate a vector $\boldsymbol{x}\in\mathbb{C}^{N}$ of $N=2000$ samples of a Gaussian random process with a power spectral density function of

S(f)=\begin{cases}10^{3}&\text{if}\ f\in[0.18,0.22]\\ 10^{9}&\text{if}\ f\in[0.28,0.32]\\ 10^{2}&\text{if}\ f\in[0.38,0.42]\\ 10^{1}&\text{if}\ f\in[0.78,0.82]\\ 10^{0}&\text{else}\end{cases}.

This simulates an antenna receiving signals from four narrowband sources with some background noise. Note that one source is significantly stronger than the other three sources. In Figure 3, we plot:

•

the periodogram of $\boldsymbol{x}$ ,
•

the multitaper spectral estimate of $\boldsymbol{x}$ with $W=\tfrac{1}{100}$ and $K=\left\lfloor 2NW\right\rfloor-1=39$ tapers,
•

the multitaper spectral estimate of $\boldsymbol{x}$ with $W=\tfrac{1}{100}$ and $K=29$ tapers (chosen so $\lambda_{K-1}\geq 1-10^{-9}$ ).

We note that all three spectral estimates yield large values in the frequency band $[0.28,0.32]$ . However, in the periodogram and the multitaper spectral estimate with $K=39$ tapers, the energy in the frequency band $[0.28,0.32]$ “leaks” into the frequency bands occupied by the smaller three sources. As a result, the smaller three sources are hard to detect using the periodogram or the multitaper spectral estimate with $K=39$ tapers. However, all four sources are clearly visible when looking at the multitaper spectral estimate with $K=29$ tapers. For frequencies $f$ not within $W$ of the edges of the frequency band, the multitaper spectral estimate is within a small constant factor of the true power spectral density.

5.2 Comparison of spectral estimation methods

Next, we demonstrate a few key points about selecting the bandwidth parameter $W$ and the number of tapers $K$ used in the multitaper spectral estimate. We compare the following eight methods to estimate the power spectrum of a Gaussian random process from a vector $\boldsymbol{x}\in\mathbb{C}^{N}$ of $N=2^{18}=262144$ samples:

1.

The classic periodogram
2.

A tapered periodogram using a single DPSS taper $\boldsymbol{s}_{0}$ with $2NW=8$ (chosen such that $\lambda_{0}\geq 1-10^{-9}$ )
3.

The exact multitaper spectral estimate $\widehat{S}^{\text{mt}}_{K}(f)$ with a small bandwidth parameter of $W=1.25\times 10^{-4}$ and $K=\left\lfloor 2NW\right\rfloor-1=64$ tapers
4.

The exact multitaper spectral estimate $\widehat{S}^{\text{mt}}_{K}(f)$ with a small bandwidth parameter of $W=1.25\times 10^{-4}$ and $K=53$ tapers (chosen such that $\lambda_{K-1}\geq 1-10^{-9}>\lambda_{K}$ )
5.

The approximate multitaper spectral estimate $\widetilde{S}^{\text{mt}}_{K}(f)$ with a larger bandwidth parameter of $W=2.0\times 10^{-3}$ , $K=\left\lfloor 2NW\right\rfloor-1=1047$ tapers, and a tolerance parameter of $\epsilon=10^{-9}$
6.

The approximate multitaper spectral estimate $\widetilde{S}^{\text{mt}}_{K}(f)$ with a larger bandwidth parameter of $W=2.0\times 10^{-3}$ , $K=1031$ tapers (chosen such that $\lambda_{K-1}\geq 1-10^{-9}>\lambda_{K}$ ), and a tolerance parameter of $\epsilon=10^{-9}$
7.

The exact multitaper spectral estimate with the adaptive weighting scheme suggested by Thomson[4] with a small bandwidth parameter of $W=1.25\times 10^{-4}$ and $K=\left\lfloor 2NW\right\rfloor-1=64$ tapers
8.

The exact multitaper spectral estimate with the adaptive weighting scheme suggested by Thomson[4] with a larger bandwidth parameter of $W=2.0\times 10^{-3}$ and $K=\left\lfloor 2NW\right\rfloor-1=1047$ tapers

The adaptive weighting scheme computes the single taper periodograms $\widehat{S}_{k}(f)$ for $k\in[K]$ , and then forms a weighted estimate

\widehat{S}^{\text{ad}}_{K}(f)=\dfrac{\sum_{k=0}^{K-1}\alpha_{k}(f)\widehat{S}_{k}(f)}{\sum_{k=0}^{K-1}\alpha_{k}(f)}

where the frequency dependent weights $\alpha_{k}(f)$ satisfy

\alpha_{k}(f)=\dfrac{\lambda_{k}\widehat{S}^{\text{ad}}_{K}(f)^{2}}{\left(\lambda_{k}\widehat{S}^{\text{ad}}_{K}(f)+(1-\lambda_{k})\sigma^{2}\right)^{2}}

where $\sigma^{2}=\tfrac{1}{N}\|\boldsymbol{x}\|_{2}^{2}$ . Of course, solving for the weights directly is difficult, so this method requires initializing the weights and alternating between updating the estimate $\widehat{S}^{\text{ad}}_{K}(f)$ and updating the weights $\alpha_{k}(f)$ . This weighting procedure is designed to keep all the weights large at frequencies where $S(f)$ is large and reduce the weights of the last few tapers at frequencies where $S(f)$ is small. Effectively, this allows the spectrum to be estimated with more tapers at frequencies where $S(f)$ is large while simultaneously reducing the spectral leakage from the last few tapers at frequencies where $S(f)$ is small. The cost is the increased computation time due to setting the weights iteratively. For more details regarding this adaptive scheme, see [8].

In Figure 4, we plot the power spectrum and the eight estimates for a single realization $\boldsymbol{x}\in\mathbb{C}^{N}$ of the Gaussian random process. Additionally, for $1000$ realizations $\boldsymbol{x}_{i}\in\mathbb{C}^{N}$ , $i=1,\ldots,1000$ of the Gaussian random process, we compute a spectral estimate using each of the above eight methods. In Figure 5, we plot the empirical mean logarithmic deviation in dB, i.e.,

\dfrac{1}{1000}\sum_{i=1}^{1000}\left|10\log_{10}\dfrac{\widehat{S}[\boldsymbol{x}_{i}](f)}{S(f)}\right|

for each of the eight methods. In Table 1, we list the average time needed to precompute the DPSS tapers, the average time to compute the spectral estimate after the tapers are computed, and the average of the empirical mean logarithmic deviation in dB across the frequency spectrum.

We make the following observations:

•

The periodogram and the single taper periodogram (methods 1 and 2) are too noisy to be useful spectral estimates.
•

Methods 3, 4, and 7 yield a noticeably noisier spectral estimate than methods 5, 6, and 8. This is due to the fact that methods 5, 6, and 8 use a larger bandwidth parameter and more tapers.
•

The spectral estimates obtained with methods 1, 3 and 5 suffer from spectral leakage, i.e., the error is large at frequencies $f$ where $S(f)$ is small, as can be seen in Figure 5. This is due to the fact that they use $K=\left\lfloor 2NW\right\rfloor-1$ tapers, and thus, include tapered periodograms $\widehat{S}_{k}(f)$ for which $\lambda_{k}$ is not very close to $1$ .
•

Methods 4 and 6 avoid using tapered periodograms $\widehat{S}_{k}(f)$ for which $\lambda_{k}<1-10^{-9}$ and methods 7 and 8 use these tapered periodograms but assign a low weight to them at frequencies where $S(f)$ is small. Hence, methods 4, 6, 7, and 8 are able to mitigate the spectral leakage phenomenon.
•

Methods 5 and 6 are slightly faster than methods 3 and 4 due to the fact that our approximate multitaper spectral estimate only needs to compute $\#\{k:\epsilon<\lambda_{k}<1-\epsilon\}=36$ tapers and $36$ tapered periodograms.
•

Method 7 takes noticeably longer than methods 3 and 4, and method 8 takes considerably longer than methods 5 and 6. This is because the iterative method for computing the adaptive weights requires many iterations to converge when the underlying power spectral density has a high dynamic range.
•

Methods 6 and 8 exhibit very similar performance. This is to be expected, as using a weighted average of $1047$ tapered periodograms is similar to using the unweighted average of the first $1031$ tapered periodograms. The empirical mean logarithmic deviation is larger at frequencies where $S(f)$ is rapidly varying and smaller at frequencies where $S(f)$ is slowly varying. This is to be expected as the local bias (caused due to the smoothing effect of the tapers) dominates the variance at these frequencies.

Method	Avg MLD (dB)	Avg Precomputation Time (sec)	Avg Run Time (sec)
1	5.83030	N/A	0.0111
2	4.41177	0.1144	0.0111
3	1.48562	4.3893	0.3771
4	0.47836	3.6274	0.3093
5	1.56164	3.3303	0.2561
6	0.12534	3.3594	0.2533
7	0.45080	4.3721	1.4566
8	0.12532	78.9168	17.2709

Table 1: Table of mean logarithmic deviations (averaged across entire frequency spectrum), precomputation times, and computation times for each of the eight spectral estimation methods.

From this, we can draw several conclusions. First, by slightly trimming the number of tapers from $K=\left\lfloor 2NW\right\rfloor-1$ to $2NW-O(\log(NW))$ , one can significantly mitigate the amount of spectral leakage in the spectral estimate. Second, using a larger bandwidth parameter and averaging over more tapered periodograms will result in a less noisy spectral estimate. Third, our fast approximate multitaper spectral estimate can yield a more accurate estimate in the same amount of time as an exact multitaper spectral estimate because our fast approximation needs to compute significantly fewer tapers and tapered periodograms. Fourth, a multitaper spectral estimate with $\left\lfloor 2NW\right\rfloor-1$ tapers and adaptive weights can yield a slightly lower error than a multitaper spectral estimate with $2NW-O(\log(NW))$ tapers and fixed weights, but as $2NW$ increases, this effect becomes minimal and not worth the increased computational cost.

5.3 Fast multitaper spectral estimation

Finally, we demonstrate that the runtime for computing the approximate multitaper spectral estimate $\widetilde{S}^{\text{mt}}_{K}(f)$ scales roughly linearly with the number of samples. We vary the signal length $N$ over several values between $2^{10}$ and $2^{20}$ . For each value of $N$ , we set the bandwidth parameter to be $W=0.08\cdot N^{-1/5}$ as this is similar to what is asymptotically optimal for a twice differentiable spectrum [36]. We then choose a number of tapers $K$ such that $\lambda_{K-1}\geq 1-10^{-3}>\lambda_{K}$ , and then compute the approximate multitaper spectral estimate $\widetilde{S}^{\text{mt}}_{K}(f)$ at $f\in[N]/N$ for each of the tolerances $\epsilon=10^{-4}$ , $10^{-8}$ , and $10^{-12}$ . If $N\leq 2^{17}$ , we also compute the exact multitaper spectral estimate $\widehat{S}^{\text{mt}}_{K}(f)$ at $f\in[N]/N$ . For $N>2^{17}$ , we did not evaluate the exact multitaper spectral estimate due to the excessive computational requirements.

We split the work needed to produce the spectral estimate into a precomputation stage and a computation stage. The precomputation stage involves computing the Slepian tapers which are required for the spectral estimate. In applications where a multitaper spectral estimate needs to be computed for several signals (using the same parameters $N,W,K$ ), computing the tapers only needs to be done once. The exact multitaper spectral estimate $\widehat{S}^{\text{mt}}_{K}(f)$ requires computing $\boldsymbol{s}_{k}$ for $k\in[K]$ , while the approximate multitaper spectral estimate $\widetilde{S}^{\text{mt}}_{K}(f)$ requires computing $\boldsymbol{s}_{k}$ for $k\in\mathcal{I}_{2}\cup\mathcal{I}_{3}=\{k:\epsilon<\lambda_{k}<1-\epsilon\}$ . The computation stage involves evaluating $\widehat{S}^{\text{mt}}_{K}(f)$ or $\widetilde{S}^{\text{mt}}_{K}(f)$ at $f\in[N]/N$ with the required tapers $\boldsymbol{s}_{k}$ already computed.

In Figures 6 and 7, we plot the precomputation and computation time respectively versus the signal length $N$ for the exact multitaper spectral estimate $\widehat{S}^{\text{mt}}_{K}(f)$ as well as the approximate multitaper spectral estimate $\widetilde{S}^{\text{mt}}_{K}(f)$ for $\epsilon=10^{-4}$ , $10^{-8}$ , and $10^{-12}$ . The times were averaged over $100$ trials of the procedure described above. The plots are on a log-log scale. We note that the precomputation and computation times for the approximate multitaper spectral estimate grow roughly linearly with $N$ while the precomputation and computation times for the exact multitaper spectral estimate grow significantly faster. Also, the precomputation and computation times for the approximate multitaper spectral estimate with the very small tolerance $\epsilon=10^{-12}$ are only two to three times more than those for the larger tolerance $\epsilon=10^{-4}$ .

6 Conclusions

Thomson’s multitaper method is a widely-used tool for spectral estimation. In this paper, we have presented a linear algebra based perspective of Thomson’s multitaper method. Specifically, the multitaper spectral estimate $\widehat{S}^{\text{mt}}_{K}(f)$ can be viewed as the $2$ -norm energy of the projection of the signal $\boldsymbol{x}\in\mathbb{C}^{N}$ onto a $K$ -dimensional subspace $\mathcal{S}_{f}\subset\mathbb{C}^{N}$ . The subspace $\mathcal{S}_{f}$ is chosen to be the optimal $K$ -dimensional subspace for representing sinusoids with frequencies between $f-W$ and $f+W$ .

We have also provided non-asymptotic bounds on the bias, variance, covariance, and probability tails of the multitaper spectral estimate. In particular, if the multitaper spectral estimate uses $K=2NW-O(\log(NW))$ tapers, the effects of spectral leakage are mitigated, and our non-asymptotic bounds are very similar to existing asymptotic results. However, when the traditional choice of $K=\left\lfloor 2NW\right\rfloor-1$ tapers is used, the non-asymptotic bounds are significantly worse, and our simulations show the multitaper spectral estimate suffers severely from spectral leakage.

Finally, we have derived an $\epsilon$ -approximation to the multitaper spectral estimate which can be evaluated at a grid of $L$ evenly spaced frequencies in $O(L\log L\log(NW)\log\tfrac{1}{\epsilon})$ operations. Since evaluating the exact multitaper spectral estimate at a grid of $L$ evenly spaced frequencies requires $O(KL\log L)$ operations, our $\epsilon$ -approximation is faster when $K\gtrsim\log(NW)\log\tfrac{1}{\epsilon}$ tapers are required. In problems involving a large number of samples, minimizing the mean-squared error of the multitaper spectral estimate requires using a large number of tapers, and thus, our $\epsilon$ -approximation can drastically reduce the required computation.

References

[1] S. Karnik, J. Romberg, and M. A. Davenport. Fast multitaper spectral estimation. In 2019 13th International conference on Sampling Theory and Applications (SampTA), pages 1–4, Bordeaux, France, July 2019. IEEE.
[2] J. H. McClellan, R. W. Schafer, and M. A. Yoder. Signal Processing First. Pearson, Upper Saddle River, NJ, 2003.
[3] D. Slepian. Prolate spheroidal wave functions, Fourier analysis, and uncertainty. V – The discrete case. Bell Sys. Tech. J., 57(5):1371–1430, 1978.
[4] D. J. Thomson. Spectrum estimation and harmonic analysis. Proc. IEEE, 70(9):1055–1096, September 1982.
[5] S. Haykin. Cognitive radio: brain-empowered wireless communications. IEEE J. Sel. Areas Commun., 23(2):201–220, February 2005.
[6] B. Farhang-Boroujeny. Filter Bank Spectrum Sensing for Cognitive Radios. IEEE Trans. Signal Process., 56(5):1801–1811, May 2008.
[7] B. Farhang-Boroujeny and R. Kempter. Multicarrier communication techniques for spectrum sensing and communication in cognitive radios. IEEE Communications Magazine, 46(4):80–85, April 2008.
[8] S. Haykin, D. J. Thomson, and J. H. Reed. Spectrum sensing for cognitive radio. Proc. IEEE, 97(5):849–877, May 2009.
[9] E. Axell, G. Leus, E. Larsson, and H. Poor. Spectrum sensing for cognitive radio : State-of-the-art and recent advances. IEEE Signal Processing Magazine, 29(3):101–116, May 2012.
[10] K. N. Hamdy, M. Ali, and A. H. Tewfik. Low bit rate high quality audio coding with combined harmonic and wavelet representations. In 1996 IEEE International Conference on Acoustics, Speech, and Signal Processing (ICASSP) Conference Proceedings, volume 2, pages 1045–1048, Atlanta, GA, USA, 1996. IEEE.
[11] T. Painter and A. Spanias. Perceptual coding of digital audio. Proc. IEEE, 88(4):451–515, April 2000.
[12] A. Delorme and S. Makeig. EEGLAB: an open source toolbox for analysis of single-trial EEG dynamics including independent component analysis. J. Neuroscience Methods, 134(1):9–21, March 2004.
[13] A. Delorme, T. Sejnowski, and S. Makeig. Enhanced detection of artifacts in EEG data using higher-order statistics and independent component analysis. NeuroImage, 34(4):1443–1449, February 2007.
[14] R. R. Llinas, U. Ribary, D. Jeanmonod, E. Kronberg, and P. P. Mitra. Thalamocortical dysrhythmia: A neurological and neuropsychiatric syndrome characterized by magnetoencephalography. Proc. Nat. Acad. Sci., 96(26):15222–15227, December 1999.
[15] B. Pesaran, J. S. Pezaris, M. Sahani, P. P. Mitra, and R. A. Andersen. Temporal structure in neuronal activity during working memory in macaque parietal cortex. Nature Neuroscience, 5(8):805–811, August 2002.
[16] J. Csicsvari, B. Jamieson, K. D. Wise, and G. Buzsáki. Mechanisms of gamma oscillations in the hippocampus of the behaving rat. Neuron, 37(2):311–322, January 2003.
[17] P. P. Mitra and B. Pesaran. Analysis of dynamic brain imaging data. Biophysical Journal, 76(2):691–708, February 1999.
[18] M. W. Jones and M. A. Wilson. Theta rhythms coordinate hippocampal–prefrontal interactions in a spatial memory task. PLoS Biology, 3(12):e402, November 2005.
[19] G. Bond, W. Showers, M. Cheseby, R. Lotti, P. Almasi, P. de Menocal, P. Priore, H. Cullen, I. Hajdas, and G. Bonani. A pervasive millennial-scale cycle in north atlantic holocene and glacial climates. Science, 278(5341):1257–1266, November 1997.
[20] M. Ghil, M. R. Allen, M. D. Dettinger, K. Ide, D. Kondrashov, M. E. Mann, A. W. Robertson, A. Saunders, Y. Tian, F. Varadi, and P. Yiou. Advanced spectral methods for climatic time series. Rev. Geophys., 40(1):3–1–3–41, 2002.
[21] R. Vautard and M. Ghil. Singular spectrum analysis in nonlinear dynamics, with applications to paleoclimatic time series. Physica D: Nonlinear Phenomena, 35(3):395–424, May 1989.
[22] M. E. Mann and J. M. Lees. Robust estimation of background noise and signal detection in climatic time series. Climatic Change, 33(3):409–445, July 1996.
[23] S. Minobe. A 50-70 year climatic oscillation over the north pacific and north america. Geophys. Res. Lett., 24(6):683–686, March 1997.
[24] J. Jouzel, N. I. Barkov, J. M. Barnola, M. Bender, J. Chappellaz, C. Genthon, V. M. Kotlyakov, V. Lipenkov, C. Lorius, J. R. Petit, D. Raynaud, G. Raisbeck, C. Ritz, T. Sowers, M. Stievenard, F. Yiou, and P. Yiou. Extending the Vostok ice-core record of palaeoclimate to the penultimate glacial period. Nature, 364(6436):407–412, July 1993.
[25] D. J. Thomson. Time series analysis of holocene climate data. Philos. Trans. Roy. Soc. London. Ser. A, Math. and Physical Sci., 330(1615):601–616, 1990.
[26] D. Brunton, A. Rodrigo, and E. Marks. Ecstatic display calls of the Adélie penguin honestly predict male condition and breeding success. Behaviour, 147(2):165–184, 2010.
[27] S. Saar and P. P. Mitra. A technique for characterizing the development of rhythms in bird song. PLoS ONE, 3(1):e1461, January 2008.
[28] M. Hansson-Sandsten, M. Tarka, J. Caissy-Martineau, B. Hansson, and D. Hasselquist. SVD-based classification of bird singing in different time-frequency domains using multitapers. In 2011 19th European Signal Processing Conference, pages 966–970, August 2011. ISSN: 2076-1465.
[29] A. Leonardo and M. Konishi. Decrystallization of adult birdsong by perturbation of auditory feedback. Nature, 399(6735):466–470, June 1999.
[30] O. Tchernichovski, F. Nottebohm, C. E. Ho, B. Pesaran, and P. P. Mitra. A procedure for an automated measurement of song similarity. Animal Behaviour, 59(6):1167–1176, June 2000.
[31] M. A. Wieczorek. Gravity and topography of the terrestrial planets. In G. Schubert, editor, Treatise on Geophysics, pages 165–206. Elsevier, Amsterdam, 2nd edition, January 2007.
[32] S. G. Claudepierre, S. R. Elkington, and M. Wiltberger. Solar wind driving of magnetospheric ULF waves: Pulsations driven by velocity shear at the magnetopause: ULF waves driven by velocity shear at the magnetopause. J. Geophys. Res.: Space Phys., 113(A5), May 2008.
[33] B. Allen and J. D. Romano. Detecting a stochastic background of gravitational radiation: Signal processing strategies and sensitivities. Phys. Rev. D, 59(10):102001, March 1999.
[34] D. B. Percival and A. T. Walden. Spectral analysis for physical applications: multitaper and conventional univariate techniques. Cambridge University Press, Cambridge ; New York, NY, USA, 1993.
[35] A. T. Walden. A unified view of multitaper multivariate spectral estimation. Biometrika, 87(4):767–788, December 2000.
[36] L. D. Abreu and J. L. Romero. MSE estimates for multitaper spectral estimation and off-grid compressive sensing. IEEE Trans. Inf. Theory, 63(12):7770–7776, December 2017.
[37] C. L. Haley and M. Anitescu. Optimal bandwidth for multitaper spectrum estimation. IEEE Sig. Proc. Lett., 24(11):1696–1700, November 2017.
[38] G. G. Stokes. On a method of detecting inequalities of unknown periods in a series of observations. Proc. Roy. Soc. London, 29:122–123, 1879.
[39] A. Schuster. On the investigation of hidden periodicities with application to a supposed 26 day period of meteorological phenomena. Terr. Magnetism, 3(1):13–41, 1898.
[40] F. J. Harris. On the use of windows for harmonic analysis with the discrete Fourier transform. Proc. IEEE, 66(1):51–83, 1978.
[41] J. M. Varah. The prolate matrix. Linear Algebra and its Applications, 187:269–278, July 1993.
[42] A. W. Bojanczyk, R. P. Brent, F. R. De Hoog, and D. R. Sweet. On the stability of the Bareiss and related Toeplitz factorization algorithms. SIAM J. Matrix Anal. Appl., 16(1):40–57, January 1995.
[43] H. Stark and J. W. Woods. Probability, random processes, and estimation theory for engineers. Prentice-Hall, Englewood Cliffs, N.J, 1986.
[44] S. Karnik, J. Romberg, and M. A. Davenport. Improved bounds for the eigenvalues of prolate spheroidal wave functions and discrete prolate spheroidal sequences. arXiv:2006.00427 [math], September 2020. arXiv: 2006.00427.
[45] Z. Zhu and M. B. Wakin. Approximating sampled sinusoids and multiband signals using multiband modulated DPSS dictionaries. J. Fourier Anal. Appl., 23(6):1263–1310, December 2017.
[46] M. Boulsane, N. Bourguiba, and A. Karoui. Discrete prolate spheroidal wave functions: Further spectral analysis and some related applications. Journal of Scientific Computing, 82(3):1–19, February 2020.
[47] S. Karnik, Z. Zhu, M. B. Wakin, J. Romberg, and M. A. Davenport. The fast Slepian transform. Applied and Computational Harmonic Analysis, 46(3):624–652, May 2019.
[48] K. S. Lii and M. Rosenblatt. Prolate spheroidal spectral estimates. Statist. & Probability. Lett., 78(11):1339–1348, August 2008.
[49] P. H. M. Janssen and P. Stoica. On the expectation of the product of four matrix-valued Gaussian random variables. IEEE Trans. Autom. Control, 33(9):867–870, September 1988.
[50] S. Janson. Tail bounds for sums of geometric and exponential variables. Statist. & Probability. Lett., 135:1–6, April 2018.
[51] D. L. Hanson and F. T. Wright. A bound on tail probabilities for quadratic forms in independent random variables. Ann. Math. Statist., 42(3):1079–1083, June 1971.
[52] F. T. Wright. A bound on tail probabilities for quadratic forms in independent random variables whose distributions are not necessarily symmetric. Ann. Probability, 1(6):1068–1070, December 1973.
[53] D. M. Gruenbacher and D. R. Hummels. A simple algorithm for generating discrete prolate spheroidal sequences. IEEE Trans. Signal Process., 42(11):3276–3278, November 1994.

Appendix A Proof of Results in Section 3

A.1 Norms of Gaussian random variables

In this subsection, we develop results on the $2$ -norm of linear combinations of Gaussian random variables, which will be critical to our proofs of the theorems in Section 3.

First, we state a result from [49] regarding the product of four jointly complex Gaussian vectors.

Lemma 1.

[49] Suppose $\boldsymbol{a},\boldsymbol{b},\boldsymbol{c},\boldsymbol{d}\in\mathbb{C}^{K}$ have a joint complex Gaussian distribution. Then

\mathbb{E}[\boldsymbol{a}^{*}\boldsymbol{b}\boldsymbol{c}^{*}\boldsymbol{d}]=\mathbb{E}[\boldsymbol{a}^{*}\boldsymbol{b}]\mathbb{E}[\boldsymbol{c}^{*}\boldsymbol{d}]+\mathbb{E}[\boldsymbol{c}^{*}\otimes\boldsymbol{a}^{*}]\mathbb{E}[\boldsymbol{d}\otimes\boldsymbol{b}]+\mathbb{E}[\boldsymbol{a}^{*}\mathbb{E}[\boldsymbol{b}\boldsymbol{c}^{*}]\boldsymbol{d}]+2\mathbb{E}[\boldsymbol{a}^{*}]\mathbb{E}[\boldsymbol{b}]\mathbb{E}[\boldsymbol{c}^{*}]\mathbb{E}[\boldsymbol{d}],

where $\otimes$ denotes the Kronecker product.

It should be noted that [49] proves a more general result for the product of four joint complex Gaussian matrices, but this is a bit harder to state. So we give the result for vectors.

We now prove a result regarding the norm-squared of linear combinations of complex Gaussians.

Lemma 2.

Let $\boldsymbol{x}\sim\mathcal{CN}(\boldsymbol{0},\boldsymbol{R})$ for some positive semidefinite $\boldsymbol{R}\in\mathbb{C}^{N\times N}$ and let $\boldsymbol{U},\boldsymbol{V}\in\mathbb{C}^{K\times N}$ . Then, we have:

\mathbb{E}\left[\left\|\boldsymbol{U}\boldsymbol{x}\right\|_{2}^{2}\right]=\operatorname{tr}[\boldsymbol{U}\boldsymbol{R}\boldsymbol{U}^{*}]\quad\text{and}\quad\operatorname{Cov}\left[\left\|\boldsymbol{U}\boldsymbol{x}\right\|_{2}^{2},\left\|\boldsymbol{V}\boldsymbol{x}\right\|_{2}^{2}\right]=\|\boldsymbol{U}\boldsymbol{R}\boldsymbol{V}^{*}\|_{F}^{2}.

Proof.

The expectation of $\|\boldsymbol{U}\boldsymbol{x}\|_{2}^{2}$ can be computed as follows

\mathbb{E}\left[\left\|\boldsymbol{U}\boldsymbol{x}\right\|_{2}^{2}\right]=\mathbb{E}\left[\operatorname{tr}\left[\boldsymbol{U}\boldsymbol{x}\boldsymbol{x}^{*}\boldsymbol{U}^{*}\right]\right]=\operatorname{tr}\left[\mathbb{E}\left[\boldsymbol{U}\boldsymbol{x}\boldsymbol{x}^{*}\boldsymbol{U}^{*}\right]\right]=\operatorname{tr}\left[\boldsymbol{U}\mathbb{E}\left[\boldsymbol{x}\boldsymbol{x}^{*}\right]\boldsymbol{U}^{*}\right]=\operatorname{tr}\left[\boldsymbol{U}\boldsymbol{R}\boldsymbol{U}^{*}\right].

Since the entries of $\boldsymbol{x}$ are jointly Gaussian with mean $0$ , the entries of $\boldsymbol{U}\boldsymbol{x}$ and $\boldsymbol{V}\boldsymbol{x}$ are also jointly Gaussian with mean $0$ . Then, by applying Lemma 1 (with $\boldsymbol{a}=\boldsymbol{b}=\boldsymbol{U}\boldsymbol{x}$ and $\boldsymbol{c}=\boldsymbol{d}=\boldsymbol{V}\boldsymbol{x}$ ), we find that

	$\displaystyle\operatorname{Cov}\left[\left\\|\boldsymbol{U}\boldsymbol{x}\right\\|_{2}^{2},\left\\|\boldsymbol{V}\boldsymbol{x}\right\\|_{2}^{2}\right]$	$\displaystyle=\mathbb{E}\left[\left\\|\boldsymbol{U}\boldsymbol{x}\right\\|_{2}^{2}\cdot\left\\|\boldsymbol{V}\boldsymbol{x}\right\\|_{2}^{2}\right]-\mathbb{E}\left[\left\\|\boldsymbol{U}\boldsymbol{x}\right\\|_{2}^{2}\right]\mathbb{E}\left[\left\\|\boldsymbol{V}\boldsymbol{x}\right\\|_{2}^{2}\right]$
		$\displaystyle=\mathbb{E}\left[\boldsymbol{x}^{}\boldsymbol{U}^{}\boldsymbol{U}\boldsymbol{x}\boldsymbol{x}^{}\boldsymbol{V}^{}\boldsymbol{V}\boldsymbol{x}\right]-\mathbb{E}\left[\boldsymbol{x}^{}\boldsymbol{U}^{}\boldsymbol{U}\boldsymbol{x}\right]\mathbb{E}\left[\boldsymbol{x}^{}\boldsymbol{V}^{}\boldsymbol{V}\boldsymbol{x}\right]$
		$\displaystyle=\mathbb{E}\left[\boldsymbol{x}^{}\boldsymbol{V}^{}\otimes\boldsymbol{x}^{}\boldsymbol{U}^{}\right]\mathbb{E}\left[\boldsymbol{V}\boldsymbol{x}\otimes\boldsymbol{U}\boldsymbol{x}\right]+\mathbb{E}\left[\boldsymbol{x}^{}\boldsymbol{U}^{}\mathbb{E}\left[\boldsymbol{U}\boldsymbol{x}\boldsymbol{x}^{}\boldsymbol{V}^{}\right]\boldsymbol{V}\boldsymbol{x}\right]$
		$\displaystyle\ \ \ \ \ +2\mathbb{E}[\boldsymbol{x}^{}\boldsymbol{U}^{}]\mathbb{E}[\boldsymbol{U}\boldsymbol{x}]\mathbb{E}[\boldsymbol{x}^{}\boldsymbol{V}^{}]\mathbb{E}[\boldsymbol{V}\boldsymbol{x}].$

We proceed to evaluate each of these three terms.

Since $\boldsymbol{x}\sim\mathcal{CN}(\boldsymbol{0},\boldsymbol{R})$ , we can write $\boldsymbol{x}=\boldsymbol{R}^{1/2}\boldsymbol{y}$ where $\boldsymbol{y}\sim\mathcal{CN}(\boldsymbol{0},{\bf I})$ and $\boldsymbol{R}^{1/2}$ is the unique positive semidefinite squareroot of $\boldsymbol{R}$ . Then, by using the identity $\boldsymbol{X}_{1}\boldsymbol{X}_{2}\otimes\boldsymbol{X}_{3}\boldsymbol{X}_{4}=(\boldsymbol{X}_{1}\otimes\boldsymbol{X}_{3})(\boldsymbol{X}_{2}\otimes\boldsymbol{X}_{4})$ for appropriately sized matrices $\boldsymbol{X}_{1},\boldsymbol{X}_{2},\boldsymbol{X}_{3},\boldsymbol{X}_{4}$ , we get

\mathbb{E}[\boldsymbol{V}\boldsymbol{x}\otimes\boldsymbol{U}\boldsymbol{x}]=\mathbb{E}\left[\boldsymbol{V}\boldsymbol{R}^{1/2}\boldsymbol{y}\otimes\boldsymbol{U}\boldsymbol{R}^{1/2}\boldsymbol{y}\right]=\mathbb{E}\left[\left(\boldsymbol{V}\boldsymbol{R}^{1/2}\otimes\boldsymbol{U}\boldsymbol{R}^{1/2}\right)(\boldsymbol{y}\otimes\boldsymbol{y})\right]=\left(\boldsymbol{V}\boldsymbol{R}^{1/2}\otimes\boldsymbol{U}\boldsymbol{R}^{1/2}\right)\mathbb{E}[\boldsymbol{y}\otimes\boldsymbol{y}]

Since the entries of $\boldsymbol{y}$ are i.i.d. $\mathcal{CN}(0,1)$ , $\mathbb{E}\left[\boldsymbol{y}[n]\boldsymbol{y}[n^{\prime}]\right]=0$ for all indices $n,n^{\prime}=0,\ldots,N-1$ . (Note that $\mathbb{E}[\boldsymbol{y}[n]^{2}]\neq 0$ if the entries of $\boldsymbol{y}$ were i.i.d. $\mathcal{N}(0,1)$ instead of $\mathcal{CN}(0,1)$ .) Hence, $\mathbb{E}[\boldsymbol{y}\otimes\boldsymbol{y}]=\boldsymbol{0}$ , and thus, $\mathbb{E}[\boldsymbol{V}\boldsymbol{x}\otimes\boldsymbol{U}\boldsymbol{x}]=(\boldsymbol{V}\boldsymbol{R}^{1/2}\otimes\boldsymbol{U}\boldsymbol{R}^{1/2})\mathbb{E}[\boldsymbol{y}\otimes\boldsymbol{y}]=\boldsymbol{0}$ . Similarly, $\mathbb{E}[\boldsymbol{x}^{*}\boldsymbol{V}^{*}\otimes\boldsymbol{x}^{*}\boldsymbol{U}^{*}]=\boldsymbol{0}^{*}$ , and so, $\mathbb{E}\left[\boldsymbol{x}^{*}\boldsymbol{V}^{*}\otimes\boldsymbol{x}^{*}\boldsymbol{U}^{*}\right]\mathbb{E}\left[\boldsymbol{V}\boldsymbol{x}\otimes\boldsymbol{U}\boldsymbol{x}\right]=0$ .

Using the cyclic property of the trace operator, linearity of the trace and expectation operators, and the fact that $\mathbb{E}[\boldsymbol{x}\boldsymbol{x}^{*}]=\boldsymbol{R}$ , we find that

	$\displaystyle E\left[\boldsymbol{x}^{}\boldsymbol{U}^{}\mathbb{E}\left[\boldsymbol{U}\boldsymbol{x}\boldsymbol{x}^{}\boldsymbol{V}^{}\right]\boldsymbol{V}\boldsymbol{x}\right]$	$\displaystyle=\mathbb{E}\left[\operatorname{tr}\left[\mathbb{E}\left[\boldsymbol{U}\boldsymbol{x}\boldsymbol{x}^{}\boldsymbol{V}^{}\right]\boldsymbol{V}\boldsymbol{x}\boldsymbol{x}^{}\boldsymbol{U}^{}\right]\right]$
		$\displaystyle=\operatorname{tr}\left[\mathbb{E}\left[\mathbb{E}\left[\boldsymbol{U}\boldsymbol{x}\boldsymbol{x}^{}\boldsymbol{V}^{}\right]\boldsymbol{V}\boldsymbol{x}\boldsymbol{x}^{}\boldsymbol{U}^{}\right]\right]$
		$\displaystyle=\operatorname{tr}\left[\boldsymbol{U}\mathbb{E}\left[\boldsymbol{x}\boldsymbol{x}^{}\right]\boldsymbol{V}^{}\boldsymbol{V}\mathbb{E}\left[\boldsymbol{x}\boldsymbol{x}^{}\right]\boldsymbol{U}^{}\right]$
		$\displaystyle=\operatorname{tr}\left[\boldsymbol{U}\boldsymbol{R}\boldsymbol{V}^{}\boldsymbol{V}\boldsymbol{R}\boldsymbol{U}^{}\right]$
		$\displaystyle=\left\\|\boldsymbol{U}\boldsymbol{R}\boldsymbol{V}^{*}\right\\|_{F}^{2}.$

Finally, since $\mathbb{E}[\boldsymbol{x}]=\boldsymbol{0}$ , we have $\mathbb{E}[\boldsymbol{U}\boldsymbol{x}]=\mathbb{E}[\boldsymbol{V}\boldsymbol{x}]=\boldsymbol{0}$ and $\mathbb{E}[\boldsymbol{x}^{*}\boldsymbol{U}^{*}]=\mathbb{E}[\boldsymbol{x}^{*}\boldsymbol{V}^{*}]=\boldsymbol{0}^{*}$ . Hence,

2\mathbb{E}[\boldsymbol{x}^{*}\boldsymbol{U}^{*}]\mathbb{E}[\boldsymbol{U}\boldsymbol{x}]\mathbb{E}[\boldsymbol{x}^{*}\boldsymbol{V}^{*}]\mathbb{E}[\boldsymbol{V}\boldsymbol{x}]=0.

Adding these three terms yields,

	$\displaystyle\operatorname{Cov}\left[\left\\|\boldsymbol{U}\boldsymbol{x}\right\\|_{2}^{2},\left\\|\boldsymbol{V}\boldsymbol{x}\right\\|_{2}^{2}\right]$	$\displaystyle=\mathbb{E}\left[\boldsymbol{x}^{}\boldsymbol{V}^{}\otimes\boldsymbol{x}^{}\boldsymbol{U}^{}\right]\mathbb{E}\left[\boldsymbol{V}\boldsymbol{x}\otimes\boldsymbol{U}\boldsymbol{x}\right]+\mathbb{E}\left[\boldsymbol{x}^{}\boldsymbol{U}^{}\mathbb{E}\left[\boldsymbol{U}\boldsymbol{x}\boldsymbol{x}^{}\boldsymbol{V}^{}\right]\boldsymbol{V}\boldsymbol{x}\right]$
		$\displaystyle\ \ \ \ \ +2\mathbb{E}[\boldsymbol{x}^{}\boldsymbol{U}^{}]\mathbb{E}[\boldsymbol{U}\boldsymbol{x}]\mathbb{E}[\boldsymbol{x}^{}\boldsymbol{V}^{}]\mathbb{E}[\boldsymbol{V}\boldsymbol{x}].$
		$\displaystyle=\left\\|\boldsymbol{U}\boldsymbol{R}\boldsymbol{V}^{*}\right\\|_{F}^{2}.$

∎

A.2 Concentration of norms of Gaussian random variables

Next, we state a result from [50] regarding concentration bounds for sums of independent exponential random variables.

Lemma 3.

[50] Let $Z_{0},\ldots,Z_{N-1}$ be independent exponential random variables with $\mathbb{E}[Z_{n}]=\mu_{n}$ . Then, the sum

Z:=\displaystyle\sum_{n=0}^{N-1}Z_{n}

satisfies

\displaystyle\mathbb{P}\left\{Z\geq\beta\mathbb{E}[Z]\right\}\leq\beta^{-1}e^{-\kappa(\beta-1-\ln\beta)}\quad\text{for}\quad\beta>1,

and

\displaystyle\mathbb{P}\left\{Z\leq\beta\mathbb{E}[Z]\right\}\leq e^{-\kappa(\beta-1-\ln\beta)}\quad\text{for}\quad 0<\beta<1,

where

\kappa=\dfrac{\sum\limits_{n=0}^{N-1}\mu_{n}}{\max\limits_{n=0,\ldots,N-1}\mu_{n}}.

We now apply this lemma to derive concentration bounds for $\|\boldsymbol{A}\boldsymbol{x}\|_{2}^{2}$ , where $\boldsymbol{x}$ is a vector of Gaussian random variables and $\boldsymbol{A}$ is a matrix.

Lemma 4.

Let $\boldsymbol{x}\sim\mathcal{CN}(\boldsymbol{0},\boldsymbol{R})$ for some positive semidefinite $\boldsymbol{R}\in\mathbb{C}^{N\times N}$ . Also, let $\boldsymbol{A}\in\mathbb{C}^{K\times N}$ . Then, the random variable $\left\|\boldsymbol{A}\boldsymbol{x}\right\|_{2}^{2}$ satisfies

\mathbb{P}\left\{\left\|\boldsymbol{A}\boldsymbol{x}\right\|_{2}^{2}\geq\beta\mathbb{E}\left[\left\|\boldsymbol{A}\boldsymbol{x}\right\|_{2}^{2}\right]\right\}\leq\beta^{-1}e^{-\kappa(\beta-1-\ln\beta)}\quad\text{for}\quad\beta>1,

and

\mathbb{P}\left\{\left\|\boldsymbol{A}\boldsymbol{x}\right\|_{2}^{2}\leq\beta\mathbb{E}\left[\left\|\boldsymbol{A}\boldsymbol{x}\right\|_{2}^{2}\right]\right\}\leq e^{-\kappa(\beta-1-\ln\beta)}\quad\text{for}\quad 0<\beta<1,

where

\kappa=\dfrac{\operatorname{tr}\left[\boldsymbol{A}\boldsymbol{R}\boldsymbol{A}^{*}\right]}{\left\|\boldsymbol{A}\boldsymbol{R}\boldsymbol{A}^{*}\right\|}.

Proof.

Since $\boldsymbol{x}\sim\mathcal{CN}(\boldsymbol{0},\boldsymbol{R})$ , we can write $\boldsymbol{x}=\boldsymbol{R}^{1/2}\boldsymbol{y}$ where $\boldsymbol{y}\sim\mathcal{CN}(\boldsymbol{0},{\bf I})$ and $\boldsymbol{R}^{1/2}$ is the unique positive semidefinite squareroot of $\boldsymbol{R}$ . Then, using eigendecomposition, we can write $\boldsymbol{R}^{1/2}\boldsymbol{A}^{*}\boldsymbol{A}\boldsymbol{R}^{1/2}=\boldsymbol{W}\boldsymbol{D}\boldsymbol{W}^{*}$ where $\boldsymbol{W}$ is unitary and $\boldsymbol{D}=\operatorname{diag}(d_{1},\ldots,\,\mathrm{d}_{n})$ . Since $\boldsymbol{W}$ is unitary, $\boldsymbol{z}:=\boldsymbol{W}^{*}\boldsymbol{y}\sim\mathcal{CN}(\boldsymbol{0},{\bf I}_{N})$ . Then, we have:

\left\|\boldsymbol{A}\boldsymbol{x}\right\|_{2}^{2}=\boldsymbol{x}^{*}\boldsymbol{A}^{*}\boldsymbol{A}\boldsymbol{x}=\boldsymbol{y}^{*}\boldsymbol{R}^{1/2}\boldsymbol{A}^{*}\boldsymbol{A}\boldsymbol{R}^{1/2}\boldsymbol{y}=\boldsymbol{y}^{*}\boldsymbol{W}\boldsymbol{D}\boldsymbol{W}^{*}\boldsymbol{y}=\boldsymbol{z}^{*}\boldsymbol{D}\boldsymbol{z}=\sum_{n=0}^{N-1}d_{n}\left|\boldsymbol{z}[n]\right|^{2}.

Since $\boldsymbol{z}\sim\mathcal{CN}(\boldsymbol{0},{\bf I}_{N})$ , we have that $\boldsymbol{z}[0],\ldots,\boldsymbol{z}[N-1]$ are i.i.d. $\mathcal{CN}(0,1)$ . Hence, $d_{n}|\boldsymbol{z}[n]|^{2}\sim\text{Exp}(d_{n})$ , and are independent. If we apply Lemma 3, the fact that the trace and operator norm are invariant under unitary similarity transforms, and the matrix identities $\operatorname{tr}[\boldsymbol{X}\boldsymbol{X}^{*}]=\operatorname{tr}[\boldsymbol{X}^{*}\boldsymbol{X}]$ and $\|\boldsymbol{X}\boldsymbol{X}^{*}\|=\|\boldsymbol{X}^{*}\boldsymbol{X}\|$ , we obtain

\displaystyle\mathbb{P}\left\{\left\|\boldsymbol{A}\boldsymbol{x}\right\|_{2}^{2}\geq\beta\mathbb{E}\left[\left\|\boldsymbol{A}\boldsymbol{x}\right\|_{2}^{2}\right]\right\}\leq\beta^{-1}e^{-\kappa(\beta-1-\ln\beta)}\quad\text{for}\quad\beta>1,

and

\displaystyle\mathbb{P}\left\{\left\|\boldsymbol{A}\boldsymbol{x}\right\|_{2}^{2}\leq\beta\mathbb{E}\left[\left\|\boldsymbol{A}\boldsymbol{x}\right\|_{2}^{2}\right]\right\}\leq e^{-\kappa(\beta-1-\ln\beta)}\quad\text{for}\quad 0<\beta<1,

where

\kappa=\dfrac{\sum\limits_{n=0}^{N-1}d_{n}}{\max\limits_{n=0,\ldots,N-1}d_{n}}=\dfrac{\operatorname{tr}[\boldsymbol{D}]}{\|\boldsymbol{D}\|}=\dfrac{\operatorname{tr}[\boldsymbol{W}\boldsymbol{D}\boldsymbol{W}^{*}]}{\|\boldsymbol{W}\boldsymbol{D}\boldsymbol{W}^{*}\|}=\dfrac{\operatorname{tr}\left[\boldsymbol{R}^{1/2}\boldsymbol{A}^{*}\boldsymbol{A}\boldsymbol{R}^{1/2}\right]}{\left\|\boldsymbol{R}^{1/2}\boldsymbol{A}^{*}\boldsymbol{A}\boldsymbol{R}^{1/2}\right\|}=\dfrac{\operatorname{tr}\left[\boldsymbol{A}\boldsymbol{R}\boldsymbol{A}^{*}\right]}{\left\|\boldsymbol{A}\boldsymbol{R}\boldsymbol{A}^{*}\right\|}.

∎

We note that similar bounds can be obtained by applying the Hanson-Wright Inequality [51, 52].

A.3 Intermediate results

We continue by presenting a lemma showing that certain matrices can be represented as an integral of a frequency dependent rank- $1$ matrix.

Lemma 5.

For any frequency $f\in\mathbb{R}$ , define a complex sinusoid $\boldsymbol{e}_{f}\in\mathbb{C}^{N}$ by $\boldsymbol{e}_{f}[n]=e^{j2\pi fn}$ for $n=0,\ldots,N-1$ . Then, we have:

\boldsymbol{B}=\int_{-W}^{W}\boldsymbol{e}_{f}\boldsymbol{e}_{f}^{*}\,df,\quad{\bf I}=\int_{-1/2}^{1/2}\boldsymbol{e}_{f}\boldsymbol{e}_{f}^{*}\,df,\quad\text{and}\quad\int_{\Omega}\boldsymbol{e}_{f}\boldsymbol{e}_{f}^{*}\,df={\bf I}-\boldsymbol{B},

where $\Omega=[-\tfrac{1}{2},\tfrac{1}{2}]\setminus[-W,W]$ . Furthermore, if $x(t)$ is a stationary, ergodic, zero-mean, Gaussian stochastic process $x(t)$ with power spectral density $S(f)$ , and $\boldsymbol{x}\in\mathbb{C}^{N}$ is a vector of equispaced samples $\boldsymbol{x}[n]=x(n)$ for $n=0,\ldots,N-1$ , then the covariance matrix of $\boldsymbol{x}$ can be written as

\boldsymbol{R}:=\mathbb{E}[\boldsymbol{x}\boldsymbol{x}^{*}]=\int_{-1/2}^{1/2}S(f)\boldsymbol{e}_{f}\boldsymbol{e}_{f}^{*}\,df.

Proof.

For any $m,n=0,\ldots,N-1$ , we have

\int_{-W}^{W}\boldsymbol{e}_{f}[m]\overline{\boldsymbol{e}_{f}[n]}\,df=\int_{-W}^{W}e^{j2\pi f(m-n)}\,df=\dfrac{\sin[2\pi W(m-n)]}{\pi(m-n)}=\boldsymbol{B}[m,n],

and

\int_{-1/2}^{1/2}\boldsymbol{e}_{f}[m]\overline{\boldsymbol{e}_{f}[n]}\,df=\int_{-1/2}^{1/2}e^{j2\pi f(m-n)}\,df=\begin{cases}1&\text{if}\ m=n\\ 0&\text{if}\ m\neq n\end{cases}={\bf I}[m,n].

We can put these into matrix form as

\int_{-W}^{W}\boldsymbol{e}_{f}\boldsymbol{e}_{f}^{*}\,df=\boldsymbol{B}\quad\text{and}\quad\int_{-1/2}^{1/2}\boldsymbol{e}_{f}\boldsymbol{e}_{f}^{*}\,df={\bf I}.

From this, it follows that

\int_{\Omega}\boldsymbol{e}_{f}\boldsymbol{e}_{f}^{*}\,df=\int_{-1/2}^{1/2}\boldsymbol{e}_{f}\boldsymbol{e}_{f}^{*}\,df-\int_{-W}^{W}\boldsymbol{e}_{f}\boldsymbol{e}_{f}^{*}\,df={\bf I}-\boldsymbol{B}.

Finally, using the definition of the power spectral density, we have

\boldsymbol{R}[m,n]=\mathbb{E}[\boldsymbol{x}[m]\overline{\boldsymbol{x}[n]}]=\int_{-1/2}^{1/2}S(f)e^{j2\pi f(m-n)}\,df=\int_{-1/2}^{1/2}S(f)\boldsymbol{e}_{f}[m]\overline{\boldsymbol{e}_{f}[n]}\,df

for $m,n=0,\ldots,N-1$ . Again, we can put this into matrix form as

\boldsymbol{R}:=\mathbb{E}[\boldsymbol{x}\boldsymbol{x}^{*}]=\int_{-1/2}^{1/2}S(f)\boldsymbol{e}_{f}\boldsymbol{e}_{f}^{*}\,df.

∎

Next, we show that the expectation of the multitaper spectral estimate is the convolution of the power spectral density $S(f)$ with the spectral window $\psi(f)$ .

Lemma 6.

The expectation of the multitaper spectral estimate can be written as

\mathbb{E}\left[\widehat{S}^{\text{mt}}_{K}(f)\right]=\int_{-1/2}^{1/2}S(f-f^{\prime})\psi(f^{\prime})\,df^{\prime}

where

\psi(f):=\displaystyle\dfrac{1}{K}\operatorname{tr}\left[\boldsymbol{S}_{K}^{*}\boldsymbol{e}_{-f}\boldsymbol{e}_{-f}^{*}\boldsymbol{S}_{K}\right]=\dfrac{1}{K}\left\|\boldsymbol{S}_{K}^{*}\boldsymbol{e}_{-f}\right\|_{2}^{2}=\dfrac{1}{K}\sum_{k=0}^{K-1}\left|\sum_{n=0}^{N-1}\boldsymbol{s}_{k}[n]e^{-j2\pi fn}\right|^{2}

is the spectral window of the multitaper estimate.

Proof.

Since $\widehat{S}^{\text{mt}}_{K}(f)=\dfrac{1}{K}\left\|\boldsymbol{S}_{K}^{*}\boldsymbol{E}_{f}^{*}\boldsymbol{x}\right\|_{2}^{2}$ where $\boldsymbol{x}\sim\mathcal{CN}(\boldsymbol{0},\boldsymbol{R})$ , by lemma 2, we have

\mathbb{E}\left[\widehat{S}^{\text{mt}}_{K}(f)\right]=\dfrac{1}{K}\operatorname{tr}\left[\boldsymbol{S}_{K}^{*}\boldsymbol{E}_{f}^{*}\boldsymbol{R}\boldsymbol{E}_{f}\boldsymbol{S}_{K}\right].

We can rewrite this expression as follows:

	$\displaystyle\mathbb{E}\left[\widehat{S}^{\text{mt}}_{K}(f)\right]$	$\displaystyle=\dfrac{1}{K}\operatorname{tr}\left[\boldsymbol{S}_{K}^{}\boldsymbol{E}_{f}^{}\boldsymbol{R}\boldsymbol{E}_{f}\boldsymbol{S}_{K}\right]$
		$\displaystyle=\dfrac{1}{K}\operatorname{tr}\left[\boldsymbol{S}_{K}^{}\boldsymbol{E}_{f}^{}\left(\int_{-1/2}^{1/2}S(f^{\prime})\boldsymbol{e}_{f^{\prime}}\boldsymbol{e}_{f^{\prime}}^{*}\,df^{\prime}\right)\boldsymbol{E}_{f}\boldsymbol{S}_{K}\right]$
		$\displaystyle=\dfrac{1}{K}\int_{-1/2}^{1/2}S(f^{\prime})\operatorname{tr}\left[\boldsymbol{S}_{K}^{}\boldsymbol{E}_{f}^{}\boldsymbol{e}_{f^{\prime}}\boldsymbol{e}_{f^{\prime}}^{*}\boldsymbol{E}_{f}\boldsymbol{S}_{K}\right]\,df^{\prime}$
		$\displaystyle=\dfrac{1}{K}\int_{-1/2}^{1/2}S(f^{\prime})\operatorname{tr}\left[\boldsymbol{S}_{K}^{}\boldsymbol{e}_{f^{\prime}-f}\boldsymbol{e}_{f^{\prime}-f}^{}\boldsymbol{S}_{K}\right]\,df^{\prime}$
		$\displaystyle=\int_{-1/2}^{1/2}S(f^{\prime})\psi(f-f^{\prime})\,df^{\prime},$
		$\displaystyle=\int_{f-1/2}^{f+1/2}S(f-f^{\prime})\psi(f^{\prime})\,df^{\prime}$
		$\displaystyle=\int_{-1/2}^{1/2}S(f-f^{\prime})\psi(f^{\prime})\,df^{\prime},$

where the last line follows since $S(f)$ is $1$ -periodic (by definition) and $\psi(f)=\dfrac{1}{K}\|\boldsymbol{S}_{K}^{*}\boldsymbol{e}_{-f}\|_{2}^{2}$ is $1$ -periodic (because $\boldsymbol{e}_{-f}[n]=e^{-2\pi fn}$ is $1$ -periodic for all $n$ ).

Finally, it’s easy to check that $\psi(f)$ can be written in the following alternate forms:

\psi(f)=\dfrac{1}{K}\operatorname{tr}\left[\boldsymbol{S}_{K}^{*}\boldsymbol{e}_{-f}\boldsymbol{e}_{-f}^{*}\boldsymbol{S}_{K}\right]=\dfrac{1}{K}\boldsymbol{e}_{-f}^{*}\boldsymbol{S}_{K}\boldsymbol{S}_{K}^{*}\boldsymbol{e}_{-f}=\dfrac{1}{K}\|\boldsymbol{S}_{K}^{*}\boldsymbol{e}_{-f}\|_{2}^{2}=\dfrac{1}{K}\sum_{k=0}^{K-1}\left|\sum_{n=0}^{N-1}\boldsymbol{s}_{k}[n]e^{-j2\pi fn}\right|^{2}.

∎

Lemma 7.

The spectral window $\psi(f)$ defined in Lemma 6, satisfies the following properties:

•

$\psi(f)=\psi(-f)$ for all $f\in\mathbb{R}$
•

$\displaystyle\int_{-W}^{W}\psi(f)\,df=1-\Sigma^{(1)}_{K}$ and $\displaystyle\int_{\Omega}\psi(f)\,df=\Sigma^{(1)}_{K}$ where $\Omega=[-\tfrac{1}{2},\tfrac{1}{2}]\setminus[-W,W]$
•

$0\leq\psi(f)\leq\dfrac{N}{K}$

Proof.

First, the spectral window is an even function since

\psi(-f)=\dfrac{1}{K}\sum_{k=0}^{K-1}\left|\sum_{n=0}^{N-1}\boldsymbol{s}_{k}[n]e^{j2\pi fn}\right|^{2}=\dfrac{1}{K}\sum_{k=0}^{K-1}\left|\sum_{n=0}^{N-1}\boldsymbol{s}_{k}[n]e^{-j2\pi fn}\right|^{2}=\psi(f)

for all $f\in\mathbb{R}$ .

As a consequence, we can set $\boldsymbol{S}_{K}=\begin{bmatrix}\boldsymbol{s}_{0}&\boldsymbol{s}_{1}&\cdots&\boldsymbol{s}_{K-1}\end{bmatrix}\in\mathbb{R}^{N\times K}$ and write $\psi(f)=\psi(-f)=\dfrac{1}{K}\operatorname{tr}\left[\boldsymbol{S}_{K}^{*}\boldsymbol{e}_{f}\boldsymbol{e}_{f}^{*}\boldsymbol{S}_{K}\right]=\dfrac{1}{K}\left\|\boldsymbol{S}_{K}^{*}\boldsymbol{e}_{f}\right\|_{2}^{2}$ . Then, the integral of the spectral window $\psi(f)$ over $[-W,W]$ is

	$\displaystyle\int_{-W}^{W}\psi(f)\,df$	$\displaystyle=\int_{-W}^{W}\dfrac{1}{K}\operatorname{tr}\left[\boldsymbol{S}_{K}^{}\boldsymbol{e}_{f}\boldsymbol{e}_{f}^{}\boldsymbol{S}_{K}\right]\,df$
		$\displaystyle=\dfrac{1}{K}\operatorname{tr}\left[\boldsymbol{S}_{K}^{}\left(\int_{-W}^{W}\boldsymbol{e}_{f}\boldsymbol{e}_{f}^{}\,df\right)\boldsymbol{S}_{K}\right]$
		$\displaystyle=\dfrac{1}{K}\operatorname{tr}\left[\boldsymbol{S}_{K}^{*}\boldsymbol{B}\boldsymbol{S}_{K}\right]$
		$\displaystyle=\dfrac{1}{K}\operatorname{tr}\left[\boldsymbol{\Lambda}_{K}\right]$
		$\displaystyle=\dfrac{1}{K}\sum_{k=0}^{K-1}\lambda_{k}$
		$\displaystyle=1-\dfrac{1}{K}\sum_{k=0}^{K-1}(1-\lambda_{k})$
		$\displaystyle=1-\Sigma_{K}^{(1)},$

where we have used the notation $\boldsymbol{\Lambda}_{K}=\operatorname{diag}(\lambda_{0},\ldots,\lambda_{K-1})$ . Similarly, the integral of the spectral window $\psi(f)$ over $\Omega=[-\tfrac{1}{2},-W)\cup(W,\tfrac{1}{2}]$ is

	$\displaystyle\int_{\Omega}\psi(f)\,df$	$\displaystyle=\int_{\Omega}\dfrac{1}{K}\operatorname{tr}\left[\boldsymbol{S}_{K}^{}\boldsymbol{e}_{f}\boldsymbol{e}_{f}^{}\boldsymbol{S}_{K}\right]\,df$
		$\displaystyle=\dfrac{1}{K}\operatorname{tr}\left[\boldsymbol{S}_{K}^{}\left(\int_{\Omega}\boldsymbol{e}_{f}\boldsymbol{e}_{f}^{}\,df\right)\boldsymbol{S}_{K}\right]$
		$\displaystyle=\dfrac{1}{K}\operatorname{tr}\left[\boldsymbol{S}_{K}^{*}({\bf I}-\boldsymbol{B})\boldsymbol{S}_{K}\right]$
		$\displaystyle=\dfrac{1}{K}\operatorname{tr}\left[{\bf I}-\boldsymbol{\Lambda}_{K}\right]$
		$\displaystyle=\dfrac{1}{K}\sum_{k=0}^{K-1}(1-\lambda_{k})$
		$\displaystyle=\Sigma_{K}^{(1)}.$

Finally, the spectral window is bounded by

0\leq\psi(f)=\dfrac{1}{K}\left\|\boldsymbol{S}_{K}^{*}\boldsymbol{e}_{f}\right\|_{2}^{2}\leq\dfrac{1}{K}\left\|\boldsymbol{S}_{K}^{*}\right\|^{2}\left\|\boldsymbol{e}_{f}\right\|_{2}^{2}=\dfrac{N}{K},

where we have used the facts that $\|\boldsymbol{S}_{K}^{*}\|=1$ (because $\boldsymbol{S}_{K}$ is orthonormal) and $\left\|\boldsymbol{e}_{f}\right\|_{2}^{2}=N$ . ∎

A.4 Proof of Theorem 1

By Lemma 6, we have

\mathbb{E}\left[\widehat{S}^{\text{mt}}_{K}(f)\right]=\int_{-1/2}^{1/2}S(f-f^{\prime})\psi(f^{\prime})\,df^{\prime}.

We now split the expression for the bias into two pieces as follows:

	$\displaystyle\operatorname{Bias}\left[\widehat{S}^{\text{mt}}_{K}(f)\right]$	$\displaystyle=\left\|\mathbb{E}\widehat{S}^{\text{mt}}_{K}(f)-S(f)\right\|$
		$\displaystyle=\left\|\int_{-1/2}^{1/2}S(f-f^{\prime})\psi(f^{\prime})\,df^{\prime}-S(f)\right\|$
		$\displaystyle=\left\|\int_{-W}^{W}S(f-f^{\prime})\psi(f^{\prime})\,df^{\prime}+\int_{\Omega}S(f-f^{\prime})\psi(f^{\prime})\,df^{\prime}-S(f)\right\|$
		$\displaystyle\leq\underbrace{\left\|\int_{-W}^{W}S(f-f^{\prime})\psi(f^{\prime})\,df^{\prime}-S(f)\right\|}_{\text{local bias}}+\underbrace{\left\|\int_{\Omega}S(f-f^{\prime})\psi(f^{\prime})\,df^{\prime}\right\|}_{\text{broadband bias}}.$

Since $S(f)$ is twice continuously differentiable, for any $f^{\prime}\in[-W,W]$ , there exists a $\xi_{f^{\prime}}$ between $f-f^{\prime}$ and $f$ such that

S(f-f^{\prime})=S(f)-S^{\prime}(f)f^{\prime}+\dfrac{1}{2}S^{\prime\prime}(\xi_{f^{\prime}})f^{\prime 2}.

Then, since $\left|S^{\prime\prime}(\xi_{f^{\prime}})\right|\leq\displaystyle\max_{\xi\in[f-W,f+W]}|S^{\prime\prime}(\xi)|=M^{\prime\prime}_{f}$ , $\int_{-W}^{W}\psi(f)\,df=1-\Sigma^{(1)}_{K}$ , and $0\leq\psi(f^{\prime})\leq\tfrac{N}{K}$ for all $f^{\prime}\in\mathbb{R}$ , we can bound the local bias as follows:

	$\displaystyle\left\|\int_{-W}^{W}S(f-f^{\prime})\psi(f^{\prime})\,df^{\prime}-S(f)\right\|$	$\displaystyle=\left\|\int_{-W}^{W}\left(S(f)-S^{\prime}(f)f^{\prime}+\dfrac{1}{2}S^{\prime\prime}(\xi_{f^{\prime}})f^{\prime 2}\right)\psi(f^{\prime})\,df^{\prime}-S(f)\right\|$
		$\displaystyle=\left\|\int_{-W}^{W}S(f)\psi(f^{\prime})\,df^{\prime}-S(f)-\int_{-W}^{W}\underbrace{S^{\prime}(f)f^{\prime}\psi(f^{\prime})}_{\text{odd w.r.t.}\ f^{\prime}}\,df^{\prime}+\dfrac{1}{2}\int_{-W}^{W}S^{\prime\prime}(\xi_{f^{\prime}})f^{\prime 2}\psi(f^{\prime})\,df^{\prime}\right\|$
		$\displaystyle=\left\|\int_{-W}^{W}S(f)\psi(f^{\prime})\,df^{\prime}-S(f)-0+\dfrac{1}{2}\int_{-W}^{W}S^{\prime\prime}(\xi_{f^{\prime}})f^{\prime 2}\psi(f^{\prime})\,df^{\prime}\right\|$
		$\displaystyle\leq\left\|\int_{-W}^{W}S(f)\psi(f^{\prime})\,df^{\prime}-S(f)\right\|+\left\|\dfrac{1}{2}\int_{-W}^{W}S^{\prime\prime}(\xi_{f^{\prime}})f^{\prime 2}\psi(f^{\prime})\,df^{\prime}\right\|$
		$\displaystyle=\left\|S(f)(1-\Sigma_{K}^{(1)})-S(f)\right\|+\left\|\dfrac{1}{2}\int_{-W}^{W}S^{\prime\prime}(\xi_{f^{\prime}})f^{\prime 2}\psi(f^{\prime})\,df^{\prime}\right\|$
		$\displaystyle\leq S(f)\Sigma_{K}^{(1)}+\dfrac{1}{2}\int_{-W}^{W}\|S^{\prime\prime}(\xi_{f^{\prime}})\|\|f^{\prime}\|^{2}\psi(f^{\prime})\,df^{\prime}$
		$\displaystyle\leq S(f)\Sigma_{K}^{(1)}+\dfrac{1}{2}\int_{-W}^{W}M^{\prime\prime}_{f}\|f^{\prime}\|^{2}\dfrac{N}{K}\,df^{\prime}$
		$\displaystyle=S(f)\Sigma_{K}^{(1)}+\dfrac{M^{\prime\prime}_{f}NW^{3}}{3K}$
		$\displaystyle\leq M_{f}\Sigma_{K}^{(1)}+\dfrac{M^{\prime\prime}_{f}NW^{3}}{3K}$

Since $S(f^{\prime})\leq\displaystyle\max_{f^{\prime}\in\mathbb{R}}S(f^{\prime})=M$ , we can bound the broadband bias as follows:

	$\displaystyle\left\|\int_{\Omega}S(f-f^{\prime})\psi(f^{\prime})\,df^{\prime}\right\|$	$\displaystyle=\int_{\Omega}S(f-f^{\prime})\psi(f^{\prime})\,df^{\prime}$
		$\displaystyle\leq\int_{\Omega}M\psi(f^{\prime})\,df^{\prime}$
		$\displaystyle=M\Sigma_{K}^{(1)}$

Combining the bounds on the local bias and broadband bias yields

	$\displaystyle\operatorname{Bias}\left[\widehat{S}^{\text{mt}}_{K}(f)\right]$	$\displaystyle\leq\underbrace{\left\|\int_{-W}^{W}S(f-f^{\prime})\psi(f^{\prime})\,df^{\prime}-S(f)\right\|}_{\text{local bias}}+\underbrace{\left\|\int_{\Omega}S(f-f^{\prime})\psi(f^{\prime})\,df^{\prime}\right\|}_{\text{broadband bias}}$
		$\displaystyle\leq M_{f}\Sigma_{K}^{(1)}+\dfrac{M^{\prime\prime}_{f}NW^{3}}{3K}+M\Sigma_{K}^{(1)}$
		$\displaystyle=\dfrac{M^{\prime\prime}_{f}NW^{3}}{3K}+(M+M_{f})\Sigma_{K}^{(1)},$

which establishes Theorem 1.

A.5 Proof of Theorem 2

Without the assumption that $S(f)$ is twice differentiable, we can still obtain a bound on the bias. Using the bounds $m_{f}=\displaystyle\min_{\xi\in[f-W,f+W]}S(\xi)\leq S(f^{\prime})\leq\max_{\xi\in[f-W,f+W]}S(\xi)=M_{f}$ and $0\leq S(f^{\prime})\leq\max_{\xi\in\mathbb{R}}S(\xi)=M$ along with the integrals $\int_{-W}^{W}\psi(f)\,df=1-\Sigma^{(1)}_{K}$ and $\int_{\Omega}\psi(f)\,df=\Sigma^{(1)}_{K}$ , we can obtain the following upper bound on $\mathbb{E}\widehat{S}^{\text{mt}}_{K}(f)-S(f)$ :

	$\displaystyle\mathbb{E}\widehat{S}^{\text{mt}}_{K}(f)-S(f)$	$\displaystyle=\int_{-1/2}^{1/2}S(f-f^{\prime})\psi(f^{\prime})\,df^{\prime}-S(f^{\prime})$
		$\displaystyle=\int_{-W}^{W}S(f-f^{\prime})\psi(f^{\prime})\,df^{\prime}+\int_{\Omega}S(f-f^{\prime})\psi(f^{\prime})\,df^{\prime}-S(f^{\prime})$
		$\displaystyle\leq\int_{-W}^{W}M_{f}\psi(f^{\prime})\,df^{\prime}+\int_{\Omega}M\psi(f^{\prime})\,df^{\prime}-m_{f}$
		$\displaystyle=M_{f}(1-\Sigma^{(1)}_{K})+M\Sigma^{(1)}_{K}-m_{f}$
		$\displaystyle=(M_{f}-m_{f})(1-\Sigma^{(1)}_{K})+(M-m_{f})\Sigma^{(1)}_{K}$
		$\displaystyle\leq(M_{f}-m_{f})(1-\Sigma^{(1)}_{K})+M\Sigma^{(1)}_{K},$

Similarly, we can obtain the following lower bound on $\mathbb{E}\widehat{S}^{\text{mt}}_{K}(f)-S(f)$ :

	$\displaystyle\mathbb{E}\widehat{S}^{\text{mt}}_{K}(f)-S(f)$	$\displaystyle=\int_{-1/2}^{1/2}S(f-f^{\prime})\psi(f^{\prime})\,df^{\prime}-S(f^{\prime})$
		$\displaystyle=\int_{-W}^{W}S(f-f^{\prime})\psi(f^{\prime})\,df^{\prime}+\int_{\Omega}S(f-f^{\prime})\psi(f^{\prime})\,df^{\prime}-S(f^{\prime})$
		$\displaystyle\geq\int_{-W}^{W}m_{f}\psi(f^{\prime})\,df^{\prime}+\int_{\Omega}0\psi(f^{\prime})\,df^{\prime}-M_{f}$
		$\displaystyle=m_{f}(1-\Sigma^{(1)}_{K})+0-M_{f}$
		$\displaystyle=-(M_{f}-m_{f})(1-\Sigma^{(1)}_{K})-M_{f}\Sigma^{(1)}_{K}$
		$\displaystyle\geq-(M_{f}-m_{f})(1-\Sigma^{(1)}_{K})-M\Sigma^{(1)}_{K}.$

From the above two bounds, we have

\operatorname{Bias}\left[\widehat{S}^{\text{mt}}_{K}(f)\right]=\left|\mathbb{E}\widehat{S}^{\text{mt}}_{K}(f)-S(f)\right|\leq(M_{f}-m_{f})(1-\Sigma_{K}^{(1)})+M\Sigma_{K}^{(1)},

which establishes Theorem 2.

A.6 Proof of Theorem 3

\operatorname{Var}\left[\widehat{S}^{\text{mt}}_{K}(f)\right]=\operatorname{Var}\left[\dfrac{1}{K}\left\|\boldsymbol{S}_{K}^{*}\boldsymbol{E}_{f}^{*}\boldsymbol{x}\right\|_{2}^{2}\right]=\dfrac{1}{K^{2}}\operatorname{Cov}\left[\left\|\boldsymbol{S}_{K}^{*}\boldsymbol{E}_{f}^{*}\boldsymbol{x}\right\|_{2}^{2},\left\|\boldsymbol{S}_{K}^{*}\boldsymbol{E}_{f}^{*}\boldsymbol{x}\right\|_{2}^{2}\right]=\dfrac{1}{K^{2}}\left\|\boldsymbol{S}_{K}^{*}\boldsymbol{E}_{f}^{*}\boldsymbol{R}\boldsymbol{E}_{f}\boldsymbol{S}_{K}\right\|_{F}^{2}.

We focus on bounding the Frobenius norm of $\boldsymbol{S}_{K}^{*}\boldsymbol{E}_{f}^{*}\boldsymbol{R}\boldsymbol{E}_{f}\boldsymbol{S}_{K}$ . To do this, we first split it into two pieces - an integral over $[-W,W]$ and an integral over $\Omega=[-\tfrac{1}{2},\tfrac{1}{2}]\setminus[-W,W]$ :

	$\displaystyle\boldsymbol{S}_{K}^{}\boldsymbol{E}_{f}^{}\boldsymbol{R}\boldsymbol{E}_{f}\boldsymbol{S}_{K}$	$\displaystyle=\boldsymbol{S}_{K}^{}\boldsymbol{E}_{f}^{}\left(\int_{-1/2}^{1/2}S(f^{\prime})\boldsymbol{e}_{f^{\prime}}\boldsymbol{e}_{f^{\prime}}^{*}\,df^{\prime}\right)\boldsymbol{E}_{f}\boldsymbol{S}_{K}$
		$\displaystyle=\boldsymbol{S}_{K}^{}\left(\int_{-1/2}^{1/2}S(f^{\prime})\boldsymbol{E}_{f}^{}\boldsymbol{e}_{f^{\prime}}\boldsymbol{e}_{f^{\prime}}^{*}\boldsymbol{E}_{f}\,df^{\prime}\right)\boldsymbol{S}_{K}$
		$\displaystyle=\boldsymbol{S}_{K}^{}\left(\int_{-1/2}^{1/2}S(f^{\prime})\boldsymbol{e}_{f^{\prime}-f}\boldsymbol{e}_{f^{\prime}-f}^{}\,df^{\prime}\right)\boldsymbol{S}_{K}$
		$\displaystyle=\boldsymbol{S}_{K}^{}\left(\int_{f-1/2}^{f+1/2}S(f+f^{\prime})\boldsymbol{e}_{f^{\prime}}\boldsymbol{e}_{f^{\prime}}^{}\,df^{\prime}\right)\boldsymbol{S}_{K}$
		$\displaystyle=\boldsymbol{S}_{K}^{}\left(\int_{-1/2}^{1/2}S(f+f^{\prime})\boldsymbol{e}_{f^{\prime}}\boldsymbol{e}_{f^{\prime}}^{}\,df^{\prime}\right)\boldsymbol{S}_{K}$
		$\displaystyle=\boldsymbol{S}_{K}^{}\left(\int_{-W}^{W}S(f+f^{\prime})\boldsymbol{e}_{f^{\prime}}\boldsymbol{e}_{f^{\prime}}^{}\,df^{\prime}+\int_{\Omega}S(f+f^{\prime})\boldsymbol{e}_{f^{\prime}}\boldsymbol{e}_{f^{\prime}}^{*}\,df^{\prime}\right)\boldsymbol{S}_{K}$
		$\displaystyle=\boldsymbol{S}_{K}^{}\left(\int_{-W}^{W}S(f+f^{\prime})\boldsymbol{e}_{f^{\prime}}\boldsymbol{e}_{f^{\prime}}^{}\,df^{\prime}\right)\boldsymbol{S}_{K}+\boldsymbol{S}_{K}^{}\left(\int_{\Omega}S(f+f^{\prime})\boldsymbol{e}_{f^{\prime}}\boldsymbol{e}_{f^{\prime}}^{}\,df^{\prime}\right)\boldsymbol{S}_{K}.$

We will proceed by bounding the Frobenius norm of the two pieces above, and then applying the triangle inequality. Since $S(f)\leq\displaystyle\max_{f\in\mathbb{R}}S(f)=M$ for all $f\in\mathbb{R}$ , we trivially have

\boldsymbol{S}_{K}^{*}\left(\int_{\Omega}S(f+f^{\prime})\boldsymbol{e}_{f^{\prime}}\boldsymbol{e}_{f^{\prime}}^{*}\,df^{\prime}\right)\boldsymbol{S}_{K}\preceq\boldsymbol{S}_{K}^{*}\left(\int_{\Omega}M\boldsymbol{e}_{f^{\prime}}\boldsymbol{e}_{f^{\prime}}^{*}\,df^{\prime}\right)\boldsymbol{S}_{K}=\boldsymbol{S}_{K}^{*}[M({\bf I}-\boldsymbol{B})]\boldsymbol{S}_{K}=M({\bf I}-\boldsymbol{\Lambda}_{K}).

Then, since $\boldsymbol{P}\preceq\boldsymbol{Q}$ implies $\|\boldsymbol{P}\|_{F}\leq\|\boldsymbol{Q}\|_{F}$ , we have

\left\|\boldsymbol{S}_{K}^{*}\left(\int_{\Omega}S(f+f^{\prime})\boldsymbol{e}_{f^{\prime}}\boldsymbol{e}_{f^{\prime}}^{*}\,df^{\prime}\right)\boldsymbol{S}_{K}\right\|_{F}\leq\left\|M({\bf I}-\boldsymbol{\Lambda}_{K})\right\|_{F}=\sqrt{\sum_{k=0}^{K-1}M^{2}(1-\lambda_{k})^{2}}=M\sqrt{K}\Sigma^{(2)}_{K}.

Obtaining a good bound on the first piece requires a more intricate argument. We define $\mathbbm{1}_{W}(f)=1$ if $f\in[-W,W]$ and $\mathbbm{1}_{W}(f)=0$ if $f\in[-\tfrac{1}{2},\tfrac{1}{2}]\setminus[-W,W]$ . For convenience, we also extend $\mathbbm{1}_{W}(f)$ to $f\in\mathbb{R}$ such that $\mathbbm{1}_{W}(f)$ is $1$ -periodic. With this definition, we can write

	$\displaystyle\int_{-W}^{W}S(f+f^{\prime})\boldsymbol{e}_{f^{\prime}}\boldsymbol{e}_{f^{\prime}}^{*}\,df^{\prime}$	$\displaystyle=\int_{-1/2}^{1/2}S(f+f^{\prime})\mathbbm{1}_{W}(f^{\prime})\boldsymbol{e}_{f^{\prime}}\boldsymbol{e}_{f^{\prime}}^{*}\,df^{\prime}$
		$\displaystyle=\int_{0}^{1}S(f+f^{\prime})\mathbbm{1}_{W}(f^{\prime})\boldsymbol{e}_{f^{\prime}}\boldsymbol{e}_{f^{\prime}}^{*}\,df^{\prime}$
		$\displaystyle=\sum_{\ell=0}^{N-1}\int_{\tfrac{\ell}{N}}^{\tfrac{\ell+1}{N}}S(f+f^{\prime})\mathbbm{1}_{W}(f^{\prime})\boldsymbol{e}_{f^{\prime}}\boldsymbol{e}_{f^{\prime}}^{*}\,df^{\prime}$
		$\displaystyle=\sum_{\ell=0}^{N-1}\int_{0}^{\tfrac{1}{N}}S(f+f^{\prime}+\tfrac{\ell}{N})\mathbbm{1}_{W}(f^{\prime}+\tfrac{\ell}{N})\boldsymbol{e}_{f^{\prime}+\tfrac{\ell}{N}}\boldsymbol{e}_{f^{\prime}+\tfrac{\ell}{N}}^{*}\,df^{\prime}$
		$\displaystyle=\int_{0}^{\tfrac{1}{N}}\sum_{\ell=0}^{N-1}S(f+f^{\prime}+\tfrac{\ell}{N})\mathbbm{1}_{W}(f^{\prime}+\tfrac{\ell}{N})\boldsymbol{e}_{f^{\prime}+\tfrac{\ell}{N}}\boldsymbol{e}_{f^{\prime}+\tfrac{\ell}{N}}^{*}\,df^{\prime},$

where the second line holds since $S(f)$ , $\mathbbm{1}_{W}(f)$ , and $\boldsymbol{e}_{f}$ are $1$ -periodic. Now, for any $f^{\prime}\in\mathbb{R}$ , the vectors $\left\{\tfrac{1}{\sqrt{N}}\boldsymbol{e}_{f^{\prime}+\tfrac{\ell}{N}}\right\}_{\ell=0}^{N-1}$ form an orthonormal basis of $\mathbb{C}^{N}$ . Hence, we have

\left\|\sum_{\ell=0}^{N-1}a_{\ell}\boldsymbol{e}_{f^{\prime}+\tfrac{\ell}{N}}\boldsymbol{e}_{f^{\prime}+\tfrac{\ell}{N}}^{*}\right\|_{F}^{2}=N^{2}\sum_{\ell=0}^{N-1}|a_{\ell}|^{2}

for any choice of coefficients $\{a_{\ell}\}_{\ell=0}^{N-1}$ and offset frequency $f^{\prime}\in\mathbb{R}$ . By applying this formula, along with the triangle inequality and the Cauchy-Shwarz Integral inequality, we obtain

	$\displaystyle\left\\|\int_{-W}^{W}S(f+f^{\prime})\boldsymbol{e}_{f^{\prime}}\boldsymbol{e}_{f^{\prime}}^{*}\,df^{\prime}\right\\|_{F}^{2}$	$\displaystyle=\left\\|\int_{0}^{\tfrac{1}{N}}\sum_{\ell=0}^{N-1}S(f+f^{\prime}+\tfrac{\ell}{N})\mathbbm{1}_{W}(f^{\prime}+\tfrac{\ell}{N})\boldsymbol{e}_{f^{\prime}+\tfrac{\ell}{N}}\boldsymbol{e}_{f^{\prime}+\tfrac{\ell}{N}}^{*}\,df^{\prime}\right\\|_{F}^{2}$
		$\displaystyle\leq\left(\int_{0}^{\tfrac{1}{N}}\left\\|\sum_{\ell=0}^{N-1}S(f+f^{\prime}+\tfrac{\ell}{N})\mathbbm{1}_{W}(f^{\prime}+\tfrac{\ell}{N})\boldsymbol{e}_{f^{\prime}+\tfrac{\ell}{N}}\boldsymbol{e}_{f^{\prime}+\tfrac{\ell}{N}}^{*}\right\\|_{F}\,df^{\prime}\right)^{2}$
		$\displaystyle=\left(\int_{0}^{\tfrac{1}{N}}N\left(\sum_{\ell=0}^{N-1}S(f+f^{\prime}+\tfrac{\ell}{N})^{2}\mathbbm{1}_{W}(f^{\prime}+\tfrac{\ell}{N})^{2}\right)^{1/2}\,df^{\prime}\right)^{2}$
		$\displaystyle\leq\left(\int_{0}^{\tfrac{1}{N}}N^{2}\,df^{\prime}\right)\left(\int_{0}^{\tfrac{1}{N}}\sum_{\ell=0}^{N-1}S(f+f^{\prime}+\tfrac{\ell}{N})^{2}\mathbbm{1}_{W}(f^{\prime}+\tfrac{\ell}{N})^{2}\,df^{\prime}\right)$
		$\displaystyle=N\sum_{\ell=0}^{N-1}\int_{0}^{\tfrac{1}{N}}S(f+f^{\prime}+\tfrac{\ell}{N})^{2}\mathbbm{1}_{W}(f^{\prime}+\tfrac{\ell}{N})^{2}\,df^{\prime}$
		$\displaystyle=N\sum_{\ell=0}^{N-1}\int_{\tfrac{\ell}{N}}^{\tfrac{\ell+1}{N}}S(f+f^{\prime})^{2}\mathbbm{1}_{W}(f^{\prime})^{2}\,df^{\prime}$
		$\displaystyle=N\int_{0}^{1}S(f+f^{\prime})^{2}\mathbbm{1}_{W}(f^{\prime})^{2}\,df^{\prime}$
		$\displaystyle=N\int_{-1/2}^{1/2}S(f+f^{\prime})^{2}\mathbbm{1}_{W}(f^{\prime})^{2}\,df^{\prime}$
		$\displaystyle=N\int_{-W}^{W}S(f+f^{\prime})^{2}\,df^{\prime}$
		$\displaystyle=2NWR_{f}^{2}$

Since $\boldsymbol{S}_{K}\in\mathbb{R}^{N\times K}$ is orthonormal, $\|\boldsymbol{S}_{K}^{*}\boldsymbol{X}\boldsymbol{S}_{K}\|_{F}\leq\|\boldsymbol{X}\|_{F}$ for any Hermitian matrix $\boldsymbol{X}\in\mathbb{C}^{N\times N}$ . Hence,

\left\|\boldsymbol{S}_{K}^{*}\int_{-W}^{W}S(f+f^{\prime})\boldsymbol{e}_{f^{\prime}}\boldsymbol{e}_{f^{\prime}}^{*}\,df^{\prime}\boldsymbol{S}_{K}\right\|_{F}\leq\left\|\int_{-W}^{W}S(f+f^{\prime})\boldsymbol{e}_{f^{\prime}}\boldsymbol{e}_{f^{\prime}}^{*}\,df^{\prime}\right\|_{F}\leq R_{f}\sqrt{2NW}.

Finally, by applying the two bounds we’ve derived, we obtain

	$\displaystyle\left\\|\boldsymbol{S}_{K}^{}\boldsymbol{E}_{f}^{}\boldsymbol{R}\boldsymbol{E}_{f}\boldsymbol{S}_{K}\right\\|_{F}$	$\displaystyle=\left\\|\boldsymbol{S}_{K}^{}\left(\int_{-W}^{W}S(f+f^{\prime})\boldsymbol{e}_{f^{\prime}}\boldsymbol{e}_{f^{\prime}}^{}\,df^{\prime}\right)\boldsymbol{S}_{K}+\boldsymbol{S}_{K}^{}\left(\int_{\Omega}S(f+f^{\prime})\boldsymbol{e}_{f^{\prime}}\boldsymbol{e}_{f^{\prime}}^{}\,df^{\prime}\right)\boldsymbol{S}_{K}\right\\|_{F}$
		$\displaystyle\leq\left\\|\boldsymbol{S}_{K}^{}\left(\int_{-W}^{W}S(f+f^{\prime})\boldsymbol{e}_{f^{\prime}}\boldsymbol{e}_{f^{\prime}}^{}\,df^{\prime}\right)\boldsymbol{S}_{K}\right\\|_{F}+\left\\|\boldsymbol{S}_{K}^{}\left(\int_{\Omega}S(f+f^{\prime})\boldsymbol{e}_{f^{\prime}}\boldsymbol{e}_{f^{\prime}}^{}\,df^{\prime}\right)\boldsymbol{S}_{K}\right\\|_{F}$
		$\displaystyle\leq R_{f}\sqrt{2NW}+M\sqrt{K}\Sigma^{(2)}_{K},$

and thus,

\operatorname{Var}\left[\widehat{S}^{\text{mt}}_{K}(f)\right]=\dfrac{1}{K^{2}}\left\|\boldsymbol{S}_{K}^{*}\boldsymbol{E}_{f}^{*}\boldsymbol{R}\boldsymbol{E}_{f}\boldsymbol{S}_{K}\right\|_{F}^{2}\leq\dfrac{1}{K}\left(R_{f}\sqrt{\dfrac{2NW}{K}}+M\Sigma^{(2)}_{K}\right)^{2}.

A.7 Proof of Theorem 4

\operatorname{Cov}\left[\widehat{S}^{\text{mt}}_{K}(f_{1}),\widehat{S}^{\text{mt}}_{K}(f_{2})\right]=\operatorname{Cov}\left[\dfrac{1}{K}\left\|\boldsymbol{S}_{K}^{*}\boldsymbol{E}_{f_{1}}^{*}\boldsymbol{x}\right\|_{2}^{2},\dfrac{1}{K}\left\|\boldsymbol{S}_{K}^{*}\boldsymbol{E}_{f_{2}}^{*}\boldsymbol{x}\right\|_{2}^{2}\right]=\dfrac{1}{K^{2}}\left\|\boldsymbol{S}_{K}^{*}\boldsymbol{E}_{f_{1}}^{*}\boldsymbol{R}\boldsymbol{E}_{f_{2}}\boldsymbol{S}_{K}\right\|_{F}^{2}.

We focus on bounding the Frobenius norm of $\boldsymbol{S}_{K}^{*}\boldsymbol{E}_{f_{1}}^{*}\boldsymbol{R}\boldsymbol{E}_{f_{2}}\boldsymbol{S}_{K}$ . To do this, we first split it into three pieces - an integral over $[f_{1}-W,f_{1}+W]$ , an integral over $[f_{2}-W,f_{2}+W]$ , and an integral over $\Omega^{\prime}=[-\tfrac{1}{2},\tfrac{1}{2}]\setminus([f_{1}-W,f_{1}+W]\cup[f_{2}-W,f_{2}+W])$ :

	$\displaystyle\boldsymbol{S}_{K}^{}\boldsymbol{E}_{f_{1}}^{}\boldsymbol{R}\boldsymbol{E}_{f_{2}}\boldsymbol{S}_{K}$	$\displaystyle=\boldsymbol{S}_{K}^{}\boldsymbol{E}_{f_{1}}^{}\left(\int_{-1/2}^{1/2}S(f)\boldsymbol{e}_{f}\boldsymbol{e}_{f}^{*}\,df\right)\boldsymbol{E}_{f_{2}}\boldsymbol{S}_{K}$
		$\displaystyle=\int_{-1/2}^{1/2}S(f)\boldsymbol{S}_{K}^{}\boldsymbol{E}_{f_{1}}^{}\boldsymbol{e}_{f}\boldsymbol{e}_{f}^{*}\boldsymbol{E}_{f_{2}}\boldsymbol{S}_{K}\,df$
		$\displaystyle=\int_{-1/2}^{1/2}S(f)\boldsymbol{S}_{K}^{}\boldsymbol{e}_{f-f_{1}}\boldsymbol{e}_{f-f_{2}}^{}\boldsymbol{S}_{K}\,df$
		$\displaystyle=\int_{f_{1}-W}^{f_{1}+W}S(f)\boldsymbol{S}_{K}^{}\boldsymbol{e}_{f-f_{1}}\boldsymbol{e}_{f-f_{2}}^{}\boldsymbol{S}_{K}\,df+\int_{f_{2}-W}^{f_{2}+W}S(f)\boldsymbol{S}_{K}^{}\boldsymbol{e}_{f-f_{1}}\boldsymbol{e}_{f-f_{2}}^{}\boldsymbol{S}_{K}\,df$
		$\displaystyle\ \ \ \ \ +\int_{\Omega^{\prime}}S(f)\boldsymbol{S}_{K}^{}\boldsymbol{e}_{f-f_{1}}\boldsymbol{e}_{f-f_{2}}^{}\boldsymbol{S}_{K}\,df.$

By using the triangle inequality, the identity $\|\boldsymbol{x}\boldsymbol{y}^{*}\|_{F}=\|\boldsymbol{x}\|_{2}\|\boldsymbol{y}\|_{2}$ for vectors $\boldsymbol{x},\boldsymbol{y}$ , the Cauchy-Shwarz Inequality, and the facts that $\psi(f)\leq\tfrac{N}{K}$ and $\int_{\Omega}\psi(f)\,df=\Sigma^{(1)}_{K}$ , we can bound the Frobenius norm of the first piece by

	$\displaystyle\left\\|\int_{f_{1}-W}^{f_{1}+W}S(f)\boldsymbol{S}_{K}^{}\boldsymbol{e}_{f-f_{1}}\boldsymbol{e}_{f-f_{2}}^{}\boldsymbol{S}_{K}\,df\right\\|_{F}^{2}$	$\displaystyle\leq\left(\int_{f_{1}-W}^{f_{1}+W}\left\\|S(f)\boldsymbol{S}_{K}^{}\boldsymbol{e}_{f-f_{1}}\boldsymbol{e}_{f-f_{2}}^{}\boldsymbol{S}_{K}\right\\|_{F}\,df\right)^{2}$
		$\displaystyle\leq\left(\int_{f_{1}-W}^{f_{1}+W}S(f)\left\\|\boldsymbol{S}_{K}^{}\boldsymbol{e}_{f-f_{1}}\right\\|_{2}\left\\|\boldsymbol{S}_{K}^{}\boldsymbol{e}_{f-f_{2}}\right\\|_{2}\,df\right)^{2}$
		$\displaystyle\leq\left(\int_{f_{1}-W}^{f_{1}+W}S(f)^{2}\left\\|\boldsymbol{S}_{K}^{}\boldsymbol{e}_{f-f_{1}}\right\\|_{2}^{2}\,df\right)\left(\int_{f_{1}-W}^{f_{1}+W}\left\\|\boldsymbol{S}_{K}^{}\boldsymbol{e}_{f-f_{2}}\right\\|_{2}^{2}\,df\right)$
		$\displaystyle=\left(\int_{-W}^{W}S(f+f_{1})^{2}\left\\|\boldsymbol{S}_{K}^{}\boldsymbol{e}_{f}\right\\|_{2}^{2}\,df\right)\left(\int_{f_{1}-f_{2}-W}^{f_{1}-f_{2}+W}\left\\|\boldsymbol{S}_{K}^{}\boldsymbol{e}_{f}\right\\|_{2}^{2}\,df\right)$
		$\displaystyle=\left(\int_{-W}^{W}S(f+f_{1})^{2}K\psi(f)\,df\right)\left(\int_{f_{1}-f_{2}-W}^{f_{1}-f_{2}+W}K\psi(f)\,df\right)$
		$\displaystyle\leq\left(\int_{-W}^{W}S(f+f_{1})^{2}\cdot N\,df\right)\left(\int_{\Omega}K\psi(f)\,df\right)$
		$\displaystyle=2NWR_{f_{1}}^{2}\cdot K\Sigma^{(1)}_{K}$
		$\displaystyle=R_{f_{1}}^{2}\cdot 2NWK\Sigma^{(1)}_{K}$

In a nearly identical manner, we can bound the second piece by

\left\|\int_{f_{2}-W}^{f_{2}+W}S(f)\boldsymbol{S}_{K}^{*}\boldsymbol{e}_{f-f_{1}}\boldsymbol{e}_{f-f_{2}}^{*}\boldsymbol{S}_{K}\,df\right\|_{F}^{2}\leq R_{f_{2}}^{2}\cdot 2NWK\Sigma^{(1)}_{K}.

The third piece can also be bounded in a similar manner, but the details are noticeably different, so we show the derivation:

	$\displaystyle\left\\|\int_{\Omega^{\prime}}S(f)\boldsymbol{S}_{K}^{}\boldsymbol{e}_{f-f_{1}}\boldsymbol{e}_{f-f_{2}}^{}\boldsymbol{S}_{K}\,df\right\\|_{F}^{2}$	$\displaystyle\leq\left(\int_{\Omega^{\prime}}\left\\|S(f)\boldsymbol{S}_{K}^{}\boldsymbol{e}_{f-f_{1}}\boldsymbol{e}_{f-f_{2}}^{}\boldsymbol{S}_{K}\right\\|_{F}\,df\right)^{2}$
		$\displaystyle\leq\left(\int_{\Omega^{\prime}}S(f)\left\\|\boldsymbol{S}_{K}^{}\boldsymbol{e}_{f-f_{1}}\right\\|_{2}\left\\|\boldsymbol{S}_{K}^{}\boldsymbol{e}_{f-f_{2}}\right\\|_{2}\,df\right)^{2}$
		$\displaystyle\leq\left(\int_{\Omega^{\prime}}M\left\\|\boldsymbol{S}_{K}^{}\boldsymbol{e}_{f-f_{1}}\right\\|_{2}\left\\|\boldsymbol{S}_{K}^{}\boldsymbol{e}_{f-f_{2}}\right\\|_{2}\,df\right)^{2}$
		$\displaystyle\leq M^{2}\left(\int_{\Omega^{\prime}}\left\\|\boldsymbol{S}_{K}^{}\boldsymbol{e}_{f-f_{1}}\right\\|_{2}^{2}\,df\right)\left(\int_{\Omega^{\prime}}\left\\|\boldsymbol{S}_{K}^{}\boldsymbol{e}_{f-f_{2}}\right\\|_{2}^{2}\,df\right)$
		$\displaystyle=M^{2}\left(\int_{\Omega^{\prime}}K\psi(f-f_{1})\,df\right)\left(\int_{\Omega^{\prime}}K\psi(f-f_{2})\,df\right)$
		$\displaystyle\leq M^{2}\left(\int_{\Omega^{\prime}_{1}}K\psi(f-f_{1})\,df\right)\left(\int_{\Omega^{\prime}_{2}}K\psi(f-f_{2})\,df\right)$
		$\displaystyle=M^{2}\left(\int_{\Omega}K\psi(f)\,df\right)\left(\int_{\Omega}K\psi(f)\,df\right)$
		$\displaystyle=M^{2}\left(K\Sigma^{(1)}_{K}\right)^{2}$

where $\Omega^{\prime}_{1}=[-\tfrac{1}{2},\tfrac{1}{2}]\setminus[f_{1}-W,f_{1}+W]$ and $\Omega^{\prime}_{2}=[-\tfrac{1}{2},\tfrac{1}{2}]\setminus[f_{2}-W,f_{2}+W]$ .

Finally, by applying the three bounds we’ve derived, we obtain

	$\displaystyle\left\\|\boldsymbol{S}_{K}^{}\boldsymbol{E}_{f_{1}}^{}\boldsymbol{R}\boldsymbol{E}_{f_{2}}\boldsymbol{S}_{K}\right\\|_{F}$	$\displaystyle\leq\left\\|\int_{f_{1}-W}^{f_{1}+W}S(f)\boldsymbol{S}_{K}^{}\boldsymbol{e}_{f-f_{1}}\boldsymbol{e}_{f-f_{2}}^{}\boldsymbol{S}_{K}\,df\right\\|_{F}+\left\\|\int_{f_{2}-W}^{f_{2}+W}S(f)\boldsymbol{S}_{K}^{}\boldsymbol{e}_{f-f_{1}}\boldsymbol{e}_{f-f_{2}}^{}\boldsymbol{S}_{K}\,df\right\\|_{F}$
		$\displaystyle\ \ \ \ \ +\left\\|\int_{\Omega^{\prime}}S(f)\boldsymbol{S}_{K}^{}\boldsymbol{e}_{f-f_{1}}\boldsymbol{e}_{f-f_{2}}^{}\boldsymbol{S}_{K}\,df\right\\|_{F}$
		$\displaystyle\leq R_{f_{1}}\sqrt{2NWK\Sigma^{(1)}_{K}}+R_{f_{2}}\sqrt{2NWK\Sigma^{(1)}_{K}}+MK\Sigma^{(1)}_{K},$

and thus,

0\leq\operatorname{Cov}\left[\widehat{S}^{\text{mt}}_{K}(f_{1}),\widehat{S}^{\text{mt}}_{K}(f_{2})\right]=\dfrac{1}{K^{2}}\left\|\boldsymbol{S}_{K}^{*}\boldsymbol{E}_{f_{1}}^{*}\boldsymbol{R}\boldsymbol{E}_{f_{2}}\boldsymbol{S}_{K}\right\|_{F}^{2}\leq\left((R_{f_{1}}+R_{f_{2}})\sqrt{\dfrac{2NW}{K}\Sigma^{(1)}_{K}}+M\Sigma^{(1)}_{K}\right)^{2}.

A.8 Proof of Theorem 5

\mathbb{P}\left\{\widehat{S}^{\text{mt}}_{K}(f)\geq\beta\mathbb{E}\widehat{S}^{\text{mt}}_{K}(f)\right\}\leq\beta^{-1}\mathrm{e}^{-\kappa_{f}(\beta-1-\ln\beta)}\quad\text{for}\quad\beta>1,

and

\mathbb{P}\left\{\widehat{S}^{\text{mt}}_{K}(f)\leq\beta\mathbb{E}\widehat{S}^{\text{mt}}_{K}(f)\right\}\leq\mathrm{e}^{-\kappa_{f}(\beta-1-\ln\beta)}\quad\text{for}\quad 0<\beta<1,

where

\kappa_{f}=\dfrac{\operatorname{tr}\left[\boldsymbol{S}_{K}^{*}\boldsymbol{E}_{f}^{*}\boldsymbol{R}\boldsymbol{E}_{f}\boldsymbol{S}_{K}\right]}{\left\|\boldsymbol{S}_{K}^{*}\boldsymbol{E}_{f}^{*}\boldsymbol{R}\boldsymbol{E}_{f}\boldsymbol{S}_{K}\right\|}.

We can get an upper bound on $\boldsymbol{S}_{K}^{*}\boldsymbol{E}_{f}^{*}\boldsymbol{R}\boldsymbol{E}_{f}\boldsymbol{S}_{K}$ in the Loewner ordering by splitting it into two pieces as done in the proof of Theorem 3, and then bounding each piece:

	$\displaystyle\boldsymbol{S}_{K}^{}\boldsymbol{E}_{f}^{}\boldsymbol{R}\boldsymbol{E}_{f}\boldsymbol{S}_{K}$	$\displaystyle=\boldsymbol{S}_{K}^{}\left(\int_{-W}^{W}S(f+f^{\prime})\boldsymbol{e}_{f^{\prime}}\boldsymbol{e}_{f^{\prime}}^{}\,df^{\prime}\right)\boldsymbol{S}_{K}+\boldsymbol{S}_{K}^{}\left(\int_{\Omega}S(f+f^{\prime})\boldsymbol{e}_{f^{\prime}}\boldsymbol{e}_{f^{\prime}}^{}\,df^{\prime}\right)\boldsymbol{S}_{K}$
		$\displaystyle\preceq\boldsymbol{S}_{K}^{}\left(\int_{-W}^{W}M_{f}\boldsymbol{e}_{f^{\prime}}\boldsymbol{e}_{f^{\prime}}^{}\,df^{\prime}\right)\boldsymbol{S}_{K}+\boldsymbol{S}_{K}^{}\left(\int_{\Omega}M\boldsymbol{e}_{f^{\prime}}\boldsymbol{e}_{f^{\prime}}^{}\,df^{\prime}\right)\boldsymbol{S}_{K}$
		$\displaystyle=\boldsymbol{S}_{K}^{}\left(M_{f}\boldsymbol{B}\right)\boldsymbol{S}_{K}+\boldsymbol{S}_{K}^{}\left(M({\bf I}-\boldsymbol{B})\right)\boldsymbol{S}_{K}$
		$\displaystyle=M_{f}\boldsymbol{\Lambda}_{K}+M({\bf I}-\boldsymbol{\Lambda}_{K})$
		$\displaystyle=M_{f}{\bf I}+(M-M_{f})({\bf I}-\boldsymbol{\Lambda}_{K}).$

Then, by using the fact that $\boldsymbol{P}\preceq\boldsymbol{Q}\implies\|\boldsymbol{P}\|\leq\|\boldsymbol{Q}\|$ for PSD matrices $\boldsymbol{P}$ and $\boldsymbol{Q}$ , we can bound,

\left\|\boldsymbol{S}_{K}^{*}\boldsymbol{E}_{f}^{*}\boldsymbol{R}\boldsymbol{E}_{f}\boldsymbol{S}_{K}\right\|\leq\left\|M_{f}{\bf I}+(M-M_{f})({\bf I}-\boldsymbol{\Lambda}_{K})\right\|=M_{f}+(M-M_{f})(1-\lambda_{K-1}).

We can also get a lower bound on $\operatorname{tr}[\boldsymbol{S}_{K}^{*}\boldsymbol{E}_{f}^{*}\boldsymbol{R}\boldsymbol{E}_{f}\boldsymbol{S}_{K}]=K\mathbb{E}\left[\widehat{S}^{\text{mt}}_{K}(f)\right]$ by using the formula for $E\left[\widehat{S}^{\text{mt}}_{K}(f)\right]$ from Lemma 6 along with the properties of the spectral window derived in Lemma 7 as follows:

	$\displaystyle\operatorname{tr}[\boldsymbol{S}_{K}^{}\boldsymbol{E}_{f}^{}\boldsymbol{R}\boldsymbol{E}_{f}\boldsymbol{S}_{K}]$	$\displaystyle=K\mathbb{E}\left[\widehat{S}^{\text{mt}}_{K}(f)\right]$
		$\displaystyle=K\int_{-1/2}^{1/2}S(f-f^{\prime})\psi(f^{\prime})\,df^{\prime}$
		$\displaystyle\geq K\int_{-W}^{W}S(f-f^{\prime})\psi(f^{\prime})\,df^{\prime}$
		$\displaystyle=K\int_{-W}^{W}S(f-f^{\prime})\left[\dfrac{N}{K}-\left(\dfrac{N}{K}-\psi(f^{\prime})\right)\right]\,df^{\prime}$
		$\displaystyle=K\int_{-W}^{W}S(f-f^{\prime})\dfrac{N}{K}\,df^{\prime}-K\int_{-W}^{W}S(f-f^{\prime})\left(\dfrac{N}{K}-\psi(f^{\prime})\right)\,df^{\prime}$
		$\displaystyle=N\int_{f-W}^{f+W}S(f^{\prime})\,df^{\prime}-\int_{-W}^{W}S(f-f^{\prime})(N-K\psi(f^{\prime}))\,df^{\prime}$
		$\displaystyle\geq N\int_{f-W}^{f+W}S(f^{\prime})\,df^{\prime}-\int_{-W}^{W}M_{f}(N-K\psi(f^{\prime}))\,df^{\prime}$
		$\displaystyle=2NWA_{f}-\left(2NW-K\left(1-\Sigma^{(1)}_{K}\right)\right)M_{f}$
		$\displaystyle=K\left(1-\Sigma^{(1)}_{K}\right)M_{f}-2NW(M_{f}-A_{f})$

Combining the upper bound on $\left\|\boldsymbol{S}_{K}^{*}\boldsymbol{E}_{f}^{*}\boldsymbol{R}\boldsymbol{E}_{f}\boldsymbol{S}_{K}\right\|$ with the lower bound on $\operatorname{tr}[\boldsymbol{S}_{K}^{*}\boldsymbol{E}_{f}^{*}\boldsymbol{R}\boldsymbol{E}_{f}\boldsymbol{S}_{K}]$ , yields

\kappa_{f}=\dfrac{\operatorname{tr}\left[\boldsymbol{S}_{K}^{*}\boldsymbol{E}_{f}^{*}\boldsymbol{R}\boldsymbol{E}_{f}\boldsymbol{S}_{K}\right]}{\left\|\boldsymbol{S}_{K}^{*}\boldsymbol{E}_{f}^{*}\boldsymbol{R}\boldsymbol{E}_{f}\boldsymbol{S}_{K}\right\|}\geq\dfrac{K\left(1-\Sigma^{(1)}_{K}\right)M_{f}-2NW(M_{f}-A_{f})}{M_{f}+\left(M-M_{f}\right)(1-\lambda_{K-1})}.

Appendix B Proof of Results in Section 4

B.1 Fast algorithm for computing $\Psi(f)$ at grid frequencies

To begin developing our fast approximations for $\widehat{S}^{\text{mt}}_{K}(f)$ , we first show that an eigenvalue weighted sum of $N$ tapered spectral estimates can be evaluated at a grid of frequencies $f\in[L]/L$ where $L\geq 2N$ in $O(L\log L)$ operations and using $O(L)$ memory.

Lemma 8.

For any vector $\boldsymbol{x}\in\mathbb{C}^{N}$ and any integer $L\geq 2N$ , the quantity

\Psi(f):=\sum_{k=0}^{N-1}\lambda_{k}\widehat{S}_{k}(f)\quad\text{where}\quad\widehat{S}_{k}(f)=\left|\sum_{n=0}^{N-1}\boldsymbol{s}_{k}[n]\boldsymbol{x}[n]e^{-j2\pi fn}\right|^{2}

can be evaluated at the grid frequencies $f\in[L]/L$ in $O(L\log L)$ operations and using $O(L)$ memory.

Proof.

Using eigendecomposition, we can write $\boldsymbol{B}=\boldsymbol{S}\boldsymbol{\Lambda}\boldsymbol{S}^{*}$ , where

\boldsymbol{S}=\begin{bmatrix}\boldsymbol{s}_{0}&\boldsymbol{s}_{1}&\cdots&\boldsymbol{s}_{N-1}\end{bmatrix}

and

\boldsymbol{\Lambda}=\operatorname{diag}(\lambda_{0},\lambda_{1},\ldots,\lambda_{N-1}).

For any $f\in\mathbb{R}$ , we let $\boldsymbol{E}_{f}\in\mathbb{C}^{N\times N}$ be a diagonal matrix with diagonal entries

\boldsymbol{E}_{f}[n,n]=e^{j2\pi fn}\quad\text{for}\quad n\in[N].

With this definition, $\Psi(f)$ can be written as

	$\displaystyle\Psi(f)$	$\displaystyle=\sum_{k=0}^{N-1}\lambda_{k}\widehat{S}_{k}(f)$
		$\displaystyle=\sum_{k=0}^{N-1}\lambda_{k}\left\|\sum_{n=0}^{N-1}\boldsymbol{s}_{k}[n]\boldsymbol{x}[n]e^{-j2\pi fn}\right\|^{2}$
		$\displaystyle=\sum_{k=0}^{N-1}\lambda_{k}\left\|\boldsymbol{s}_{k}^{}\boldsymbol{E}_{f}^{}\boldsymbol{x}\right\|^{2}$
		$\displaystyle=\sum_{k=0}^{N-1}\boldsymbol{\Lambda}[k,k]\left\|(\boldsymbol{S}^{}\boldsymbol{E}_{f}^{}\boldsymbol{x})[k]\right\|^{2}$
		$\displaystyle=\boldsymbol{x}^{}\boldsymbol{E}_{f}\boldsymbol{S}\boldsymbol{\Lambda}\boldsymbol{S}^{}\boldsymbol{E}_{f}^{*}\boldsymbol{x}$
		$\displaystyle=\boldsymbol{x}^{}\boldsymbol{E}_{f}\boldsymbol{B}\boldsymbol{E}_{f}^{}\boldsymbol{x}$

This gives us a a formula for $\Psi(f)=\sum_{k=0}^{N-1}\lambda_{k}\widehat{S}_{k}(f)$ which does not require computing any of the Slepian tapers. We will now use the fact that $\boldsymbol{B}$ is a Toeplitz matrix to efficiently compute $\Psi(\tfrac{\ell}{L})$ for all $\ell\in[L]/L$ .

First, note that we can “extend” $\boldsymbol{B}$ to a larger circulant matrix, which is diagonalized by a Fourier Transform matrix. Specifically, define a vector of sinc samples $\boldsymbol{b}\in\mathbb{R}^{L}$ by

\boldsymbol{b}[\ell]=\begin{cases}\dfrac{\sin[2\pi W\ell]}{\pi\ell}&\text{if}\ \ell\in\{0,\ldots,N-1\}\\ 0&\text{if}\ \ell\in\{N,\ldots,L-N\}\\ \dfrac{\sin[2\pi W(L-\ell)]}{\pi(L-\ell)}&\text{if}\ \ell\in\{L-N+1,\ldots,L-1\}\end{cases},

and let $\boldsymbol{B}_{\text{ext}}\in\mathbb{R}^{L\times L}$ be defined by

\boldsymbol{B}_{\text{ext}}[m,n]=\boldsymbol{b}[m-n\pmod{L}]\ \text{for}\ m,n\in[L].

It is easy to check that $\boldsymbol{B}_{\text{ext}}[m,n]=\boldsymbol{B}[m,n]\ \text{for all}\ m,n\in[N]$ , i.e., the upper-left $N\times N$ submatrix of $\boldsymbol{B}_{\text{ext}}$ is $\boldsymbol{B}$ . Hence, we can write

\boldsymbol{B}=\boldsymbol{Z}^{*}\boldsymbol{B}_{\text{ext}}\boldsymbol{Z},

where

\boldsymbol{Z}=\begin{bmatrix}{\bf I}_{N\times N}\\ \boldsymbol{0}_{(L-N)\times N}\end{bmatrix}\in\mathbb{R}^{L\times N}

is a zeropadding matrix. Since $\boldsymbol{B}_{\text{ext}}$ is a circulant matrix whose first column is $\boldsymbol{b}$ , we can write

\boldsymbol{B}_{\text{ext}}=\boldsymbol{F}^{-1}\operatorname{diag}(\boldsymbol{F}\boldsymbol{b})\boldsymbol{F}

where $\boldsymbol{F}\in\mathbb{C}^{L\times L}$ is an FFT matrix, i.e.,

\boldsymbol{F}[m,n]=e^{-j2\pi mn/L}\ \text{for}\ m,n\in[L].

Note that with this normalization, the inverse FFT satisfies

\boldsymbol{F}^{-1}=\dfrac{1}{L}\boldsymbol{F}^{*},

as well as the conjugation identity

\boldsymbol{F}^{-1}\boldsymbol{y}=\dfrac{1}{L}\overline{\boldsymbol{F}\overline{\boldsymbol{y}}}\quad\text{for all}\quad\boldsymbol{y}\in\mathbb{C}^{N}.

Next, for any $f\in\mathbb{R}$ , let $\boldsymbol{D}_{f}\in\mathbb{C}^{L\times L}$ be a diagonal matrix with diagonal entries

\boldsymbol{D}_{f}[m,m]=e^{j2\pi fm}\quad\text{for}\quad m\in[L],

and for each $\ell\in[L]$ , let $\boldsymbol{C}_{\ell}\in\mathbb{C}^{L\times L}$ be a cyclic shift matrix, i.e.,

\boldsymbol{C}_{\ell}[m,n]=\begin{cases}1&\text{if}\ n-m\equiv\ell\pmod{L}\\ 0&\text{otherwise}\end{cases}.

Since $\boldsymbol{E}_{f}$ and $\boldsymbol{D}_{f}$ are both diagonal matrices and $\boldsymbol{E}_{f}[n,n]=\boldsymbol{D}_{f}[n,n]$ for $n\in[N]$ , we have

\boldsymbol{Z}\boldsymbol{E}_{f}^{*}=\boldsymbol{D}_{f}^{*}\boldsymbol{Z}\quad\text{and}\quad\boldsymbol{E}_{f}\boldsymbol{Z}^{*}=\boldsymbol{Z}^{*}\boldsymbol{D}_{f}

for all $f\in\mathbb{R}$ . Also, it is easy to check that cyclically shifting each column of $\boldsymbol{F}$ by $\ell$ indices is equivalent to modulating the rows of $\boldsymbol{F}$ , or more specifically

\boldsymbol{F}\boldsymbol{D}_{\ell/L}^{*}=\boldsymbol{C}_{\ell}\boldsymbol{F}\quad\text{and}\quad\boldsymbol{D}_{\ell/L}\boldsymbol{F}^{*}=\boldsymbol{F}^{*}\boldsymbol{C}_{\ell}^{*}

for all $\ell\in[L]$ . Additionally, for any vectors $\boldsymbol{p},\boldsymbol{q}\in\mathbb{C}^{L}$ , we will denote $\boldsymbol{p}\circ\boldsymbol{q}\in\mathbb{C}^{L}$ to be the pointwise multiplication of $\boldsymbol{p}$ and $\boldsymbol{q}$ , i.e.,

(\boldsymbol{p}\circ\boldsymbol{q})[\ell]=\boldsymbol{p}[\ell]\boldsymbol{q}[\ell]\quad\text{for}\quad\ell\in[L],

and $\boldsymbol{p}\circledast\boldsymbol{q}\in\mathbb{C}^{L}$ to be the circluar cross-correlation of $\boldsymbol{p}$ and $\boldsymbol{q}$ , i.e.,

(\boldsymbol{p}\circledast\boldsymbol{q})[\ell]=\sum_{\ell^{\prime}=0}^{L-1}\overline{\boldsymbol{p}[\ell^{\prime}]}\boldsymbol{q}[\ell^{\prime}+\ell\pmod{L}]\quad\text{for}\quad\ell\in[L].

Note that the circular cross-correlation of $\boldsymbol{p}$ and $\boldsymbol{q}$ can be computed using FFTs via the formula

\boldsymbol{p}\circledast\boldsymbol{q}=\boldsymbol{F}^{-1}(\overline{\boldsymbol{F}\boldsymbol{p}}\circ\boldsymbol{F}\boldsymbol{q}).

We will also use the notation $|\boldsymbol{p}|^{2}=\overline{\boldsymbol{p}}\circ\boldsymbol{p}$ for convenience.

We can now manipulate our formula for $\Psi(\tfrac{\ell}{L})$ as follows

	$\displaystyle\Psi(\tfrac{\ell}{L})$	$\displaystyle:=\boldsymbol{x}^{}\boldsymbol{E}_{\ell/L}\boldsymbol{B}\boldsymbol{E}_{\ell/L}^{}\boldsymbol{x}$
		$\displaystyle=\boldsymbol{x}^{}\boldsymbol{E}_{\ell/L}\boldsymbol{Z}^{}\boldsymbol{B}_{\text{ext}}\boldsymbol{Z}\boldsymbol{E}_{\ell/L}^{*}\boldsymbol{x}$
		$\displaystyle=\boldsymbol{x}^{}\boldsymbol{E}_{\ell/L}\boldsymbol{Z}^{}\boldsymbol{F}^{-1}\operatorname{diag}(\boldsymbol{F}\boldsymbol{b})\boldsymbol{F}\boldsymbol{Z}\boldsymbol{E}_{\ell/L}^{*}\boldsymbol{x}$
		$\displaystyle=\dfrac{1}{L}\boldsymbol{x}^{}\boldsymbol{E}_{\ell/L}\boldsymbol{Z}^{}\boldsymbol{F}^{}\operatorname{diag}(\boldsymbol{F}\boldsymbol{b})\boldsymbol{F}\boldsymbol{Z}\boldsymbol{E}_{\ell/L}^{}\boldsymbol{x}$
		$\displaystyle=\dfrac{1}{L}\boldsymbol{x}^{}\boldsymbol{Z}^{}\boldsymbol{D}_{\ell/L}\boldsymbol{F}^{}\operatorname{diag}(\boldsymbol{F}\boldsymbol{b})\boldsymbol{F}\boldsymbol{D}_{\ell/L}^{}\boldsymbol{Z}\boldsymbol{x}$
		$\displaystyle=\dfrac{1}{L}\boldsymbol{x}^{}\boldsymbol{Z}^{}\boldsymbol{F}^{}\boldsymbol{C}_{\ell}^{}\operatorname{diag}(\boldsymbol{F}\boldsymbol{b})\boldsymbol{C}_{\ell}\boldsymbol{F}\boldsymbol{Z}\boldsymbol{x}$
		$\displaystyle=\dfrac{1}{L}\sum_{\ell^{\prime}=0}^{L-1}(\boldsymbol{F}\boldsymbol{b})[\ell^{\prime}]\cdot\left\|\left(\boldsymbol{C}_{\ell}\boldsymbol{F}\boldsymbol{Z}\boldsymbol{x}\right)[\ell^{\prime}]\right\|^{2}$
		$\displaystyle=\dfrac{1}{L}\sum_{\ell^{\prime}=0}^{L-1}(\boldsymbol{F}\boldsymbol{b})[\ell^{\prime}]\cdot\left\|\left(\boldsymbol{F}\boldsymbol{Z}\boldsymbol{x}\right)[\ell^{\prime}+\ell\pmod{L}]\right\|^{2}$
		$\displaystyle=\left(\dfrac{1}{L}\overline{\boldsymbol{F}\boldsymbol{b}}\circledast\left\|\boldsymbol{F}\boldsymbol{Z}\boldsymbol{x}\right\|^{2}\right)[\ell].$
		$\displaystyle=\left(\boldsymbol{F}^{-1}\left(\dfrac{1}{L}\overline{\boldsymbol{F}\overline{\boldsymbol{F}\boldsymbol{b}}}\circ\boldsymbol{F}\left\|\boldsymbol{F}\boldsymbol{Z}\boldsymbol{x}\right\|^{2}\right)\right)[\ell]$
		$\displaystyle=\left(\boldsymbol{F}^{-1}\left(\boldsymbol{F}^{-1}\boldsymbol{F}\boldsymbol{b}\circ\boldsymbol{F}\left\|\boldsymbol{F}\boldsymbol{Z}\boldsymbol{x}\right\|^{2}\right)\right)[\ell]$
		$\displaystyle=\left(\boldsymbol{F}^{-1}\left(\boldsymbol{b}\circ\boldsymbol{F}\left\|\boldsymbol{F}\boldsymbol{Z}\boldsymbol{x}\right\|^{2}\right)\right)[\ell].$

Therefore,

\begin{bmatrix}\Psi(\tfrac{0}{L})&\Psi(\tfrac{1}{L})&\cdots&\Psi(\tfrac{L-2}{L})&\Psi(\tfrac{L-1}{L})\end{bmatrix}^{T}=\boldsymbol{F}^{-1}\left(\boldsymbol{b}\circ\boldsymbol{F}\left|\boldsymbol{F}\boldsymbol{Z}\boldsymbol{x}\right|^{2}\right).

So to compute $\Psi(f)$ at all grid frequencies $f\in[L]/L$ , we only need to compute $\boldsymbol{F}^{-1}\left(\boldsymbol{b}\circ\boldsymbol{F}\left|\boldsymbol{F}\boldsymbol{Z}\boldsymbol{x}\right|^{2}\right)$ , which can be done in $O(L\log L)$ operations using $O(L)$ memory via three length- $L$ FFTs/IFFTs and a few pointwise operations on vectors of length $L$ . ∎

B.2 Approximations for general multitaper spectral estimates

Next, we present a lemma regarding approximations to spectral estimates which use orthonormal tapers.

Lemma 9.

Let $\boldsymbol{x}\in\mathbb{C}^{N}$ be a vector of $N$ equispaced samples, and let $\{\boldsymbol{v}_{k}\}_{k=0}^{N-1}$ be any orthonormal set of tapers in $\mathbb{C}^{N}$ . For each $k\in[N]$ , define a tapered spectral estimate

V_{k}(f)=\left|\sum_{n=0}^{N-1}\boldsymbol{v}_{k}[n]\boldsymbol{x}[n]e^{-j2\pi fn}\right|^{2}.

Also, let $\left\{\gamma_{k}\right\}_{k=0}^{N-1}$ and $\left\{\widetilde{\gamma}_{k}\right\}_{k=0}^{N-1}$ be real coefficients, and then define a multitaper spectral estimate $\widehat{V}(f)$ and an approximation $\widetilde{V}(f)$ by

\widehat{V}(f)=\sum_{k=0}^{N-1}\gamma_{k}V_{k}(f)\quad\text{and}\quad\widetilde{V}(f)=\sum_{k=0}^{N-1}\widetilde{\gamma}_{k}V_{k}(f).

Then, for any frequency $f\in\mathbb{R}$ , we have

\left|\widehat{V}(f)-\widetilde{V}(f)\right|\leq\left(\max_{k}\left|\gamma_{k}-\widetilde{\gamma}_{k}\right|\right)\|x\|_{2}^{2}.

Proof.

Let $\boldsymbol{V}=\begin{bmatrix}\boldsymbol{v}_{0}&\cdots&\boldsymbol{v}_{N-1}\end{bmatrix}$ , and let $\boldsymbol{\Gamma}$ , $\widetilde{\boldsymbol{\Gamma}}\in\mathbb{R}^{N\times N}$ , and $\boldsymbol{E}_{f}\in\mathbb{C}^{N\times N}$ be diagonal matrices whose diagonal entries are $\boldsymbol{\Gamma}[n,n]=\gamma_{n}$ , $\widetilde{\boldsymbol{\Gamma}}[n,n]=\widetilde{\gamma}_{n}$ , and $\boldsymbol{E}_{f}[n,n]=e^{j2\pi fn}$ for $n\in[N]$ . Then,

	$\displaystyle\widehat{V}(f)$	$\displaystyle=\sum_{k=0}^{N-1}\gamma_{k}V_{k}(f)$
		$\displaystyle=\sum_{k=0}^{N-1}\gamma_{k}\left\|\sum_{n=0}^{N-1}\boldsymbol{v}_{k}[n]\boldsymbol{x}[n]e^{-j2\pi fn}\right\|^{2}$
		$\displaystyle=\sum_{k=0}^{N-1}\gamma_{k}\left\|\boldsymbol{v}_{k}^{}\boldsymbol{E}_{f}^{}\boldsymbol{x}\right\|^{2}$
		$\displaystyle=\sum_{k=0}^{N-1}\boldsymbol{\Gamma}[k,k]\left\|(\boldsymbol{V}^{}\boldsymbol{E}_{f}^{}\boldsymbol{x})[k]\right\|^{2}$
		$\displaystyle=\boldsymbol{x}^{}\boldsymbol{E}_{f}\boldsymbol{V}\boldsymbol{\Gamma}\boldsymbol{V}^{}\boldsymbol{E}_{f}^{*}\boldsymbol{x}.$

In a nearly identical manner,

\widetilde{V}(f)=\boldsymbol{x}^{*}\boldsymbol{E}_{f}\boldsymbol{V}\widetilde{\boldsymbol{\Gamma}}\boldsymbol{V}^{*}\boldsymbol{E}_{f}^{*}\boldsymbol{x}.

Since $\boldsymbol{V}$ is orthonormal, $\|\boldsymbol{V}\|=\|\boldsymbol{V}^{*}\|=1$ . Since $\boldsymbol{E}_{f}$ is diagonal, and all the diagonal entries have modulus $1$ , $\|\boldsymbol{E}_{f}\|=\|\boldsymbol{E}_{f}^{*}\|=1$ . Hence, for any $f\in\mathbb{R}$ , we can bound

	$\displaystyle\left\|\widehat{V}(f)-\widetilde{V}(f)\right\|$	$\displaystyle=\left\|\boldsymbol{x}^{}\boldsymbol{E}_{f}\boldsymbol{V}\left(\boldsymbol{\Gamma}-\widetilde{\boldsymbol{\Gamma}}\right)\boldsymbol{V}^{}\boldsymbol{E}_{f}^{*}\boldsymbol{x}\right\|$
		$\displaystyle\leq\\|\boldsymbol{x}\\|_{2}\\|\boldsymbol{E}_{f}\\|\\|\boldsymbol{V}\\|\\|\boldsymbol{\Gamma}-\widetilde{\boldsymbol{\Gamma}}\\|\\|\boldsymbol{V}^{}\\|\\|\boldsymbol{E}_{f}^{}\\|\\|\boldsymbol{x}\\|_{2}$
		$\displaystyle=\\|\boldsymbol{\Gamma}-\widetilde{\boldsymbol{\Gamma}}\\|\\|\boldsymbol{x}\\|_{2}^{2}$
		$\displaystyle=\left(\max_{k}\|\gamma_{k}-\widetilde{\gamma}_{k}\|\right)\\|x\\|_{2}^{2},$

as desired. ∎

B.3 Proof of Theorem 6

Recall that the indices $[N]$ are partitioned as follows:

	$\displaystyle\mathcal{I}_{1}$	$\displaystyle=\{k\in[K]:\lambda_{k}\geq 1-\epsilon\}$
	$\displaystyle\mathcal{I}_{2}$	$\displaystyle=\{k\in[K]:\epsilon<\lambda_{k}<1-\epsilon\}$
	$\displaystyle\mathcal{I}_{3}$	$\displaystyle=\{k\in[N]\setminus[K]:\epsilon<\lambda_{k}<1-\epsilon\}$
	$\displaystyle\mathcal{I}_{4}$	$\displaystyle=\{k\in[N]\setminus[K]:\lambda_{k}\leq\epsilon\}.$

Using the partitioning above, the unweighted multitaper spectral estimate $\widehat{S}^{\text{mt}}_{K}(f)$ can be written as

	$\displaystyle\widehat{S}^{\text{mt}}_{K}(f)$	$\displaystyle=\dfrac{1}{K}\sum_{k=0}^{K-1}\widehat{S}_{k}(f)$
		$\displaystyle=\sum_{k\in\mathcal{I}_{1}\cup\mathcal{I}_{2}}\dfrac{1}{K}\widehat{S}_{k}(f),$

and the approximate estimator $\widetilde{S}^{\text{mt}}_{K}(f)$ can be written as

	$\displaystyle\widetilde{S}^{\text{mt}}_{K}(f)$	$\displaystyle=\dfrac{1}{K}\Psi(f)+\dfrac{1}{K}\sum_{k\in\mathcal{I}_{2}}(1-\lambda_{k})\widehat{S}_{k}(f)-\dfrac{1}{K}\sum_{k\in\mathcal{I}_{3}}\lambda_{k}\widehat{S}_{k}(f)$
		$\displaystyle=\sum_{k=0}^{N-1}\dfrac{\lambda_{k}}{K}\widehat{S}_{k}(f)+\sum_{k\in\mathcal{I}_{2}}\dfrac{1-\lambda_{k}}{K}\widehat{S}_{k}(f)-\sum_{k\in\mathcal{I}_{3}}\dfrac{\lambda_{k}}{K}\widehat{S}_{k}(f)$
		$\displaystyle=\sum_{k\in\mathcal{I}_{1}\cup\mathcal{I}_{4}}\dfrac{\lambda_{k}}{K}\widehat{S}_{k}(f)+\sum_{k\in\mathcal{I}_{2}}\dfrac{1}{K}\widehat{S}_{k}(f)$

Thus, $\widehat{S}^{\text{mt}}_{K}(f)$ and $\widetilde{S}^{\text{mt}}_{K}(f)$ can be written as

\widehat{S}^{\text{mt}}_{K}(f)=\sum_{k=0}^{N-1}\gamma_{k}\widehat{S}_{k}(f)\quad\text{and}\quad\widetilde{S}^{\text{mt}}_{K}(f)=\sum_{k=0}^{N-1}\widetilde{\gamma}_{k}\widehat{S}_{k}(f)

where

\gamma_{k}=\begin{cases}1/K&k\in\mathcal{I}_{1}\cup\mathcal{I}_{2},\\ 0&k\in\mathcal{I}_{3}\cup\mathcal{I}_{4},\end{cases}\quad\text{and}\quad\widetilde{\gamma}_{k}=\begin{cases}\lambda_{k}/K&k\in\mathcal{I}_{1}\cup\mathcal{I}_{4},\\ 1/K&k\in\mathcal{I}_{2},\\ 0&k\in\mathcal{I}_{3}.\end{cases}

We now bound $\left|\gamma_{k}-\widetilde{\gamma}_{k}\right|$ for all $k\in[N]$ . For $k\in\mathcal{I}_{1}$ , we have $\lambda_{k}\geq 1-\epsilon$ , and thus,

\left|\gamma_{k}-\widetilde{\gamma}_{k}\right|=\left|\dfrac{1}{K}-\dfrac{\lambda_{k}}{K}\right|=\dfrac{1-\lambda_{k}}{K}\leq\dfrac{\epsilon}{K}.

For $k\in\mathcal{I}_{2}\cup\mathcal{I}_{3}$ we have $\gamma_{k}=\widetilde{\gamma}_{k}$ , i.e., $\left|\gamma_{k}-\widetilde{\gamma}_{k}\right|=0$ . For $k\in\mathcal{I}_{4}$ , we have $\lambda_{k}\leq\epsilon$ , and thus,

\left|\gamma_{k}-\widetilde{\gamma}_{k}\right|=\left|0-\dfrac{\lambda_{k}}{K}\right|=\dfrac{\lambda_{k}}{K}\leq\dfrac{\epsilon}{K}.

Hence, $\left|\gamma_{k}-\widetilde{\gamma}_{k}\right|\leq\frac{\epsilon}{K}$ for all $k\in[N]$ , and thus by Lemma 9,

\left|\widehat{S}^{\text{mt}}_{K}(f)-\widetilde{S}^{\text{mt}}_{K}(f)\right|\leq\dfrac{\epsilon}{K}\|\boldsymbol{x}\|_{2}^{2}

for all frequencies $f\in\mathbb{R}$ .

B.4 Proof of Theorem 7

To evaluate the approximate multitaper estimate

\widetilde{S}^{\text{mt}}_{K}(f)=\dfrac{1}{K}\Psi(f)+\dfrac{1}{K}\sum_{k\in\mathcal{I}_{2}}(1-\lambda_{k})\widehat{S}_{k}(f)-\dfrac{1}{K}\sum_{k\in\mathcal{I}_{3}}\lambda_{k}\widehat{S}_{k}(f)

at the $L$ grid frequencies $f\in[L]/L$ one needs to:

•

For each $k\in\mathcal{I}_{2}\cup\mathcal{I}_{3}$ , precompute the Slepian basis vectors $\boldsymbol{s}_{k}$ and eigenvalues $\lambda_{k}$ .

Computing the Slepian basis vector $\boldsymbol{s}_{k}$ and the corresponding eigenvalue $\lambda_{k}$ for a single index $k$ can be done in $O(N\log N)$ operations and $O(N)$ memory via the method described in [53]. This needs to be done for $\#(\mathcal{I}_{2}\cup\mathcal{I}_{3})=\#\{k:\epsilon<\lambda_{k}<1-\epsilon\}=O(\log(NW)\log\tfrac{1}{\epsilon})$ values of $k$ , so the total cost of this step is $O(N\log(N)\log(NW)\log\tfrac{1}{\epsilon})$ operations and $O(N\log(NW)\log\tfrac{1}{\epsilon})$ memory.
•

For $\ell\in[L]$ , evaluate $\Psi(\tfrac{\ell}{L})$ .

If $L\geq 2N$ , then evaluating $\Psi(\tfrac{\ell}{L})$ for $\ell\in[L]$ can be done in $O(L\log L)$ operations and $O(L)$ memory as shown in Lemma 8. If $N\leq L<2N$ , then $2L\geq 2N$ , so by Lemma 8, we can evaluate $\Psi(\tfrac{\ell}{2L})$ for $\ell\in[2L]$ in $O(2L\log 2L)=O(L\log L)$ operations and $O(2L)=O(L)$ memory and then simply downsample the result to obtain $\Psi(\tfrac{\ell}{L})$ for $\ell\in[L]$ .
•

For each $k\in\mathcal{I}_{2}\cup\mathcal{I}_{3}$ and each $\ell\in[L]$ , evaluate $\widehat{S}_{k}(\tfrac{\ell}{L})$ .

Evaluating $\widehat{S}_{k}(\tfrac{\ell}{L})=\left|\sum_{n=0}^{N-1}\boldsymbol{s}_{k}[n]\boldsymbol{x}[n]e^{-j2\pi n\ell/L}\right|^{2}$ for all $\ell\in[L]$ can be done by pointwise multiplying $\boldsymbol{s}_{k}$ and $\boldsymbol{x}$ , zeropadding this vector to length $L$ , computing a length- $L$ FFT, and then computing the squared magnitude of each FFT coefficient. This takes $O(L\log L)$ operations and $O(L)$ memory. This needs to be done for $\#(\mathcal{I}_{2}\cup\mathcal{I}_{3})=O(\log(NW)\log\tfrac{1}{\epsilon})$ values of $k$ , so the total cost of this step is $O(L\log L\log(NW)\log\tfrac{1}{\epsilon})$ operations and $O(L\log(NW)\log\tfrac{1}{\epsilon})$ memory.
•

For each $\ell\in[L]$ , evaluate the weighted sum above for $\widetilde{S}^{\text{mt}}_{K}(\tfrac{\ell}{L})$ .

Once $\Psi(\tfrac{\ell}{L})$ and $\widehat{S}_{k}(\tfrac{\ell}{L})$ for $k\in\mathcal{I}_{2}\cup\mathcal{I}_{3}$ are computed, evaluating $\widetilde{S}^{\text{mt}}_{K}(\tfrac{\ell}{L})$ requires $O(\#(\mathcal{I}_{2}\cup\mathcal{I}_{3}))=O(\log(NW)\log\tfrac{1}{\epsilon})$ multiplications and additions. This has to be done for each $\ell\in[L]$ , so the total cost is $O(L\log(NW)\log\tfrac{1}{\epsilon})$ operations.

Since $L\geq N$ , the total cost of evaluating the approximate multitaper estimate $\widetilde{S}^{\text{mt}}_{K}(f)$ at the $L$ grid frequencies $f\in[L]/L$ is $O(L\log L\log(NW)\log\tfrac{1}{\epsilon})$ operations and $O(L\log(NW)\log\tfrac{1}{\epsilon})$ memory.

	$\displaystyle\left\\|\text{proj}_{\mathcal{S}_{f}}(\boldsymbol{x})\right\\|_{2}^{2}$	$\displaystyle=\left\\|\boldsymbol{E}_{f}\boldsymbol{S}_{K}\boldsymbol{S}_{K}^{}\boldsymbol{E}_{f}^{}\boldsymbol{x}\right\\|_{2}^{2}$
		$\displaystyle=\left\\|\boldsymbol{S}_{K}^{}\boldsymbol{E}_{f}^{}\boldsymbol{x}\right\\|_{2}^{2}$
		$\displaystyle=\sum_{k=0}^{K-1}\left\|(\boldsymbol{S}_{K}^{}\boldsymbol{E}_{f}^{}\boldsymbol{x})[k]\right\|^{2}$
		$\displaystyle=\sum_{k=0}^{K-1}\left\|\boldsymbol{s}_{k}^{*}\boldsymbol{E}_{f}\boldsymbol{x}\right\|^{2}$
		$\displaystyle=\sum_{k=0}^{K-1}\left\|\sum_{n=0}^{N-1}\boldsymbol{s}_{k}[n]\boldsymbol{x}[n]e^{-j2\pi fn}\right\|^{2}.$

	$\displaystyle\operatorname{Cov}\left[\left\\|\boldsymbol{U}\boldsymbol{x}\right\\|_{2}^{2},\left\\|\boldsymbol{V}\boldsymbol{x}\right\\|_{2}^{2}\right]$	$\displaystyle=\mathbb{E}\left[\left\\|\boldsymbol{U}\boldsymbol{x}\right\\|_{2}^{2}\cdot\left\\|\boldsymbol{V}\boldsymbol{x}\right\\|_{2}^{2}\right]-\mathbb{E}\left[\left\\|\boldsymbol{U}\boldsymbol{x}\right\\|_{2}^{2}\right]\mathbb{E}\left[\left\\|\boldsymbol{V}\boldsymbol{x}\right\\|_{2}^{2}\right]$
		$\displaystyle=\mathbb{E}\left[\boldsymbol{x}^{}\boldsymbol{U}^{}\boldsymbol{U}\boldsymbol{x}\boldsymbol{x}^{}\boldsymbol{V}^{}\boldsymbol{V}\boldsymbol{x}\right]-\mathbb{E}\left[\boldsymbol{x}^{}\boldsymbol{U}^{}\boldsymbol{U}\boldsymbol{x}\right]\mathbb{E}\left[\boldsymbol{x}^{}\boldsymbol{V}^{}\boldsymbol{V}\boldsymbol{x}\right]$
		$\displaystyle=\mathbb{E}\left[\boldsymbol{x}^{}\boldsymbol{V}^{}\otimes\boldsymbol{x}^{}\boldsymbol{U}^{}\right]\mathbb{E}\left[\boldsymbol{V}\boldsymbol{x}\otimes\boldsymbol{U}\boldsymbol{x}\right]+\mathbb{E}\left[\boldsymbol{x}^{}\boldsymbol{U}^{}\mathbb{E}\left[\boldsymbol{U}\boldsymbol{x}\boldsymbol{x}^{}\boldsymbol{V}^{}\right]\boldsymbol{V}\boldsymbol{x}\right]$
		$\displaystyle\ \ \ \ \ +2\mathbb{E}[\boldsymbol{x}^{}\boldsymbol{U}^{}]\mathbb{E}[\boldsymbol{U}\boldsymbol{x}]\mathbb{E}[\boldsymbol{x}^{}\boldsymbol{V}^{}]\mathbb{E}[\boldsymbol{V}\boldsymbol{x}].$

	$\displaystyle E\left[\boldsymbol{x}^{}\boldsymbol{U}^{}\mathbb{E}\left[\boldsymbol{U}\boldsymbol{x}\boldsymbol{x}^{}\boldsymbol{V}^{}\right]\boldsymbol{V}\boldsymbol{x}\right]$	$\displaystyle=\mathbb{E}\left[\operatorname{tr}\left[\mathbb{E}\left[\boldsymbol{U}\boldsymbol{x}\boldsymbol{x}^{}\boldsymbol{V}^{}\right]\boldsymbol{V}\boldsymbol{x}\boldsymbol{x}^{}\boldsymbol{U}^{}\right]\right]$
		$\displaystyle=\operatorname{tr}\left[\mathbb{E}\left[\mathbb{E}\left[\boldsymbol{U}\boldsymbol{x}\boldsymbol{x}^{}\boldsymbol{V}^{}\right]\boldsymbol{V}\boldsymbol{x}\boldsymbol{x}^{}\boldsymbol{U}^{}\right]\right]$
		$\displaystyle=\operatorname{tr}\left[\boldsymbol{U}\mathbb{E}\left[\boldsymbol{x}\boldsymbol{x}^{}\right]\boldsymbol{V}^{}\boldsymbol{V}\mathbb{E}\left[\boldsymbol{x}\boldsymbol{x}^{}\right]\boldsymbol{U}^{}\right]$
		$\displaystyle=\operatorname{tr}\left[\boldsymbol{U}\boldsymbol{R}\boldsymbol{V}^{}\boldsymbol{V}\boldsymbol{R}\boldsymbol{U}^{}\right]$
		$\displaystyle=\left\\|\boldsymbol{U}\boldsymbol{R}\boldsymbol{V}^{*}\right\\|_{F}^{2}.$

	$\displaystyle\operatorname{Cov}\left[\left\\|\boldsymbol{U}\boldsymbol{x}\right\\|_{2}^{2},\left\\|\boldsymbol{V}\boldsymbol{x}\right\\|_{2}^{2}\right]$	$\displaystyle=\mathbb{E}\left[\boldsymbol{x}^{}\boldsymbol{V}^{}\otimes\boldsymbol{x}^{}\boldsymbol{U}^{}\right]\mathbb{E}\left[\boldsymbol{V}\boldsymbol{x}\otimes\boldsymbol{U}\boldsymbol{x}\right]+\mathbb{E}\left[\boldsymbol{x}^{}\boldsymbol{U}^{}\mathbb{E}\left[\boldsymbol{U}\boldsymbol{x}\boldsymbol{x}^{}\boldsymbol{V}^{}\right]\boldsymbol{V}\boldsymbol{x}\right]$
		$\displaystyle\ \ \ \ \ +2\mathbb{E}[\boldsymbol{x}^{}\boldsymbol{U}^{}]\mathbb{E}[\boldsymbol{U}\boldsymbol{x}]\mathbb{E}[\boldsymbol{x}^{}\boldsymbol{V}^{}]\mathbb{E}[\boldsymbol{V}\boldsymbol{x}].$
		$\displaystyle=\left\\|\boldsymbol{U}\boldsymbol{R}\boldsymbol{V}^{*}\right\\|_{F}^{2}.$

	$\displaystyle\operatorname{Bias}\left[\widehat{S}^{\text{mt}}_{K}(f)\right]$	$\displaystyle=\left\|\mathbb{E}\widehat{S}^{\text{mt}}_{K}(f)-S(f)\right\|$
		$\displaystyle=\left\|\int_{-1/2}^{1/2}S(f-f^{\prime})\psi(f^{\prime})\,df^{\prime}-S(f)\right\|$
		$\displaystyle=\left\|\int_{-W}^{W}S(f-f^{\prime})\psi(f^{\prime})\,df^{\prime}+\int_{\Omega}S(f-f^{\prime})\psi(f^{\prime})\,df^{\prime}-S(f)\right\|$
		$\displaystyle\leq\underbrace{\left\|\int_{-W}^{W}S(f-f^{\prime})\psi(f^{\prime})\,df^{\prime}-S(f)\right\|}_{\text{local bias}}+\underbrace{\left\|\int_{\Omega}S(f-f^{\prime})\psi(f^{\prime})\,df^{\prime}\right\|}_{\text{broadband bias}}.$

Thomson’s Multitaper Method Revisited

Abstract

1 Introduction

2 Thomson’s Multitaper Method

2.1 Traditional view

2.2 Linear algebra view

2.3 Slepian basis eigenvalues

3 Statistical Properties and Spectral Leakage

Theorem 1.

Theorem 2.

Theorem 3.

Theorem 4.

Theorem 5.

4 Fast Algorithms

Theorem 6.

Theorem 7.

5 Simulations

5.1 Spectral leakage

5.2 Comparison of spectral estimation methods

5.3 Fast multitaper spectral estimation

6 Conclusions

References

Appendix A Proof of Results in Section 3

A.1 Norms of Gaussian random variables

Lemma 1.

Lemma 2.

Proof.

A.2 Concentration of norms of Gaussian random variables

Lemma 3.

Lemma 4.

Proof.

A.3 Intermediate results

Lemma 5.

Proof.

Lemma 6.

Proof.

Lemma 7.

Proof.

A.4 Proof of Theorem 1

A.5 Proof of Theorem 2

A.6 Proof of Theorem 3

A.7 Proof of Theorem 4

A.8 Proof of Theorem 5

Appendix B Proof of Results in Section 4

B.1 Fast algorithm for computing Ψ​(f)\Psi(f) at grid frequencies

Lemma 8.

Proof.

B.2 Approximations for general multitaper spectral estimates

Lemma 9.

Proof.

B.3 Proof of Theorem 6

B.4 Proof of Theorem 7

B.1 Fast algorithm for computing $\Psi(f)$ at grid frequencies