[med-svn] [python-mne] 335/376: ENH : adding multi tapper PSD computation
Yaroslav Halchenko
debian at onerussian.com
Fri Nov 27 17:23:17 UTC 2015
This is an automated email from the git hooks/post-receive script.
yoh pushed a commit to annotated tag v0.1
in repository python-mne.
commit 5b3d1ae49d2ef439c1641d57ffbeb331000559fe
Author: Alexandre Gramfort <alexandre.gramfort at inria.fr>
Date: Wed Aug 10 13:58:01 2011 -0400
ENH : adding multi tapper PSD computation
---
.../plot_compute_raw_data_spectrum.py | 55 +++++++++++++++++
mne/time_frequency/__init__.py | 1 +
mne/time_frequency/psd.py | 70 ++++++++++++++++++++++
mne/time_frequency/tests/test_psd.py | 33 ++++++++++
4 files changed, 159 insertions(+)
diff --git a/examples/time_frequency/plot_compute_raw_data_spectrum.py b/examples/time_frequency/plot_compute_raw_data_spectrum.py
new file mode 100644
index 0000000..2ad542b
--- /dev/null
+++ b/examples/time_frequency/plot_compute_raw_data_spectrum.py
@@ -0,0 +1,55 @@
+"""
+==================================================
+Compute the spectrum of raw data using multi-taper
+==================================================
+
+This script shows how to compute the power spectral density (PSD)
+of measurements on a raw dataset.
+
+"""
+# Authors: Alexandre Gramfort <gramfort at nmr.mgh.harvard.edu>
+#
+# License: BSD (3-clause)
+
+print __doc__
+
+import numpy as np
+
+from mne import fiff
+from mne.time_frequency import compute_raw_psd
+from mne.datasets import sample
+
+###############################################################################
+# Set parameters
+data_path = sample.data_path('..')
+raw_fname = data_path + '/MEG/sample/sample_audvis_raw.fif'
+
+# Setup for reading the raw data
+raw = fiff.Raw(raw_fname)
+exclude = raw.info['bads'] + ['MEG 2443', 'EEG 053'] # bads + 2 more
+
+# picks MEG gradiometers
+picks = fiff.pick_types(raw.info, meg='grad', eeg=False, eog=False,
+ stim=False, exclude=exclude)
+
+tmin, tmax = 0, 60 # use the first 60s of data
+fmin, fmax = 2, 70 # look at frequencies between 5 and 70Hz
+NFFT = 2048 # the FFT size (NFFT). Ideally a power of 2
+psds, freqs = compute_raw_psd(raw, tmin=tmin, tmax=tmax, picks=picks,
+ fmin=fmin, fmax=fmax, NFFT=NFFT, n_jobs=1)
+
+###############################################################################
+# Compute mean and standard deviation accross channels and then plot
+psd_mean = np.mean(psds, axis=0)
+psd_std = np.std(psds, axis=0)
+
+hyp_limits = (psd_mean - psd_std, psd_mean + psd_std)
+
+import pylab as pl
+pl.close('all')
+pl.plot(freqs, psd_mean)
+pl.fill_between(freqs, hyp_limits[0], y2=hyp_limits[1], color=(1, 0, 0, .3),
+ alpha=0.5)
+pl.xlabel('Freq (Hz)')
+pl.ylabel('PSD')
+pl.show()
diff --git a/mne/time_frequency/__init__.py b/mne/time_frequency/__init__.py
index 837f961..34395c5 100644
--- a/mne/time_frequency/__init__.py
+++ b/mne/time_frequency/__init__.py
@@ -2,3 +2,4 @@
"""
from .tfr import induced_power, single_trial_power
+from .psd import compute_raw_psd
diff --git a/mne/time_frequency/psd.py b/mne/time_frequency/psd.py
new file mode 100644
index 0000000..c1ffdae
--- /dev/null
+++ b/mne/time_frequency/psd.py
@@ -0,0 +1,70 @@
+# Author : Alexandre Gramfort, gramfort at nmr.mgh.harvard.edu (2011)
+# License : BSD 3-clause
+
+import numpy as np
+import pylab as pl
+from ..parallel import parallel_func
+
+
+def compute_raw_psd(raw, tmin=0, tmax=np.inf, picks=None,
+ fmin=0, fmax=np.inf, NFFT=2048, n_jobs=1):
+ """Compute power spectral density with multi-taper
+
+ Parameters
+ ----------
+ raw: instance of Raw
+ The raw data.
+
+ tmin: float
+ Min time instant to consider
+
+ tmax: float
+ Max time instant to consider
+
+ picks: None or array of integers
+ The selection of channels to include in the computation.
+ If None, take all channels.
+
+ fmin: float
+ Min frequency of interest
+
+ fmax: float
+ Max frequency of interest
+
+ NFFT: int
+ The length of the tappers ie. the windows. The smaller
+ it is the smoother are the PSDs.
+
+ n_jobs: int
+ Number of CPUs to use in the computation.
+
+ Return
+ ------
+ psd: array of float
+ The PSD for all channels
+
+ freqs: array of float
+ The frequencies
+ """
+ start, stop = raw.time_to_index(tmin, tmax)
+ if picks is not None:
+ data, times = raw[picks, start:(stop + 1)]
+ else:
+ data, times = raw[:, start:(stop + 1)]
+
+ NFFT = int(NFFT)
+ Fs = raw.info['sfreq']
+
+ print "Effective window size : %s (s)" % (NFFT * Fs)
+
+ parallel, my_psd, n_jobs = parallel_func(pl.psd, n_jobs, verbose=0)
+ out = parallel(my_psd(d, Fs=Fs, NFFT=NFFT) for d in data)
+
+ freqs = out[0][1]
+ psd = np.array(zip(*out)[0])
+
+ mask = (freqs >= fmin) & (freqs <= fmax)
+ freqs = freqs[mask]
+ psd = psd[:, mask]
+
+ return psd, freqs
diff --git a/mne/time_frequency/tests/test_psd.py b/mne/time_frequency/tests/test_psd.py
new file mode 100644
index 0000000..88b971e
--- /dev/null
+++ b/mne/time_frequency/tests/test_psd.py
@@ -0,0 +1,33 @@
+import numpy as np
+import os.path as op
+from numpy.testing import assert_array_almost_equal
+
+from ... import fiff, Epochs, read_events
+from ...time_frequency import compute_raw_psd
+
+
+raw_fname = op.join(op.dirname(__file__), '..', '..', 'fiff', 'tests', 'data',
+ 'test_raw.fif')
+
+def test_psd():
+ """Test PSD estimation
+ """
+ raw = fiff.Raw(raw_fname)
+
+ exclude = raw.info['bads'] + ['MEG 2443', 'EEG 053'] # bads + 2 more
+
+ # picks MEG gradiometers
+ picks = fiff.pick_types(raw.info, meg='mag', eeg=False,
+ stim=False, exclude=exclude)
+
+ picks = picks[:2]
+
+ tmin, tmax = 0, 10 # use the first 60s of data
+ fmin, fmax = 2, 70 # look at frequencies between 5 and 70Hz
+ NFFT = 124 # the FFT size (NFFT). Ideally a power of 2
+ psds, freqs = compute_raw_psd(raw, tmin=tmin, tmax=tmax, picks=picks,
+ fmin=fmin, fmax=fmax, NFFT=NFFT, n_jobs=1)
+
+ assert psds.shape == (len(picks), len(freqs))
+ assert np.sum(freqs < 0) == 0
+ assert np.sum(psds < 0) == 0
--
Alioth's /usr/local/bin/git-commit-notice on /srv/git.debian.org/git/debian-med/python-mne.git
More information about the debian-med-commit
mailing list