[Scipy-svn] r6973 - in trunk: doc/release doc/source scipy/signal scipy/signal/tests
scipy-svn@scip...
scipy-svn@scip...
Sun Nov 28 17:26:36 CST 2010
Author: warren.weckesser
Date: 2010-11-28 17:26:36 -0600 (Sun, 28 Nov 2010)
New Revision: 6973
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_fir_filter_design.py
Log:
ENH: signal: added the function firwin2 for constructing FIR filters with arbitrary frequency response.
Modified: trunk/doc/release/0.9.0-notes.rst
===================================================================
--- trunk/doc/release/0.9.0-notes.rst 2010-11-28 23:26:29 UTC (rev 6972)
+++ trunk/doc/release/0.9.0-notes.rst 2010-11-28 23:26:36 UTC (rev 6973)
@@ -114,6 +114,10 @@
The function ``scipy.signal.firwin`` was enhanced to allow the
design of highpass, bandpass, bandstop and multi-band FIR filters.
+The function ``scipy.signal.firwin2`` was added. This function
+uses the window method to create a linear phase FIR filter with
+an arbitrary frequency response.
+
The functions ``scipy.signal.kaiser_atten`` and ``scipy.signal.kaiser_beta``
were added.
@@ -131,7 +135,7 @@
Removed features
================
-The deprecated modules ``helpdmod``, ``pexec`` and ``ppimport`` were removed
+The deprecated modules ``helpmod``, ``pexec`` and ``ppimport`` were removed
from ``scipy.misc``.
The ``save`` method of the ``spmatrix`` class in ``scipy.sparse``, which has
Modified: trunk/doc/source/signal.rst
===================================================================
--- trunk/doc/source/signal.rst 2010-11-28 23:26:29 UTC (rev 6972)
+++ trunk/doc/source/signal.rst 2010-11-28 23:26:36 UTC (rev 6973)
@@ -64,6 +64,7 @@
bilinear
firwin
+ firwin2
freqs
freqz
iirdesign
Modified: trunk/scipy/signal/fir_filter_design.py
===================================================================
--- trunk/scipy/signal/fir_filter_design.py 2010-11-28 23:26:29 UTC (rev 6972)
+++ trunk/scipy/signal/fir_filter_design.py 2010-11-28 23:26:36 UTC (rev 6973)
@@ -1,7 +1,8 @@
"""Functions for FIR filter design."""
-import numpy
-from numpy import pi, ceil, cos, asarray
+from math import ceil, log
+import numpy as np
+from numpy.fft import irfft
from scipy.special import sinc
import sigtools
@@ -47,7 +48,7 @@
return beta
-def kaiser_atten(N, width):
+def kaiser_atten(numtaps, width):
"""Compute the attenuation of a Kaiser FIR filter.
Given the number of taps `N` and the transition width `width`, compute the
@@ -72,7 +73,7 @@
--------
kaiserord, kaiser_beta
"""
- a = 2.285 * (N - 1) * pi * width + 7.95
+ a = 2.285 * (numtaps - 1) * np.pi * width + 7.95
return a
@@ -90,7 +91,7 @@
Returns
-------
- N : int
+ numtaps : int
The length of the kaiser window.
beta :
The beta parameter for the kaiser window.
@@ -99,9 +100,9 @@
-----
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.kaiser(numtaps, beta, sym=0)
+ signal.get_window(beta, numtaps)
+ signal.get_window(('kaiser', beta), numtaps)
The empirical equations discovered by Kaiser are used.
@@ -123,28 +124,28 @@
# 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
+ numtaps = (A - 7.95) / 2.285 / (np.pi * width) + 1
- return int(ceil(N)), beta
+ return int(ceil(numtaps)), beta
-def firwin(N, cutoff, width=None, window='hamming', pass_zero=True, scale=True):
+def firwin(numtaps, cutoff, width=None, window='hamming', pass_zero=True, scale=True):
"""
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.
+ The filter will have linear phase; it will be Type I if `numtaps` is odd and
+ Type II if `numtaps` 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
+ ValueError exception is raised if firwin is called with `numtaps` even and
having a passband whose right end is at the Nyquist rate.
Parameters
----------
- N : int
+ numtaps : int
Length of the filter (number of coefficients, i.e. the filter
- order + 1). `N` must be even if a passband includes the Nyquist
+ order + 1). `numtaps` must be even if a passband includes the Nyquist
frequency.
cutoff : float or 1D array_like
@@ -181,14 +182,14 @@
Returns
-------
h : 1D ndarray
- Coefficients of length N FIR filter.
+ Coefficients of length `numtaps` 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
+ monotonically increasing, or if `numtaps` is even but a passband
includes the Nyquist frequency.
Examples
@@ -196,38 +197,42 @@
Low-pass from 0 to f::
- >>> firwin(N, f)
+ >>> firwin(numtaps, f)
Use a specific window function::
- >>> firwin(N, f, window='nuttall')
+ >>> firwin(numtaps, f, window='nuttall')
High-pass ('stop' from 0 to f)::
- >>> firwin(N, f, pass_zero=False)
+ >>> firwin(numtaps, f, pass_zero=False)
Band-pass::
- >>> firwin(N, [f1, f2], pass_zero=False)
+ >>> firwin(numtaps, [f1, f2], pass_zero=False)
Band-stop::
- >>> firwin(N, [f1, f2])
+ >>> firwin(numtaps, [f1, f2])
Multi-band (passbands are [0, f1], [f2, f3] and [f4, 1])::
- >>>firwin(N, [f1, f2, f3, f4])
+ >>>firwin(numtaps, [f1, f2, f3, f4])
Multi-band (passbands are [f1, f2] and [f3,f4])::
- >>> firwin(N, [f1, f2, f3, f4], pass_zero=False)
+ >>> firwin(numtaps, [f1, f2, f3, f4], pass_zero=False)
+ See also
+ --------
+ scipy.signal.firwin2
+
"""
# The major enhancements to this function added in November 2010 were
# developed by Tom Krauss (see ticket #902).
- cutoff = numpy.atleast_1d(cutoff)
+ cutoff = np.atleast_1d(cutoff)
# Check for invalid input.
if cutoff.ndim > 1:
@@ -236,7 +241,7 @@
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):
+ if np.any(np.diff(cutoff) <= 0):
raise ValueError("Invalid cutoff frequencies: the frequencies must be strictly increasing.")
if width is not None:
@@ -244,27 +249,27 @@
# of the Kaiser window, and set `window`. This overrides the value of
# `window` passed in.
if isinstance(width, float):
- atten = kaiser_atten(N, width)
+ atten = kaiser_atten(numtaps, 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:
+ if pass_nyquist and numtaps % 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))
+ cutoff = np.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
+ alpha = 0.5 * (numtaps-1)
+ m = np.arange(0, numtaps) - alpha
h = 0
for left, right in bands:
h += right * sinc(right * m)
@@ -272,7 +277,7 @@
# Get and apply the window function.
from signaltools import get_window
- win = get_window(window, N, fftbins=False)
+ win = get_window(window, numtaps, fftbins=False)
h *= win
# Now handle scaling if desired.
@@ -285,13 +290,160 @@
scale_frequency = 1.0
else:
scale_frequency = 0.5 * (left + right)
- c = cos(pi * m * scale_frequency)
- s = numpy.sum(h * c)
+ c = np.cos(np.pi * m * scale_frequency)
+ s = np.sum(h * c)
h /= s
return h
+# Original version of firwin2 from scipy ticket #457, submitted by "tash".
+#
+# Rewritten by Warren Weckesser, 2010.
+
+def firwin2(numtaps, freq, gain, nfreqs=None, window='hamming', nyq=1.0):
+ """FIR filter design using the window method.
+
+ From the given frequencies `freq` and corresponding gains `gain`,
+ this function constructs an FIR filter with linear phase and
+ (approximately) the given frequency response.
+
+ Parameters
+ ----------
+ numtaps : int
+ The number of taps in the FIR filter. `numtaps` must be less than
+ `nfreqs`. If the gain at the Nyquist rate, `gain[-1]`, is not 0,
+ then `numtaps` must be odd.
+
+ freq : array-like, 1D
+ The frequency sampling points. Typically 0.0 to 1.0 with 1.0 being
+ Nyquist. The Nyquist frequency can be redefined with the argument
+ `nyq`.
+
+ The values in `freq` must be nondecreasing. A value can be repeated
+ once to implement a discontinuity. The first value in `freq` must
+ be 0, and the last value must be `nyq`.
+
+ gain : array-like
+ The filter gains at the frequency sampling points.
+
+ nfreqs : int, optional
+ The size of the interpolation mesh used to construct the filter.
+ For most efficient behavior, this should be a power of 2 plus 1
+ (e.g, 129, 257, etc). The default is one more than the smallest
+ power of 2 that is not less than `numtaps`. `nfreqs` must be greater
+ than `numtaps`.
+
+ window : string or (string, float) or float, or None, optional
+ Window function to use. Default is "hamming". See
+ `scipy.signal.get_window` for the complete list of possible values.
+ If None, no window function is applied.
+
+ nyq : float
+ Nyquist frequency. Each frequency in `freq` must be between 0 and
+ `nyq` (inclusive).
+
+ Returns
+ -------
+ taps : numpy 1D array of length `numtaps`
+ The filter coefficients of the FIR filter.
+
+ Example
+ -------
+ A lowpass FIR filter with a response that is 1 on [0.0, 0.5], and
+ that decreases linearly on [0.5, 1.0] from 1 to 0:
+
+ >>> taps = firwin2(150, [0.0, 0.5, 1.0], [1.0, 1.0, 0.0])
+ >>> print(taps[72:78])
+ [-0.02286961 -0.06362756 0.57310236 0.57310236 -0.06362756 -0.02286961]
+
+ See also
+ --------
+ scipy.signal.firwin
+
+ Notes
+ -----
+
+ From the given set of frequencies and gains, the desired response is
+ constructed in the frequency domain. The inverse FFT is applied to the
+ desired response to create the associated convolution kernel, and the
+ first `numtaps` coefficients of this kernel, scaled by `window`, are
+ returned.
+
+ The FIR filter will have linear phase. The filter is Type I if `numtaps`
+ is odd and Type II if `numtaps` is even. Because Type II filters always
+ have a zero at the Nyquist frequency, `numtaps` must be odd if `gain[-1]`
+ is not zero.
+
+ .. versionadded:: 0.9.0
+
+ References
+ ----------
+ .. [1] Oppenheim, A. V. and Schafer, R. W., "Discrete-Time Signal
+ Processing", Prentice-Hall, Englewood Cliffs, New Jersey (1989).
+ (See, for example, Section 7.4.)
+
+ .. [2] Smith, Steven W., "The Scientist and Engineer's Guide to Digital
+ Signal Processing", Ch. 17. http://www.dspguide.com/ch17/1.htm
+
+ """
+
+ if len(freq) != len(gain):
+ raise ValueError('freq and gain must be of same length.')
+
+ if nfreqs is not None and numtaps >= nfreqs:
+ raise ValueError('ntaps must be less than nfreqs, but firwin2 was '
+ 'called with ntaps=%d and nfreqs=%s' % (numtaps, nfreqs))
+
+ if freq[0] != 0 or freq[-1] != nyq:
+ raise ValueError('freq must start with 0 and end with `nyq`.')
+ d = np.diff(freq)
+ if (d < 0).any():
+ raise ValueError('The values in freq must be nondecreasing.')
+ d2 = d[:-1] + d[1:]
+ if (d2 == 0).any():
+ raise ValueError('A value in freq must not occur more than twice.')
+
+ if numtaps % 2 == 0 and gain[-1] != 0.0:
+ raise ValueError("A filter with an even number of coefficients must "
+ "have zero gain at the Nyquist rate.")
+
+ if nfreqs is None:
+ nfreqs = 1 + 2 ** int(ceil(log(numtaps,2)))
+
+ # Tweak any repeated values in freq so that interp works.
+ eps = np.finfo(float).eps
+ for k in range(len(freq)):
+ if k < len(freq)-1 and freq[k] == freq[k+1]:
+ freq[k] = freq[k] - eps
+ freq[k+1] = freq[k+1] + eps
+
+ # Linearly interpolate the desired response on a uniform mesh `x`.
+ x = np.linspace(0.0, nyq, nfreqs)
+ fx = np.interp(x, freq, gain)
+
+ # Adjust the phases of the coefficients so that the first `ntaps` of the
+ # inverse FFT are the desired filter coefficients.
+ shift = np.exp(-(numtaps-1)/2. * 1.j * np.pi * x / nyq)
+ fx2 = fx * shift
+
+ # Use irfft to compute the inverse FFT.
+ out_full = irfft(fx2)
+
+ if window is not None:
+ # Create the window to apply to the filter coefficients.
+ from signaltools import get_window
+ wind = get_window(window, numtaps, fftbins=False)
+ else:
+ wind = 1
+
+ # Keep only the first `numtaps` coefficients in `out`, and multiply by
+ # the window.
+ out = out_full[:numtaps] * wind
+
+ return out
+
+
def remez(numtaps, bands, desired, weight=None, Hz=1, type='bandpass',
maxiter=25, grid_density=16):
"""
@@ -385,6 +537,6 @@
if weight is None:
weight = [1] * len(desired)
- bands = asarray(bands).copy()
+ bands = np.asarray(bands).copy()
return sigtools._remez(numtaps, bands, desired, weight, tnum, Hz,
maxiter, grid_density)
Modified: trunk/scipy/signal/info.py
===================================================================
--- trunk/scipy/signal/info.py 2010-11-28 23:26:29 UTC (rev 6972)
+++ trunk/scipy/signal/info.py 2010-11-28 23:26:36 UTC (rev 6973)
@@ -71,7 +71,9 @@
bilinear:
Return a digital filter from an analog filter using the bilinear transform.
firwin:
- Windowed FIR filter design.
+ Windowed FIR filter design, with frequency response defined as pass and stop bands.
+ firwin2:
+ Windowed FIR filter design, with arbitrary frequency response.
freqs:
Analog filter frequency response.
freqz:
Modified: trunk/scipy/signal/tests/test_fir_filter_design.py
===================================================================
--- trunk/scipy/signal/tests/test_fir_filter_design.py 2010-11-28 23:26:29 UTC (rev 6972)
+++ trunk/scipy/signal/tests/test_fir_filter_design.py 2010-11-28 23:26:36 UTC (rev 6973)
@@ -3,7 +3,7 @@
from numpy.testing import TestCase, run_module_suite, assert_raises, \
assert_array_almost_equal, assert_
-from scipy.signal import firwin, kaiserord, freqz, remez
+from scipy.signal import firwin, firwin2, kaiserord, freqz, remez
class TestFirwin(TestCase):
@@ -189,6 +189,71 @@
assert_raises(ValueError, firwin, 40, [.25, 0.5])
+
+class TestFirwin2(TestCase):
+
+ def test_invalid_args(self):
+ # `freq` and `gain` have different lengths.
+ assert_raises(ValueError, firwin2, 50, [0, 0.5, 1], [0.0, 1.0])
+ # `nfreqs` is less than `ntaps`.
+ assert_raises(ValueError, firwin2, 50, [0, 0.5, 1], [0.0, 1.0, 1.0], nfreqs=33)
+ # Decreasing value in `freq`
+ assert_raises(ValueError, firwin2, 50, [0, 0.5, 0.4, 1.0], [0, .25, .5, 1.0])
+ # Value in `freq` repeated more than once.
+ assert_raises(ValueError, firwin2, 50, [ 0, .1, .1, .1, 1.0],
+ [0.0, 0.5, 0.75, 1.0, 1.0])
+ # `freq` does not start at 0.0.
+ assert_raises(ValueError, firwin2, 50, [0.5, 1.0], [0.0, 1.0])
+
+ def test01(self):
+ width = 0.04
+ beta = 12.0
+ ntaps = 400
+ # Filter is 1 from w=0 to w=0.5, then decreases linearly from 1 to 0 as w
+ # increases from w=0.5 to w=1 (w=1 is the Nyquist frequency).
+ freq = [0.0, 0.5, 1.0]
+ gain = [1.0, 1.0, 0.0]
+ taps = firwin2(ntaps, freq, gain, window=('kaiser', beta))
+ freq_samples = np.array([0.0, 0.25, 0.5-width/2, 0.5+width/2,
+ 0.75, 1.0-width/2])
+ freqs, response = freqz(taps, worN=np.pi*freq_samples)
+ assert_array_almost_equal(np.abs(response),
+ [1.0, 1.0, 1.0, 1.0-width, 0.5, width], decimal=5)
+
+ def test02(self):
+ width = 0.04
+ beta = 12.0
+ # ntaps must be odd for positive gain at Nyquist.
+ ntaps = 401
+ # An ideal highpass filter.
+ freq = [0.0, 0.5, 0.5, 1.0]
+ gain = [0.0, 0.0, 1.0, 1.0]
+ taps = firwin2(ntaps, freq, gain, window=('kaiser', beta))
+ freq_samples = np.array([0.0, 0.25, 0.5-width, 0.5+width, 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 test03(self):
+ width = 0.02
+ ntaps, beta = kaiserord(120, width)
+ # ntaps must be odd for positive gain at Nyquist.
+ ntaps = int(ntaps) | 1
+ freq = [0.0, 0.4, 0.4, 0.5, 0.5, 1.0]
+ gain = [1.0, 1.0, 0.0, 0.0, 1.0, 1.0]
+ taps = firwin2(ntaps, freq, gain, window=('kaiser', beta))
+ freq_samples = np.array([0.0, 0.4-width, 0.4+width, 0.45,
+ 0.5-width, 0.5+width, 0.75, 1.0])
+ freqs, response = freqz(taps, worN=np.pi*freq_samples)
+ assert_array_almost_equal(np.abs(response),
+ [1.0, 1.0, 0.0, 0.0, 0.0, 1.0, 1.0, 1.0], decimal=5)
+
+ def test_nyq(self):
+ taps1 = firwin2(80, [0.0, 0.5, 1.0], [1.0, 1.0, 0.0])
+ taps2 = firwin2(80, [0.0, 30.0, 60.0], [1.0, 1.0, 0.0], nyq=60.0)
+ assert_array_almost_equal(taps1, taps2)
+
+
class TestRemez(TestCase):
def test_hilbert(self):
More information about the Scipy-svn
mailing list