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 740f72ba1131fdc9c30a73034cb3b58d98df8c3e Author: Alexandre Gramfort <[email protected]> Date: Mon Mar 7 11:38:06 2011 -0500 new example of TF cluster --- examples/stats/plot_cluster_stats_evoked.py | 21 ++-- .../stats/plot_cluster_stats_time_frequency.py | 117 +++++++++++++++++++++ mne/tfr.py | 5 +- 3 files changed, 131 insertions(+), 12 deletions(-) diff --git a/examples/stats/plot_cluster_stats_evoked.py b/examples/stats/plot_cluster_stats_evoked.py index 53dbe64..ff8e390 100644 --- a/examples/stats/plot_cluster_stats_evoked.py +++ b/examples/stats/plot_cluster_stats_evoked.py @@ -15,8 +15,6 @@ with cluster level permutation test. print __doc__ -import numpy as np - import mne from mne import fiff from mne.stats import permutation_cluster_test @@ -42,14 +40,17 @@ include = [channel] # Read epochs for the channel of interest picks = fiff.pick_types(raw['info'], meg=False, include=include) event_id = 1 -data1, times, channel_names = mne.read_epochs(raw, events, event_id, +epochs1 = mne.Epochs(raw, events, event_id, tmin, tmax, picks=picks, baseline=(None, 0)) -condition1 = np.squeeze(np.array([d['epoch'] for d in data1])) # as 3D matrix +condition1 = epochs1.get_data() # as 3D matrix event_id = 2 -data2, times, channel_names = mne.read_epochs(raw, events, event_id, +epochs2 = mne.Epochs(raw, events, event_id, tmin, tmax, picks=picks, baseline=(None, 0)) -condition2 = np.squeeze(np.array([d['epoch'] for d in data2])) # as 3D matrix +condition2 = epochs2.get_data() # as 3D matrix + +condition1 = condition1[:,0,:] # take only one channel to get a 2D array +condition2 = condition2[:,0,:] # take only one channel to get a 2D array ############################################################################### # Compute statistic @@ -60,6 +61,7 @@ T_obs, clusters, cluster_p_values, H0 = \ ############################################################################### # Plot +times = epochs1.times import pylab as pl pl.close('all') pl.subplot(211) @@ -69,11 +71,12 @@ pl.plot(times, condition1.mean(axis=0) - condition2.mean(axis=0), pl.ylabel("MEG (T / m)") pl.legend() pl.subplot(212) -for i_c, (start, stop) in enumerate(clusters): +for i_c, c in enumerate(clusters): + c = c[0] if cluster_p_values[i_c] <= 0.05: - h = pl.axvspan(times[start], times[stop-1], color='r', alpha=0.3) + h = pl.axvspan(times[c.start], times[c.stop-1], color='r', alpha=0.3) else: - pl.axvspan(times[start], times[stop-1], color=(0.3, 0.3, 0.3), + pl.axvspan(times[c.start], times[c.stop-1], color=(0.3, 0.3, 0.3), alpha=0.3) hf = pl.plot(times, T_obs, 'g') pl.legend((h, ), ('cluster p-value < 0.05', )) diff --git a/examples/stats/plot_cluster_stats_time_frequency.py b/examples/stats/plot_cluster_stats_time_frequency.py new file mode 100644 index 0000000..a5793b0 --- /dev/null +++ b/examples/stats/plot_cluster_stats_time_frequency.py @@ -0,0 +1,117 @@ +""" +====================================================== +Non-parametric cluster statistic on single trial power +====================================================== + +This script shows how to estimate significant clusters +in time-frequency power estimates. It uses a non-parametric +statistical procedure based on permutations and cluster +level statistics. + +""" +# Authors: Alexandre Gramfort <[email protected]> +# +# License: BSD (3-clause) + +print __doc__ + +import numpy as np + +import mne +from mne import fiff +from mne.tfr import single_trial_power +from mne.stats import permutation_cluster_test +from mne.datasets import sample + +############################################################################### +# Set parameters +data_path = sample.data_path('..') +raw_fname = data_path + '/MEG/sample/sample_audvis_raw.fif' +event_fname = data_path + '/MEG/sample/sample_audvis_raw-eve.fif' +event_id = 1 +tmin = -0.2 +tmax = 0.5 + +# Setup for reading the raw data +raw = fiff.setup_read_raw(raw_fname) +events = mne.read_events(event_fname) + +include = [] +exclude = raw['info']['bads'] + ['MEG 2443', 'EEG 053'] # bads + 2 more + +# picks MEG gradiometers +picks = fiff.pick_types(raw['info'], meg='grad', eeg=False, + stim=False, include=include, exclude=exclude) + +picks = [picks[97]] +ch_name = raw['info']['ch_names'][picks[0]] + +# Load condition 1 +event_id = 1 +epochs_condition_1 = mne.Epochs(raw, events, event_id, + tmin, tmax, picks=picks, baseline=(None, 0)) +data_condition_1 = epochs_condition_1.get_data() # as 3D matrix +data_condition_1 *= 1e13 # change unit to fT / cm + +# Load condition 1 +event_id = 2 +epochs_condition_2 = mne.Epochs(raw, events, event_id, + tmin, tmax, picks=picks, baseline=(None, 0)) +data_condition_2 = epochs_condition_2.get_data() # as 3D matrix +data_condition_2 *= 1e13 # change unit to fT / cm + +# Time vector +times = 1e3 * epochs_condition_1.times # change unit to ms + +frequencies = np.arange(7, 30, 3) # define frequencies of interest +Fs = raw['info']['sfreq'] # sampling in Hz +epochs_coefs_1 = single_trial_power(data_condition_1, Fs=Fs, + frequencies=frequencies, + n_cycles=2, use_fft=False) + +epochs_coefs_2 = single_trial_power(data_condition_2, Fs=Fs, + frequencies=frequencies, + n_cycles=2, use_fft=False) + +epochs_coefs_1 = epochs_coefs_1[:,0,:,:] # only 1 channel to get a 3D matrix +epochs_coefs_2 = epochs_coefs_2[:,0,:,:] # only 1 channel to get a 3D matrix + ############################################################################### +# Compute statistic +threshold = 6.0 +T_obs, clusters, cluster_p_values, H0 = \ + permutation_cluster_test([epochs_coefs_1, epochs_coefs_2], + n_permutations=100, threshold=threshold, tail=0) + +############################################################################### +# View time-frequency plots +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_contrast = np.mean(data_condition_1, 0) - np.mean(data_condition_2, 0) +pl.plot(times, evoked_contrast.T) +pl.title('Contrast of evoked response (%s)' % ch_name) +pl.xlabel('time (ms)') +pl.ylabel('Magnetic Field (fT/cm)') +pl.xlim(times[0], times[-1]) +pl.ylim(-100, 200) + +pl.subplot(2, 1, 2) + +# Create new stats image with only significant clusters +T_obs_plot = np.nan * np.ones_like(T_obs) +for c, p_val in zip(clusters, cluster_p_values): + if p_val <= 0.05: + T_obs_plot[c] = T_obs[c] + +pl.imshow(T_obs, cmap=pl.cm.gray, extent=[times[0], times[-1], + frequencies[0], frequencies[-1]], + aspect='auto', origin='lower') +pl.imshow(T_obs_plot, cmap=pl.cm.jet, extent=[times[0], times[-1], + frequencies[0], frequencies[-1]], + aspect='auto', origin='lower') + +pl.xlabel('time (ms)') +pl.ylabel('Frequency (Hz)') +pl.title('Induced power (%s)' % ch_name) +pl.show() diff --git a/mne/tfr.py b/mne/tfr.py index 2ac0b55..986f45b 100644 --- a/mne/tfr.py +++ b/mne/tfr.py @@ -189,8 +189,7 @@ def _time_frequency(X, Ws, use_fft): return psd, plf -def single_trial_power(epochs, Fs, frequencies, use_fft=True, n_cycles=25, - n_jobs=1): +def single_trial_power(epochs, Fs, frequencies, use_fft=True, n_cycles=7): """Compute time-frequency power on single epochs """ n_frequencies = len(frequencies) @@ -210,7 +209,7 @@ def single_trial_power(epochs, Fs, frequencies, use_fft=True, n_cycles=25, for k, e in enumerate(epochs): mode = 'same' - power[k] = np.abs(_cwt(e, Ws, mode))**2 + power[k] = np.abs(list(_cwt(e, Ws, mode)))**2 return power -- 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
