[med-svn] [python-mne] 140/376: better induced power stat demo + addine baseline option in single_trial_power
Yaroslav Halchenko
debian at onerussian.com
Fri Nov 27 17:22:22 UTC 2015
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 0f489abe6143d100083d23b4100b2f75e72d4153
Author: Alexandre Gramfort <alexandre.gramfort at inria.fr>
Date: Wed Mar 16 10:44:11 2011 -0400
better induced power stat demo + addine baseline option in single_trial_power
---
.../plot_cluster_1samp_test_time_frequency.py | 44 ++++++++++++++--------
examples/time_frequency/plot_time_frequency.py | 4 +-
mne/time_frequency/tfr.py | 44 +++++++++++++++++++++-
3 files changed, 74 insertions(+), 18 deletions(-)
diff --git a/examples/stats/plot_cluster_1samp_test_time_frequency.py b/examples/stats/plot_cluster_1samp_test_time_frequency.py
index b562dc8..6e8b57a 100644
--- a/examples/stats/plot_cluster_1samp_test_time_frequency.py
+++ b/examples/stats/plot_cluster_1samp_test_time_frequency.py
@@ -58,28 +58,37 @@ epochs = mne.Epochs(raw, events, event_id,
tmin, tmax, picks=picks, baseline=(None, 0))
data = epochs.get_data() # as 3D matrix
data *= 1e13 # change unit to fT / cm
-
# Time vector
times = 1e3 * epochs.times # change unit to ms
-frequencies = np.arange(7, 30, 3) # define frequencies of interest
+evoked_data = np.mean(data, 0)
+
+# data -= evoked_data[None,:,:] # remove evoked component
+# evoked_data = np.mean(data, 0)
+
+frequencies = np.arange(8, 40, 2) # define frequencies of interest
Fs = raw.info['sfreq'] # sampling in Hz
-epochs_power = single_trial_power(data, Fs=Fs,
- frequencies=frequencies,
- n_cycles=3, use_fft=False, n_jobs=1)
+epochs_power = single_trial_power(data, Fs=Fs, frequencies=frequencies,
+ n_cycles=4, use_fft=False, n_jobs=1,
+ baseline=(-100, 0), times=times,
+ baseline_mode='ratio')
-epochs_power = epochs_power[:,0,:,:] # only 1 channel to get a 3D matrix
-# do ratio with baseline power:
-epochs_power /= np.mean(epochs_power[:,:,times < 0], axis=2)[:,:,None]
-# null hypothesis is that the ratio is 1 (set it to 0 for stats below)
-epochs_power -= 1.0
+# Crop in time to keep only what is between 0 and 400 ms
+time_mask = (times > 0) & (times < 400)
+epochs_power = epochs_power[:,:,:,time_mask]
+evoked_data = evoked_data[:,time_mask]
+times = times[time_mask]
+
+epochs_power = epochs_power[:,0,:,:]
+epochs_power = np.log(epochs_power) # take log of ratio
+# under the null hypothesis epochs_power should be now be 0
###############################################################################
# Compute statistic
-threshold = 5.0
+threshold = 3
T_obs, clusters, cluster_p_values, H0 = \
permutation_cluster_t_test(epochs_power,
- n_permutations=100, threshold=threshold, tail=1)
+ n_permutations=100, threshold=threshold, tail=0)
###############################################################################
# View time-frequency plots
@@ -87,7 +96,6 @@ import pylab as pl
pl.clf()
pl.subplots_adjust(0.12, 0.08, 0.96, 0.94, 0.2, 0.43)
pl.subplot(2, 1, 1)
-evoked_data = np.mean(data, 0)
pl.plot(times, evoked_data.T)
pl.title('Evoked response (%s)' % ch_name)
pl.xlabel('time (ms)')
@@ -103,13 +111,17 @@ for c, p_val in zip(clusters, cluster_p_values):
if p_val <= 0.05:
T_obs_plot[c] = T_obs[c]
+vmax = np.max(np.abs(T_obs))
+vmin = -vmax
pl.imshow(T_obs, cmap=pl.cm.gray, extent=[times[0], times[-1],
frequencies[0], frequencies[-1]],
- aspect='auto', origin='lower')
+ aspect='auto', origin='lower',
+ vmin=vmin, vmax=vmax)
pl.imshow(T_obs_plot, cmap=pl.cm.jet, extent=[times[0], times[-1],
frequencies[0], frequencies[-1]],
- aspect='auto', origin='lower')
-
+ aspect='auto', origin='lower',
+ vmin=vmin, vmax=vmax)
+pl.colorbar()
pl.xlabel('time (ms)')
pl.ylabel('Frequency (Hz)')
pl.title('Induced power (%s)' % ch_name)
diff --git a/examples/time_frequency/plot_time_frequency.py b/examples/time_frequency/plot_time_frequency.py
index 0333976..2ee160f 100644
--- a/examples/time_frequency/plot_time_frequency.py
+++ b/examples/time_frequency/plot_time_frequency.py
@@ -55,6 +55,8 @@ Fs = raw.info['sfreq'] # sampling in Hz
power, phase_lock = induced_power(data, Fs=Fs, frequencies=frequencies,
n_cycles=2, n_jobs=1, use_fft=False)
+power /= np.mean(power[:,:,times<0], axis=2)[:,:,None] # baseline ratio
+
###############################################################################
# View time-frequency plots
import pylab as pl
@@ -82,7 +84,7 @@ pl.imshow(phase_lock[0], extent=[times[0], times[-1],
frequencies[0], frequencies[-1]],
aspect='auto', origin='lower')
pl.xlabel('Time (s)')
-pl.ylabel('PLF')
+pl.ylabel('Frequency (Hz)')
pl.title('Phase-lock (%s)' % raw.info['ch_names'][picks[0]])
pl.colorbar()
pl.show()
diff --git a/mne/time_frequency/tfr.py b/mne/time_frequency/tfr.py
index 1678211..b5b4ee5 100644
--- a/mne/time_frequency/tfr.py
+++ b/mne/time_frequency/tfr.py
@@ -228,12 +228,13 @@ def _time_frequency(X, Ws, use_fft):
def single_trial_power(epochs, Fs, frequencies, use_fft=True, n_cycles=7,
+ baseline=None, baseline_mode='ratio', times=None,
n_jobs=1):
"""Compute time-frequency power on single epochs
Parameters
----------
- epochs : instance of Epochs | 3D array
+ epochs : instance of Epochs | array of shape [n_epochs x n_channels x n_times]
The epochs
Fs : float
Sampling rate
@@ -243,6 +244,21 @@ def single_trial_power(epochs, Fs, frequencies, use_fft=True, n_cycles=7,
Use the FFT for convolutions or not.
n_cycles : float
The number of cycles in the Morlet wavelet
+ baseline: None (default) or tuple of length 2
+ The time interval to apply baseline correction.
+ If None do not apply it. If baseline is (a, b)
+ the interval is between "a (s)" and "b (s)".
+ If a is None the beginning of the data is used
+ and if b is None then b is set to the end of the interval.
+ If baseline is equal ot (None, None) all the time
+ interval is used.
+ baseline_mode : None | 'ratio' | 'zscore'
+ Do baseline correction with ratio (power is divided by mean
+ power during baseline) or zscore (power is divided by standard
+ deviatio of power during baseline after substracting the mean,
+ power = [power - mean(power_baseline)] / std(power_baseline))
+ times : array
+ Required to define baseline
n_jobs : int
The number of epochs to process at the same time
@@ -268,6 +284,8 @@ def single_trial_power(epochs, Fs, frequencies, use_fft=True, n_cycles=7,
my_cwt = cwt
parallel = list
+ print "Computing time-frequency power on single epochs..."
+
power = np.empty((n_epochs, n_channels, n_frequencies, n_times),
dtype=np.float)
if n_jobs == 1:
@@ -279,6 +297,30 @@ def single_trial_power(epochs, Fs, frequencies, use_fft=True, n_cycles=7,
for k, tfr in enumerate(tfrs):
power[k] = np.abs(tfr)**2
+ # Run baseline correction
+ if baseline is not None:
+ if times is None:
+ raise ValueError('times parameter is required to define baseline')
+ print "Applying baseline correction ..."
+ bmin = baseline[0]
+ bmax = baseline[1]
+ if bmin is None:
+ imin = 0
+ else:
+ imin = int(np.where(times >= bmin)[0][0])
+ if bmax is None:
+ imax = len(times)
+ else:
+ imax = int(np.where(times <= bmax)[0][-1]) + 1
+ mean_baseline_power = np.mean(power[:,:,:,imin:imax], axis=3)
+ if baseline_mode is 'ratio':
+ power /= mean_baseline_power[:,:,:,None]
+ elif baseline_mode is 'zscore':
+ power -= mean_baseline_power[:,:,:,None]
+ power /= np.std(power[:,:,:,imin:imax], axis=3)[:,:,:,None]
+ else:
+ print "No baseline correction applied..."
+
return power
--
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