[Scipy-svn] r6973 - in trunk: doc/release doc/source scipy/signal scipy/signal/tests

scipy-svn@scip... scipy-svn@scip...
Sun Nov 28 17:26:36 CST 2010


Author: warren.weckesser
Date: 2010-11-28 17:26:36 -0600 (Sun, 28 Nov 2010)
New Revision: 6973

Modified:
   trunk/doc/release/0.9.0-notes.rst
   trunk/doc/source/signal.rst
   trunk/scipy/signal/fir_filter_design.py
   trunk/scipy/signal/info.py
   trunk/scipy/signal/tests/test_fir_filter_design.py
Log:
ENH: signal: added the function firwin2 for constructing FIR filters with arbitrary frequency response.

Modified: trunk/doc/release/0.9.0-notes.rst
===================================================================
--- trunk/doc/release/0.9.0-notes.rst	2010-11-28 23:26:29 UTC (rev 6972)
+++ trunk/doc/release/0.9.0-notes.rst	2010-11-28 23:26:36 UTC (rev 6973)
@@ -114,6 +114,10 @@
 The function ``scipy.signal.firwin`` was enhanced to allow the
 design of highpass, bandpass, bandstop and multi-band FIR filters.
 
+The function ``scipy.signal.firwin2`` was added.  This function
+uses the window method to create a linear phase FIR filter with
+an arbitrary frequency response.
+
 The functions ``scipy.signal.kaiser_atten`` and ``scipy.signal.kaiser_beta``
 were added. 
 
@@ -131,7 +135,7 @@
 Removed features
 ================
 
-The deprecated modules ``helpdmod``, ``pexec`` and ``ppimport`` were removed
+The deprecated modules ``helpmod``, ``pexec`` and ``ppimport`` were removed
 from ``scipy.misc``.
 
 The ``save`` method of the ``spmatrix`` class in ``scipy.sparse``, which has

Modified: trunk/doc/source/signal.rst
===================================================================
--- trunk/doc/source/signal.rst	2010-11-28 23:26:29 UTC (rev 6972)
+++ trunk/doc/source/signal.rst	2010-11-28 23:26:36 UTC (rev 6973)
@@ -64,6 +64,7 @@
 
    bilinear
    firwin
+   firwin2
    freqs
    freqz
    iirdesign

Modified: trunk/scipy/signal/fir_filter_design.py
===================================================================
--- trunk/scipy/signal/fir_filter_design.py	2010-11-28 23:26:29 UTC (rev 6972)
+++ trunk/scipy/signal/fir_filter_design.py	2010-11-28 23:26:36 UTC (rev 6973)
@@ -1,7 +1,8 @@
 """Functions for FIR filter design."""
 
-import numpy
-from numpy import pi, ceil, cos, asarray
+from math import ceil, log
+import numpy as np
+from numpy.fft import irfft
 from scipy.special import sinc
 import sigtools
 
@@ -47,7 +48,7 @@
     return beta
 
 
-def kaiser_atten(N, width):
+def kaiser_atten(numtaps, width):
     """Compute the attenuation of a Kaiser FIR filter.
     
     Given the number of taps `N` and the transition width `width`, compute the
@@ -72,7 +73,7 @@
     --------
     kaiserord, kaiser_beta
     """
-    a = 2.285 * (N - 1) * pi * width + 7.95
+    a = 2.285 * (numtaps - 1) * np.pi * width + 7.95
     return a
 
 
@@ -90,7 +91,7 @@
 
     Returns
     -------
-    N : int
+    numtaps : int
         The length of the kaiser window.
     beta :
         The beta parameter for the kaiser window.
@@ -99,9 +100,9 @@
     -----
     There are several ways to obtain the Kaiser window:
 
-      signal.kaiser(N, beta, sym=0)
-      signal.get_window(beta, N)
-      signal.get_window(('kaiser', beta), N)
+      signal.kaiser(numtaps, beta, sym=0)
+      signal.get_window(beta, numtaps)
+      signal.get_window(('kaiser', beta), numtaps)
 
     The empirical equations discovered by Kaiser are used.
 
