[med-svn] [python-mne] 319/353: ENH : add function to compute source PSD from raw file
Yaroslav Halchenko
debian at onerussian.com
Fri Nov 27 17:25:24 UTC 2015
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 <alexandre.gramfort at inria.fr>
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 <gramfort at nmr.mgh.harvard.edu>
+#
+# 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
More information about the debian-med-commit
mailing list