[Scipy-svn] r6873 - in trunk/scipy/signal: . tests
scipy-svn at scipy.org
scipy-svn at scipy.org
Sat Nov 13 14:28:08 EST 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.
+
+ See Also
+ --------
+ 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 = fig.add_subplot(111)
+ >>> 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.
-
- See Also
- --------
- 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 = fig.add_subplot(111)
- >>> 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:
warnings.simplefilter("always", BadCoefficients)
-
-
-
-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()
More information about the Scipy-svn
mailing list