This is an automated email from the git hooks/post-receive script. yoh pushed a commit to tag 0.4 in repository python-mne.
commit 6c72874b57da1df8527538b8eec2d795bec6b2d4 Author: Alexandre Gramfort <[email protected]> Date: Sat Jul 28 12:09:33 2012 +0200 ENH : add function to compute source PSD from raw file --- .../time_frequency/plot_source_power_spectrum.py | 55 +++++++++ mne/fiff/raw.py | 3 +- mne/minimum_norm/__init__.py | 3 +- mne/minimum_norm/tests/test_time_frequency.py | 21 +++- mne/minimum_norm/time_frequency.py | 129 ++++++++++++++++++++- 5 files changed, 207 insertions(+), 4 deletions(-) diff --git a/examples/time_frequency/plot_source_power_spectrum.py b/examples/time_frequency/plot_source_power_spectrum.py new file mode 100644 index 0000000..918169d --- /dev/null +++ b/examples/time_frequency/plot_source_power_spectrum.py @@ -0,0 +1,55 @@ +""" +========================================================= +Compute power spectrum densities of the sources with dSPM +========================================================= + +Returns an STC file containing the PSD (in dB) of each of the sources. + +""" +# Authors: Alexandre Gramfort <[email protected]> +# +# License: BSD (3-clause) + +print __doc__ + +import mne +from mne import fiff +from mne.datasets import sample +from mne.minimum_norm import read_inverse_operator, compute_source_psd + +############################################################################### +# Set parameters +data_path = sample.data_path('..') +raw_fname = data_path + '/MEG/sample/sample_audvis_raw.fif' +fname_inv = data_path + '/MEG/sample/sample_audvis-meg-oct-6-meg-inv.fif' +fname_label = data_path + '/MEG/sample/labels/Aud-lh.label' + +# Setup for reading the raw data +raw = fiff.Raw(raw_fname, verbose=False) +events = mne.find_events(raw) +inverse_operator = read_inverse_operator(fname_inv) +raw.info['bads'] = ['MEG 2443', 'EEG 053'] + +# picks MEG gradiometers +picks = fiff.pick_types(raw.info, meg=True, eeg=False, eog=True, + stim=False, include=[], exclude=raw.info['bads']) + +tmin, tmax = 0, 120 # use the first 120s of data +fmin, fmax = 4, 100 # look at frequencies between 4 and 100Hz +NFFT = 2048 # the FFT size (NFFT). Ideally a power of 2 +label = mne.read_label(fname_label) + +stc = compute_source_psd(raw, inverse_operator, lambda2=1. / 9., method="dSPM", + tmin=tmin, tmax=tmax, fmin=fmin, fmax=fmax, + pick_normal=True, NFFT=NFFT, label=label) + +stc.save('psd_dSPM') + +############################################################################### +# View PSD of sources in label +import pylab as pl +pl.plot(1e3 * stc.times, stc.data.T) +pl.xlabel('Frequency (Hz)') +pl.ylabel('PSD (dB)') +pl.title('Source Power Spectrum (PSD)') +pl.show() diff --git a/mne/fiff/raw.py b/mne/fiff/raw.py index 6617552..68d4ea4 100644 --- a/mne/fiff/raw.py +++ b/mne/fiff/raw.py @@ -832,7 +832,8 @@ def read_raw_segment(raw, start=0, stop=None, sel=None, data_buffer=None, # Done? if this['last'] >= stop - 1: - print ' [done]' + if verbose: + print ' [done]' break times = (np.arange(start, stop) - raw.first_samp) / raw.info['sfreq'] diff --git a/mne/minimum_norm/__init__.py b/mne/minimum_norm/__init__.py index 06a815a..eaf5dc3 100644 --- a/mne/minimum_norm/__init__.py +++ b/mne/minimum_norm/__init__.py @@ -3,4 +3,5 @@ from .inverse import read_inverse_operator, apply_inverse, \ apply_inverse_raw, make_inverse_operator, \ apply_inverse_epochs, write_inverse_operator -from .time_frequency import source_band_induced_power, source_induced_power +from .time_frequency import source_band_induced_power, source_induced_power, \ + compute_source_psd diff --git a/mne/minimum_norm/tests/test_time_frequency.py b/mne/minimum_norm/tests/test_time_frequency.py index 8c8c6f0..08a71ca 100644 --- a/mne/minimum_norm/tests/test_time_frequency.py +++ b/mne/minimum_norm/tests/test_time_frequency.py @@ -8,7 +8,8 @@ from ...datasets import sample from ... import fiff, find_events, Epochs from ...label import read_label from ..inverse import read_inverse_operator -from ..time_frequency import source_band_induced_power, source_induced_power +from ..time_frequency import source_band_induced_power, source_induced_power, \ + compute_source_psd examples_folder = op.join(op.dirname(__file__), '..', '..', '..', 'examples') @@ -73,3 +74,21 @@ def test_tfr_with_inverse_operator(): assert_true(np.all(phase_lock > 0)) assert_true(np.all(phase_lock <= 1)) assert_true(np.max(power) > 10) + + +def test_source_psd(): + """Test source PSD computation in label""" + raw = fiff.Raw(fname_data) + inverse_operator = read_inverse_operator(fname_inv) + label = read_label(fname_label) + tmin, tmax = 0, 20 # seconds + fmin, fmax = 55, 65 # Hz + NFFT = 2048 + stc = compute_source_psd(raw, inverse_operator, lambda2=1. / 9., method="dSPM", + tmin=tmin, tmax=tmax, fmin=fmin, fmax=fmax, + pick_normal=True, NFFT=NFFT, label=label, overlap=0.1) + assert_true(stc.times[0] >= fmin * 1e-3) + assert_true(stc.times[-1] <= fmax * 1e-3) + # Time max at line frequency (60 Hz in US) + assert_true(59e-3 <= stc.times[np.argmax(np.sum(stc.data, axis=0))] <= 61e-3) + \ No newline at end of file diff --git a/mne/minimum_norm/time_frequency.py b/mne/minimum_norm/time_frequency.py index 48a903d..5264abf 100644 --- a/mne/minimum_norm/time_frequency.py +++ b/mne/minimum_norm/time_frequency.py @@ -3,7 +3,7 @@ # License: BSD (3-clause) import numpy as np -from scipy import linalg +from scipy import linalg, signal, fftpack from ..fiff.constants import FIFF from ..time_frequency.tfr import cwt, morlet @@ -299,3 +299,130 @@ def source_induced_power(epochs, inverse_operator, frequencies, label=None, verbose=True, copy=False) return power, plv + + +def compute_source_psd(raw, inverse_operator, lambda2=1. / 9., method="dSPM", + tmin=None, tmax=None, fmin=0., fmax=200., + NFFT=2048, overlap=0.5, pick_normal=False, label=None, + nave=1, pca=True): + """Compute source power spectrum density (PSD) + + Parameters + ---------- + raw : instance of Raw + The raw data + inverse_operator : dict + The inverse operator + lambda2: float + The regularization parameter + method: "MNE" | "dSPM" | "sLORETA" + Use mininum norm, dSPM or sLORETA + tmin : float | None + The beginning of the time interval of interest (in seconds). If None + start from the beginning of the file. + tmax : float | None + The end of the time interval of interest (in seconds). If None + stop at the end of the file. + fmin : float + The lower frequency of interest + fmax : float + The upper frequency of interest + NFFT: int + Window size for the FFT. Should be a power of 2. + overlap: float + The overlap fraction between windows. Should be between 0 and 1. + 0 means no overlap. + pick_normal : bool + If True, rather than pooling the orientations by taking the norm, + only the radial component is kept. This is only implemented + when working with loose orientations. + label: Label + Restricts the source estimates to a given label + nave : int + The number of averages used to scale the noise covariance matrix. + pca: bool + If True, the true dimension of data is estimated before running + the time frequency transforms. It reduces the computation times + e.g. with a dataset that was maxfiltered (true dim is 64) + + Returns + ------- + stc : SourceEstimate + The PSD (in dB) of each of the sources. + """ + + print 'Considering frequencies %g ... %g Hz' % (fmin, fmax) + + inv = prepare_inverse_operator(inverse_operator, nave, lambda2, method) + is_free_ori = inverse_operator['source_ori'] == FIFF.FIFFV_MNE_FREE_ORI + + # + # Pick the correct channels from the data + # + sel = _pick_channels_inverse_operator(raw.ch_names, inv) + print 'Picked %d channels from the data' % len(sel) + print 'Computing inverse...', + # + # Simple matrix multiplication followed by combination of the + # three current components + # + # This does all the data transformations to compute the weights for the + # eigenleads + # + K, noise_norm, vertno = _assemble_kernel(inv, label, method, pick_normal) + + if pca: + U, s, Vh = linalg.svd(K, full_matrices=False) + rank = np.sum(s > 1e-8 * s[0]) + K = s[:rank] * U[:, :rank] + Vh = Vh[:rank] + print 'Reducing data rank to %d' % rank + else: + Vh = None + + start, stop = 0, raw.last_samp + 1 - raw.first_samp + if tmin is not None: + start = raw.time_to_index(tmin)[0] + if tmax is not None: + stop = raw.time_to_index(tmax)[0] + 1 + NFFT = int(NFFT) + Fs = raw.info['sfreq'] + window = signal.hanning(NFFT) + freqs = fftpack.fftfreq(NFFT, 1. / Fs) + freqs_mask = (freqs >= 0) & (freqs >= fmin) & (freqs <= fmax) + freqs = freqs[freqs_mask] + fstep = np.mean(np.diff(freqs)) + psd = np.zeros((noise_norm.size, np.sum(freqs_mask))) + n_windows = 0 + + for this_start in np.arange(start, stop, int(NFFT * (1. - overlap))): + data, _ = raw[sel, this_start:this_start + NFFT] + if data.shape[1] < NFFT: + print "Skipping last buffer" + break + + if Vh is not None: + data = np.dot(Vh, data) # reducing data rank + + data *= window[None, :] + + data_fft = fftpack.fft(data)[:, freqs_mask] + sol = np.dot(K, data_fft) + + if is_free_ori and not pick_normal: + sol = combine_xyz(sol, square=True) + else: + sol = np.abs(sol) ** 2 + + if method != "MNE": + sol *= noise_norm ** 2 + + psd += sol + n_windows += 1 + + psd /= n_windows + + psd = 10 * np.log10(psd) + + stc = _make_stc(psd, fmin * 1e-3, fstep * 1e-3, vertno) + return stc -- 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
