0
# Time-Frequency Analysis
1
2
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.
3
4
## Capabilities
5
6
### Power Spectral Density (PSD) Estimation
7
8
Compute power spectral density using various methods for frequency domain analysis.
9
10
```python { .api }
11
def psd_welch(inst: Union[Raw, Epochs], fmin: float = 0, fmax: float = np.inf,
12
tmin: Optional[float] = None, tmax: Optional[float] = None,
13
picks: Optional[Union[str, List]] = None, proj: bool = False,
14
n_fft: int = 2048, n_overlap: int = 0, n_per_seg: Optional[int] = None,
15
average: str = 'mean', window: str = 'hamming',
16
verbose: Optional[Union[bool, str, int]] = None, **kwargs) -> Spectrum:
17
"""
18
Compute power spectral density using Welch's method.
19
20
Parameters:
21
- inst: Raw or Epochs instance
22
- fmin: Minimum frequency
23
- fmax: Maximum frequency
24
- tmin: Start time (Epochs only)
25
- tmax: End time (Epochs only)
26
- picks: Channel selection
27
- proj: Apply projections
28
- n_fft: FFT length
29
- n_overlap: Overlap between segments
30
- n_per_seg: Length of each segment
31
- average: Averaging method ('mean' or 'median')
32
- window: Window function
33
- verbose: Verbosity level
34
35
Returns:
36
Spectrum object containing PSD
37
"""
38
39
def psd_multitaper(inst: Union[Raw, Epochs], fmin: float = 0, fmax: float = np.inf,
40
tmin: Optional[float] = None, tmax: Optional[float] = None,
41
picks: Optional[Union[str, List]] = None, proj: bool = False,
42
bandwidth: Optional[float] = None, adaptive: bool = False,
43
low_bias: bool = True, normalization: str = 'length',
44
verbose: Optional[Union[bool, str, int]] = None, **kwargs) -> Spectrum:
45
"""
46
Compute power spectral density using multitaper method.
47
48
Parameters:
49
- inst: Raw or Epochs instance
50
- fmin: Minimum frequency
51
- fmax: Maximum frequency
52
- tmin: Start time (Epochs only)
53
- tmax: End time (Epochs only)
54
- picks: Channel selection
55
- proj: Apply projections
56
- bandwidth: Multitaper bandwidth
57
- adaptive: Use adaptive weighting
58
- low_bias: Reduce spectral bias
59
- normalization: Normalization method
60
- verbose: Verbosity level
61
62
Returns:
63
Spectrum object containing PSD
64
"""
65
66
def psd_array_welch(x: ArrayLike, sfreq: float, fmin: float = 0, fmax: float = np.inf,
67
n_fft: int = 2048, n_overlap: int = 0, n_per_seg: Optional[int] = None,
68
average: str = 'mean', window: str = 'hamming') -> Tuple[ArrayLike, ArrayLike]:
69
"""
70
Compute PSD of array using Welch's method.
71
72
Parameters:
73
- x: Data array (n_channels, n_times)
74
- sfreq: Sampling frequency
75
- fmin: Minimum frequency
76
- fmax: Maximum frequency
77
- n_fft: FFT length
78
- n_overlap: Overlap between segments
79
- n_per_seg: Length of each segment
80
- average: Averaging method
81
- window: Window function
82
83
Returns:
84
Tuple of (psd_array, frequencies)
85
"""
86
```
87
88
### Time-Frequency Decomposition
89
90
Analyze temporal evolution of spectral content using various time-frequency methods.
91
92
```python { .api }
93
def tfr_morlet(epochs: Epochs, freqs: ArrayLike, n_cycles: Union[float, ArrayLike],
94
use_fft: bool = True, return_itc: bool = True, decim: int = 1,
95
n_jobs: int = 1, picks: Optional[Union[str, List]] = None,
96
zero_mean: bool = True, average: bool = True, output: str = 'power',
97
verbose: Optional[Union[bool, str, int]] = None) -> Union[AverageTFR, EpochsTFR]:
98
"""
99
Compute time-frequency representation using Morlet wavelets.
100
101
Parameters:
102
- epochs: Epochs data
103
- freqs: Frequencies of interest
104
- n_cycles: Number of cycles for wavelets
105
- use_fft: Use FFT-based convolution
106
- return_itc: Return inter-trial coherence
107
- decim: Decimation factor
108
- n_jobs: Number of parallel jobs
109
- picks: Channel selection
110
- zero_mean: Zero-mean wavelets
111
- average: Average across epochs
112
- output: Output type ('power', 'phase', 'complex')
113
- verbose: Verbosity level
114
115
Returns:
116
AverageTFR or EpochsTFR object
117
"""
118
119
def tfr_multitaper(epochs: Epochs, freqs: ArrayLike, n_cycles: Union[float, ArrayLike],
120
time_bandwidth: float = 4.0, use_fft: bool = True, return_itc: bool = True,
121
decim: int = 1, n_jobs: int = 1, picks: Optional[Union[str, List]] = None,
122
zero_mean: bool = True, average: bool = True, output: str = 'power',
123
verbose: Optional[Union[bool, str, int]] = None) -> Union[AverageTFR, EpochsTFR]:
124
"""
125
Compute time-frequency representation using multitaper method.
126
127
Parameters:
128
- epochs: Epochs data
129
- freqs: Frequencies of interest
130
- n_cycles: Number of cycles
131
- time_bandwidth: Time-bandwidth product
132
- use_fft: Use FFT-based convolution
133
- return_itc: Return inter-trial coherence
134
- decim: Decimation factor
135
- n_jobs: Number of parallel jobs
136
- picks: Channel selection
137
- zero_mean: Zero-mean tapers
138
- average: Average across epochs
139
- output: Output type
140
- verbose: Verbosity level
141
142
Returns:
143
AverageTFR or EpochsTFR object
144
"""
145
146
def tfr_stockwell(epochs: Epochs, fmin: Optional[float] = None, fmax: Optional[float] = None,
147
n_fft: Optional[int] = None, width: float = 1.0, decim: int = 1,
148
return_itc: bool = False, n_jobs: int = 1,
149
verbose: Optional[Union[bool, str, int]] = None) -> AverageTFR:
150
"""
151
Compute time-frequency representation using Stockwell transform.
152
153
Parameters:
154
- epochs: Epochs data
155
- fmin: Minimum frequency
156
- fmax: Maximum frequency
157
- n_fft: FFT length
158
- width: Width parameter
159
- decim: Decimation factor
160
- return_itc: Return inter-trial coherence
161
- n_jobs: Number of parallel jobs
162
- verbose: Verbosity level
163
164
Returns:
165
AverageTFR object
166
"""
167
168
def stft(x: ArrayLike, wsize: int, tstep: int = 1, verbose: Optional[Union[bool, str, int]] = None) -> ArrayLike:
169
"""
170
Compute short-time Fourier transform.
171
172
Parameters:
173
- x: Signal array
174
- wsize: Window size
175
- tstep: Time step
176
- verbose: Verbosity level
177
178
Returns:
179
STFT coefficients
180
"""
181
182
def istft(X: ArrayLike, wsize: int, tstep: int = 1, verbose: Optional[Union[bool, str, int]] = None) -> ArrayLike:
183
"""
184
Compute inverse short-time Fourier transform.
185
186
Parameters:
187
- X: STFT coefficients
188
- wsize: Window size
189
- tstep: Time step
190
- verbose: Verbosity level
191
192
Returns:
193
Reconstructed signal
194
"""
195
```
196
197
### Cross-Spectral Density (CSD)
198
199
Compute cross-spectral density for connectivity and coherence analysis.
200
201
```python { .api }
202
def csd_fourier(epochs: Epochs, fmin: float = 0, fmax: float = np.inf,
203
tmin: Optional[float] = None, tmax: Optional[float] = None,
204
picks: Optional[Union[str, List]] = None, n_fft: int = 2048,
205
n_overlap: int = 0, n_per_seg: Optional[int] = None,
206
verbose: Optional[Union[bool, str, int]] = None) -> CrossSpectralDensity:
207
"""
208
Compute cross-spectral density using Fourier method.
209
210
Parameters:
211
- epochs: Epochs data
212
- fmin: Minimum frequency
213
- fmax: Maximum frequency
214
- tmin: Start time
215
- tmax: End time
216
- picks: Channel selection
217
- n_fft: FFT length
218
- n_overlap: Overlap between segments
219
- n_per_seg: Length of each segment
220
- verbose: Verbosity level
221
222
Returns:
223
CrossSpectralDensity object
224
"""
225
226
def csd_morlet(epochs: Epochs, frequencies: ArrayLike, n_cycles: Union[float, ArrayLike] = 7,
227
use_fft: bool = True, decim: int = 1, n_jobs: int = 1,
228
zero_mean: bool = True, verbose: Optional[Union[bool, str, int]] = None) -> CrossSpectralDensity:
229
"""
230
Compute cross-spectral density using Morlet wavelets.
231
232
Parameters:
233
- epochs: Epochs data
234
- frequencies: Frequencies of interest
235
- n_cycles: Number of cycles for wavelets
236
- use_fft: Use FFT-based convolution
237
- decim: Decimation factor
238
- n_jobs: Number of parallel jobs
239
- zero_mean: Zero-mean wavelets
240
- verbose: Verbosity level
241
242
Returns:
243
CrossSpectralDensity object
244
"""
245
246
def csd_multitaper(epochs: Epochs, fmin: float, fmax: float, tmin: Optional[float] = None,
247
tmax: Optional[float] = None, picks: Optional[Union[str, List]] = None,
248
n_fft: int = 2048, bandwidth: Optional[float] = None,
249
adaptive: bool = False, low_bias: bool = True,
250
normalization: str = 'length', n_jobs: int = 1,
251
verbose: Optional[Union[bool, str, int]] = None) -> CrossSpectralDensity:
252
"""
253
Compute cross-spectral density using multitaper method.
254
255
Parameters:
256
- epochs: Epochs data
257
- fmin: Minimum frequency
258
- fmax: Maximum frequency
259
- tmin: Start time
260
- tmax: End time
261
- picks: Channel selection
262
- n_fft: FFT length
263
- bandwidth: Multitaper bandwidth
264
- adaptive: Use adaptive weighting
265
- low_bias: Reduce spectral bias
266
- normalization: Normalization method
267
- n_jobs: Number of parallel jobs
268
- verbose: Verbosity level
269
270
Returns:
271
CrossSpectralDensity object
272
"""
273
```
274
275
### Wavelets and Window Functions
276
277
Generate wavelets and window functions for time-frequency analysis.
278
279
```python { .api }
280
def morlet(sfreq: float, freqs: ArrayLike, n_cycles: Union[float, ArrayLike] = 7,
281
sigma: Optional[float] = None, zero_mean: bool = True) -> ArrayLike:
282
"""
283
Generate Morlet wavelets.
284
285
Parameters:
286
- sfreq: Sampling frequency
287
- freqs: Frequencies for wavelets
288
- n_cycles: Number of cycles
289
- sigma: Standard deviation parameter
290
- zero_mean: Make wavelets zero-mean
291
292
Returns:
293
Array of Morlet wavelets
294
"""
295
296
def dpss_windows(N: int, half_nbw: float, Kmax: int, sym: bool = True,
297
norm: Optional[str] = None, return_ratios: bool = False) -> Union[ArrayLike, Tuple]:
298
"""
299
Generate DPSS (Slepian) windows for multitaper analysis.
300
301
Parameters:
302
- N: Window length
303
- half_nbw: Half bandwidth
304
- Kmax: Maximum number of windows
305
- sym: Symmetric windows
306
- norm: Normalization method
307
- return_ratios: Return eigenvalue ratios
308
309
Returns:
310
DPSS windows array or tuple with ratios
311
"""
312
313
def stftfreq(wsize: int, sfreq: float = 2 * np.pi) -> ArrayLike:
314
"""
315
Get frequencies for STFT.
316
317
Parameters:
318
- wsize: Window size
319
- sfreq: Sampling frequency
320
321
Returns:
322
Frequency array
323
"""
324
```
325
326
### Time-Frequency Data Classes
327
328
Container classes for time-frequency representations with analysis and visualization methods.
329
330
```python { .api }
331
class Spectrum:
332
"""Power spectral density container."""
333
334
def __init__(self, data: ArrayLike, freqs: ArrayLike, info: Info,
335
verbose: Optional[Union[bool, str, int]] = None):
336
"""
337
Initialize Spectrum object.
338
339
Parameters:
340
- data: PSD data array
341
- freqs: Frequency array
342
- info: Measurement info
343
- verbose: Verbosity level
344
"""
345
346
def plot(self, picks: Optional[Union[str, List]] = None, exclude: str = 'bads',
347
amplitude: bool = True, xscale: str = 'linear', area_mode: Optional[str] = None,
348
area_alpha: float = 0.33, color: str = 'black', average: bool = False,
349
spatial_colors: Union[bool, str] = True, sphere: Optional[float] = None,
350
show: bool = True) -> Figure:
351
"""
352
Plot power spectral density.
353
354
Parameters:
355
- picks: Channel selection
356
- exclude: Channels to exclude
357
- amplitude: Plot amplitude instead of power
358
- xscale: X-axis scale ('linear' or 'log')
359
- area_mode: Area plotting mode
360
- area_alpha: Transparency for area
361
- color: Line color
362
- average: Average across channels
363
- spatial_colors: Use spatial colors
364
- sphere: Sphere radius for channel positions
365
- show: Show plot immediately
366
367
Returns:
368
Figure object
369
"""
370
371
def plot_topomap(self, bands: Optional[Dict[str, Tuple[float, float]]] = None,
372
ch_type: Optional[str] = None, normalize: bool = False,
373
agg_fun: Optional[callable] = None, dB: bool = False,
374
sensors: Union[bool, str] = True, show_names: Union[bool, callable] = False,
375
mask: Optional[ArrayLike] = None, mask_params: Optional[Dict] = None,
376
contours: int = 6, outlines: str = 'head', sphere: Optional[float] = None,
377
image_interp: str = 'bilinear', extrapolate: str = 'auto',
378
border: str = 'mean', res: int = 64, size: int = 1,
379
cmap: Optional[str] = None, vlim: Tuple[Optional[float], Optional[float]] = (None, None),
380
cnorm: Optional[str] = None, colorbar: bool = True, cbar_fmt: str = '%3.1f',
381
units: Optional[str] = None, axes: Optional[plt.Axes] = None,
382
show: bool = True) -> Figure:
383
"""
384
Plot topographic maps of spectral power in frequency bands.
385
386
Parameters:
387
- bands: Frequency bands dictionary
388
- ch_type: Channel type
389
- normalize: Normalize by total power
390
- agg_fun: Aggregation function
391
- dB: Plot in decibels
392
- sensors: Show sensor locations
393
- show_names: Show channel names
394
- mask: Mask for statistical significance
395
- mask_params: Mask plotting parameters
396
- contours: Number of contour lines
397
- outlines: Outline style
398
- sphere: Sphere radius
399
- image_interp: Image interpolation
400
- extrapolate: Extrapolation method
401
- border: Border handling
402
- res: Resolution
403
- size: Figure size multiplier
404
- cmap: Colormap
405
- vlim: Value limits
406
- cnorm: Color normalization
407
- colorbar: Show colorbar
408
- cbar_fmt: Colorbar format
409
- units: Units for colorbar
410
- axes: Matplotlib axes
411
- show: Show plot immediately
412
413
Returns:
414
Figure object
415
"""
416
417
class AverageTFR:
418
"""Averaged time-frequency representation."""
419
420
def plot(self, picks: Optional[Union[str, List]] = None, baseline: Optional[Tuple[float, float]] = None,
421
mode: str = 'mean', tmin: Optional[float] = None, tmax: Optional[float] = None,
422
fmin: Optional[float] = None, fmax: Optional[float] = None,
423
vmin: Optional[float] = None, vmax: Optional[float] = None,
424
cmap: Optional[str] = None, dB: bool = False, colorbar: bool = True,
425
show: bool = True, title: Optional[str] = None, axes: Optional[plt.Axes] = None,
426
verbose: Optional[Union[bool, str, int]] = None) -> Figure:
427
"""
428
Plot time-frequency representation.
429
430
Parameters:
431
- picks: Channel selection
432
- baseline: Baseline period for correction
433
- mode: Baseline correction mode
434
- tmin: Start time
435
- tmax: End time
436
- fmin: Minimum frequency
437
- fmax: Maximum frequency
438
- vmin: Minimum value for colormap
439
- vmax: Maximum value for colormap
440
- cmap: Colormap
441
- dB: Plot in decibels
442
- colorbar: Show colorbar
443
- show: Show plot immediately
444
- title: Plot title
445
- axes: Matplotlib axes
446
- verbose: Verbosity level
447
448
Returns:
449
Figure object
450
"""
451
452
def plot_topomap(self, tmin: Optional[float] = None, tmax: Optional[float] = None,
453
fmin: Optional[float] = None, fmax: Optional[float] = None,
454
ch_type: Optional[str] = None, baseline: Optional[Tuple[float, float]] = None,
455
mode: str = 'mean', sensors: Union[bool, str] = True,
456
show_names: Union[bool, callable] = False, mask: Optional[ArrayLike] = None,
457
mask_params: Optional[Dict] = None, contours: int = 6,
458
outlines: str = 'head', sphere: Optional[float] = None,
459
image_interp: str = 'bilinear', extrapolate: str = 'auto',
460
border: str = 'mean', res: int = 64, size: int = 2,
461
cmap: Optional[str] = None, vlim: Tuple[Optional[float], Optional[float]] = (None, None),
462
cnorm: Optional[str] = None, colorbar: bool = True, cbar_fmt: str = '%1.1e',
463
units: Optional[str] = None, axes: Optional[plt.Axes] = None,
464
show: bool = True, verbose: Optional[Union[bool, str, int]] = None) -> Figure:
465
"""
466
Plot topographic maps of time-frequency data.
467
468
Returns:
469
Figure object
470
"""
471
472
class CrossSpectralDensity:
473
"""Cross-spectral density container."""
474
475
def __init__(self, data: ArrayLike, ch_names: List[str], frequencies: ArrayLike,
476
n_fft: int, tmin: Optional[float] = None, tmax: Optional[float] = None,
477
projs: Optional[List] = None, bads: List[str] = (),
478
verbose: Optional[Union[bool, str, int]] = None):
479
"""
480
Initialize CrossSpectralDensity object.
481
482
Parameters:
483
- data: CSD data array
484
- ch_names: Channel names
485
- frequencies: Frequency array
486
- n_fft: FFT length used
487
- tmin: Start time
488
- tmax: End time
489
- projs: Projections applied
490
- bads: Bad channels
491
- verbose: Verbosity level
492
"""
493
494
def mean(self, fmin: Optional[float] = None, fmax: Optional[float] = None) -> 'CrossSpectralDensity':
495
"""
496
Average CSD across frequency range.
497
498
Parameters:
499
- fmin: Minimum frequency
500
- fmax: Maximum frequency
501
502
Returns:
503
Averaged CrossSpectralDensity object
504
"""
505
506
def sum(self, fmin: Optional[float] = None, fmax: Optional[float] = None) -> 'CrossSpectralDensity':
507
"""
508
Sum CSD across frequency range.
509
510
Parameters:
511
- fmin: Minimum frequency
512
- fmax: Maximum frequency
513
514
Returns:
515
Summed CrossSpectralDensity object
516
"""
517
```
518
519
## Usage Examples
520
521
### Power Spectral Density Analysis
522
523
```python
524
import mne
525
import numpy as np
526
527
# Load epochs data
528
epochs = mne.read_epochs('sample-epo.fif')
529
530
# Compute PSD using Welch method
531
spectrum = epochs.compute_psd(method='welch', fmin=1, fmax=40, n_fft=2048)
532
533
# Plot PSD
534
spectrum.plot(picks='eeg', average=True)
535
536
# Plot topographic maps of power in frequency bands
537
bands = {'alpha': (8, 12), 'beta': (13, 30)}
538
spectrum.plot_topomap(bands=bands, normalize=True)
539
540
# Compute PSD using multitaper method
541
spectrum_mt = epochs.compute_psd(method='multitaper', fmin=1, fmax=40,
542
bandwidth=4.0, adaptive=True)
543
```
544
545
### Time-Frequency Analysis
546
547
```python
548
import mne
549
import numpy as np
550
551
# Load epochs data
552
epochs = mne.read_epochs('sample-epo.fiv')
553
554
# Define frequencies and cycles
555
freqs = np.logspace(np.log10(1), np.log10(40), 30) # 1-40 Hz
556
n_cycles = freqs / 2 # Adaptive cycles
557
558
# Compute time-frequency representation using Morlet wavelets
559
power, itc = tfr_morlet(epochs, freqs=freqs, n_cycles=n_cycles,
560
use_fft=True, return_itc=True, decim=3)
561
562
# Plot power and ITC
563
power.plot(['MEG 1332'], baseline=(-.5, 0), mode='logratio')
564
itc.plot(['MEG 1332'], baseline=(-.5, 0), mode='logratio')
565
566
# Plot topographic maps at specific time-frequency points
567
power.plot_topomap(ch_type='grad', tmin=0.1, tmax=0.2,
568
fmin=8, fmax=12, baseline=(-.5, 0))
569
570
# Compute TFR using multitaper method
571
power_mt = tfr_multitaper(epochs, freqs=freqs, n_cycles=n_cycles,
572
time_bandwidth=4.0, use_fft=True, decim=3)
573
```
574
575
### Cross-Spectral Density for Connectivity
576
577
```python
578
import mne
579
from mne.connectivity import spectral_connectivity_epochs
580
581
# Load epochs data
582
epochs = mne.read_epochs('sample-epo.fif')
583
584
# Compute cross-spectral density using Morlet wavelets
585
csd = csd_morlet(epochs, frequencies=[10, 20], n_cycles=7)
586
587
# Use CSD for connectivity analysis
588
connectivity = spectral_connectivity_epochs(
589
epochs, method='coh', mode='multitaper',
590
sfreq=epochs.info['sfreq'], fmin=8, fmax=12
591
)
592
593
# Compute CSD using Fourier method for broadband analysis
594
csd_fourier = csd_fourier(epochs, fmin=1, fmax=40, n_fft=2048)
595
```
596
597
## Types
598
599
```python { .api }
600
import numpy as np
601
from typing import Union, Optional, List, Dict, Tuple, Any
602
603
ArrayLike = Union[np.ndarray, List, Tuple]
604
Figure = Any # matplotlib.figure.Figure
605
```