Multitaper

In signal processing, multitaper analysis is a spectral density estimation technique developed by David J. Thomson. It can estimate the power spectrum SX of a stationary ergodic finite-variance random process X, given a finite contiguous realization of X as data.

Motivation The multitaper method overcomes some of the limitations of non-parametric Fourier analysis. When applying the Fourier transform to extract spectral information from a signal, we assume that each Fourier coefficient is a reliable representation of the amplitude and relative phase of the corresponding component frequency. This assumption, however, is not generally valid for empirical data. For instance, a single trial represents only one noisy realization of the underlying process of interest. A comparable situation arises in statistics when estimating measures of central tendency i.e., it is bad practice to estimate qualities of a population using individuals or very small samples. Likewise, a single sample of a process does not necessarily provide a reliable estimate of its spectral properties. Moreover, the naive power spectral density obtained from the signal's raw Fourier transform is a biased estimate of the true spectral content.

These problems are often overcome by averaging over many realizations of the same event after applying a taper to each trial. However, this method is unreliable with small data sets and undesirable when one does not wish to attenuate signal components that vary across trials. Furthermore, even when many trials are available the untapered periodogram is generally biased (with the exception of white noise) and the bias depends upon the length of each realization, not the number of realizations recorded. Applying a single taper reduces bias but at the cost of increased estimator variance due to attenuation of activity at the start and end of each recorded segment of the signal. The multitaper method partially obviates these problems by obtaining multiple independent estimates from the same sample. Each data taper is multiplied element-wise by the signal to provide a windowed trial from which one estimates the power at each component frequency. As each taper is pairwise orthogonal to all other tapers, the window functions are uncorrelated with one another. The final spectrum is obtained by averaging over all the tapered spectra thus recovering some of the information that is lost due to partial attenuation of the signal that results from applying individual tapers. This method is especially useful when a small number of trials is available as it reduces the estimator variance beyond what is possible with single taper methods. Moreover, even when many trials are available the multitaper approach is useful as it permits more rigorous control of the trade-off between bias and variance than what is possible in the single taper case. Thomson chose the Slepian functions or discrete prolate spheroidal sequences as tapers since these vectors are mutually orthogonal and possess desirable spectral concentration properties (see the section on Slepian sequences). In practice, a weighted average is often used to compensate for increased energy loss at higher order tapers.

Formulation Consider a p-dimensional zero mean stationary stochastic process

X ( t ) = [ X ( 1 , t ) , X ( 2 , t ) , … , X ( p , t ) ] T {\displaystyle \mathbf {X} (t)={\lbrack X(1,t),X(2,t),\dots ,X(p,t)\rbrack }^{T}}

Here T denotes the matrix transposition. In neurophysiology for example, p refers to the total number of channels and hence X ( t ) {\displaystyle \mathbf {X} (t)} can represent simultaneous measurement of electrical activity of those p channels. Let the sampling interval between observations be Δ t {\displaystyle \Delta t} , so that the Nyquist frequency is f N = 1 / ( 2 Δ t ) {\displaystyle f_{N}=1/(2\Delta t)} . The multitaper spectral estimator utilizes several different data tapers which are orthogonal to each other. The multitaper cross-spectral estimator between channel l and m is the average of K direct cross-spectral estimators between the same pair of channels (l and m) and hence takes the form

S ^ l m ( f ) = 1 K ∑ k = 0 K − 1 S ^ k l m ( f ) . {\displaystyle {\hat {S}}^{lm}(f)={\frac {1}{K}}\sum _{k=0}^{K-1}{\hat {S}}_{k}^{lm}(f).}

Here, S ^ k l m ( f ) {\displaystyle {\hat {S}}_{k}^{lm}(f)} (for 0 ≤ k ≤ K − 1 {\displaystyle 0\leq k\leq K-1} ) is the kth direct cross spectral estimator between channel l and m and is given by

… excerpt ends here. Continue reading the full article.