ENH: signal: add the 'nyq' keyword argument to firwin(), to be consis…

…tent with the API of firwin2().

scipy/signal/fir_filter_design.py

@@ -129,7 +129,8 @@ def kaiserord(ripple, width):
     return int(ceil(numtaps)), beta
 
 
-def firwin(numtaps, cutoff, width=None, window='hamming', pass_zero=True, scale=True):
+def firwin(numtaps, cutoff, width=None, window='hamming', pass_zero=True,
+                                                        scale=True, nyq=1.0):
     """
     FIR filter design using the window method.
     
@@ -145,25 +146,25 @@ def firwin(numtaps, cutoff, width=None, window='hamming', pass_zero=True, scale=
     ----------
     numtaps : int
         Length of the filter (number of coefficients, i.e. the filter
-        order + 1).  `numtaps` must be even if a passband includes the Nyquist
-        frequency.
+        order + 1).  `numtaps` 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`.
+        Cutoff frequency of filter (expressed in the same units as `nyq`)
+        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 `nyq`.  The values 0 and
+        `nyq` 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)
+        If `width` is not None, then assume it is the approximate width
+        of the transition region (expressed in the same units as `nyq`)
         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.
+        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.
@@ -175,10 +176,14 @@ def firwin(numtaps, cutoff, width=None, window='hamming', pass_zero=True, scale=
         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);
+            `nyq` (the Nyquist rate) if the first passband ends at
+                `nyq` (i.e the filter is a single band highpass filter);
             center of first passband otherwise.
 
+    nyq : float
+        Nyquist frequency.  Each frequency in `cutoff` must be between 0
+        and `nyq`.
+
     Returns
     -------
     h : 1D ndarray
@@ -187,8 +192,8 @@ def firwin(numtaps, cutoff, width=None, window='hamming', pass_zero=True, scale=
     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
+        If any value in `cutoff` is less than or equal to 0 or greater
+        than or equal to `nyq`, if the values in `cutoff` are not strictly
         monotonically increasing, or if `numtaps` is even but a passband
         includes the Nyquist frequency.
 
@@ -232,28 +237,24 @@ def firwin(numtaps, cutoff, width=None, window='hamming', pass_zero=True, scale=
     # The major enhancements to this function added in November 2010 were
     # developed by Tom Krauss (see ticket #902).
 
-    cutoff = np.atleast_1d(cutoff)
-    
+    cutoff = np.atleast_1d(cutoff) / float(nyq)
+
     # 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.")
+        raise ValueError("Invalid cutoff frequency: frequencies must be greater than 0 and less than nyq.")
     if np.any(np.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(numtaps, width)
-            beta = kaiser_beta(atten)
-            window = ('kaiser', beta)
-        else:
-            raise ValueError("Invalid value for width: %s", width)
+        # A width was given.  Find the beta parameter of the Kaiser window
+        # and set `window`.  This overrides the value of `window` passed in.
+        atten = kaiser_atten(numtaps, float(width)/nyq)
+        beta = kaiser_beta(atten)
+        window = ('kaiser', beta)
 
     pass_nyquist = bool(cutoff.size & 1) ^ pass_zero
     if pass_nyquist and numtaps % 2 == 0:

scipy/signal/tests/test_fir_filter_design.py

@@ -165,6 +165,25 @@ def test_multi(self):
                 [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_nyq(self):
+        """Test the nyq keyword."""
+        nyquist = 1000
+        width = 40.0
+        relative_width = width/nyquist
+        ntaps, beta = kaiserord(120, relative_width)
+        taps = firwin(ntaps, cutoff=[300, 700], window=('kaiser', beta),
+                        pass_zero=False, scale=False, nyq=nyquist)
+
+        # 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, 200, 300-width/2, 300+width/2, 500,
+                                700-width/2, 700+width/2, 800, 1000])
+        freqs, response = freqz(taps, worN=np.pi*freq_samples/nyquist)
+        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_bad_cutoff(self):
         """Test that invalid cutoff argument raises ValueError."""
         # cutoff values must be greater than 0 and less than 1.
@@ -180,7 +199,9 @@ def test_bad_cutoff(self):
         assert_raises(ValueError, firwin, 99, [])
         # 2D array not allowed.
         assert_raises(ValueError, firwin, 99, [[0.1, 0.2],[0.3, 0.4]])               
-
+        # cutoff values must be less than nyq.
+        assert_raises(ValueError, firwin, 99, 50.0, nyq=40)
+        assert_raises(ValueError, firwin, 99, [10, 20, 30], nyq=25)
 
     def test_even_highpass_raises_value_error(self):
         """Test that attempt to create a highpass filter with an even number
@@ -190,6 +211,7 @@ def test_even_highpass_raises_value_error(self):
 
 
 
+
 class TestFirwin2(TestCase):
 
     def test_invalid_args(self):