[med-svn] [python-mne] 274/353: ENH : cleanup and remove duplicated gabor atoms code
Yaroslav Halchenko
debian at onerussian.com
Fri Nov 27 17:25:14 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 0a4c44ba414a65da877e4e36ba53022b0cbd9bb3
Author: Alexandre Gramfort <alexandre.gramfort at inria.fr>
Date: Mon Jul 16 16:14:13 2012 +0200
ENH : cleanup and remove duplicated gabor atoms code
---
examples/plot_simulate_evoked_data.py | 41 +++++++-------
mne/simulation/__init__.py | 3 +-
mne/simulation/sim_evoked.py | 101 ++++++++--------------------------
mne/time_frequency/__init__.py | 2 +-
4 files changed, 47 insertions(+), 100 deletions(-)
diff --git a/examples/plot_simulate_evoked_data.py b/examples/plot_simulate_evoked_data.py
index cf0f3ee..788a12a 100644
--- a/examples/plot_simulate_evoked_data.py
+++ b/examples/plot_simulate_evoked_data.py
@@ -16,9 +16,10 @@ import mne
from mne.fiff.pick import pick_types_evoked, pick_types_forward
from mne.forward import apply_forward
from mne.datasets import sample
-from mne.time_frequency import fir_filter_raw
+from mne.time_frequency import fir_filter_raw, morlet
from mne.viz import plot_evoked, plot_sparse_source_estimates
-from mne.simulation.sim_evoked import source_signal, generate_stc, generate_noise_evoked, add_noise
+from mne.simulation import generate_stc, generate_noise_evoked, \
+ add_noise_evoked
###############################################################################
# Load real data as templates
@@ -42,26 +43,29 @@ evoked_template = mne.fiff.read_evoked(ave_fname, setno=0, baseline=None)
evoked_template = pick_types_evoked(evoked_template, meg=True, eeg=True,
exclude=raw.info['bads'])
-tmin = -0.1
-sfreq = 1000 # Hz
-tstep = 1. / sfreq
-n_samples = 300
-timesamples = np.linspace(tmin, tmin + n_samples * tstep, n_samples)
-
label_names = ['Aud-lh', 'Aud-rh']
labels = [mne.read_label(data_path + '/MEG/sample/labels/%s.label' % ln)
for ln in label_names]
-mus = [[0.030, 0.060, 0.120], [0.040, 0.060, 0.140]]
-sigmas = [[0.01, 0.02, 0.03], [0.01, 0.02, 0.03]]
-amps = [[40 * 1e-9, 40 * 1e-9, 30 * 1e-9], [30 * 1e-9, 40 * 1e-9, 40 * 1e-9]]
-freqs = [[0, 0, 0], [0, 0, 0]]
-phis = [[0, 0, 0], [0, 0, 0]]
+###############################################################################
+# Generate source time courses and the correspond evoked data
+snr = 6 # dB
+tmin = -0.1
+sfreq = 1000. # Hz
+tstep = 1. / sfreq
+n_samples = 600
+times = np.linspace(tmin, tmin + n_samples * tstep, n_samples)
+
+# Generate times series from 2 Morlet wavelets
+stc_data = np.zeros((len(labels), len(times)))
+Ws = morlet(sfreq, [3, 10], n_cycles=[1, 1.5])
+stc_data[0][:len(Ws[0])] = np.real(Ws[0])
+stc_data[1][:len(Ws[1])] = np.real(Ws[1])
+stc_data *= 100 * 1e-9 # use nAm as unit
-SNR = 6
-dB = True
+# time translation
+stc_data[1] = np.roll(stc_data[1], 80)
-stc_data = source_signal(mus, sigmas, amps, freqs, phis, timesamples)
stc = generate_stc(fwd, labels, stc_data, tmin, tstep, random_state=0)
evoked = apply_forward(fwd, stc, evoked_template)
@@ -70,8 +74,7 @@ evoked = apply_forward(fwd, stc, evoked_template)
picks = mne.fiff.pick_types(raw.info, meg=True)
fir_filter = fir_filter_raw(raw, order=5, picks=picks, tmin=60, tmax=180)
noise = generate_noise_evoked(evoked, noise_cov, n_samples, fir_filter)
-
-evoked_noise = add_noise(evoked, noise, SNR, timesamples, tmin=0.0, tmax=0.2, dB=dB)
+evoked_noise = add_noise_evoked(evoked, noise, snr, times, tmin=0.0, tmax=0.2)
###############################################################################
# Plot
@@ -82,4 +85,4 @@ pl.figure()
pl.psd(evoked_noise.data[0])
pl.figure()
-plot_evoked(evoked)
+plot_evoked(evoked_noise)
diff --git a/mne/simulation/__init__.py b/mne/simulation/__init__.py
index 918184b..64fdc38 100644
--- a/mne/simulation/__init__.py
+++ b/mne/simulation/__init__.py
@@ -1,4 +1,5 @@
"""Data simulation code
"""
-from .sim_evoked import select_source_in_label, generate_stc
\ No newline at end of file
+from .sim_evoked import select_source_in_label, generate_stc, \
+ generate_noise_evoked, add_noise_evoked
diff --git a/mne/simulation/sim_evoked.py b/mne/simulation/sim_evoked.py
index 97e4b72..4924517 100644
--- a/mne/simulation/sim_evoked.py
+++ b/mne/simulation/sim_evoked.py
@@ -13,66 +13,11 @@ from ..minimum_norm.inverse import _make_stc
from ..utils import check_random_state
-def gaboratomr(timesamples, sigma, mu, k, phase):
- """Computes a real-valued Gabor atom
+def generate_noise_evoked(evoked, noise_cov, n_samples, fir_filter=None,
+ random_state=None):
+ """Creates noise as a multivariate Gaussian
- Parameters
- ----------
- timesamples : array
- samples in seconds
- sigma : float
- the variance of the gauss function.
- mu : float
- the mean of the gauss function.
- mu : float
- number of modulation of the cosine function.
- phase : float
- the phase of the modulated cosine function.
-
- Returns
- -------
- gnorm : array
- real_valued gabor atom with amplitude = 1
- """
- N = len(timesamples)
- g = 1.0 / (sigma * np.sqrt(2 * np.pi)) * np.exp(-0.5 * ((timesamples - mu) / sigma) ** 2) *\
- np.cos(2 * np.pi * k / N * np.arange(0, N) + phase)
- gnorm = g / np.max(np.abs(g))
- return gnorm
-
-
-def source_signal(mus, sigmas, amps, freqs, phis, timesamples):
- """Simulates source signal as sum of Gabor atoms
-
- Parameters
- ----------
- mu : list
- the means of the gauss functions.
- sigma : list
- the variances of the gauss functions.
- amps : list
- amplitudes of the Gabor atoms.
- freqs : list
- numbers of modulation of the cosine function.
- phase : list
- the phases of the modulated cosine function.
- timesamples : array
- samples in seconds
-
- Returns
- -------
- signal : array
- simulated source signal
- """
- data = np.zeros((len(mus), len(timesamples)))
- for k in range(len(mus)):
- for m, s, a, f, p in zip(mus[k], sigmas[k], amps[k], freqs[k], phis[k]):
- data[k] += gaboratomr(timesamples, s, m, f, p) * a
- return data
-
-
-def generate_noise_evoked(evoked, noise_cov, n_samples, fir_filter=None, random_state=None):
- """Creates noise as a multivariate random process with specified cov matrix.
+ The spatial covariance of the noise is given from the cov matrix.
Parameters
----------
@@ -96,15 +41,17 @@ def generate_noise_evoked(evoked, noise_cov, n_samples, fir_filter=None, random_
noise_cov = pick_channels_cov(noise_cov, include=noise.info['ch_names'])
rng = check_random_state(random_state)
n_channels = np.zeros(noise.info['nchan'])
- noise.data = rng.multivariate_normal(n_channels, noise_cov.data, n_samples).T
+ noise.data = rng.multivariate_normal(n_channels, noise_cov.data,
+ n_samples).T
if fir_filter is not None:
noise.data = signal.lfilter([1], fir_filter, noise.data, axis=-1)
return noise
-def add_noise(evoked, noise, SNR, timesamples, tmin=None, tmax=None, dB=False):
- """Adds noise to evoked object with specified SNR. SNR is computed in the
- interval from tmin to tmax. No deepcopy of evoked applied.
+def add_noise_evoked(evoked, noise, snr, times, tmin=None, tmax=None):
+ """Adds noise to evoked object with specified SNR.
+
+ SNR is computed in the interval from tmin to tmax.
Parameters
----------
@@ -112,35 +59,31 @@ def add_noise(evoked, noise, SNR, timesamples, tmin=None, tmax=None, dB=False):
an instance of evoked with signal
noise : evoked object
an instance of evoked with noise
- SNR : float
- signal to noise ratio
+ snr : float
+ signal to noise ratio in dB. It corresponds to
+ 10 * log10( var(signal) / var(noise) )
timesamples : array
samples in seconds
tmin : float
start time before event
tmax : float
end time after event
- dB : bool
- SNR in dB or not
Returns
-------
evoked : evoked object
- an instance of evoked
+ An instance of evoked corrupted by noise
"""
+ evoked = copy.deepcopy(evoked)
if tmin is None:
- tmin = np.min(timesamples)
+ tmin = np.min(times)
if tmax is None:
- tmax = np.max(timesamples)
- tmask = (timesamples >= tmin) & (timesamples <= tmax)
- if dB:
- SNRtemp = 20 * np.log10(np.sqrt(np.mean((evoked.data[:, tmask] ** 2).ravel()) / \
- np.mean((noise.data ** 2).ravel())))
- noise.data = 10 ** ((SNRtemp - float(SNR)) / 20) * noise.data
- else:
- SNRtemp = np.sqrt(np.mean((evoked.data[:, tmask] ** 2).ravel()) / \
- np.mean((noise.data ** 2).ravel()))
- noise.data = SNRtemp / SNR * noise.data
+ tmax = np.max(times)
+ tmask = (times >= tmin) & (times <= tmax)
+ tmp = np.mean((evoked.data[:, tmask] ** 2).ravel()) / \
+ np.mean((noise.data ** 2).ravel())
+ tmp = 10 * np.log10(tmp)
+ noise.data = 10 ** ((tmp - float(snr)) / 20) * noise.data
evoked.data += noise.data
return evoked
diff --git a/mne/time_frequency/__init__.py b/mne/time_frequency/__init__.py
index be88845..cee6d90 100644
--- a/mne/time_frequency/__init__.py
+++ b/mne/time_frequency/__init__.py
@@ -1,6 +1,6 @@
"""Time frequency analysis tools
"""
-from .tfr import induced_power, single_trial_power
+from .tfr import induced_power, single_trial_power, morlet
from .psd import compute_raw_psd
from .ar import yule_walker, ar_raw, fir_filter_raw
--
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