spectral_connectivity.transforms.Multitaper#

class Multitaper(time_series: ndarray[tuple[int, ...], dtype[floating]], sampling_frequency: float = 1000, time_halfbandwidth_product: float = 3, detrend_type: str | None = 'constant', time_window_duration: float | None = None, time_window_step: float | None = None, n_tapers: int | None = None, tapers: ndarray[tuple[int, ...], dtype[floating]] | None = None, start_time: float | ndarray[tuple[int, ...], dtype[floating]] = 0, n_fft_samples: int | None = None, n_time_samples_per_window: int | None = None, n_time_samples_per_step: int | None = None, is_low_bias: bool = True)[source]#

Bases: object

Multitaper spectral analysis for robust power spectral density estimation.

Transforms time-domain signals to frequency domain using multiple orthogonal tapering windows (Slepian sequences). This approach reduces spectral leakage and provides better spectral estimates than single-taper methods.

Parameters:

time_series (NDArray[floating], shape (n_time_samples, n_trials, n_signals)) –
Input time series data. Must be 3D array. - n_time_samples: number of time points - n_trials: number of trials (use 1 for single trial) - n_signals: number of signals/channels Multiple trials are averaged in the spectral domain.

Important: If your data is 1D or 2D, use prepare_time_series() helper function to convert it to the required 3D format.
sampling_frequency (float, default=1000) – Sampling rate in Hz of the time series data.
time_halfbandwidth_product (float, default=3) –
Time-bandwidth product (often denoted as NW) controlling the trade-off between frequency resolution and variance reduction.

Effect on analysis: - Determines frequency resolution: Δf = 2·NW / T (where T is window duration) - Determines number of tapers: n_tapers = floor(2·NW) - 1 - Higher values → More spectral smoothing, better variance reduction - Lower values → Better frequency resolution, less averaging

Typical values: - NW = 2: Minimal smoothing, 3 tapers, best frequency resolution - NW = 3: Balanced trade-off (recommended default), 5 tapers - NW = 4: More smoothing, 7 tapers, strong variance reduction - NW = 5+: Heavy smoothing, 9+ tapers, very strong variance reduction

Examples: For 1 Hz frequency resolution with 1 second windows: use NW ≤ 0.5 For 5 Hz frequency resolution with 1 second windows: use NW ≤ 2.5 For 10 Hz frequency resolution with 1 second windows: use NW ≤ 5

Use estimate_frequency_resolution() to calculate the resolution for different parameter combinations, or use suggest_parameters() to get recommendations for your specific data and analysis goals.
detrend_type ({"constant", "linear", None}, default="constant") – Type of detrending applied to each time window: - “constant”: remove DC component - “linear”: remove linear trend - None: no detrending
time_window_duration (float, optional) – Duration in seconds of sliding time windows. If None, analyzes entire time series (no time resolution).
time_window_step (float, optional) – Step size in seconds between consecutive time windows. If None, uses non-overlapping windows (step = window duration).
n_tapers (int, optional) – Number of DPSS tapers to use. If None, computed as floor(2 * time_halfbandwidth_product) - 1.
tapers (NDArray[floating], shape (n_time_samples_per_window, n_tapers), optional) – Pre-computed tapering windows. If None, DPSS tapers are computed automatically.
start_time (float or NDArray[floating], default=0) – Start time in seconds of the time series data.
n_fft_samples (int, optional) – Length of FFT. If None, uses next power of 2 >= n_time_samples_per_window.
n_time_samples_per_window (int, optional) – Number of samples per time window. Computed from time_window_duration if not provided.
n_time_samples_per_step (int, optional) – Number of samples to advance between windows. Computed from time_window_step if not provided.
is_low_bias (bool, default=True) – If True, exclude tapers with eigenvalues < MIN_EIGENVALUE_THRESHOLD (0.9) to reduce bias.

fft[source]#

