[med-svn] [python-mne] 212/353: FIX : Backporting scipy.signal.firwin2
Yaroslav Halchenko
debian at onerussian.com
Fri Nov 27 17:24:59 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 a11e228fb496a9ce6b5f0338e5d3ca6c4a43175b
Author: Alexandre Gramfort <alexandre.gramfort at inria.fr>
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
More information about the debian-med-commit
mailing list