MNE-Python provides comprehensive tools for analyzing MEG, EEG, and other neuroimaging data with advanced source estimation and connectivity analysis.
—
Advanced spectral analysis tools including power spectral density estimation, time-frequency decomposition using wavelets and multitaper methods, and cross-spectral density computation for connectivity analysis.
Compute power spectral density using various methods for frequency domain analysis.
def psd_welch(inst: Union[Raw, Epochs], fmin: float = 0, fmax: float = np.inf,
tmin: Optional[float] = None, tmax: Optional[float] = None,
picks: Optional[Union[str, List]] = None, proj: bool = False,
n_fft: int = 2048, n_overlap: int = 0, n_per_seg: Optional[int] = None,
average: str = 'mean', window: str = 'hamming',
verbose: Optional[Union[bool, str, int]] = None, **kwargs) -> Spectrum:
"""
Compute power spectral density using Welch's method.
Parameters:
- inst: Raw or Epochs instance
- fmin: Minimum frequency
- fmax: Maximum frequency
- tmin: Start time (Epochs only)
- tmax: End time (Epochs only)
- picks: Channel selection
- proj: Apply projections
- n_fft: FFT length
- n_overlap: Overlap between segments
- n_per_seg: Length of each segment
- average: Averaging method ('mean' or 'median')
- window: Window function
- verbose: Verbosity level
Returns:
Spectrum object containing PSD
"""
def psd_multitaper(inst: Union[Raw, Epochs], fmin: float = 0, fmax: float = np.inf,
tmin: Optional[float] = None, tmax: Optional[float] = None,
picks: Optional[Union[str, List]] = None, proj: bool = False,
bandwidth: Optional[float] = None, adaptive: bool = False,
low_bias: bool = True, normalization: str = 'length',
verbose: Optional[Union[bool, str, int]] = None, **kwargs) -> Spectrum:
"""
Compute power spectral density using multitaper method.
Parameters:
- inst: Raw or Epochs instance
- fmin: Minimum frequency
- fmax: Maximum frequency
- tmin: Start time (Epochs only)
- tmax: End time (Epochs only)
- picks: Channel selection
- proj: Apply projections
- bandwidth: Multitaper bandwidth
- adaptive: Use adaptive weighting
- low_bias: Reduce spectral bias
- normalization: Normalization method
- verbose: Verbosity level
Returns:
Spectrum object containing PSD
"""
def psd_array_welch(x: ArrayLike, sfreq: float, fmin: float = 0, fmax: float = np.inf,
n_fft: int = 2048, n_overlap: int = 0, n_per_seg: Optional[int] = None,
average: str = 'mean', window: str = 'hamming') -> Tuple[ArrayLike, ArrayLike]:
"""
Compute PSD of array using Welch's method.
Parameters:
- x: Data array (n_channels, n_times)
- sfreq: Sampling frequency
- fmin: Minimum frequency
- fmax: Maximum frequency
- n_fft: FFT length
- n_overlap: Overlap between segments
- n_per_seg: Length of each segment
- average: Averaging method
- window: Window function
Returns:
Tuple of (psd_array, frequencies)
"""Analyze temporal evolution of spectral content using various time-frequency methods.
def tfr_morlet(epochs: Epochs, freqs: ArrayLike, n_cycles: Union[float, ArrayLike],
use_fft: bool = True, return_itc: bool = True, decim: int = 1,
n_jobs: int = 1, picks: Optional[Union[str, List]] = None,
zero_mean: bool = True, average: bool = True, output: str = 'power',
verbose: Optional[Union[bool, str, int]] = None) -> Union[AverageTFR, EpochsTFR]:
"""
Compute time-frequency representation using Morlet wavelets.
Parameters:
- epochs: Epochs data
- freqs: Frequencies of interest
- n_cycles: Number of cycles for wavelets
- use_fft: Use FFT-based convolution
- return_itc: Return inter-trial coherence
- decim: Decimation factor
- n_jobs: Number of parallel jobs
- picks: Channel selection
- zero_mean: Zero-mean wavelets
- average: Average across epochs
- output: Output type ('power', 'phase', 'complex')
- verbose: Verbosity level
Returns:
AverageTFR or EpochsTFR object
"""
def tfr_multitaper(epochs: Epochs, freqs: ArrayLike, n_cycles: Union[float, ArrayLike],
time_bandwidth: float = 4.0, use_fft: bool = True, return_itc: bool = True,
decim: int = 1, n_jobs: int = 1, picks: Optional[Union[str, List]] = None,
zero_mean: bool = True, average: bool = True, output: str = 'power',
verbose: Optional[Union[bool, str, int]] = None) -> Union[AverageTFR, EpochsTFR]:
"""
Compute time-frequency representation using multitaper method.
Parameters:
- epochs: Epochs data
- freqs: Frequencies of interest
- n_cycles: Number of cycles
- time_bandwidth: Time-bandwidth product
- use_fft: Use FFT-based convolution
- return_itc: Return inter-trial coherence
- decim: Decimation factor
- n_jobs: Number of parallel jobs
- picks: Channel selection
- zero_mean: Zero-mean tapers
- average: Average across epochs
- output: Output type
- verbose: Verbosity level
Returns:
AverageTFR or EpochsTFR object
"""
def tfr_stockwell(epochs: Epochs, fmin: Optional[float] = None, fmax: Optional[float] = None,
n_fft: Optional[int] = None, width: float = 1.0, decim: int = 1,
return_itc: bool = False, n_jobs: int = 1,
verbose: Optional[Union[bool, str, int]] = None) -> AverageTFR:
"""
Compute time-frequency representation using Stockwell transform.
Parameters:
- epochs: Epochs data
- fmin: Minimum frequency
- fmax: Maximum frequency
- n_fft: FFT length
- width: Width parameter
- decim: Decimation factor
- return_itc: Return inter-trial coherence
- n_jobs: Number of parallel jobs
- verbose: Verbosity level
Returns:
AverageTFR object
"""
def stft(x: ArrayLike, wsize: int, tstep: int = 1, verbose: Optional[Union[bool, str, int]] = None) -> ArrayLike:
"""
Compute short-time Fourier transform.
Parameters:
- x: Signal array
- wsize: Window size
- tstep: Time step
- verbose: Verbosity level
Returns:
STFT coefficients
"""
def istft(X: ArrayLike, wsize: int, tstep: int = 1, verbose: Optional[Union[bool, str, int]] = None) -> ArrayLike:
"""
Compute inverse short-time Fourier transform.
Parameters:
- X: STFT coefficients
- wsize: Window size
- tstep: Time step
- verbose: Verbosity level
Returns:
Reconstructed signal
"""Compute cross-spectral density for connectivity and coherence analysis.
def csd_fourier(epochs: Epochs, fmin: float = 0, fmax: float = np.inf,
tmin: Optional[float] = None, tmax: Optional[float] = None,
picks: Optional[Union[str, List]] = None, n_fft: int = 2048,
n_overlap: int = 0, n_per_seg: Optional[int] = None,
verbose: Optional[Union[bool, str, int]] = None) -> CrossSpectralDensity:
"""
Compute cross-spectral density using Fourier method.
Parameters:
- epochs: Epochs data
- fmin: Minimum frequency
- fmax: Maximum frequency
- tmin: Start time
- tmax: End time
- picks: Channel selection
- n_fft: FFT length
- n_overlap: Overlap between segments
- n_per_seg: Length of each segment
- verbose: Verbosity level
Returns:
CrossSpectralDensity object
"""
def csd_morlet(epochs: Epochs, frequencies: ArrayLike, n_cycles: Union[float, ArrayLike] = 7,
use_fft: bool = True, decim: int = 1, n_jobs: int = 1,
zero_mean: bool = True, verbose: Optional[Union[bool, str, int]] = None) -> CrossSpectralDensity:
"""
Compute cross-spectral density using Morlet wavelets.
Parameters:
- epochs: Epochs data
- frequencies: Frequencies of interest
- n_cycles: Number of cycles for wavelets
- use_fft: Use FFT-based convolution
- decim: Decimation factor
- n_jobs: Number of parallel jobs
- zero_mean: Zero-mean wavelets
- verbose: Verbosity level
Returns:
CrossSpectralDensity object
"""
def csd_multitaper(epochs: Epochs, fmin: float, fmax: float, tmin: Optional[float] = None,
tmax: Optional[float] = None, picks: Optional[Union[str, List]] = None,
n_fft: int = 2048, bandwidth: Optional[float] = None,
adaptive: bool = False, low_bias: bool = True,
normalization: str = 'length', n_jobs: int = 1,
verbose: Optional[Union[bool, str, int]] = None) -> CrossSpectralDensity:
"""
Compute cross-spectral density using multitaper method.
Parameters:
- epochs: Epochs data
- fmin: Minimum frequency
- fmax: Maximum frequency
- tmin: Start time
- tmax: End time
- picks: Channel selection
- n_fft: FFT length
- bandwidth: Multitaper bandwidth
- adaptive: Use adaptive weighting
- low_bias: Reduce spectral bias
- normalization: Normalization method
- n_jobs: Number of parallel jobs
- verbose: Verbosity level
Returns:
CrossSpectralDensity object
"""Generate wavelets and window functions for time-frequency analysis.
def morlet(sfreq: float, freqs: ArrayLike, n_cycles: Union[float, ArrayLike] = 7,
sigma: Optional[float] = None, zero_mean: bool = True) -> ArrayLike:
"""
Generate Morlet wavelets.
Parameters:
- sfreq: Sampling frequency
- freqs: Frequencies for wavelets
- n_cycles: Number of cycles
- sigma: Standard deviation parameter
- zero_mean: Make wavelets zero-mean
Returns:
Array of Morlet wavelets
"""
def dpss_windows(N: int, half_nbw: float, Kmax: int, sym: bool = True,
norm: Optional[str] = None, return_ratios: bool = False) -> Union[ArrayLike, Tuple]:
"""
Generate DPSS (Slepian) windows for multitaper analysis.
Parameters:
- N: Window length
- half_nbw: Half bandwidth
- Kmax: Maximum number of windows
- sym: Symmetric windows
- norm: Normalization method
- return_ratios: Return eigenvalue ratios
Returns:
DPSS windows array or tuple with ratios
"""
def stftfreq(wsize: int, sfreq: float = 2 * np.pi) -> ArrayLike:
"""
Get frequencies for STFT.
Parameters:
- wsize: Window size
- sfreq: Sampling frequency
Returns:
Frequency array
"""Container classes for time-frequency representations with analysis and visualization methods.
class Spectrum:
"""Power spectral density container."""
def __init__(self, data: ArrayLike, freqs: ArrayLike, info: Info,
verbose: Optional[Union[bool, str, int]] = None):
"""
Initialize Spectrum object.
Parameters:
- data: PSD data array
- freqs: Frequency array
- info: Measurement info
- verbose: Verbosity level
"""
def plot(self, picks: Optional[Union[str, List]] = None, exclude: str = 'bads',
amplitude: bool = True, xscale: str = 'linear', area_mode: Optional[str] = None,
area_alpha: float = 0.33, color: str = 'black', average: bool = False,
spatial_colors: Union[bool, str] = True, sphere: Optional[float] = None,
show: bool = True) -> Figure:
"""
Plot power spectral density.
Parameters:
- picks: Channel selection
- exclude: Channels to exclude
- amplitude: Plot amplitude instead of power
- xscale: X-axis scale ('linear' or 'log')
- area_mode: Area plotting mode
- area_alpha: Transparency for area
- color: Line color
- average: Average across channels
- spatial_colors: Use spatial colors
- sphere: Sphere radius for channel positions
- show: Show plot immediately
Returns:
Figure object
"""
def plot_topomap(self, bands: Optional[Dict[str, Tuple[float, float]]] = None,
ch_type: Optional[str] = None, normalize: bool = False,
agg_fun: Optional[callable] = None, dB: bool = False,
sensors: Union[bool, str] = True, show_names: Union[bool, callable] = False,
mask: Optional[ArrayLike] = None, mask_params: Optional[Dict] = None,
contours: int = 6, outlines: str = 'head', sphere: Optional[float] = None,
image_interp: str = 'bilinear', extrapolate: str = 'auto',
border: str = 'mean', res: int = 64, size: int = 1,
cmap: Optional[str] = None, vlim: Tuple[Optional[float], Optional[float]] = (None, None),
cnorm: Optional[str] = None, colorbar: bool = True, cbar_fmt: str = '%3.1f',
units: Optional[str] = None, axes: Optional[plt.Axes] = None,
show: bool = True) -> Figure:
"""
Plot topographic maps of spectral power in frequency bands.
Parameters:
- bands: Frequency bands dictionary
- ch_type: Channel type
- normalize: Normalize by total power
- agg_fun: Aggregation function
- dB: Plot in decibels
- sensors: Show sensor locations
- show_names: Show channel names
- mask: Mask for statistical significance
- mask_params: Mask plotting parameters
- contours: Number of contour lines
- outlines: Outline style
- sphere: Sphere radius
- image_interp: Image interpolation
- extrapolate: Extrapolation method
- border: Border handling
- res: Resolution
- size: Figure size multiplier
- cmap: Colormap
- vlim: Value limits
- cnorm: Color normalization
- colorbar: Show colorbar
- cbar_fmt: Colorbar format
- units: Units for colorbar
- axes: Matplotlib axes
- show: Show plot immediately
Returns:
Figure object
"""
class AverageTFR:
"""Averaged time-frequency representation."""
def plot(self, picks: Optional[Union[str, List]] = None, baseline: Optional[Tuple[float, float]] = None,
mode: str = 'mean', tmin: Optional[float] = None, tmax: Optional[float] = None,
fmin: Optional[float] = None, fmax: Optional[float] = None,
vmin: Optional[float] = None, vmax: Optional[float] = None,
cmap: Optional[str] = None, dB: bool = False, colorbar: bool = True,
show: bool = True, title: Optional[str] = None, axes: Optional[plt.Axes] = None,
verbose: Optional[Union[bool, str, int]] = None) -> Figure:
"""
Plot time-frequency representation.
Parameters:
- picks: Channel selection
- baseline: Baseline period for correction
- mode: Baseline correction mode
- tmin: Start time
- tmax: End time
- fmin: Minimum frequency
- fmax: Maximum frequency
- vmin: Minimum value for colormap
- vmax: Maximum value for colormap
- cmap: Colormap
- dB: Plot in decibels
- colorbar: Show colorbar
- show: Show plot immediately
- title: Plot title
- axes: Matplotlib axes
- verbose: Verbosity level
Returns:
Figure object
"""
def plot_topomap(self, tmin: Optional[float] = None, tmax: Optional[float] = None,
fmin: Optional[float] = None, fmax: Optional[float] = None,
ch_type: Optional[str] = None, baseline: Optional[Tuple[float, float]] = None,
mode: str = 'mean', sensors: Union[bool, str] = True,
show_names: Union[bool, callable] = False, mask: Optional[ArrayLike] = None,
mask_params: Optional[Dict] = None, contours: int = 6,
outlines: str = 'head', sphere: Optional[float] = None,
image_interp: str = 'bilinear', extrapolate: str = 'auto',
border: str = 'mean', res: int = 64, size: int = 2,
cmap: Optional[str] = None, vlim: Tuple[Optional[float], Optional[float]] = (None, None),
cnorm: Optional[str] = None, colorbar: bool = True, cbar_fmt: str = '%1.1e',
units: Optional[str] = None, axes: Optional[plt.Axes] = None,
show: bool = True, verbose: Optional[Union[bool, str, int]] = None) -> Figure:
"""
Plot topographic maps of time-frequency data.
Returns:
Figure object
"""
class CrossSpectralDensity:
"""Cross-spectral density container."""
def __init__(self, data: ArrayLike, ch_names: List[str], frequencies: ArrayLike,
n_fft: int, tmin: Optional[float] = None, tmax: Optional[float] = None,
projs: Optional[List] = None, bads: List[str] = (),
verbose: Optional[Union[bool, str, int]] = None):
"""
Initialize CrossSpectralDensity object.
Parameters:
- data: CSD data array
- ch_names: Channel names
- frequencies: Frequency array
- n_fft: FFT length used
- tmin: Start time
- tmax: End time
- projs: Projections applied
- bads: Bad channels
- verbose: Verbosity level
"""
def mean(self, fmin: Optional[float] = None, fmax: Optional[float] = None) -> 'CrossSpectralDensity':
"""
Average CSD across frequency range.
Parameters:
- fmin: Minimum frequency
- fmax: Maximum frequency
Returns:
Averaged CrossSpectralDensity object
"""
def sum(self, fmin: Optional[float] = None, fmax: Optional[float] = None) -> 'CrossSpectralDensity':
"""
Sum CSD across frequency range.
Parameters:
- fmin: Minimum frequency
- fmax: Maximum frequency
Returns:
Summed CrossSpectralDensity object
"""import mne
import numpy as np
# Load epochs data
epochs = mne.read_epochs('sample-epo.fif')
# Compute PSD using Welch method
spectrum = epochs.compute_psd(method='welch', fmin=1, fmax=40, n_fft=2048)
# Plot PSD
spectrum.plot(picks='eeg', average=True)
# Plot topographic maps of power in frequency bands
bands = {'alpha': (8, 12), 'beta': (13, 30)}
spectrum.plot_topomap(bands=bands, normalize=True)
# Compute PSD using multitaper method
spectrum_mt = epochs.compute_psd(method='multitaper', fmin=1, fmax=40,
bandwidth=4.0, adaptive=True)import mne
import numpy as np
# Load epochs data
epochs = mne.read_epochs('sample-epo.fiv')
# Define frequencies and cycles
freqs = np.logspace(np.log10(1), np.log10(40), 30) # 1-40 Hz
n_cycles = freqs / 2 # Adaptive cycles
# Compute time-frequency representation using Morlet wavelets
power, itc = tfr_morlet(epochs, freqs=freqs, n_cycles=n_cycles,
use_fft=True, return_itc=True, decim=3)
# Plot power and ITC
power.plot(['MEG 1332'], baseline=(-.5, 0), mode='logratio')
itc.plot(['MEG 1332'], baseline=(-.5, 0), mode='logratio')
# Plot topographic maps at specific time-frequency points
power.plot_topomap(ch_type='grad', tmin=0.1, tmax=0.2,
fmin=8, fmax=12, baseline=(-.5, 0))
# Compute TFR using multitaper method
power_mt = tfr_multitaper(epochs, freqs=freqs, n_cycles=n_cycles,
time_bandwidth=4.0, use_fft=True, decim=3)import mne
from mne.connectivity import spectral_connectivity_epochs
# Load epochs data
epochs = mne.read_epochs('sample-epo.fif')
# Compute cross-spectral density using Morlet wavelets
csd = csd_morlet(epochs, frequencies=[10, 20], n_cycles=7)
# Use CSD for connectivity analysis
connectivity = spectral_connectivity_epochs(
epochs, method='coh', mode='multitaper',
sfreq=epochs.info['sfreq'], fmin=8, fmax=12
)
# Compute CSD using Fourier method for broadband analysis
csd_fourier = csd_fourier(epochs, fmin=1, fmax=40, n_fft=2048)import numpy as np
from typing import Union, Optional, List, Dict, Tuple, Any
ArrayLike = Union[np.ndarray, List, Tuple]
Figure = Any # matplotlib.figure.FigureInstall with Tessl CLI
npx tessl i tessl/pypi-mne