[Scipy-svn] r6872 - in trunk: doc/release doc/source scipy/signal scipy/signal/tests
scipy-svn at scipy.org
scipy-svn at scipy.org
Sat Nov 13 12:19:03 EST 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