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