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