spectral_connectivity.transforms.Multitaper#

class Multitaper(time_series: ndarray[tuple[int, ...], dtype[floating]], sampling_frequency: float, 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 = 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, fft_workers: int | None = None, taper_weighting: Literal['uniform', 'eigen', 'adaptive'] = 'uniform', adaptive_max_iterations: int = 50, adaptive_tolerance: float = 1e-08)[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) – Samples per second; required because it labels the frequency axis and scales power.
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, rounded to the nearest sample (the time_window_duration attribute reports the rounded duration). If None, analyzes entire time series (no time resolution).
time_window_step (float, optional) – Step size in seconds between consecutive time windows, rounded to the nearest sample (the time_window_step attribute reports the rounded step). 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, default=0) – Time in seconds of the first sample. Must be a scalar: time is one-dimensional, so per-trial start times are not supported.
n_fft_samples (int, optional) – Length of FFT. If None, uses a value >= n_time_samples_per_window chosen to be fast for the FFT algorithm. The value is determined by scipy.fft.next_fast_len (or cupy.fft.next_fast_len when GPU is enabled).
n_time_samples_per_window (int, optional) – Number of samples per time window. Computed (to the nearest sample) from time_window_duration if not provided; when both are given they must resolve to the same count, otherwise a ValueError is raised.
n_time_samples_per_step (int, optional) – Number of samples to advance between windows. Computed (to the nearest sample) from time_window_step if not provided; when both are given they must resolve to the same count, otherwise a ValueError is raised.
is_low_bias (bool, default=True) – If True, exclude tapers with eigenvalues < MIN_EIGENVALUE_THRESHOLD (0.9) to reduce bias.
fft_workers (int, optional) – Number of parallel worker threads for the CPU FFT (forwarded to scipy.fft.fft; -1 uses all cores). None (the default) keeps SciPy’s single-threaded default, which avoids oversubscribing CPUs when the analysis is already parallelized at a higher level (e.g. across trials or subjects). This is a CPU-only option; it is ignored on the GPU backend, whose FFT is already parallel.
taper_weighting ({"uniform", "eigen", "adaptive"}, default="uniform") –
How DPSS eigencoefficients contribute to spectra. "uniform" keeps historical equal weighting (a good default). "eigen" down-weights the less concentrated tapers by their eigenvalues. "adaptive" iteratively estimates frequency- and signal-specific Thomson weights; reach for it when high spectral dynamic range or line noise makes some tapers systematically leakier, at the cost of extra iterations for a typically modest bias reduction. Non-uniform modes require internally generated DPSS tapers.

The adaptive weights are applied to the coefficients themselves: each signal’s weights d_k(f) are RMS-normalized over tapers and multiplied into its eigencoefficients, so the downstream taper mean of |coefficient|**2 is exactly Thomson’s adaptive auto-spectrum sum_k d_k^2 |Y_k|^2 / sum_k d_k^2. A cross-spectrum between signals x and y, however, becomes sum_k d_xk d_yk Y_xk conj(Y_yk) / sqrt(sum_k d_xk^2 sum_k d_yk^2) rather than Thomson’s joint ... / sum_k d_xk d_yk. By the Cauchy-Schwarz inequality the two differ by the cosine similarity of the two weight vectors, so adaptive coherence (and every normalized cross-signal measure) is shrunk by that factor wherever the signals’ spectra, and hence their weights, differ: a perfectly coherent pair with different spectra does not reach |coherence| = 1. The shrinkage vanishes when the weights agree (white spectra, or the "eigen" mode, whose weights are signal-independent).
adaptive_max_iterations (int, default=50) – Iteration ceiling for Thomson adaptive weighting.
adaptive_tolerance (float, default=1e-8) – Relative convergence tolerance for the adaptive spectrum estimate.

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,)

See also

spectral_connectivity.connectivity.Connectivity: Compute connectivity measures from a Multitaper transform.
spectral_connectivity.transforms.prepare_time_series: Reshape 1-D/2-D input to the required 3-D format.

Notes

A Multitaper represents one immutable transform configuration. Constructor arrays are copied, and array properties return detached read-only copies. To change the data or a parameter, create a new instance.

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.

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.default_rng(0).standard_normal((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}")
FFT shape: (1, 1, 7, 5000, 64)
>>> print(f"Frequencies: {len(mt.frequencies)} bins, max = {mt.frequencies.max():.1f} Hz")
Frequencies: 5000 bins, max = 499.8 Hz

Manual reshaping with np.newaxis (advanced):

>>> # Generate test signal: 50Hz + noise
>>> fs = 1000  # 1 kHz sampling
>>> t = np.arange(0, 1, 1/fs)
>>> rng = np.random.default_rng(0)
>>> signal = np.sin(2*np.pi*50*t) + 0.1*rng.standard_normal(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.default_rng(0).standard_normal((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

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.
observations_are_independent: Whether the trial/taper observations may be counted as independent.
start_time: Independent read-only copy of the transform’s start-time coordinate.
taper_eigenvalues: DPSS spectral-concentration ratios used for taper weighting.
tapers: Return the tapers used for the multitaper function.
time: Return time of each time bin.
time_bins_are_independent: Whether the time windows may be counted as independent observations.
time_series: Independent read-only copy of the input time-series snapshot.
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.
`is_one_sided`	Multitaper returns the full two-sided FFT spectrum (both positive and negative frequencies), so consumers must not assume a one-sided layout.
`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.
`observations_are_independent`	Whether the trial/taper observations may be counted as independent.
`start_time`	Independent read-only copy of the transform's start-time coordinate.
`taper_eigenvalues`	DPSS spectral-concentration ratios used for taper weighting.
`tapers`	Return the tapers used for the multitaper function.
`time`	Return time of each time bin.
`time_bins_are_independent`	Whether the time windows may be counted as independent observations.
`time_series`	Independent read-only copy of the input time-series snapshot.
`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

is_one_sided = False#: Multitaper returns the full two-sided FFT spectrum (both positive and negative frequencies), so consumers must not assume a one-sided layout.

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.

An explicit n_time_samples_per_step is used as given. Otherwise, if time_window_step is set, the step is the nearest whole number of samples in that duration. If neither is set, the step defaults to the window length (non-overlapping windows).

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

property observations_are_independent: bool#

Whether the trial/taper observations may be counted as independent.

DPSS tapers are orthogonal, so for a process that is smooth across the taper bandwidth the eigencoefficients of one window are approximately uncorrelated (Thomson 1982; Percival & Walden 1993, ch. 7), and trials are independent realizations. Connectivity reads this flag to decide whether n_observations (trials x tapers) counts effectively independent samples: the jackknife, the debiased measures (pairwise phase consistency, debiased squared PLI/WPLI) and the zero-coherence significance test (Beta(1, n_observations - 1) null) all assume it does. Always True for this transform. The flag describes the observations within one window; correlation between overlapping windows, which matters only for expectations that also average over time, is reported by time_bins_are_independent.

Return type:: bool

property start_time: Any#: Independent read-only copy of the transform’s start-time coordinate.

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.default_rng(0).standard_normal((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 Hz
Time-halfbandwidth product:    3
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:          1000