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

scipy-svn@scip... scipy-svn@scip...
Sat Nov 13 11:19:03 CST 2010


Author: warren.weckesser
Date: 2010-11-13 11:19:02 -0600 (Sat, 13 Nov 2010)
New Revision: 6872

Added:
   trunk/scipy/signal/tests/test_fir_filter_design.py
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_filter_design.py
Log:
ENH: signal: firwin() can now create highpass, bandpass, bandstop and multi-band FIR filters.

Modified: trunk/doc/release/0.9.0-notes.rst
===================================================================
--- trunk/doc/release/0.9.0-notes.rst	2010-11-12 00:57:30 UTC (rev 6871)
+++ trunk/doc/release/0.9.0-notes.rst	2010-11-13 17:19:02 UTC (rev 6872)
@@ -108,5 +108,15 @@
 equation systems (``scipy.linalg.solve_triangular``).
 
 
+Improved FIR filter design functions (``scipy.signal``)
+-------------------------------------------------------
+
+The function ``scipy.signal.firwin`` was enhanced to allow the
+design of highpass, bandpass, bandstop and multi-band FIR filters.
+
+The functions ``scipy.signal.kaiser_atten`` and ``scipy.signal.kaiser_beta``
+were added. 
+
+
 Removed features
 ================

Modified: trunk/doc/source/signal.rst
===================================================================
--- trunk/doc/source/signal.rst	2010-11-12 00:57:30 UTC (rev 6871)
+++ trunk/doc/source/signal.rst	2010-11-13 17:19:02 UTC (rev 6872)
@@ -68,6 +68,8 @@
    freqz
    iirdesign
    iirfilter
+   kaiser_atten
+   kaiser_beta
    kaiserord
    remez
 

Modified: trunk/scipy/signal/fir_filter_design.py
===================================================================
--- trunk/scipy/signal/fir_filter_design.py	2010-11-12 00:57:30 UTC (rev 6871)
+++ trunk/scipy/signal/fir_filter_design.py	2010-11-13 17:19:02 UTC (rev 6872)
@@ -1,10 +1,80 @@
 """Functions for FIR filter design."""
 
 import numpy
-from numpy import pi, ceil
-from scipy import special
+from numpy import pi, ceil, cos
+from scipy.special import sinc
 
+# Some notes on function parameters:
+#
+# `cutoff` and `width` are given as a numbers between 0 and 1.  These
+# are relative frequencies, expressed as a fraction of the Nyquist rate.
+# For example, if the Nyquist rate is 2KHz, then width=0.15 is a width
+# of 300 Hz.
+#
+# The `order` of a FIR filter is one less than the number of taps.
+# This is a potential source of confusion, so in the following code,
+# we will always use the number of taps as the parameterization of
+# the 'size' of the filter. The "number of taps" means the number
+# of coefficients, which is the same as the length of the impulse
+# response of the filter.
 
+
+def kaiser_beta(a):
+    """Compute the Kaiser parameter `beta`, given the attenuation `a`.
+    
+    Parameters
+    ----------
+    a : float
+        The desired attenuation in the stopband and maximum ripple in
+        the passband, in dB.  This should be a *positive* number.
+
+    Returns
+    -------
+    beta : float
+        The `beta` parameter to be used in the formula for a Kaiser window.
+    
+    References
+    ----------
+    Oppenheim, Schafer, "Discrete-Time Signal Processing", p.475-476.
+    """
+    if  a > 50:
+        beta = 0.1102 * (a - 8.7)
+    elif a > 21:
+        beta = 0.5842 * (a - 21)**0.4 + 0.07886 * (a - 21)
+    else:
+        beta = 0.0
+    return beta
+
+
+def kaiser_atten(N, width):
+    """Compute the attenuation of a Kaiser FIR filter.
+    
+    Given the number of taps `N` and the transition width `width`, compute the
+    attenuation `a` in dB, given by Kaiser's formula:
+
+        a = 2.285 * (N - 1) * pi * width + 7.95
+
+    Parameters
+    ----------
+    N : int
+        The number of taps in the FIR filter.
+    width : float
+        The desired width of the transition region between passband and stopband
+        (or, in general, at any discontinuity) for the filter.
+
+    Returns
+    -------
+    a : float
+        The attenuation of the ripple, in dB.
+        
+    See Also
+    --------
+    kaiserord, kaiser_beta
+    """
+    a = 2.285 * (N - 1) * pi * width + 7.95
+    return a
+
+
 def kaiserord(ripple, width):
     """Design a Kaiser window to limit ripple and width of transition region.
 
@@ -20,7 +90,7 @@
     Returns
     -------
     N : int
-        The order parameter for the kaiser window.
+        The length of the kaiser window.
     beta :
         The beta parameter for the kaiser window.
 
@@ -29,62 +99,193 @@
     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.get_window(beta, N)
+      signal.get_window(('kaiser', beta), N)
 
     The empirical equations discovered by Kaiser are used.
 
+    See Also
+    --------
+    kaiser_beta, kaiser_atten
+
     References
     ----------
     Oppenheim, Schafer, "Discrete-Time Signal Processing", p.475-476.
 
     """
     A = abs(ripple)  # in case somebody is confused as to what's meant
