# Multitaper

> Mediated Wiki article. Canonical URL: https://mediated.wiki/source/Multitaper
> Markdown URL: https://mediated.wiki/source/Multitaper.md
> Source: https://en.wikipedia.org/wiki/Multitaper
> Source revision: 1331444903
> License: Creative Commons Attribution-ShareAlike 4.0 International (https://creativecommons.org/licenses/by-sa/4.0/)

{{Short description|Spectral density estimation technique}}
[[Image:MultiTaper Plot.JPG|thumb|Comparison of [periodogram](/source/periodogram) (black) and multitaper estimate (red) of a single trial local field potential measurement.  This estimate used 9 tapers.]]

In [signal processing](/source/signal_processing), '''multitaper''' analysis is a [spectral density estimation](/source/spectral_density_estimation) technique developed by [David J. Thomson](/source/David_J._Thomson).<ref>Thomson, D. J. (1982) ''Spectrum estimation and harmonic analysis''. ''Proceedings of the IEEE'', 70, 1055&ndash;1096</ref><ref>*{{Citation | last1=Press | first1=WH | last2=Teukolsky | first2=SA | last3=Vetterling | first3=WT | last4=Flannery | first4=BP | year=2007 | title=Numerical Recipes: The Art of Scientific Computing | edition=3rd | publisher=Cambridge University Press |  publication-place=New York | isbn=978-0-521-88068-8 | chapter=Section 13.4.3. ''Multitaper methods and Slepian functions''}}</ref> It can [estimate](/source/Estimation_theory) the [power spectrum](/source/power_spectrum) ''S''<sub>''X''</sub> of a [stationary](/source/stationary_process) [ergodic](/source/ergodic) finite-variance [random process](/source/random_process) ''X'', given a finite contiguous [realization](/source/realization_(probability)) of ''X'' as data.

==Motivation==

The multitaper method overcomes some of the limitations of non-parametric [Fourier analysis](/source/Fourier_analysis).  When applying the [Fourier transform](/source/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](/source/mean) 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](/source/power_spectral_density) obtained from the signal's raw Fourier transform is a [biased](/source/estimator_bias) estimate of the true spectral content.

thumb|The importance of averaging in (cross-)spectral density estimation.<ref>{{Cite journal | last1 = Simons | first1 = F. J. | last2 = Korenaga | first2 = J. | last3 = Zuber | first3 = M. T. | title = Isostatic response of the Australian lithosphere: Estimation of effective elastic thickness and anisotropy using multitaper analysis | journal =   | volume = 105 | issue = B8 | pages=19163–19184 | year = 2000}}</ref> (a) Synthetically generated noisy signal with two coherent frequencies at 0.03 and 0.6 Hz. (b) Multitaper (MT) spectral density estimates. (c) Coherence squared estimates using Slepian multitaper analysis (thick line, unshaded) and Welch overlapping segment analysis (WOSA) (thin line, shaded area). (d) Estimate of the phase of the cross-spectral density estimate using MT (solid) and WOSA (dashed). At 0.03 Hz the signals are in phase, while at 0.6 Hz they are <math>-\pi/4</math> out of phase.

These problems are often overcome by averaging over many realizations of the same event after applying a [taper](/source/Window_function) 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](/source/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](/source/Window_function) 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<ref>{{Cite book | last1 = Simons | first1 = F. J. | last2 = Plattner | first2 = A. | chapter = Scalar and Vector Slepian Functions, Spherical Signal Estimation and Spectral Analysis | title = Handbook of Geomathematics | year = 2015 | pages = 2563–2608 | doi = 10.1007/978-3-642-54551-1_30| isbn = 978-3-642-54550-4 }}
</ref> or discrete prolate spheroidal sequences as tapers since these vectors are mutually orthogonal and possess desirable [spectral concentration](/source/Spectral_concentration_problem) properties (see the section on Slepian sequences).  In practice, a [weighted average](/source/weighted_average) is often used to compensate for increased energy loss at higher order tapers.<ref>Percival, D. B., and A. T. Walden. ''Spectral Analysis for Physical Applications: Multitaper and Conventional Univariate Techniques''. Cambridge: Cambridge University Press, 1993.</ref>

==Formulation==
Consider a p-dimensional zero mean [stationary stochastic process](/source/Stationary_process)

: <math>\mathbf{X}(t) = {\lbrack X(1,t), X(2,t), \dots , X(p,t)
\rbrack}^T</math>

Here ''T'' denotes the matrix transposition.  In [neurophysiology](/source/neurophysiology) for example, ''p'' refers to the total number of channels and
hence <math>\mathbf{X}(t)</math> can represent simultaneous measurement of
electrical activity of those ''p'' channels. Let the sampling interval
between observations be <math> \Delta t</math>, so that the [Nyquist frequency](/source/Nyquist_frequency) is <math>f_N=1/(2 \Delta t)</math>.

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

: <math>\hat{S}^{lm} (f)= \frac{1}{K} \sum_{k=0}^{K-1} \hat{S}_k^{lm}(f).</math>

Here, <math>\hat{S}_{k}^{lm}(f)</math> (for <math>0 \leq k \leq K-1</math>) is the ''k''<sup>th</sup> direct cross spectral estimator between channel ''l'' and ''m'' and is given by

: <math> \hat{S}_{k}^{lm}(f) = \frac{1}{N\Delta t}
         {\lbrack J_{k}^{l}(f) \rbrack}^{*} {\lbrack J_{k}^{m}(f)
         \rbrack}, 
</math>

where

: <math>J_k^l(f) = \sum_{t=1}^N h_{t,k}X(l,t) e^{-i 2\pi ft\Delta t}.
</math>

