# [Scipy-svn] r6873 - in trunk/scipy/signal: . tests

scipy-svn@scip... scipy-svn@scip...
Sat Nov 13 13:28:08 CST 2010

```Author: warren.weckesser
Date: 2010-11-13 13:28:08 -0600 (Sat, 13 Nov 2010)
New Revision: 6873

Modified:
trunk/scipy/signal/fir_filter_design.py
trunk/scipy/signal/signaltools.py
trunk/scipy/signal/tests/test_filter_design.py
trunk/scipy/signal/tests/test_fir_filter_design.py
Log:
REF: signal: move remez from signaltools.py to fir_filter_design.py

Modified: trunk/scipy/signal/fir_filter_design.py
===================================================================
--- trunk/scipy/signal/fir_filter_design.py	2010-11-13 17:19:02 UTC (rev 6872)
+++ trunk/scipy/signal/fir_filter_design.py	2010-11-13 19:28:08 UTC (rev 6873)
@@ -1,8 +1,9 @@
"""Functions for FIR filter design."""

import numpy
-from numpy import pi, ceil, cos
+from numpy import pi, ceil, cos, asarray
from scipy.special import sinc
+import sigtools

# Some notes on function parameters:
#
@@ -289,3 +290,101 @@
h /= s

return h
+
+
+def remez(numtaps, bands, desired, weight=None, Hz=1, type='bandpass',
+          maxiter=25, grid_density=16):
+    """
+    Calculate the minimax optimal filter using the Remez exchange algorithm.
+
+    Calculate the filter-coefficients for the finite impulse response
+    (FIR) filter whose transfer function minimizes the maximum error
+    between the desired gain and the realized gain in the specified
+    frequency bands using the Remez exchange algorithm.
+
+    Parameters
+    ----------
+    numtaps : int
+        The desired number of taps in the filter. The number of taps is
+        the number of terms in the filter, or the filter order plus one.
+    bands : array_like
+        A monotonic sequence containing the band edges in Hz.
+        All elements must be non-negative and less than half the sampling
+        frequency as given by `Hz`.
+    desired : array_like
+        A sequence half the size of bands containing the desired gain
+        in each of the specified bands.
+    weight : array_like, optional
+        A relative weighting to give to each band region. The length of
+        `weight` has to be half the length of `bands`.
+    Hz : scalar, optional
+        The sampling frequency in Hz. Default is 1.
+    type : {'bandpass', 'differentiator', 'hilbert'}, optional
+        The type of filter:
+
+          'bandpass' : flat response in bands. This is the default.
+
+          'differentiator' : frequency proportional response in bands.
+
+          'hilbert' : filter with odd symmetry, that is, type III
+                      (for even order) or type IV (for odd order)
+                      linear phase filters.
+
+    maxiter : int, optional
+        Maximum number of iterations of the algorithm. Default is 25.
+    grid_density : int, optional
+        Grid density. The dense grid used in `remez` is of size
+        ``(numtaps + 1) * grid_density``. Default is 16.
+
+    Returns
+    -------
+    out : ndarray
+        A rank-1 array containing the coefficients of the optimal
+        (in a minimax sense) filter.
+
+    --------
+    freqz : Compute the frequency response of a digital filter.
+
+    References
+    ----------
+    .. [1] J. H. McClellan and T. W. Parks, "A unified approach to the
+           design of optimum FIR linear phase digital filters",
+           IEEE Trans. Circuit Theory, vol. CT-20, pp. 697-701, 1973.
+    .. [2] J. H. McClellan, T. W. Parks and L. R. Rabiner, "A Computer
+           Program for Designing Optimum FIR Linear Phase Digital
+           Filters", IEEE Trans. Audio Electroacoust., vol. AU-21,
+           pp. 506-525, 1973.
+
+    Examples
+    --------
+    We want to construct a filter with a passband at 0.2-0.4 Hz, and
+    stop bands at 0-0.1 Hz and 0.45-0.5 Hz. Note that this means that the
+    behavior in the frequency ranges between those bands is unspecified and
+    may overshoot.
+
+    >>> bpass = sp.signal.remez(72, [0, 0.1, 0.2, 0.4, 0.45, 0.5], [0, 1, 0])
+    >>> freq, response = sp.signal.freqz(bpass)
+    >>> ampl = np.abs(response)
+
+    >>> import matplotlib.pyplot as plt
+    >>> fig = plt.figure()
+    >>> ax1.semilogy(freq/(2*np.pi), ampl, 'b-') # freq in Hz
+    [<matplotlib.lines.Line2D object at 0xf486790>]
+    >>> plt.show()
+
+    """
+    # Convert type
+    try:
+        tnum = {'bandpass':1, 'differentiator':2, 'hilbert':3}[type]
+    except KeyError:
+        raise ValueError("Type must be 'bandpass', 'differentiator', or 'hilbert'")
+
+    # Convert weight
+    if weight is None:
+        weight = [1] * len(desired)
+
+    bands = asarray(bands).copy()
+    return sigtools._remez(numtaps, bands, desired, weight, tnum, Hz,
+                           maxiter, grid_density)

Modified: trunk/scipy/signal/signaltools.py
===================================================================
--- trunk/scipy/signal/signaltools.py	2010-11-13 17:19:02 UTC (rev 6872)
+++ trunk/scipy/signal/signaltools.py	2010-11-13 19:28:08 UTC (rev 6873)
@@ -540,103 +540,8 @@

return sigtools._medfilt2d(image, kernel_size)

-def remez(numtaps, bands, desired, weight=None, Hz=1, type='bandpass',
-          maxiter=25, grid_density=16):
-    """
-    Calculate the minimax optimal filter using the Remez exchange algorithm.

-    Calculate the filter-coefficients for the finite impulse response
-    (FIR) filter whose transfer function minimizes the maximum error
-    between the desired gain and the realized gain in the specified
-    frequency bands using the Remez exchange algorithm.

-    Parameters
-    ----------
-    numtaps : int
-        The desired number of taps in the filter. The number of taps is
-        the number of terms in the filter, or the filter order plus one.
-    bands : array_like
-        A monotonic sequence containing the band edges in Hz.
-        All elements must be non-negative and less than half the sampling
-        frequency as given by `Hz`.
-    desired : array_like
-        A sequence half the size of bands containing the desired gain
-        in each of the specified bands.
-    weight : array_like, optional
-        A relative weighting to give to each band region. The length of
-        `weight` has to be half the length of `bands`.
-    Hz : scalar, optional
-        The sampling frequency in Hz. Default is 1.
-    type : {'bandpass', 'differentiator', 'hilbert'}, optional
-        The type of filter:
-
-          'bandpass' : flat response in bands. This is the default.
-
-          'differentiator' : frequency proportional response in bands.
-
-          'hilbert' : filter with odd symmetry, that is, type III
-                      (for even order) or type IV (for odd order)
-                      linear phase filters.
-
-    maxiter : int, optional
-        Maximum number of iterations of the algorithm. Default is 25.
-    grid_density : int, optional
-        Grid density. The dense grid used in `remez` is of size
-        ``(numtaps + 1) * grid_density``. Default is 16.
-
-    Returns
-    -------
-    out : ndarray
-        A rank-1 array containing the coefficients of the optimal
-        (in a minimax sense) filter.
-
-    --------
-    freqz : Compute the frequency response of a digital filter.
-
-    References
-    ----------
-    .. [1] J. H. McClellan and T. W. Parks, "A unified approach to the
-           design of optimum FIR linear phase digital filters",
-           IEEE Trans. Circuit Theory, vol. CT-20, pp. 697-701, 1973.
-    .. [2] J. H. McClellan, T. W. Parks and L. R. Rabiner, "A Computer
-           Program for Designing Optimum FIR Linear Phase Digital
-           Filters", IEEE Trans. Audio Electroacoust., vol. AU-21,
-           pp. 506-525, 1973.
-
-    Examples
-    --------
-    We want to construct a filter with a passband at 0.2-0.4 Hz, and
-    stop bands at 0-0.1 Hz and 0.45-0.5 Hz. Note that this means that the
-    behavior in the frequency ranges between those bands is unspecified and
-    may overshoot.
-
-    >>> bpass = sp.signal.remez(72, [0, 0.1, 0.2, 0.4, 0.45, 0.5], [0, 1, 0])
-    >>> freq, response = sp.signal.freqz(bpass)
-    >>> ampl = np.abs(response)
-
-    >>> import matplotlib.pyplot as plt
-    >>> fig = plt.figure()
-    >>> ax1.semilogy(freq/(2*np.pi), ampl, 'b-') # freq in Hz
-    [<matplotlib.lines.Line2D object at 0xf486790>]
-    >>> plt.show()
-
-    """
-    # Convert type
-    try:
-        tnum = {'bandpass':1, 'differentiator':2, 'hilbert':3}[type]
-    except KeyError:
-        raise ValueError("Type must be 'bandpass', 'differentiator', or 'hilbert'")
-
-    # Convert weight
-    if weight is None:
-        weight = [1] * len(desired)
-
-    bands = asarray(bands).copy()
-    return sigtools._remez(numtaps, bands, desired, weight, tnum, Hz,
-                           maxiter, grid_density)
-
def lfilter(b, a, x, axis=-1, zi=None):
"""
Filter data along one-dimension with an IIR or FIR filter.

Modified: trunk/scipy/signal/tests/test_filter_design.py
===================================================================
--- trunk/scipy/signal/tests/test_filter_design.py	2010-11-13 17:19:02 UTC (rev 6872)
+++ trunk/scipy/signal/tests/test_filter_design.py	2010-11-13 19:28:08 UTC (rev 6873)
@@ -1,9 +1,9 @@
import warnings

import numpy as np
-from numpy.testing import TestCase, assert_array_almost_equal, assert_
+from numpy.testing import TestCase, assert_array_almost_equal

-from scipy.signal import tf2zpk, bessel, BadCoefficients, freqz, remez
+from scipy.signal import tf2zpk, bessel, BadCoefficients

class TestTf2zpk(TestCase):
@@ -37,36 +37,3 @@
pass
finally:
-
-
-
-class TestRemez(TestCase):
-
-    def test_hilbert(self):
-        N = 11 # number of taps in the filter
-        a = 0.1 # width of the transition band
-
-        # design an unity gain hilbert bandpass filter from w to 0.5-w
-        h = remez(11, [ a, 0.5-a ], [ 1 ], type='hilbert')
-
-        # make sure the filter has correct # of taps
-        assert_(len(h) == N, "Number of Taps")
-
-        # make sure it is type III (anti-symmtric tap coefficients)
-        assert_array_almost_equal(h[:(N-1)/2], -h[:-(N-1)/2-1:-1])
-
-        # Since the requested response is symmetric, all even coeffcients
-        # should be zero (or in this case really small)
-        assert_((abs(h[1::2]) < 1e-15).all(), "Even Coefficients Equal Zero")
-
-        # now check the frequency response
-        w, H = freqz(h, 1)
-        f = w/2/np.pi
-        Hmag = abs(H)
-
-        # should have a zero at 0 and pi (in this case close to zero)
-        assert_((Hmag[ [0,-1] ] < 0.02).all(), "Zero at zero and pi")
-
-        # check that the pass band is close to unity
-        idx = (f > a) * (f < 0.5-a)
-        assert_((abs(Hmag[idx] - 1) < 0.015).all(), "Pass Band Close To Unity")

Modified: trunk/scipy/signal/tests/test_fir_filter_design.py
===================================================================
--- trunk/scipy/signal/tests/test_fir_filter_design.py	2010-11-13 17:19:02 UTC (rev 6872)
+++ trunk/scipy/signal/tests/test_fir_filter_design.py	2010-11-13 19:28:08 UTC (rev 6873)
@@ -1,9 +1,9 @@

import numpy as np
from numpy.testing import TestCase, run_module_suite, assert_raises, \
-        assert_array_almost_equal
+        assert_array_almost_equal, assert_

-from scipy.signal import firwin, kaiserord, freqz
+from scipy.signal import firwin, kaiserord, freqz, remez

class TestFirwin(TestCase):
@@ -189,5 +189,37 @@
assert_raises(ValueError, firwin, 40, [.25, 0.5])

+class TestRemez(TestCase):
+
+    def test_hilbert(self):
+        N = 11 # number of taps in the filter
+        a = 0.1 # width of the transition band
+
+        # design an unity gain hilbert bandpass filter from w to 0.5-w
+        h = remez(11, [ a, 0.5-a ], [ 1 ], type='hilbert')
+
+        # make sure the filter has correct # of taps
+        assert_(len(h) == N, "Number of Taps")
+
+        # make sure it is type III (anti-symmtric tap coefficients)
+        assert_array_almost_equal(h[:(N-1)/2], -h[:-(N-1)/2-1:-1])
+
+        # Since the requested response is symmetric, all even coeffcients
+        # should be zero (or in this case really small)
+        assert_((abs(h[1::2]) < 1e-15).all(), "Even Coefficients Equal Zero")
+
+        # now check the frequency response
+        w, H = freqz(h, 1)
+        f = w/2/np.pi
+        Hmag = abs(H)
+
+        # should have a zero at 0 and pi (in this case close to zero)
+        assert_((Hmag[ [0,-1] ] < 0.02).all(), "Zero at zero and pi")
+
+        # check that the pass band is close to unity
+        idx = (f > a) * (f < 0.5-a)
+        assert_((abs(Hmag[idx] - 1) < 0.015).all(), "Pass Band Close To Unity")
+
+
if __name__ == "__main__":
run_module_suite()

```