-    if (A>50):
-        beta = 0.1102*(A-8.7)
-    elif (A>21):
-        beta = 0.5842*(A-21)**0.4 + 0.07886*(A-21)
-    else:
-        beta = 0.0
-    N = (A-8)/2.285/(pi*width)
-    return ceil(N), beta
+    if A < 8:
+        # Formula for N is not valid in this range.
+        raise ValueError("Requested maximum ripple attentuation %f is too "
+                            "small for the Kaiser formula." % A)
+    beta = kaiser_beta(A)
 
-def firwin(N, cutoff, width=None, window='hamming'):
+    # 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
+
+    return int(ceil(N)), beta
+
+
+def firwin(N, cutoff, width=None, window='hamming', pass_zero=True, scale=True):
     """
-    FIR Filter Design using windowed ideal filter method.
+    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.
+    
+    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
+    having a passband whose right end is at the Nyquist rate.
 
     Parameters
     ----------
     N : int
-        Order of filter (number of taps).
-    cutoff : float
-        Cutoff frequency of filter (normalized so that 1 corresponds to Nyquist
-        or pi radians / sample)
-    width : float
-        If `width` is not None, then assume it is the approximate width of the
-        transition region (normalized so that 1 corresonds to pi) for use in
-        kaiser FIR filter design.
-    window : str. optional
-        Desired window to use. See `get_window` for a list of windows and
-        required parameters. Default is 'hamming'.
+        Length of the filter (number of coefficients, i.e. the filter
+        order + 1).  `N` must be even if a passband includes the Nyquist
+        frequency.
 
+    cutoff : float or 1D array_like
+        Cutoff frequency of filter (normalized so that 1 corresponds to
+        Nyquist or pi radians / sample) OR an array of cutoff frequencies
+        (that is, band edges). In the latter case, the frequencies in 
+        `cutoff` should be positive and monotonically increasing between
+        0 and 1.  The values 0 and 1 must not be included in `cutoff`.
+
+    width : float or None
+        If `width` is not None, then assume it is the approximate width of
+        the transition region (normalized so that 1 corresponds to pi)
+        for use in Kaiser FIR filter design.  In this case, the `window`
+        argument is ignored.
+
+    window : string or tuple of string and parameter values
+        Desired window to use. See `scipy.signal.get_window` for a list of
+        windows and required parameters.
+
+    pass_zero : bool
+        If True, the gain at the frequency 0 (i.e. the "DC gain") is 1.
+        Otherwise the DC gain is 0.
+
+    scale : bool
+        Set to True to scale the coefficients so that the frequency
+        response is exactly unity at a certain frequency.
+        That frequency is either:
+            0 (DC) if the first passband starts at 0 (i.e. pass_zero
+                is True);
+            1 (Nyquist) if the first passband ends at 1 (i.e the
+                filter is a single band highpass filter);
+            center of first passband otherwise.
+
     Returns
     -------
-    h : ndarray
+    h : 1D ndarray
         Coefficients of length N 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
+        includes the Nyquist frequency.
+
+    Examples
+    --------
+    
+    Low-pass from 0 to f::
+    
+    >>> firwin(N, f)
+    
+    Use a specific window function::
+    
+    >>> firwin(N, f, window='nuttall')
+    
+    High-pass ('stop' from 0 to f)::
+     
+    >>> firwin(N, f, pass_zero=False)
+
+    Band-pass::
+    
+    >>> firwin(N, [f1, f2], pass_zero=False)
+    
+    Band-stop::
+    
+    >>> firwin(N, [f1, f2]) 
+
+    Multi-band (passbands are [0, f1], [f2, f3] and [f4, 1])::
+
+    >>>firwin(N, [f1, f2, f3, f4])
+    
+    Multi-band (passbands are [f1, f2] and [f3,f4])::
+    
+    >>> firwin(N, [f1, f2, f3, f4], pass_zero=False)
+
     """
 