@@ -123,28 +124,28 @@
 
     # Kaiser's formula (as given in Oppenheim and Schafer) is for the filter
     # order, so we have to add 1 to get the number of taps.  
-    N = (A - 7.95) / 2.285 / (pi * width) + 1
+    numtaps = (A - 7.95) / 2.285 / (np.pi * width) + 1
 
-    return int(ceil(N)), beta
+    return int(ceil(numtaps)), beta
 
 
-def firwin(N, cutoff, width=None, window='hamming', pass_zero=True, scale=True):
+def firwin(numtaps, cutoff, width=None, window='hamming', pass_zero=True, scale=True):
     """
     FIR filter design using the window method.
     
     This function computes the coefficients of a finite impulse response filter.
-    The filter will have linear phase; it will be Type I if `N` is odd and
-    Type II if `N` is even.
+    The filter will have linear phase; it will be Type I if `numtaps` is odd and
+    Type II if `numtaps` is even.
     
     Type II filters always have zero response at the Nyquist rate, so a
-    ValueError exception is raised if firwin is called with `N` even and
+    ValueError exception is raised if firwin is called with `numtaps` even and
     having a passband whose right end is at the Nyquist rate.
 
     Parameters
     ----------
-    N : int
+    numtaps : int
         Length of the filter (number of coefficients, i.e. the filter
-        order + 1).  `N` must be even if a passband includes the Nyquist
+        order + 1).  `numtaps` must be even if a passband includes the Nyquist
         frequency.
 
     cutoff : float or 1D array_like
@@ -181,14 +182,14 @@
     Returns
     -------
     h : 1D ndarray
-        Coefficients of length N FIR filter.
+        Coefficients of length `numtaps` FIR filter.
 
     Raises
     ------
     ValueError
         If any value in cutoff is less than or equal to 0 or greater
         than or equal to 1, if the values in cutoff are not strictly
-        monotonically increasing, or if `N` is even but a passband
+        monotonically increasing, or if `numtaps` is even but a passband
         includes the Nyquist frequency.
 
     Examples
@@ -196,38 +197,42 @@
     
     Low-pass from 0 to f::
     
-    >>> firwin(N, f)
+    >>> firwin(numtaps, f)
     
     Use a specific window function::
     
-    >>> firwin(N, f, window='nuttall')
+    >>> firwin(numtaps, f, window='nuttall')
     
     High-pass ('stop' from 0 to f)::
      
-    >>> firwin(N, f, pass_zero=False)
+    >>> firwin(numtaps, f, pass_zero=False)
 
     Band-pass::
     
-    >>> firwin(N, [f1, f2], pass_zero=False)
+    >>> firwin(numtaps, [f1, f2], pass_zero=False)
     
     Band-stop::
     
-    >>> firwin(N, [f1, f2]) 
+    >>> firwin(numtaps, [f1, f2]) 
 
     Multi-band (passbands are [0, f1], [f2, f3] and [f4, 1])::
 
-    >>>firwin(N, [f1, f2, f3, f4])
+    >>>firwin(numtaps, [f1, f2, f3, f4])
     
     Multi-band (passbands are [f1, f2] and [f3,f4])::
     
-    >>> firwin(N, [f1, f2, f3, f4], pass_zero=False)
+    >>> firwin(numtaps, [f1, f2, f3, f4], pass_zero=False)
 
+    See also
+    --------
+    scipy.signal.firwin2
+
     """
 
     # The major enhancements to this function added in November 2010 were
     # developed by Tom Krauss (see ticket #902).
 
-    cutoff = numpy.atleast_1d(cutoff)
+    cutoff = np.atleast_1d(cutoff)
     
     # Check for invalid input.
     if cutoff.ndim > 1:
@@ -236,7 +241,7 @@
         raise ValueError("At least one cutoff frequency must be given.")
     if cutoff.min() <= 0 or cutoff.max() >= 1:
         raise ValueError("Invalid cutoff frequency: frequencies must be greater than 0 and less than 1.")
-    if numpy.any(numpy.diff(cutoff) <= 0):
+    if np.any(np.diff(cutoff) <= 0):
         raise ValueError("Invalid cutoff frequencies: the frequencies must be strictly increasing.")
 
     if width is not None:
