braintools.analysis.spectral.multitaper_psd

Contents

braintools.analysis.spectral.multitaper_psd#

braintools.analysis.spectral.multitaper_psd(x, dt, nw=4.0, n_tapers=None, *, remove_dc=True)#

Estimate the power spectral density with DPSS (Slepian) multitapers.

The signal is multiplied by n_tapers discrete prolate spheroidal sequences [1] with time-half-bandwidth product nw. The one-sided periodograms of the tapered signals are averaged with weights equal to the taper concentration ratios lambda_k (non-adaptive weighting, [2]):

\[S(f) = \frac{2}{f_s} \frac{\sum_k \lambda_k |X_k(f)|^2}{\sum_k \lambda_k},\]

where the factor 2 is not applied at DC and (even T) Nyquist.

Parameters:

  • x (Array | ndarray | bool | number | bool | int | float | complex | Quantity) – Real signal of shape (T, ...); time is axis 0. A brainunit.Quantity is accepted.

  • dt (float | Quantity) – Sampling interval (seconds if a plain number). Static.

  • nw (float) – Time-half-bandwidth product NW; the frequency resolution (full bandwidth) is 2 * nw / (T * dt) Hz. Static, 0.5 <= nw < T / 2.

  • n_tapers (int | None) – Number of tapers K, 1 <= K <= T. Static. Defaults to int(2 * nw) - 1 (the tapers with good spectral concentration).

  • remove_dc (bool) – Subtract the mean of each channel before tapering (as MNE does).

Return type:

Tuple[Array, Array | Quantity]

Returns:

  • freqs (jax.Array) – Frequencies in Hz, shape (T // 2 + 1,).

  • psd (jax.Array or brainunit.Quantity) – Power spectral density, shape (T // 2 + 1, ...). If x has unit U, the result has unit U**2 * second (U**2/Hz).

Raises:

ValueError – If nw or n_tapers is out of range, or dt is not positive.

See also

spectrogram

Short-time Fourier power spectrogram.

band_power

Welch PSD integrated over a frequency band.

Notes

The tapers are periodic (sym=False) L2-normalised DPSS windows from scipy.signal.windows.dpss, designed on the host. The result equals

mne.time_frequency.psd_array_multitaper(
    x.T, sfreq=fs, bandwidth=2 * nw * fs / T, adaptive=False,
    low_bias=False, normalization='full', remove_dc=True)

when n_tapers = int(2 * nw): with low_bias=False MNE keeps int(2 * nw) tapers, whereas the default here drops the last, poorly concentrated one.

The function is jittable and vmappable and differentiable with respect to x.

References

Examples

>>> import jax.numpy as jnp
>>> from braintools.analysis.spectral import multitaper_psd
>>> t = jnp.arange(1000) * 1e-3
>>> x = jnp.sin(2 * jnp.pi * 40 * t)
>>> freqs, psd = multitaper_psd(x, dt=1e-3, nw=3.0)
>>> psd.shape
(501,)
>>> float(freqs[jnp.argmax(psd)])
40.0
Contents