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 <[email protected]> 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 <[email protected]> +# +# 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, [email protected] (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 _______________________________________________ debian-med-commit mailing list [email protected] http://lists.alioth.debian.org/cgi-bin/mailman/listinfo/debian-med-commit
