[med-svn] [python-mne] 225/353: ENH : add support for n_cycles varying with frequency
Yaroslav Halchenko
debian at onerussian.com
Fri Nov 27 17:25:02 UTC 2015
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 <alexandre.gramfort at inria.fr>
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
More information about the debian-med-commit
mailing list