+    # The major enhancements to this function added in November 2010 were
+    # developed by Tom Krauss (see ticket #902).
+
+    cutoff = numpy.atleast_1d(cutoff)
+    
+    # Check for invalid input.
+    if cutoff.ndim > 1:
+        raise ValueError("The cutoff argument must be at most one-dimensional.")
+    if cutoff.size == 0:
+        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):
+        raise ValueError("Invalid cutoff frequencies: the frequencies must be strictly increasing.")
+
+    if width is not None:
+        # A width was given.  Verify that it is a float, find the beta parameter
+        # of the Kaiser window, and set `window`.  This overrides the value of
+        # `window` passed in.
+        if isinstance(width, float):
+            atten = kaiser_atten(N, 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:
+        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))
+
+    # `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
+    h = 0
+    for left, right in bands:
+        h += right * sinc(right * m)
+        h -= left * sinc(left * m)
+
+    # Get and apply the window function.
     from signaltools import get_window
-    if isinstance(width,float):
-        A = 2.285*N*width + 8
-        if (A < 21): beta = 0.0
-        elif (A <= 50): beta = 0.5842*(A-21)**0.4 + 0.07886*(A-21)
-        else: beta = 0.1102*(A-8.7)
-        window=('kaiser',beta)
-
-    win = get_window(window,N,fftbins=1)
-    alpha = N//2
-    m = numpy.arange(0,N)
-    h = win*special.sinc(cutoff*(m-alpha))
-    return h / numpy.sum(h,axis=0)
+    win = get_window(window, N, fftbins=False)
+    h *= win
+        
+    # Now handle scaling if desired.
+    if scale:
+        # Get the first passband.
+        left, right = bands[0]
+        if left == 0:
+            scale_frequency = 0.0
+        elif right == 1:
+            scale_frequency = 1.0
+        else:
+            scale_frequency = 0.5 * (left + right)
+        c = cos(pi * m * scale_frequency)
+        s = numpy.sum(h * c)
+        h /= s
+ 
+    return h

Modified: trunk/scipy/signal/info.py
===================================================================
--- trunk/scipy/signal/info.py	2010-11-12 00:57:30 UTC (rev 6871)
+++ trunk/scipy/signal/info.py	2010-11-13 17:19:02 UTC (rev 6872)
@@ -82,6 +82,11 @@
         IIR filter design given order and critical frequencies.
     invres:
         Inverse partial fraction expansion.
+    kaiser_beta:
+        Compute the Kaiser parameter beta, given the desired FIR filter attenuation.
+    kaiser_atten:
+        Compute the attenuation of a Kaiser FIR filter, given the number of taps
+        and the transition width at discontinuities in the frequency response.
     kaiserord:
         Design a Kaiser window to limit ripple and width of transition region.
     remez:

Modified: trunk/scipy/signal/tests/test_filter_design.py
===================================================================
--- trunk/scipy/signal/tests/test_filter_design.py	2010-11-12 00:57:30 UTC (rev 6871)
+++ trunk/scipy/signal/tests/test_filter_design.py	2010-11-13 17:19:02 UTC (rev 6872)
@@ -3,7 +3,7 @@
 import numpy as np
 from numpy.testing import TestCase, assert_array_almost_equal, assert_
 
-from scipy.signal import tf2zpk, bessel, BadCoefficients, kaiserord, firwin, freqz, remez
+from scipy.signal import tf2zpk, bessel, BadCoefficients, freqz, remez
 
 
 class TestTf2zpk(TestCase):
@@ -39,17 +39,7 @@
             warnings.simplefilter("always", BadCoefficients)
 
 
-class TestFirWin(TestCase):
 
-    def test_lowpass(self):
-        width = 0.04
-        ntaps, beta = kaiserord(120, width)
-        taps = firwin(ntaps, cutoff=0.5, window=('kaiser', beta))
-        freq_samples = np.array([0.0, 0.25, 0.5-width/2, 0.5+width/2, 0.75, 1.0])
-        freqs, response = freqz(taps, worN=np.pi*freq_samples)
-        assert_array_almost_equal(np.abs(response),
-                                    [1.0, 1.0, 1.0, 0.0, 0.0, 0.0], decimal=5)
-
 class TestRemez(TestCase):
 
     def test_hilbert(self):

