[Scipysvn] r6873  in trunk/scipy/signal: . tests
scipysvn@scip...
scipysvn@scip...
Sat Nov 13 13:28:08 CST 2010
Author: warren.weckesser
Date: 20101113 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 20101113 17:19:02 UTC (rev 6872)
+++ trunk/scipy/signal/fir_filter_design.py 20101113 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 filtercoefficients 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 nonnegative 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 rank1 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. CT20, pp. 697701, 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. AU21,
+ pp. 506525, 1973.
+
+ Examples
+ 
+ We want to construct a filter with a passband at 0.20.4 Hz, and
+ stop bands at 00.1 Hz and 0.450.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 20101113 17:19:02 UTC (rev 6872)
+++ trunk/scipy/signal/signaltools.py 20101113 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 filtercoefficients 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 nonnegative 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 rank1 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. CT20, pp. 697701, 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. AU21,
 pp. 506525, 1973.

 Examples
 
 We want to construct a filter with a passband at 0.20.4 Hz, and
 stop bands at 00.1 Hz and 0.450.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 onedimension with an IIR or FIR filter.
Modified: trunk/scipy/signal/tests/test_filter_design.py
===================================================================
 trunk/scipy/signal/tests/test_filter_design.py 20101113 17:19:02 UTC (rev 6872)
+++ trunk/scipy/signal/tests/test_filter_design.py 20101113 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.5w
 h = remez(11, [ a, 0.5a ], [ 1 ], type='hilbert')

 # make sure the filter has correct # of taps
 assert_(len(h) == N, "Number of Taps")

 # make sure it is type III (antisymmtric tap coefficients)
 assert_array_almost_equal(h[:(N1)/2], h[:(N1)/21:1])

 # Since the requested response is symmetric, all even coeffcients
 # should be zero (or in this case really small)
 assert_((abs(h[1::2]) < 1e15).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.5a)
 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 20101113 17:19:02 UTC (rev 6872)
+++ trunk/scipy/signal/tests/test_fir_filter_design.py 20101113 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.5w
+ h = remez(11, [ a, 0.5a ], [ 1 ], type='hilbert')
+
+ # make sure the filter has correct # of taps
+ assert_(len(h) == N, "Number of Taps")
+
+ # make sure it is type III (antisymmtric tap coefficients)
+ assert_array_almost_equal(h[:(N1)/2], h[:(N1)/21:1])
+
+ # Since the requested response is symmetric, all even coeffcients
+ # should be zero (or in this case really small)
+ assert_((abs(h[1::2]) < 1e15).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.5a)
+ assert_((abs(Hmag[idx]  1) < 0.015).all(), "Pass Band Close To Unity")
+
+
if __name__ == "__main__":
run_module_suite()
More information about the Scipysvn
mailing list