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 46ba5143a4a7f6fddb5979923c7e9f2935db5034 Author: Alexandre Gramfort <[email protected]> Date: Sun Jun 24 16:35:29 2012 +0200 ENH : add support for n_cycles varying with frequency --- .../plot_source_label_time_frequency.py | 3 +- examples/time_frequency/plot_time_frequency.py | 3 +- mne/minimum_norm/time_frequency.py | 8 +-- mne/time_frequency/tests/test_tfr.py | 3 +- mne/time_frequency/tfr.py | 62 +++++++++++++--------- 5 files changed, 48 insertions(+), 31 deletions(-) diff --git a/examples/time_frequency/plot_source_label_time_frequency.py b/examples/time_frequency/plot_source_label_time_frequency.py index c60e66a..3907f13 100644 --- a/examples/time_frequency/plot_source_label_time_frequency.py +++ b/examples/time_frequency/plot_source_label_time_frequency.py @@ -52,9 +52,10 @@ epochs = mne.Epochs(raw, events, event_id, tmin, tmax, picks=picks, # Compute a source estimate per frequency band frequencies = np.arange(7, 30, 2) # define frequencies of interest label = mne.read_label(fname_label) +n_cycles = frequencies / float(7) # different number of cycle per frequency power, phase_lock = source_induced_power(epochs, inverse_operator, frequencies, label, baseline=(-0.1, 0), baseline_mode='percent', - n_cycles=2, n_jobs=1) + n_cycles=n_cycles, n_jobs=1) power = np.mean(power, axis=0) # average over sources phase_lock = np.mean(phase_lock, axis=0) # average over sources diff --git a/examples/time_frequency/plot_time_frequency.py b/examples/time_frequency/plot_time_frequency.py index fdecfa4..66a26bd 100644 --- a/examples/time_frequency/plot_time_frequency.py +++ b/examples/time_frequency/plot_time_frequency.py @@ -52,10 +52,11 @@ data = data[:, 97:98, :] evoked_data = evoked_data[97:98, :] frequencies = np.arange(7, 30, 3) # define frequencies of interest +n_cycles = frequencies / float(7) # different number of cycle per frequency Fs = raw.info['sfreq'] # sampling in Hz decim = 3 power, phase_lock = induced_power(data, Fs=Fs, frequencies=frequencies, - n_cycles=2, n_jobs=1, use_fft=False, + n_cycles=n_cycles, n_jobs=1, use_fft=False, decim=decim) # baseline corrections with ratio diff --git a/mne/minimum_norm/time_frequency.py b/mne/minimum_norm/time_frequency.py index ffc421c..18e89d5 100644 --- a/mne/minimum_norm/time_frequency.py +++ b/mne/minimum_norm/time_frequency.py @@ -34,8 +34,8 @@ def source_band_induced_power(epochs, inverse_operator, bands, label=None, The regularization parameter of the minimum norm method: "MNE" | "dSPM" | "sLORETA" Use mininum norm, dSPM or sLORETA - n_cycles: int - Number of cycles + n_cycles: float | array of float + Number of cycles. Fixed number or one per frequency. df: float delta frequency within bands decim: int @@ -253,8 +253,8 @@ def source_induced_power(epochs, inverse_operator, frequencies, label=None, The regularization parameter of the minimum norm method: "MNE" | "dSPM" | "sLORETA" Use mininum norm, dSPM or sLORETA - n_cycles: int - Number of cycles + n_cycles: float | array of float + Number of cycles. Fixed number or one per frequency. decim: int Temporal decimation factor use_fft: bool diff --git a/mne/time_frequency/tests/test_tfr.py b/mne/time_frequency/tests/test_tfr.py index b94b0c3..ca4b34f 100644 --- a/mne/time_frequency/tests/test_tfr.py +++ b/mne/time_frequency/tests/test_tfr.py @@ -39,8 +39,9 @@ def test_time_frequency(): frequencies = np.arange(6, 20, 5) # define frequencies of interest Fs = raw.info['sfreq'] # sampling in Hz + n_cycles = frequencies / float(4) power, phase_lock = induced_power(data, Fs=Fs, frequencies=frequencies, - n_cycles=2, use_fft=True) + n_cycles=n_cycles, use_fft=True) assert_true(power.shape == (len(picks), len(frequencies), len(times))) assert_true(power.shape == phase_lock.shape) diff --git a/mne/time_frequency/tfr.py b/mne/time_frequency/tfr.py index 94df8b1..3587a0c 100644 --- a/mne/time_frequency/tfr.py +++ b/mne/time_frequency/tfr.py @@ -26,8 +26,8 @@ def morlet(Fs, freqs, n_cycles=7, sigma=None): freqs : array frequency range of interest (1 x Frequencies) - n_cycles : float - Number of oscillations if wavelet + n_cycles: float | array of float + Number of cycles. Fixed number or one per frequency. sigma : float, (optional) It controls the width of the wavelet ie its temporal @@ -44,12 +44,20 @@ def morlet(Fs, freqs, n_cycles=7, sigma=None): Wavelets time series """ Ws = list() - for f in freqs: + n_cycles = np.atleast_1d(n_cycles) + if (n_cycles.size != 1) and (n_cycles.size != len(freqs)): + raise ValueError("n_cycles should be fixed or defined for " + "each frequency.") + for k, f in enumerate(freqs): + if len(n_cycles) != 1: + this_n_cycles = n_cycles[k] + else: + this_n_cycles = n_cycles[0] # fixed or scale-dependent window if sigma is None: - sigma_t = n_cycles / (2.0 * np.pi * f) + sigma_t = this_n_cycles / (2.0 * np.pi * f) else: - sigma_t = n_cycles / (2.0 * np.pi * sigma) + sigma_t = this_n_cycles / (2.0 * np.pi * sigma) # this scaling factor is proportional to (Tallon-Baudry 98): # (sigma_t*sqrt(pi))^(-1/2); t = np.arange(0, 5 * sigma_t, 1.0 / Fs) @@ -89,6 +97,9 @@ def _cwt_fft(X, Ws, mode="same"): # precompute FFTs of Ws fft_Ws = np.empty((n_freqs, fsize), dtype=np.complex128) for i, W in enumerate(Ws): + if len(W) > n_times: + raise ValueError('Wavelet is too long for such a short signal. ' + 'Reduce the number of cycles.') fft_Ws[i] = fftn(W, [fsize]) for k, x in enumerate(X): @@ -123,6 +134,9 @@ def _cwt_convolve(X, Ws, mode='same'): tfr = np.zeros((n_freqs, n_times), dtype=np.complex128) for i, W in enumerate(Ws): ret = np.convolve(x, W, mode=mode) + if len(W) > len(x): + raise ValueError('Wavelet is too long for such a short signal. ' + 'Reduce the number of cycles.') if mode == "valid": sz = abs(W.size - n_times) + 1 offset = (n_times - sz) / 2 @@ -139,12 +153,14 @@ def cwt_morlet(X, Fs, freqs, use_fft=True, n_cycles=7.0): ---------- X : array of shape [n_signals, n_times] signals (one per line) - Fs : float sampling Frequency - freqs : array Array of frequencies of interest + use_fft : bool + Compute convolution with FFT or temoral convolution. + n_cycles: float | array of float + Number of cycles. Fixed number or one per frequency. Returns ------- @@ -178,13 +194,10 @@ def cwt(X, Ws, use_fft=True, mode='same'): ---------- X : array of shape [n_signals, n_times] signals (one per line) - Ws : list of array Wavelets time series - use_fft : bool Use FFT for convolutions - mode : 'same' | 'valid' | 'full' Convention for convolution @@ -230,14 +243,14 @@ def _time_frequency(X, Ws, use_fft): return psd, plf -def single_trial_power(epochs, Fs, frequencies, use_fft=True, n_cycles=7, +def single_trial_power(data, 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 Epochs | array of shape [n_epochs, n_channels, n_times] + data : array of shape [n_epochs, n_channels, n_times] The epochs Fs : float Sampling rate @@ -245,8 +258,9 @@ def single_trial_power(epochs, Fs, frequencies, use_fft=True, n_cycles=7, The frequencies use_fft : bool Use the FFT for convolutions or not. - n_cycles : float - The number of cycles in the Morlet wavelet + n_cycles: float | array of float + Number of cycles in the Morlet wavelet. Fixed number + or one per frequency. 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) @@ -272,7 +286,7 @@ def single_trial_power(epochs, Fs, frequencies, use_fft=True, n_cycles=7, """ mode = 'same' n_frequencies = len(frequencies) - n_epochs, n_channels, n_times = epochs.shape + n_epochs, n_channels, n_times = data.shape # Precompute wavelets for given frequency range to save time Ws = morlet(Fs, frequencies, n_cycles=n_cycles) @@ -284,11 +298,11 @@ def single_trial_power(epochs, Fs, frequencies, use_fft=True, n_cycles=7, power = np.empty((n_epochs, n_channels, n_frequencies, n_times), dtype=np.float) if n_jobs == 1: - for k, e in enumerate(epochs): + for k, e in enumerate(data): power[k] = np.abs(cwt(e, Ws, mode)) ** 2 else: # Precompute tf decompositions in parallel - tfrs = parallel(my_cwt(e, Ws, use_fft, mode) for e in epochs) + tfrs = parallel(my_cwt(e, Ws, use_fft, mode) for e in data) for k, tfr in enumerate(tfrs): power[k] = np.abs(tfr) ** 2 @@ -298,7 +312,7 @@ def single_trial_power(epochs, Fs, frequencies, use_fft=True, n_cycles=7, return power -def induced_power(epochs, Fs, frequencies, use_fft=True, n_cycles=7, +def induced_power(data, Fs, frequencies, use_fft=True, n_cycles=7, decim=1, n_jobs=1): """Compute time induced power and inter-trial phase-locking factor @@ -306,7 +320,7 @@ def induced_power(epochs, Fs, frequencies, use_fft=True, n_cycles=7, Parameters ---------- - epochs : array + data : array 3D array of shape [n_epochs, n_channels, n_times] Fs : float sampling Frequency @@ -315,8 +329,8 @@ def induced_power(epochs, Fs, frequencies, use_fft=True, n_cycles=7, use_fft : bool Compute transform with fft based convolutions or temporal convolutions. - n_cycles : int - The number of cycles in the wavelet + n_cycles: float | array of float + Number of cycles. Fixed number or one per frequency. decim: int Temporal decimation factor n_jobs : int @@ -332,7 +346,7 @@ def induced_power(epochs, Fs, frequencies, use_fft=True, n_cycles=7, Phase locking factor in [0, 1] (Channels x Frequencies x Timepoints) """ n_frequencies = len(frequencies) - n_epochs, n_channels, n_times = epochs[:, :, ::decim].shape + n_epochs, n_channels, n_times = data[:, :, ::decim].shape # Precompute wavelets for given frequency range to save time Ws = morlet(Fs, frequencies, n_cycles=n_cycles) @@ -342,13 +356,13 @@ def induced_power(epochs, Fs, frequencies, use_fft=True, n_cycles=7, plf = np.empty((n_channels, n_frequencies, n_times), dtype=np.complex) for c in range(n_channels): - X = np.squeeze(epochs[:, c, :]) + X = np.squeeze(data[:, c, :]) this_psd, this_plf = _time_frequency(X, Ws, use_fft) psd[c], plf[c] = this_psd[:, ::decim], this_plf[:, ::decim] else: parallel, my_time_frequency, _ = parallel_func(_time_frequency, n_jobs) - psd_plf = parallel(my_time_frequency(np.squeeze(epochs[:, c, :]), + psd_plf = parallel(my_time_frequency(np.squeeze(data[:, c, :]), Ws, use_fft) for c in range(n_channels)) -- 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
