Multitaper

Last updated
Comparison of periodogram (black) and multitaper estimate (red) of a single trial local field potential measurement. This estimate used 9 tapers. MultiTaper Plot.JPG
Comparison of periodogram (black) and multitaper estimate (red) of a single trial local field potential measurement. This estimate used 9 tapers.

In signal processing, multitaper is a spectral density estimation technique developed by David J. Thomson. [1] 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.

Contents

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 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. [2]

Formulation

Consider a p-dimensional zero mean stationary stochastic process

Here T denotes the matrix transposition. In neurophysiology for example, p refers to the total number of channels and hence can represent simultaneous measurement of electrical activity of those p channels. Let the sampling interval between observations be , so that the Nyquist frequency is .

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

Here, (for ) is the kth direct cross spectral estimator between channel l and m and is given by

where

The three leading Slepian sequences for T=1000 and 2WT=6. Note that each higher order sequence has an extra zero crossing. DPSS figure.svg
The three leading Slepian sequences for T=1000 and 2WT=6. Note that each higher order sequence has an extra zero crossing.

The Slepian sequences

The sequence is the data taper for the kth direct cross-spectral estimator and is chosen as follows:

We choose a set of K orthogonal data tapers such that each one provides a good protection against leakage. These are given by the Slepian sequences, [3] after David Slepian (also known in literature as discrete prolate spheroidal sequences or DPSS for short) with parameter W and orders k = 0 to K  1. The maximum order K is chosen to be less than the Shannon number . The quantity 2W defines the resolution bandwidth for the spectral concentration problem and . When l = m, we get the multitaper estimator for the auto-spectrum of the lth channel. In recent years, a dictionary based on modulated DPSS was proposed as an overcomplete alternative to DPSS. [4]

See also Window function:DPSS or Slepian window

Applications

This technique is currently used in the spectral analysis toolkit of Chronux. An extensive treatment about the application of this method to analyze multi-trial, multi-channel data generated in neuroscience experiments, biomedical engineering and others can be found here. Not limited to time series, the multitaper method can be reformulated for spectral estimation on the sphere using Slepian functions constructed from spherical harmonics [5] for applications in geophysics and cosmology [6] [7] among others.

See also

Related Research Articles

<span class="mw-page-title-main">Autocorrelation</span> Correlation of a signal with a time-shifted copy of itself, as a function of shift

Autocorrelation, sometimes known as serial correlation in the discrete time case, is the correlation of a signal with a delayed copy of itself as a function of delay. Informally, it is the similarity between observations of a random variable as a function of the time lag between them. The analysis of autocorrelation is a mathematical tool for finding repeating patterns, such as the presence of a periodic signal obscured by noise, or identifying the missing fundamental frequency in a signal implied by its harmonic frequencies. It is often used in signal processing for analyzing functions or series of values, such as time domain signals.

<span class="mw-page-title-main">Least squares</span> Approximation method in statistics

The method of least squares is a standard approach in regression analysis to approximate the solution of overdetermined systems by minimizing the sum of the squares of the residuals made in the results of each individual equation.

The Whittaker–Shannon interpolation formula or sinc interpolation is a method to construct a continuous-time bandlimited function from a sequence of real numbers. The formula dates back to the works of E. Borel in 1898, and E. T. Whittaker in 1915, and was cited from works of J. M. Whittaker in 1935, and in the formulation of the Nyquist–Shannon sampling theorem by Claude Shannon in 1949. It is also commonly called Shannon's interpolation formula and Whittaker's interpolation formula. E. T. Whittaker, who published it in 1915, called it the Cardinal series.

In statistics, point estimation involves the use of sample data to calculate a single value which is to serve as a "best guess" or "best estimate" of an unknown population parameter. More formally, it is the application of a point estimator to the data to obtain a point estimate.

<span class="mw-page-title-main">Spectral density</span> Relative importance of certain frequencies in a composite signal

The power spectrum of a time series describes the distribution of power into frequency components composing that signal. According to Fourier analysis, any physical signal can be decomposed into a number of discrete frequencies, or a spectrum of frequencies over a continuous range. The statistical average of a certain signal or sort of signal as analyzed in terms of its frequency content, is called its spectrum.

<span class="mw-page-title-main">Window function</span> Function used in signal processing

In signal processing and statistics, a window function is a mathematical function that is zero-valued outside of some chosen interval, normally symmetric around the middle of the interval, usually near a maximum in the middle, and usually tapering away from the middle. Mathematically, when another function or waveform/data-sequence is "multiplied" by a window function, the product is also zero-valued outside the interval: all that is left is the part where they overlap, the "view through the window". Equivalently, and in actual practice, the segment of data within the window is first isolated, and then only that data is multiplied by the window function values. Thus, tapering, not segmentation, is the main purpose of window functions.

In signal processing, a periodogram is an estimate of the spectral density of a signal. The term was coined by Arthur Schuster in 1898. Today, the periodogram is a component of more sophisticated methods. It is the most common tool for examining the amplitude vs frequency characteristics of FIR filters and window functions. FFT spectrum analyzers are also implemented as a time-sequence of periodograms.

<span class="mw-page-title-main">Cramér–Rao bound</span> Lower bound on variance of an estimator

In estimation theory and statistics, the Cramér–Rao bound (CRB) relates to estimation of a deterministic parameter. The result is named in honor of Harald Cramér and C. R. Rao, but has also been derived independently by Maurice Fréchet, Georges Darmois, and by Alexander Aitken and Harold Silverstone. It states that the precision of any unbiased estimator is at most the Fisher information; or (equivalently) the inverse of the Fisher information is a lower bound on its variance.

