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 a11e228fb496a9ce6b5f0338e5d3ca6c4a43175b Author: Alexandre Gramfort <[email protected]> Date: Wed Jun 20 11:45:39 2012 +0300 FIX : Backporting scipy.signal.firwin2 --- doc/source/whats_new.rst | 2 + mne/filter.py | 6 +- mne/utils.py | 159 ++++++++++++++++++++++++++++++++++++++++++++++- 3 files changed, 163 insertions(+), 4 deletions(-) diff --git a/doc/source/whats_new.rst b/doc/source/whats_new.rst index 5913ccc..0985f76 100644 --- a/doc/source/whats_new.rst +++ b/doc/source/whats_new.rst @@ -7,6 +7,8 @@ Current Changelog ~~~~~~~~~ + - Backporting scipy.signa.firwin2 so filtering works with old scipy by `Alex Gramfort`_. + - Add Raw.filter method to more easily band pass data by `Alex Gramfort`_. - Add tmin + tmax parameters in mne.compute_covariance to estimate noise covariance in epochs baseline without creating new epochs by `Alex Gramfort`_. diff --git a/mne/filter.py b/mne/filter.py index 074c05e..99eb8cc 100644 --- a/mne/filter.py +++ b/mne/filter.py @@ -1,8 +1,8 @@ import warnings import numpy as np -from scipy import signal from scipy.fftpack import fft, ifft +from .utils import firwin2 # back port for old scipy def is_power2(num): """Test if number is a power of 2 @@ -178,7 +178,7 @@ def _filter(x, Fs, freq, gain, filter_length=None): N = len(x) - H = signal.firwin2(N, freq, gain) + H = firwin2(N, freq, gain) # Make zero-phase filter function B = np.abs(fft(H)) @@ -195,7 +195,7 @@ def _filter(x, Fs, freq, gain, filter_length=None): # Gain at Nyquist freq: 1: make N EVEN, 0: make N ODD N += 1 - H = signal.firwin2(N, freq, gain) + H = firwin2(N, freq, gain) xf = _overlap_add_filter(x, H, zero_phase=True) return xf diff --git a/mne/utils.py b/mne/utils.py index 31f826a..20f3832 100644 --- a/mne/utils.py +++ b/mne/utils.py @@ -5,7 +5,9 @@ # License: BSD (3-clause) import warnings - +import numpy as np +from math import ceil, log +from numpy.fft import irfft # Following deprecated class copied from scikit-learn @@ -90,3 +92,158 @@ class deprecated(object): if olddoc: newdoc = "%s\n\n%s" % (newdoc, olddoc) return newdoc + + +############################################################################### +# Back porting firwin2 for older scipy + +# Original version of firwin2 from scipy ticket #457, submitted by "tash". +# +# Rewritten by Warren Weckesser, 2010. + +def _firwin2(numtaps, freq, gain, nfreqs=None, window='hamming', nyq=1.0): + """FIR filter design using the window method. + + From the given frequencies `freq` and corresponding gains `gain`, + this function constructs an FIR filter with linear phase and + (approximately) the given frequency response. + + Parameters + ---------- + numtaps : int + The number of taps in the FIR filter. `numtaps` must be less than + `nfreqs`. If the gain at the Nyquist rate, `gain[-1]`, is not 0, + then `numtaps` must be odd. + + freq : array-like, 1D + The frequency sampling points. Typically 0.0 to 1.0 with 1.0 being + Nyquist. The Nyquist frequency can be redefined with the argument + `nyq`. + + The values in `freq` must be nondecreasing. A value can be repeated + once to implement a discontinuity. The first value in `freq` must + be 0, and the last value must be `nyq`. + + gain : array-like + The filter gains at the frequency sampling points. + + nfreqs : int, optional + The size of the interpolation mesh used to construct the filter. + For most efficient behavior, this should be a power of 2 plus 1 + (e.g, 129, 257, etc). The default is one more than the smallest + power of 2 that is not less than `numtaps`. `nfreqs` must be greater + than `numtaps`. + + window : string or (string, float) or float, or None, optional + Window function to use. Default is "hamming". See + `scipy.signal.get_window` for the complete list of possible values. + If None, no window function is applied. + + nyq : float + Nyquist frequency. Each frequency in `freq` must be between 0 and + `nyq` (inclusive). + + Returns + ------- + taps : numpy 1D array of length `numtaps` + The filter coefficients of the FIR filter. + + Examples + -------- + A lowpass FIR filter with a response that is 1 on [0.0, 0.5], and + that decreases linearly on [0.5, 1.0] from 1 to 0: + + >>> taps = firwin2(150, [0.0, 0.5, 1.0], [1.0, 1.0, 0.0]) + >>> print(taps[72:78]) + [-0.02286961 -0.06362756 0.57310236 0.57310236 -0.06362756 -0.02286961] + + See also + -------- + scipy.signal.firwin + + Notes + ----- + + From the given set of frequencies and gains, the desired response is + constructed in the frequency domain. The inverse FFT is applied to the + desired response to create the associated convolution kernel, and the + first `numtaps` coefficients of this kernel, scaled by `window`, are + returned. + + The FIR filter will have linear phase. The filter is Type I if `numtaps` + is odd and Type II if `numtaps` is even. Because Type II filters always + have a zero at the Nyquist frequency, `numtaps` must be odd if `gain[-1]` + is not zero. + + .. versionadded:: 0.9.0 + + References + ---------- + .. [1] Oppenheim, A. V. and Schafer, R. W., "Discrete-Time Signal + Processing", Prentice-Hall, Englewood Cliffs, New Jersey (1989). + (See, for example, Section 7.4.) + + .. [2] Smith, Steven W., "The Scientist and Engineer's Guide to Digital + Signal Processing", Ch. 17. http://www.dspguide.com/ch17/1.htm + + """ + + if len(freq) != len(gain): + raise ValueError('freq and gain must be of same length.') + + if nfreqs is not None and numtaps >= nfreqs: + raise ValueError('ntaps must be less than nfreqs, but firwin2 was ' + 'called with ntaps=%d and nfreqs=%s' % (numtaps, nfreqs)) + + if freq[0] != 0 or freq[-1] != nyq: + raise ValueError('freq must start with 0 and end with `nyq`.') + d = np.diff(freq) + if (d < 0).any(): + raise ValueError('The values in freq must be nondecreasing.') + d2 = d[:-1] + d[1:] + if (d2 == 0).any(): + raise ValueError('A value in freq must not occur more than twice.') + + if numtaps % 2 == 0 and gain[-1] != 0.0: + raise ValueError("A filter with an even number of coefficients must " + "have zero gain at the Nyquist rate.") + + if nfreqs is None: + nfreqs = 1 + 2 ** int(ceil(log(numtaps,2))) + + # Tweak any repeated values in freq so that interp works. + eps = np.finfo(float).eps + for k in range(len(freq)): + if k < len(freq)-1 and freq[k] == freq[k+1]: + freq[k] = freq[k] - eps + freq[k+1] = freq[k+1] + eps + + # Linearly interpolate the desired response on a uniform mesh `x`. + x = np.linspace(0.0, nyq, nfreqs) + fx = np.interp(x, freq, gain) + + # Adjust the phases of the coefficients so that the first `ntaps` of the + # inverse FFT are the desired filter coefficients. + shift = np.exp(-(numtaps-1)/2. * 1.j * np.pi * x / nyq) + fx2 = fx * shift + + # Use irfft to compute the inverse FFT. + out_full = irfft(fx2) + + if window is not None: + # Create the window to apply to the filter coefficients. + from scipy.signal.signaltools import get_window + wind = get_window(window, numtaps, fftbins=False) + else: + wind = 1 + + # Keep only the first `numtaps` coefficients in `out`, and multiply by + # the window. + out = out_full[:numtaps] * wind + + return out + +try: + from scipy.signal import firwin2 +except ImportError: + firwin2 = _firwin2 -- 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
