Skip to main content

MultitaperPsd

Namespace: RTPROCESSINGLIB  ·  Library: DSP Library

Python equivalent

mne.time_frequency.psd_array_multitaper in MNE-Python.

#include <dsp/multitaper_psd.h>

class UTILSLIB::MultitaperPsd

Multitaper power spectral density estimator using DPSS (Slepian) tapers.

Applies multiple orthogonal DPSS tapers to the data, computes the FFT of each tapered segment, and averages the resulting periodograms weighted by the taper eigenvalues. This provides a PSD estimate with reduced variance compared to a single-taper (periodogram) approach.

src/examples/ex_dsp_spectral/main.cpp
const MultitaperPsdResult psd = MultitaperPsd::compute(signal, kSFreq, 2.0); // NW = 2, i.e. +-1 Hz resolution
const double binWidth = psd.vecFreqs(1) - psd.vecFreqs(0);
const double power = psd.matPsd.row(0).sum() * binWidth; // one-sided PSD integrates to the variance

Static Methods​

compute(matData, sfreq, halfBandwidth, nTapers, nFft)​

Compute multitaper PSD for every channel of a data matrix.

Parameters:

  • matData : const Eigen::MatrixXd & Data matrix (n_channels × n_times).

  • sfreq : double Sampling frequency in Hz.

  • halfBandwidth : double Half-bandwidth parameter (NW); default 4.0.

  • nTapers : int Number of DPSS tapers; -1 → as MNE-Python: those of the first floor(2*halfBandwidth) with concentration > 0.9.

  • nFft : int FFT length; -1 → n_times.

Returns:

Example​

Used in examples​

ex_dsp_spectral​

Spectral and filtering tools of the dsp library on signals with known content.

src/examples/ex_dsp_spectral/main.cpp, line 134
const MultitaperPsdResult psd = MultitaperPsd::compute(signal, kSFreq, 2.0); // NW = 2, i.e. +-1 Hz resolution
const double binWidth = psd.vecFreqs(1) - psd.vecFreqs(0);
const double power = psd.matPsd.row(0).sum() * binWidth; // one-sided PSD integrates to the variance
// mne.time_frequency.psd_array_multitaper(bandwidth=4 Hz, adaptive=False): 20 Hz 0.31659, 55 Hz 0.079135
const double psd20 = psd.matPsd(0, 40);
const double psd55 = psd.matPsd(0, 110);
ok &= expect(frequencyOfMax(psd.matPsd.row(0), psd.vecFreqs) == 20.0 && std::fabs(psd20 - 0.316593412495819) < 1e-3 && std::fabs(psd55 - 0.07913494397010905) < 1e-3 && std::fabs(power - 0.625) < 0.01,
             QString("multitaper PSD at 20/55 Hz %1/%2 = mne, integral %3").arg(psd20).arg(psd55).arg(power));

Authors of this file​