@@ -244,27 +249,27 @@
         # of the Kaiser window, and set `window`.  This overrides the value of
         # `window` passed in.
         if isinstance(width, float):
-            atten = kaiser_atten(N, width)
+            atten = kaiser_atten(numtaps, width)
             beta = kaiser_beta(atten)
             window = ('kaiser', beta)
         else:
             raise ValueError("Invalid value for width: %s", width)
 
     pass_nyquist = bool(cutoff.size & 1) ^ pass_zero
-    if pass_nyquist and N % 2 == 0:
+    if pass_nyquist and numtaps % 2 == 0:
         raise ValueError("A filter with an even number of coefficients must "
                             "have zero response at the Nyquist rate.")
 
     # Insert 0 and/or 1 at the ends of cutoff so that the length of cutoff is even,
     # and each pair in cutoff corresponds to passband.
-    cutoff = numpy.hstack(([0.0]*pass_zero, cutoff, [1.0]*pass_nyquist))
+    cutoff = np.hstack(([0.0]*pass_zero, cutoff, [1.0]*pass_nyquist))
 
     # `bands` is a 2D array; each row gives the left and right edges of a passband.
     bands = cutoff.reshape(-1,2)
 
     # Build up the coefficients.
-    alpha = 0.5 * (N-1)
-    m = numpy.arange(0, N) - alpha
+    alpha = 0.5 * (numtaps-1)
+    m = np.arange(0, numtaps) - alpha
     h = 0
     for left, right in bands:
         h += right * sinc(right * m)
@@ -272,7 +277,7 @@
 
     # Get and apply the window function.
     from signaltools import get_window
-    win = get_window(window, N, fftbins=False)
+    win = get_window(window, numtaps, fftbins=False)
     h *= win
         
     # Now handle scaling if desired.
@@ -285,13 +290,160 @@
             scale_frequency = 1.0
         else:
             scale_frequency = 0.5 * (left + right)
-        c = cos(pi * m * scale_frequency)
-        s = numpy.sum(h * c)
+        c = np.cos(np.pi * m * scale_frequency)
+        s = np.sum(h * c)
         h /= s
  
     return h
 
 
