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 e915ad04a348c03ec1a7379454e142398d60e336 Author: Alexandre Gramfort <[email protected]> Date: Mon Jun 25 15:10:03 2012 +0200 ENH : add basic support for temporal whitening with AR model --- doc/source/whats_new.rst | 2 + examples/time_frequency/plot_temporal_whitening.py | 63 ++++++++++++ mne/time_frequency/__init__.py | 1 + mne/time_frequency/ar.py | 111 +++++++++++++++++++++ mne/time_frequency/tests/test_ar.py | 42 ++++++++ 5 files changed, 219 insertions(+) diff --git a/doc/source/whats_new.rst b/doc/source/whats_new.rst index 6dc2680..1a6d0aa 100644 --- a/doc/source/whats_new.rst +++ b/doc/source/whats_new.rst @@ -7,6 +7,8 @@ Current Changelog ~~~~~~~~~ + - Fit AR models to raw data for temporal whitening by `Alex Gramfort`_. + - speedup + reduce memory of mne.morph_data by `Alex Gramfort`_. - Backporting scipy.signa.firwin2 so filtering works with old scipy by `Alex Gramfort`_. diff --git a/examples/time_frequency/plot_temporal_whitening.py b/examples/time_frequency/plot_temporal_whitening.py new file mode 100644 index 0000000..eeebbfd --- /dev/null +++ b/examples/time_frequency/plot_temporal_whitening.py @@ -0,0 +1,63 @@ +""" +================================ +Temporal whitening with AR model +================================ + +This script shows how to fit an AR model to data and use it +to temporally whiten the signals. + +""" +# Authors: Alexandre Gramfort <[email protected]> +# +# License: BSD (3-clause) + +print __doc__ + +import numpy as np +from scipy import signal +import pylab as pl + +import mne +from mne.time_frequency import ar_raw +from mne.datasets import sample +data_path = sample.data_path('..') + +raw_fname = data_path + '/MEG/sample/sample_audvis_raw.fif' +proj_fname = data_path + '/MEG/sample/sample_audvis_ecg_proj.fif' + +raw = mne.fiff.Raw(raw_fname) +proj = mne.read_proj(proj_fname) +raw.info['projs'] += proj +raw.info['bads'] = ['MEG 2443', 'EEG 053'] # mark bad channels + +# Set up pick list: Gradiometers - bad channels +picks = mne.fiff.pick_types(raw.info, meg='grad', exclude=raw.info['bads']) + +order = 5 # define model order +picks = picks[:5] + +# Estimate AR models on raw data +coefs = ar_raw(raw, order=order, picks=picks, tmin=60, tmax=180) +mean_coefs = np.mean(coefs, axis=0) # mean model accross channels + +filt = np.r_[1, -mean_coefs] # filter coefficient +d, times = raw[0, 1e4:2e4] # look at one channel from now on +d = d.ravel() # make flat vector +innovation = signal.convolve(d, filt, 'valid') +d_ = signal.lfilter([1], filt, innovation) # regenerate the signal +d_ = np.r_[d_[0] * np.ones(order), d_] # dummy samples to keep signal length + +############################################################################### +# Plot the different time series and PSDs +pl.close('all') +pl.figure() +pl.plot(d[:100], label='signal') +pl.plot(d_[:100], label='regenerated signal') +pl.legend() + +pl.figure() +pl.psd(d, Fs=raw.info['sfreq'], NFFT=2048) +pl.psd(innovation, Fs=raw.info['sfreq'], NFFT=2048) +pl.psd(d_, Fs=raw.info['sfreq'], NFFT=2048, linestyle='--') +pl.legend(('Signal', 'Innovation', 'Regenerated signal')) +pl.show() diff --git a/mne/time_frequency/__init__.py b/mne/time_frequency/__init__.py index 34395c5..4826123 100644 --- a/mne/time_frequency/__init__.py +++ b/mne/time_frequency/__init__.py @@ -3,3 +3,4 @@ from .tfr import induced_power, single_trial_power from .psd import compute_raw_psd +from .ar import yule_walker, ar_raw diff --git a/mne/time_frequency/ar.py b/mne/time_frequency/ar.py new file mode 100644 index 0000000..3daafc1 --- /dev/null +++ b/mne/time_frequency/ar.py @@ -0,0 +1,111 @@ +# Authors: Alexandre Gramfort <[email protected]> +# The statsmodels folks for AR yule_walker +# +# License: BSD (3-clause) + +import numpy as np +from scipy.linalg import toeplitz + + +# XXX : Back ported from statsmodels + +def yule_walker(X, order=1, method="unbiased", df=None, inv=False, demean=True): + """ + Estimate AR(p) parameters from a sequence X using Yule-Walker equation. + + Unbiased or maximum-likelihood estimator (mle) + + See, for example: + + http://en.wikipedia.org/wiki/Autoregressive_moving_average_model + + Parameters + ---------- + X : array-like + 1d array + order : integer, optional + The order of the autoregressive process. Default is 1. + method : string, optional + Method can be "unbiased" or "mle" and this determines denominator in + estimate of autocorrelation function (ACF) at lag k. If "mle", the + denominator is n=X.shape[0], if "unbiased" the denominator is n-k. + The default is unbiased. + df : integer, optional + Specifies the degrees of freedom. If `df` is supplied, then it is assumed + the X has `df` degrees of freedom rather than `n`. Default is None. + inv : bool + If inv is True the inverse of R is also returned. Default is False. + demean : bool + True, the mean is subtracted from `X` before estimation. + + Returns + ------- + rho + The autoregressive coefficients + sigma + TODO + + """ +#TODO: define R better, look back at notes and technical notes on YW. +#First link here is useful +#http://www-stat.wharton.upenn.edu/~steele/Courses/956/ResourceDetails/YuleWalkerAndMore.htm + method = str(method).lower() + if method not in ["unbiased", "mle"]: + raise ValueError("ACF estimation method must be 'unbiased' or 'MLE'") + X = np.array(X) + if demean: + X -= X.mean() # automatically demean's X + n = df or X.shape[0] + + if method == "unbiased": # this is df_resid ie., n - p + denom = lambda k: n - k + else: + denom = lambda k: n + if X.ndim > 1 and X.shape[1] != 1: + raise ValueError("expecting a vector to estimate AR parameters") + r = np.zeros(order+1, np.float64) + r[0] = (X**2).sum() / denom(0) + for k in range(1,order+1): + r[k] = (X[0:-k]*X[k:]).sum() / denom(k) + R = toeplitz(r[:-1]) + + rho = np.linalg.solve(R, r[1:]) + sigmasq = r[0] - (r[1:]*rho).sum() + if inv == True: + return rho, np.sqrt(sigmasq), np.linalg.inv(R) + else: + return rho, np.sqrt(sigmasq) + + +def ar_raw(raw, order, picks, tmin=None, tmax=None): + """Fit AR model on raw data + + Fit AR models for each channels and returns the models + coefficients for each of them. + + Parameters + ---------- + raw : Raw instance + The raw data + order : int + The AR model order + picks : array of int + The channels indices to include + tmin : float + The beginning of time interval in seconds. + tmax : float + The end of time interval in seconds. + + Returns + ------- + coefs : array + Sets of coefficients for each channel + """ + start, stop = raw.time_to_index(tmin, tmax) + data, times = raw[picks, start:(stop + 1)] + + coefs = np.empty((len(data), order)) + for k, d in enumerate(data): + this_coefs, _ = yule_walker(d, order=order) + coefs[k, :] = this_coefs + return coefs diff --git a/mne/time_frequency/tests/test_ar.py b/mne/time_frequency/tests/test_ar.py new file mode 100644 index 0000000..dfd6819 --- /dev/null +++ b/mne/time_frequency/tests/test_ar.py @@ -0,0 +1,42 @@ +import os.path as op +import numpy as np +from numpy.testing import assert_array_almost_equal +from nose.tools import assert_true +import nose + +from ... import fiff +from .. import yule_walker, ar_raw + +raw_fname = op.join(op.dirname(__file__), '..', '..', 'fiff', 'tests', 'data', + 'test_raw.fif') + + +def test_yule_walker(): + """Test Yule-Walker against statsmodels + """ + try: + from statsmodels.regression.linear_model import yule_walker as sm_yw + d = np.random.randn(100) + sm_rho, sm_sigma = sm_yw(d, order=2) + rho, sigma = yule_walker(d, order=2) + assert_array_almost_equal(sm_sigma, sigma) + assert_array_almost_equal(sm_rho, rho) + except ImportError: + raise nose.SkipTest("XFailed Test") + + +def test_ar_raw(): + raw = fiff.Raw(raw_fname) + + # picks MEG gradiometers + picks = fiff.pick_types(raw.info, meg='grad') + + picks = picks[:2] + + tmin, tmax = 0, 10 # use the first s of data + order = 2 + coefs = ar_raw(raw, picks=picks, order=order, tmin=tmin, tmax=tmax) + mean_coefs = np.mean(coefs, axis=0) + + assert_true(coefs.shape == (len(picks), order)) + assert_true(0.9 < mean_coefs[0] < 1.1) -- 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