Added: trunk/scipy/signal/tests/test_fir_filter_design.py
===================================================================
--- trunk/scipy/signal/tests/test_fir_filter_design.py	                        (rev 0)
+++ trunk/scipy/signal/tests/test_fir_filter_design.py	2010-11-13 17:19:02 UTC (rev 6872)
@@ -0,0 +1,193 @@
+
+import numpy as np
+from numpy.testing import TestCase, run_module_suite, assert_raises, \
+        assert_array_almost_equal
+
+from scipy.signal import firwin, kaiserord, freqz
+
+
+class TestFirwin(TestCase):
+    
+    def check_response(self, h, expected_response, tol=.05):
+        N = len(h)
+        alpha = 0.5 * (N-1)
+        m = np.arange(0,N) - alpha   # time indices of taps
+        for freq, expected in expected_response:
+            actual = abs(np.sum(h*np.exp(-1.j*np.pi*m*freq)))
+            mse = abs(actual-expected)**2
+            self.assertTrue(mse < tol, 'response not as expected, mse=%g > %g'\
+               %(mse, tol))
+    
+    def test_response(self):
+        N = 51
+        f = .5
+        # increase length just to try even/odd
+        h = firwin(N, f) # low-pass from 0 to f 
+        self.check_response(h, [(.25,1), (.75,0)])
+
+        h = firwin(N+1, f, window='nuttall') # specific window 
+        self.check_response(h, [(.25,1), (.75,0)])
+
+        h = firwin(N+2, f, pass_zero=False) # stop from 0 to f --> high-pass
+        self.check_response(h, [(.25,0), (.75,1)])
+        
+        f1, f2, f3, f4 = .2, .4, .6, .8
+        h = firwin(N+3, [f1, f2], pass_zero=False) # band-pass filter
+        self.check_response(h, [(.1,0), (.3,1), (.5,0)])
+
+        h = firwin(N+4, [f1, f2]) # band-stop filter 
+        self.check_response(h, [(.1,1), (.3,0), (.5,1)])
+        
+        h = firwin(N+5, [f1, f2, f3, f4], pass_zero=False, scale=False) 
+        self.check_response(h, [(.1,0), (.3,1), (.5,0), (.7,1), (.9,0)])
+
+        h = firwin(N+6, [f1, f2, f3, f4])  # multiband filter
+        self.check_response(h, [(.1,1), (.3,0), (.5,1), (.7,0), (.9,1)])
+
+        h = firwin(N+7, 0.1, width=.03) # low-pass 
+        self.check_response(h, [(.05,1), (.75,0)])
+
+        h = firwin(N+8, 0.1, pass_zero=False) # high-pass 
+        self.check_response(h, [(.05,0), (.75,1)])
+
+    def mse(self, h, bands):
+        """Compute mean squared error versus ideal response across frequency 
+        band.
+          h -- coefficients
+          bands -- list of (left, right) tuples relative to 1==Nyquist of 
+            passbands
+        """
+        w, H = freqz(h, worN=1024)
+        f = w/np.pi
+        passIndicator = np.zeros(len(w), bool)
+        for left, right in bands:
+            passIndicator |= (f>=left) & (f<right)
+        Hideal = np.where(passIndicator, 1, 0)
+        mse = np.mean(abs(abs(H)-Hideal)**2)
+        return mse
+        
+    def test_scaling(self):
+        """
+        For one lowpass, bandpass, and highpass example filter, this test 
+        checks two things:
+          - the mean squared error over the frequency domain of the unscaled
+            filter is smaller than the scaled filter (true for rectangular
+            window)
+          - the response of the scaled filter is exactly unity at the center
+            of the first passband
+        """
+        N = 11
+        cases = [
+            ([.5],      True,   (0, 1)),
+            ([0.2, .6], False,  (.4, 1)),
+            ([.5],      False,  (1, 1)),
+        ]
+        for cutoff, pass_zero, expected_response in cases:
+            h = firwin(N, cutoff, scale=False, pass_zero=pass_zero, window='ones')
+            hs = firwin(N, cutoff, scale=True, pass_zero=pass_zero, window='ones')
+            if len(cutoff) == 1:
+                if pass_zero:
+                    cutoff = [0] + cutoff
+                else:
+                    cutoff = cutoff + [1]
+            self.assertTrue(self.mse(h, [cutoff]) < self.mse(hs, [cutoff]), 
+                'least squares violation')
+            self.check_response(hs, [expected_response], 1e-12)
+
+
+class TestFirWinMore(TestCase):
+    """Different author, different style, different tests..."""
+
+    def test_lowpass(self):
+        width = 0.04
+        ntaps, beta = kaiserord(120, width)
+        taps = firwin(ntaps, cutoff=0.5, window=('kaiser', beta), scale=False)
+
+        # Check the symmetry of taps.
+        assert_array_almost_equal(taps[:ntaps/2], taps[ntaps:ntaps-ntaps/2-1:-1])
+
+        # Check the gain at a few samples where we know it should be approximately 0 or 1.
+        freq_samples = np.array([0.0, 0.25, 0.5-width/2, 0.5+width/2, 0.75, 1.0])
+        freqs, response = freqz(taps, worN=np.pi*freq_samples)
+        assert_array_almost_equal(np.abs(response),
+                                    [1.0, 1.0, 1.0, 0.0, 0.0, 0.0], decimal=5)
+
+    def test_highpass(self):
+        width = 0.04
+        ntaps, beta = kaiserord(120, width)
+
+        # Ensure that ntaps is odd.
+        ntaps |= 1
+
+        taps = firwin(ntaps, cutoff=0.5, window=('kaiser', beta),
+                        pass_zero=False, scale=False)
+
+        # Check the symmetry of taps.
+        assert_array_almost_equal(taps[:ntaps/2], taps[ntaps:ntaps-ntaps/2-1:-1])
+
+        # Check the gain at a few samples where we know it should be approximately 0 or 1.
+        freq_samples = np.array([0.0, 0.25, 0.5-width/2, 0.5+width/2, 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 test_bandpass(self):
+        width = 0.04
+        ntaps, beta = kaiserord(120, width)
+        taps = firwin(ntaps, cutoff=[0.3, 0.7], window=('kaiser', beta),
+                        pass_zero=False, scale=False)
+
+        # Check the symmetry of taps.
+        assert_array_almost_equal(taps[:ntaps/2], taps[ntaps:ntaps-ntaps/2-1:-1])
+
+        # Check the gain at a few samples where we know it should be approximately 0 or 1.
+        freq_samples = np.array([0.0, 0.2, 0.3-width/2, 0.3+width/2, 0.5,
+                                0.7-width/2, 0.7+width/2, 0.8, 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, 0.0, 0.0, 0.0], decimal=5)
+
+    def test_multi(self):
+        width = 0.04
+        ntaps, beta = kaiserord(120, width)
+        taps = firwin(ntaps, cutoff=[0.2, 0.5, 0.8], window=('kaiser', beta),
+                        pass_zero=True, scale=False)
+
+        # Check the symmetry of taps.
+        assert_array_almost_equal(taps[:ntaps/2], taps[ntaps:ntaps-ntaps/2-1:-1])
+
+        # Check the gain at a few samples where we know it should be approximately 0 or 1.
+        freq_samples = np.array([0.0, 0.1, 0.2-width/2, 0.2+width/2, 0.35,
+                                0.5-width/2, 0.5+width/2, 0.65,
+                                0.8-width/2, 0.8+width/2, 0.9, 1.0])
+        freqs, response = freqz(taps, worN=np.pi*freq_samples)
+        assert_array_almost_equal(np.abs(response),
+                [1.0, 1.0, 1.0, 0.0, 0.0, 0.0, 1.0, 1.0, 1.0, 0.0, 0.0, 0.0],
+                decimal=5)
+
+    def test_bad_cutoff(self):
+        """Test that invalid cutoff argument raises ValueError."""
+        # cutoff values must be greater than 0 and less than 1.
+        assert_raises(ValueError, firwin, 99, -0.5)
+        assert_raises(ValueError, firwin, 99, 1.5)
+        # Don't allow 0 or 1 in cutoff.
+        assert_raises(ValueError, firwin, 99, [0, 0.5])
+        assert_raises(ValueError, firwin, 99, [0.5, 1])
+        # cutoff values must be strictly increasing.
+        assert_raises(ValueError, firwin, 99, [0.1, 0.5, 0.2])
+        assert_raises(ValueError, firwin, 99, [0.1, 0.5, 0.5])
+        # Must have at least one cutoff value.
+        assert_raises(ValueError, firwin, 99, [])
+        # 2D array not allowed.
+        assert_raises(ValueError, firwin, 99, [[0.1, 0.2],[0.3, 0.4]])               
+
+
+    def test_even_highpass_raises_value_error(self):
+        """Test that attempt to create a highpass filter with an even number
+        of taps raises a ValueError exception."""
+        assert_raises(ValueError, firwin, 40, 0.5, pass_zero=False)
+        assert_raises(ValueError, firwin, 40, [.25, 0.5])
+
+
+if __name__ == "__main__":
+    run_module_suite()



More information about the Scipy-svn mailing list