+# 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.
+
+    Example
+    -------
+    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 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
+
+
 def remez(numtaps, bands, desired, weight=None, Hz=1, type='bandpass',
           maxiter=25, grid_density=16):
     """
@@ -385,6 +537,6 @@
     if weight is None:
         weight = [1] * len(desired)
 
-    bands = asarray(bands).copy()
+    bands = np.asarray(bands).copy()
     return sigtools._remez(numtaps, bands, desired, weight, tnum, Hz,
                            maxiter, grid_density)

Modified: trunk/scipy/signal/info.py
===================================================================
--- trunk/scipy/signal/info.py	2010-11-28 23:26:29 UTC (rev 6972)
+++ trunk/scipy/signal/info.py	2010-11-28 23:26:36 UTC (rev 6973)
@@ -71,7 +71,9 @@
     bilinear:
         Return a digital filter from an analog filter using the bilinear transform.
     firwin:
-        Windowed FIR filter design.
+        Windowed FIR filter design, with frequency response defined as pass and stop bands.
+    firwin2:
+        Windowed FIR filter design, with arbitrary frequency response.
     freqs:
         Analog filter frequency response.
     freqz:

Modified: trunk/scipy/signal/tests/test_fir_filter_design.py
===================================================================
--- trunk/scipy/signal/tests/test_fir_filter_design.py	2010-11-28 23:26:29 UTC (rev 6972)
+++ trunk/scipy/signal/tests/test_fir_filter_design.py	2010-11-28 23:26:36 UTC (rev 6973)
@@ -3,7 +3,7 @@
 from numpy.testing import TestCase, run_module_suite, assert_raises, \
         assert_array_almost_equal, assert_
 
-from scipy.signal import firwin, kaiserord, freqz, remez
+from scipy.signal import firwin, firwin2, kaiserord, freqz, remez
 
 
 class TestFirwin(TestCase):
@@ -189,6 +189,71 @@
         assert_raises(ValueError, firwin, 40, [.25, 0.5])
 
 
+
+class TestFirwin2(TestCase):
+
+    def test_invalid_args(self):
+        # `freq` and `gain` have different lengths.
+        assert_raises(ValueError, firwin2, 50, [0, 0.5, 1], [0.0, 1.0])
+        # `nfreqs` is less than `ntaps`.
+        assert_raises(ValueError, firwin2, 50, [0, 0.5, 1], [0.0, 1.0, 1.0], nfreqs=33)
+        # Decreasing value in `freq`
+        assert_raises(ValueError, firwin2, 50, [0, 0.5, 0.4, 1.0], [0, .25, .5, 1.0])
+        # Value in `freq` repeated more than once.
+        assert_raises(ValueError, firwin2, 50, [  0,   .1, .1,   .1, 1.0],
+                                             [0.0, 0.5, 0.75, 1.0, 1.0])
+        # `freq` does not start at 0.0.
+        assert_raises(ValueError, firwin2, 50, [0.5, 1.0], [0.0, 1.0])
+
+    def test01(self):
+        width = 0.04
+        beta = 12.0
+        ntaps = 400
+        # Filter is 1 from w=0 to w=0.5, then decreases linearly from 1 to 0 as w
+        # increases from w=0.5 to w=1  (w=1 is the Nyquist frequency). 
+        freq = [0.0, 0.5, 1.0]
+        gain = [1.0, 1.0, 0.0] 
+        taps = firwin2(ntaps, freq, gain, window=('kaiser', beta))
+        freq_samples = np.array([0.0, 0.25, 0.5-width/2, 0.5+width/2,
+                                                        0.75, 1.0-width/2])
+        freqs, response = freqz(taps, worN=np.pi*freq_samples)
+        assert_array_almost_equal(np.abs(response),
+                        [1.0, 1.0, 1.0, 1.0-width, 0.5, width], decimal=5)
+
+    def test02(self):
+        width = 0.04
+        beta = 12.0
+        # ntaps must be odd for positive gain at Nyquist.
+        ntaps = 401
+        # An ideal highpass filter.
+        freq = [0.0, 0.5, 0.5, 1.0]
+        gain = [0.0, 0.0, 1.0, 1.0]
+        taps = firwin2(ntaps, freq, gain, window=('kaiser', beta))
+        freq_samples = np.array([0.0, 0.25, 0.5-width, 0.5+width, 0.75, 1.0])
+        freqs, response = freqz(taps, worN=np.pi*freq_samples)
+        assert_array_almost_equal(np.abs(response),
+                                [0.0, 0.0, 0.0, 1.0, 1.0, 1.0], decimal=5)
+
+    def test03(self):
+        width = 0.02
+        ntaps, beta = kaiserord(120, width)
+        # ntaps must be odd for positive gain at Nyquist.
+        ntaps = int(ntaps) | 1
+        freq = [0.0, 0.4, 0.4, 0.5, 0.5, 1.0]
+        gain = [1.0, 1.0, 0.0, 0.0, 1.0, 1.0]
+        taps = firwin2(ntaps, freq, gain, window=('kaiser', beta))
+        freq_samples = np.array([0.0, 0.4-width, 0.4+width, 0.45, 
+                                    0.5-width, 0.5+width, 0.75, 1.0])
+        freqs, response = freqz(taps, worN=np.pi*freq_samples)
+        assert_array_almost_equal(np.abs(response),
+                    [1.0, 1.0, 0.0, 0.0, 0.0, 1.0, 1.0, 1.0], decimal=5)
+
+    def test_nyq(self):
+        taps1 = firwin2(80, [0.0, 0.5, 1.0], [1.0, 1.0, 0.0])
+        taps2 = firwin2(80, [0.0, 30.0, 60.0], [1.0, 1.0, 0.0], nyq=60.0)
+        assert_array_almost_equal(taps1, taps2)
+
+
 class TestRemez(TestCase):
 
     def test_hilbert(self):



More information about the Scipy-svn mailing list