thumb|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 <math>\lbrace h_{t,k} \rbrace </math> is the data taper for the
''k''<sup>th</sup> direct cross-spectral estimator <math>\hat{S}_k^{lm}(f)</math> 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 function](/source/Slepian_function)s or [discrete prolate spheroidal sequences](/source/Window_function),<ref>Slepian, D. (1978) "Prolate spheroidal wave functions, Fourier analysis, and uncertainty &ndash; V: The discrete case." ''Bell System Technical Journal'', 57, 1371&ndash;1430</ref> after [David Slepian](/source/David_Slepian) (also known in literature as discrete prolate spheroidal sequences or DPSS for short) with parameter ''W'' and orders ''k''&nbsp;=&nbsp;0 to ''K''&nbsp;&minus;&nbsp;1.  The maximum order ''K'' is chosen to be less than the [Shannon number](/source/Shannon_number) <math>2NW\Delta t</math>. The quantity 2''W'' defines the resolution bandwidth for the [spectral concentration problem](/source/spectral_concentration_problem) and <math>W \in
(0,f_{N})</math>. When ''l''&nbsp;=&nbsp;''m'', we get the multitaper estimator for the auto-spectrum of the ''l''<sup>th</sup> channel. In recent years, a dictionary based on modulated DPSS was proposed as an overcomplete alternative to DPSS.<ref>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.</ref>

==Applications==
Not limited to time series, the multitaper method is easily extensible to multiple Cartesian dimensions using custom Slepian functions,<ref>{{Cite journal |  last1= Simons | first1 = F. J. | last2 = Wang  | first2= D. V. | title = Spatiospectral concentration in the Cartesian plane | journal =    GEM: International Journal on Geomathematics| year =  2011 | volume = 2 | pages = 1–36 | doi =10.1007/s13137-011-0016-z| arxiv = 1007.5226 }}.</ref> and can be reformulated for spectral estimation on the sphere using Slepian functions constructed from [spherical harmonics](/source/spherical_harmonics)<ref>{{Cite journal | last1 = Simons | first1 = F. J. | last2 = Dahlen | first2 = F. A. | last3 = Wieczorek | first3 = M. A. | doi = 10.1137/S0036144504445765 | title = Spatiospectral Concentration on a Sphere | journal = SIAM Review | volume = 48 | issue = 3 | pages = 504–536 | year = 2006 | arxiv = math/0408424 | bibcode = 2006SIAMR..48..504S }}</ref> for applications in [geophysics](/source/geophysics) and [cosmology](/source/cosmology)<ref>{{Cite journal | last1 = Wieczorek | first1 = M. A. | last2 = Simons | first2 = F. J. | doi = 10.1007/s00041-006-6904-1 | title = Minimum-variance multitaper spectral estimation on the sphere | journal = Journal of Fourier Analysis and Applications | volume = 13 | issue = 6 | pages = 665 | year = 2007  | arxiv = 1306.3254 }}</ref><ref>{{Cite journal | last1 = Dahlen | first1 = F. A. | last2 = Simons | first2 = F. J. | doi = 10.1111/j.1365-246X.2008.03854.x | title = Spectral estimation on a sphere in geophysics and cosmology | journal = Geophysical Journal International | volume = 174 | issue = 3 | pages = 774 | year = 2008 | doi-access = free | arxiv = 0705.3083 | bibcode = 2008GeoJI.174..774D }}</ref> among others. An extensive treatment about the application of this method to analyze multi-trial, multi-channel data generated in [neuroscience](/source/neuroscience), [biomedical engineering](/source/biomedical_engineering) and elsewhere can be found [http://global.oup.com/academic/product/observed-brain-dynamics-9780195178081?q=9780195178081&lang=en&cc=us here]. This technique is currently used in the [spectral analysis](/source/frequency_spectrum) toolkit of [Chronux](/source/Chronux).

==See also==
*[Periodogram](/source/Periodogram)

== References ==
{{Reflist}}

== External links ==
* [https://github.com/pantheras/mtpsd] C++/Octave libraries for the multitaper method, including adaptive weighting (hosted on GitHub)
* [http://www.atmos.ucla.edu/tcd/ssa/guide/mann/mann.html] Documentation on the multitaper method from the SSA-MTM Toolkit implementation
* [https://www.mit.edu/~gprieto/mwlib.html] Fortran 90 library with additional multivariate applications
* [https://github.com/csdms-contrib/slepian_alpha] [MATLAB](/source/MATLAB) code base to generate ''spherical scalar'' Slepian functions
* [https://github.com/csdms-contrib/slepian_charlie] [MATLAB](/source/MATLAB) code base to perform ''spherical multitaper'' analysis
* [https://github.com/csdms-contrib/slepian_foxtrot] [MATLAB](/source/MATLAB) code base to generate ''Cartesian'' Slepian functions
* [https://github.com/csdms-contrib/slepian_golf] [MATLAB](/source/MATLAB) code base to generate ''spherical vector'' Slepian functions
* [https://code.google.com/p/pymutt/] Python module
* [https://cran.r-project.org/web/packages/multitaper/index.html] [R (programming language)](/source/R_(programming_language)) multitaper Package
* [https://github.com/krahim/dpss-splus-script] [S-Plus](/source/S-Plus) script to generate Slepian sequences (dpss)

Category:Frequency-domain analysis
Category:Signal processing
Category:Time–frequency analysis
Category:Signal estimation

---
Adapted from the Wikipedia article [Multitaper](https://en.wikipedia.org/wiki/Multitaper) by Wikipedia contributors ([contributor history](https://en.wikipedia.org/wiki/Multitaper?action=history)). Available under [Creative Commons Attribution-ShareAlike 4.0 International](https://creativecommons.org/licenses/by-sa/4.0/). Changes may have been made.