Complex-valued FFT coefficients with shape (n_time_windows, n_trials, n_tapers, n_frequencies, n_signals).

Type:: NDArray[complex128]

frequencies#

Frequency values in Hz corresponding to FFT bins.

Type:: NDArray[float64], shape (n_frequencies,)

time#

Time values in seconds for center of each time window.

Type:: NDArray[float64], shape (n_time_windows,)

Examples

Using the helper function (recommended for 2D data):

>>> import numpy as np
>>> from spectral_connectivity.transforms import Multitaper, prepare_time_series
>>> # EEG recording: 5 seconds at 1000 Hz, 64 channels
>>> eeg_data = np.random.randn(5000, 64)  # Shape: (n_time, n_channels)
>>> eeg_3d = prepare_time_series(eeg_data, axis='signals')
>>> mt = Multitaper(eeg_3d, sampling_frequency=1000, time_halfbandwidth_product=4)
>>> print(f"FFT shape: {mt.fft().shape}")
>>> print(f"Frequencies: {len(mt.frequencies)} bins, max = {mt.frequencies[-1]:.1f} Hz")

Manual reshaping with np.newaxis (advanced):

>>> # Generate test signal: 50Hz + noise
>>> fs = 1000  # 1 kHz sampling
>>> t = np.arange(0, 1, 1/fs)
>>> signal = np.sin(2*np.pi*50*t) + 0.1*np.random.randn(len(t))
>>> # Manually reshape to 3D: (n_time, n_trials, n_signals)
>>> data = signal[:, np.newaxis, np.newaxis]  # Shape: (1000, 1, 1)
>>> mt = Multitaper(data, sampling_frequency=fs, time_halfbandwidth_product=4)

Multiple trials (already 3D):

>>> # Epoched data: 100 trials, 5 channels, 1 second each at 1000 Hz
>>> epoched_data = np.random.randn(1000, 100, 5)  # (n_time, n_trials, n_signals)
>>> mt = Multitaper(epoched_data, sampling_frequency=1000)
>>> print(f"Trials: {mt.n_trials}, Signals: {mt.n_signals}")  # Trials: 100, Signals: 5

Notes

The multitaper method uses discrete prolate spheroidal sequences (DPSS) as tapers, which are optimal for spectral analysis in the sense of minimizing spectral leakage while maximizing energy concentration in the frequency band of interest.

References

[1]

Thomson, D. J. (1982). Spectrum estimation and harmonic analysis. Proceedings of the IEEE, 70(9), 1055-1096.

[2]

Percival, D. B., & Walden, A. T. (1993). Spectral Analysis for Physical Applications. Cambridge University Press.

Attributes:

frequencies: Return frequency of each frequency bin.
frequency_resolution: Return range of frequencies the transform is able to resolve.
n_fft_samples: Return number of frequency bins.
n_signals: Return number of signals computed.
n_tapers: Return number of desired tapers.
n_time_samples_per_step: Return number of samples to step between windows.
n_time_samples_per_window: Return number of samples per time bin.
n_trials: Return number of trials computed.
nyquist_frequency: Return maximum resolvable frequency.
tapers: Return the tapers used for the multitaper function.
time: Return time of each time bin.
time_window_duration: Return duration of each time bin.
time_window_step: Return how much each time window slides.

Methods

`fft`()	Compute the fast Fourier transform using the multitaper method.
`summarize_parameters`()	Generate a human-readable summary of the multitaper analysis parameters.

Methods

`fft`	Compute the fast Fourier transform using the multitaper method.
`summarize_parameters`	Generate a human-readable summary of the multitaper analysis parameters.

Attributes