Importance sampling is a Monte Carlo method for evaluating properties of a particular distribution, while only having samples generated from a different distribution than the distribution of interest. Its introduction in statistics is generally attributed to a paper by Teun Kloek and Herman K. van Dijk in 1978, but its precursors can be found in statistical physics as early as 1949. Importance sampling is also related to umbrella sampling in computational physics. Depending on the application, the term may refer to the process of sampling from this alternative distribution, the process of inference, or both.

<span class="mw-page-title-main">Sensor array</span> Group of sensors used to increase gain or dimensionality over a single sensor

A sensor array is a group of sensors, usually deployed in a certain geometry pattern, used for collecting and processing electromagnetic or acoustic signals. The advantage of using a sensor array over using a single sensor lies in the fact that an array adds new dimensions to the observation, helping to estimate more parameters and improve the estimation performance. For example an array of radio antenna elements used for beamforming can increase antenna gain in the direction of the signal while decreasing the gain in other directions, i.e., increasing signal-to-noise ratio (SNR) by amplifying the signal coherently. Another example of sensor array application is to estimate the direction of arrival of impinging electromagnetic waves. The related processing method is called array signal processing. A third examples includes chemical sensor arrays, which utilize multiple chemical sensors for fingerprint detection in complex mixtures or sensing environments. Application examples of array signal processing include radar/sonar, wireless communications, seismology, machine condition monitoring, astronomical observations fault diagnosis, etc.

In mathematics, the discrete-time Fourier transform (DTFT), also called the finite Fourier transform, is a form of Fourier analysis that is applicable to a sequence of values.

Estimation theory is a branch of statistics that deals with estimating the values of parameters based on measured empirical data that has a random component. The parameters describe an underlying physical setting in such a way that their value affects the distribution of the measured data. An estimator attempts to approximate the unknown parameters using the measurements. In estimation theory, two approaches are generally considered:

<span class="mw-page-title-main">David Slepian</span> American mathematician

David S. Slepian was an American mathematician. He is best known for his work with algebraic coding theory, probability theory, and distributed source coding. He was colleagues with Claude Shannon and Richard Hamming at Bell Labs.

<span class="mw-page-title-main">Wavelet transform</span> Mathematical technique used in data compression and analysis

In mathematics, a wavelet series is a representation of a square-integrable function by a certain orthonormal series generated by a wavelet. This article provides a formal, mathematical definition of an orthonormal wavelet and of the integral wavelet transform.

In estimation theory and decision theory, a Bayes estimator or a Bayes action is an estimator or decision rule that minimizes the posterior expected value of a loss function. Equivalently, it maximizes the posterior expectation of a utility function. An alternative way of formulating an estimator within Bayesian statistics is maximum a posteriori estimation.

In statistical signal processing, the goal of spectral density estimation (SDE) or simply spectral estimation is to estimate the spectral density of a signal from a sequence of time samples of the signal. Intuitively speaking, the spectral density characterizes the frequency content of the signal. One purpose of estimating the spectral density is to detect any periodicities in the data, by observing peaks at the frequencies corresponding to these periodicities.

The prolate spheroidal wave functions are eigenfunctions of the Laplacian in prolate spheroidal coordinates, adapted to boundary conditions on certain ellipsoids of revolution. Related are the oblate spheroidal wave functions.

<span class="mw-page-title-main">Spectral concentration problem</span>

The spectral concentration problem in Fourier analysis refers to finding a time sequence of a given length whose discrete Fourier transform is maximally localized on a given frequency interval, as measured by the spectral concentration.

Multidimension spectral estimation is a generalization of spectral estimation, normally formulated for one-dimensional signals, to multidimensional signals or multivariate data, such as wave vectors.

In statistics, Whittle likelihood is an approximation to the likelihood function of a stationary Gaussian time series. It is named after the mathematician and statistician Peter Whittle, who introduced it in his PhD thesis in 1951. It is commonly used in time series analysis and signal processing for parameter estimation and signal detection.

References

  1. Thomson, D. J. (1982) "Spectrum estimation and harmonic analysis." Proceedings of the IEEE, 70, 10551096
  2. Percival, D. B., and A. T. Walden. Spectral Analysis for Physical Applications: Multitaper and Conventional Univariate Techniques. Cambridge: Cambridge University Press, 1993.
  3. Slepian, D. (1978) "Prolate spheroidal wave functions, Fourier analysis, and uncertainty V: The discrete case." Bell System Technical Journal, 57, 13711430
  4. E. Sejdić, M. Luccini, S. Primak, K. Baddour, T. Willink, “Channel estimation using modulated discrete prolate spheroidal sequences based frames,” in Proc. of IEEE International Conference on Acoustics, Speech, and Signal Processing (ICASSP 2008), Las Vegas, Nevada, USA, March 31-April 04, 2008, pp. 2849-2852.
  5. Simons, F. J.; Dahlen, F. A.; Wieczorek, M. A. (2006). "Spatiospectral Concentration on a Sphere". SIAM Review. 48 (3): 504–536. arXiv: math/0408424 . Bibcode:2006SIAMR..48..504S. doi:10.1137/S0036144504445765. S2CID   27519592.
  6. Wieczorek, M. A.; Simons, F. J. (2007). "Minimum-Variance Multitaper Spectral Estimation on the Sphere". Journal of Fourier Analysis and Applications. 13 (6): 665. arXiv: 1306.3254 . doi:10.1007/s00041-006-6904-1. S2CID   7664891.
  7. Dahlen, F. A.; Simons, F. J. (2008). "Spectral estimation on a sphere in geophysics and cosmology". Geophysical Journal International. 174 (3): 774. arXiv: 0705.3083 . Bibcode:2008GeoJI.174..774D. doi:10.1111/j.1365-246X.2008.03854.x. S2CID   119253186.