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 55b831b6b741cb1da16d983e7798e31d3e0680b2 Author: Alexandre Gramfort <[email protected]> Date: Thu Feb 3 13:52:40 2011 -0500 adding routines to compute simple time frequency analysis with morlet wavelet --- examples/plot_time_frequency.py | 85 +++++++++++++++++ mne/__init__.py | 1 + mne/tfr.py | 204 ++++++++++++++++++++++++++++++++++++++++ 3 files changed, 290 insertions(+) diff --git a/examples/plot_time_frequency.py b/examples/plot_time_frequency.py new file mode 100644 index 0000000..b53eb4a --- /dev/null +++ b/examples/plot_time_frequency.py @@ -0,0 +1,85 @@ +""" +========================================================= +Time frequency : Induced power and inter-trial phase-lock +========================================================= + +This script shows how to compute induced power and inter-trial +phase-lock for a list of epochs read in a raw file given +a list of events. + +""" +# Authors: Alexandre Gramfort <[email protected]> +# +# License: BSD (3-clause) + +print __doc__ + +import os +import numpy as np + +import mne +from mne import fiff +from mne import time_frequency + +############################################################################### +# Set parameters +raw_fname = os.environ['MNE_SAMPLE_DATASET_PATH'] +raw_fname += '/MEG/sample/sample_audvis_raw.fif' +event_fname = os.environ['MNE_SAMPLE_DATASET_PATH'] +event_fname += '/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]] +data, times, channel_names = mne.read_epochs(raw, events, event_id, + tmin, tmax, picks=picks, baseline=(None, 0)) +epochs = np.array([d['epoch'] for d in data]) # as 3D matrix +evoked_data = np.mean(epochs, axis=0) # compute evoked fields + +frequencies = np.arange(4, 30, 3) # define frequencies of interest +Fs = raw['info']['sfreq'] # sampling in Hz +power, phase_lock = time_frequency(epochs, Fs=Fs, frequencies=frequencies, + n_cycles=2) + +############################################################################### +# View time-frequency plots +import pylab as pl +pl.clf() +pl.subplots_adjust(0.1, 0.08, 0.98, 0.94, 0.2, 0.43) +pl.subplot(3, 1, 1) +pl.plot(times, evoked_data.T) +pl.title('Evoked response (%s)' % raw['info']['ch_names'][picks[0]]) +pl.xlabel('time (s)') +pl.ylabel('T / m') +pl.xlim(times[0], times[-1]) + +pl.subplot(3, 1, 2) +pl.imshow(20*np.log10(power[0]), extent=[times[0], times[-1], + frequencies[0], frequencies[-1]], + aspect='auto', origin='lower') +pl.xlabel('Time (s)') +pl.ylabel('Frequency (Hz)') +pl.title('Induced power (%s)' % raw['info']['ch_names'][picks[0]]) +pl.colorbar() + +pl.subplot(3, 1, 3) +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.title('Phase-lock (%s)' % raw['info']['ch_names'][picks[0]]) +pl.colorbar() +pl.show() diff --git a/mne/__init__.py b/mne/__init__.py index 95ae1fb..f6e87f4 100644 --- a/mne/__init__.py +++ b/mne/__init__.py @@ -7,4 +7,5 @@ from .stc import read_stc, write_stc from .bem_surfaces import read_bem_surfaces from .inverse import read_inverse_operator, compute_inverse from .epochs import read_epochs +from .tfr import time_frequency import fiff diff --git a/mne/tfr.py b/mne/tfr.py new file mode 100644 index 0000000..0fce274 --- /dev/null +++ b/mne/tfr.py @@ -0,0 +1,204 @@ +"""A module which implements the continuous wavelet transform +with complex Morlet wavelets. + +Author : Alexandre Gramfort, [email protected] (2011) +License : BSD 3-clause + +inspired by Matlab code from Sheraz Khan & Brainstorm & SPM +""" + +from math import sqrt +import numpy as np +from scipy import linalg +from scipy.fftpack import fftn, ifftn + + +def morlet(Fs, freqs, n_cycles=7, sigma=None): + """Compute Wavelets for the given frequency range + + Parameters + ---------- + Fs : float + Sampling Frequency + + freqs : array + frequency range of interest (1 x Frequencies) + + n_cycles : float + Number of oscillations if wavelet + + sigma : float, (optional) + It controls the width of the wavelet ie its temporal + resolution. If sigma is None the temporal resolution + is adapted with the frequency like for all wavelet transform. + The higher the frequency the shorter is the wavelet. + If sigma is fixed the temporal resolution is fixed + like for the short time Fourier transform and the number + of oscillations increases with the frequency. + + Returns + ------- + Ws : list of array + Wavelets time series + """ + Ws = list() + for f in freqs: + # fixed or scale-dependent window + if sigma is None: + sigma_t = n_cycles / (2.0 * np.pi * f) + else: + sigma_t = 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) + t = np.r_[-t[::-1], t[1:]] + W = np.exp(2.0 * 1j * np.pi * f *t) + W *= np.exp(-t**2 / (2.0 * sigma_t**2)) + W /= sqrt(0.5) * linalg.norm(W.ravel()) + Ws.append(W) + return Ws + + +def _centered(arr, newsize): + # Return the center newsize portion of the array. + newsize = np.asarray(newsize) + currsize = np.array(arr.shape) + startind = (currsize - newsize) / 2 + endind = startind + newsize + myslice = [slice(startind[k], endind[k]) for k in range(len(endind))] + return arr[tuple(myslice)] + + +def _cwt_morlet_fft(x, Fs, freqs, mode="same", Ws=None): + """Compute cwt with fft based convolutions + """ + x = np.asarray(x) + freqs = np.asarray(freqs) + + # Precompute wavelets for given frequency range to save time + n_samples = x.size + n_freqs = freqs.size + + if Ws is None: + Ws = morlet(Fs, freqs) + + Ws_max_size = max(W.size for W in Ws) + size = n_samples + Ws_max_size - 1 + # Always use 2**n-sized FFT + fsize = 2**np.ceil(np.log2(size)) + fft_x = fftn(x, [fsize]) + + if mode == "full": + tfr = np.zeros((n_freqs, fsize), dtype=np.complex128) + elif mode == "same" or mode == "valid": + tfr = np.zeros((n_freqs, n_samples), dtype=np.complex128) + + for i, W in enumerate(Ws): + ret = ifftn(fft_x * fftn(W, [fsize]))[:n_samples + W.size - 1] + if mode == "valid": + sz = abs(W.size - n_samples) + 1 + offset = (n_samples - sz) / 2 + tfr[i, offset:(offset + sz)] = _centered(ret, sz) + else: + tfr[i] = _centered(ret, n_samples) + return tfr + + +def _cwt_morlet_convolve(x, Fs, freqs, mode='same', Ws=None): + """Compute time freq decomposition with temporal convolutions + """ + x = np.asarray(x) + freqs = np.asarray(freqs) + + if Ws is None: + Ws = morlet(Fs, freqs) + + n_samples = x.size + # Compute convolutions + tfr = np.zeros((freqs.size, len(x)), dtype=np.complex128) + for i, W in enumerate(Ws): + ret = np.convolve(x, W, mode=mode) + if mode == "valid": + sz = abs(W.size - n_samples) + 1 + offset = (n_samples - sz) / 2 + tfr[i, offset:(offset + sz)] = ret + else: + tfr[i] = ret + return tfr + + +def cwt_morlet(x, Fs, freqs, use_fft=True, n_cycles=7.0): + """Compute time freq decomposition with Morlet wavelets + + Parameters + ---------- + x : array + signal + + Fs : float + sampling Frequency + + freqs : array + Array of frequencies of interest + + Returns + ------- + tfr : 2D array + Time Frequency Decomposition (Frequencies x Timepoints) + """ + mode = 'same' + # mode = "valid" + + # Precompute wavelets for given frequency range to save time + Ws = morlet(Fs, freqs, n_cycles=n_cycles) + + if use_fft: + return _cwt_morlet_fft(x, Fs, freqs, mode, Ws) + else: + return _cwt_morlet_convolve(x, Fs, freqs, mode, Ws) + + +def time_frequency(epochs, Fs, frequencies, use_fft=True, n_cycles=25): + """Compute time induced power and inter-trial phase-locking factor + + The time frequency decomposition is done with Morlet wavelets + + Parameters + ---------- + epochs : array + 3D array of shape [n_epochs, n_channels, n_times] + + Fs : float + sampling Frequency + + frequencies : array + Array of frequencies of interest + + use_fft : bool + Compute transform with fft based convolutions or temporal + convolutions. + + n_cycles : int + The number of cycles in the wavelet + + Returns + ------- + power : 2D array + Induced power (Channels x Frequencies x Timepoints) + phase_lock : 2D array + Phase locking factor in [0, 1] (Channels x Frequencies x Timepoints) + """ + n_frequencies = len(frequencies) + n_epochs, n_channels, n_times = epochs.shape + psd = np.zeros((n_channels, n_frequencies, n_times)) # PSD + plf = np.zeros((n_channels, n_frequencies, n_times)) # phase lock + for c in range(n_channels): + for e in range(n_epochs): + tfr = cwt_morlet(epochs[e, c, :].ravel(), Fs, frequencies, + use_fft=use_fft, n_cycles=n_cycles) + psd[c,:,:] += np.abs(tfr) + plf[c,:,:] += tfr / psd[c,:,:] + + psd /= n_epochs + plf = np.abs(plf) / n_epochs + return psd, plf -- 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