`frequencies`	Return frequency of each frequency bin.
`frequency_resolution`	Return range of frequencies the transform is able to resolve.
`n_fft_samples`	Return number of frequency bins.
`n_signals`	Return number of signals computed.
`n_tapers`	Return number of desired tapers.
`n_time_samples_per_step`	Return number of samples to step between windows.
`n_time_samples_per_window`	Return number of samples per time bin.
`n_trials`	Return number of trials computed.
`nyquist_frequency`	Return maximum resolvable frequency.
`tapers`	Return the tapers used for the multitaper function.
`time`	Return time of each time bin.
`time_window_duration`	Return duration of each time bin.
`time_window_step`	Return how much each time window slides.

fft() → ndarray[tuple[int, ...], dtype[complexfloating]][source]#

Compute the fast Fourier transform using the multitaper method.

Returns:: fourier_coefficients – Shape (n_time_windows, n_trials, n_tapers, n_fft_samples, n_signals). Complex-valued Fourier coefficients.
Return type:: array

property frequencies: ndarray[tuple[int, ...], dtype[floating]]#

Return frequency of each frequency bin.

Returns:: Frequency values in Hz corresponding to FFT bins.
Return type:: NDArray[float64], shape (n_frequencies,)

property frequency_resolution: float#

Return range of frequencies the transform is able to resolve.

Given the time-frequency tradeoff.

Returns:: Frequency resolution in Hz.
Return type:: float

property n_fft_samples: int#

Return number of frequency bins.

Returns:: Number of FFT samples.
Return type:: int

property n_signals: int#

Return number of signals computed.

Returns:: Number of signals in the time series.
Return type:: int

property n_tapers: int#

Return number of desired tapers.

Note that the number of tapers may be less than this number if the bias of the tapers is too high (eigenvalues > MIN_EIGENVALUE_THRESHOLD = 0.9).

Returns:: Number of tapers to use.
Return type:: int

property n_time_samples_per_step: int#

Return number of samples to step between windows.

If time_window_step is set, then calculate the n_time_samples_per_step based on the time window duration. If time_window_step and n_time_samples_per_step are both not set, default the window step size to the same size as the window.

Returns:: Number of samples to advance between windows.
Return type:: int

property n_time_samples_per_window: int#

Return number of samples per time bin.

Returns:: Number of time samples in each window.
Return type:: int
Raises:: ValueError – If neither n_time_samples_per_window nor time_window_duration is set.

property n_trials: int#

Return number of trials computed.

Returns:: Number of trials in the time series.
Return type:: int

property nyquist_frequency: float#

Return maximum resolvable frequency.

Returns:: Nyquist frequency in Hz.
Return type:: float

summarize_parameters() → str[source]#

Generate a human-readable summary of the multitaper analysis parameters.

This method displays key parameters and their implications for your analysis, making it easier to understand and communicate your spectral analysis settings.

Returns:: summary – A formatted string containing: - Input parameters (sampling frequency, time-halfbandwidth product) - Derived parameters (n_tapers, frequency resolution) - Data dimensions (n_signals, n_trials, n_time_samples) - Frequency range (0 to Nyquist)
Return type:: str

Examples

>>> import numpy as np
>>> from spectral_connectivity.transforms import Multitaper
>>> data = np.random.randn(5000, 1, 64)  # 5s, 64 EEG channels
>>> mt = Multitaper(
...     data,
...     sampling_frequency=1000,
...     time_window_duration=1.0,
...     time_halfbandwidth_product=3,
... )
>>> print(mt.summarize_parameters())
Multitaper Spectral Analysis Configuration
===========================================

Data Shape
----------
Time samples:    5000 (5.00 seconds)
Signals:         64
Trials:          1

Spectral Parameters
-------------------
Sampling frequency:            1000.0 Hz
Time-halfbandwidth product:    3.0
Number of tapers:              5

Time Windowing
--------------
Window duration:  1.000 s (1000 samples)
Window step:      1.000 s (non-overlapping)
Number of windows: 5

Frequency Analysis
------------------
Frequency resolution: 6.0 Hz
Nyquist frequency:    500.0 Hz
Frequency range:      0.0 - 500.0 Hz
FFT samples:          1024