[SciPy-user] firwin upgrades
Tom K.
tpk at kraussfamily.org
Wed May 6 21:33:55 EDT 2009
OK, here's a stab at the new functionality.
First the new routine (excerpt of filter_design.py), then the new test file
(entire test_filter_design.py).
What do y'all think about adding this functionality into the scipy
repository? Can I check it in somewhere, or should I submit a patch (to
where)?
The only thing I can think of in addition is that the multiband API we came
up with is a little bit non-intuitive (since it creates a passband starting
from nyquist at the end of the cutoff list, then alternates between
stopbands and passbands iterating backwards through the cutoff list). I
would propose accepting a list of (left, right) tuples as a clearer, easier
to use alternative. Any thoughts?
class PassBand(object):
def __init__(self, left, right):
self.left = left
self.right = right
def __repr__(self):
return "<PassBand %(left)g, %(right)g>"%self.__dict__
def firwin(N, cutoff, width=None, window='hamming', btype='pass',
scale=True):
"""FIR Filter Design using windowed ideal filter method.
Examples (preferred syntax):
firwin(N, [0, f], btype='pass') # low-pass from 0 to f
firwin(N, [0, f], btype='pass', window='nuttall') # specific window
firwin(N, [0, f], btype='stop') # stop from 0 to f --> high-pass
firwin(N, [f1, f2], btype='pass') # band-pass filter
firwin(N, [f1, f2], btype='stop') # band-stop filter
firwin(N, [f1, f2, f3, f4], btype='pass') # multiband filter:
starting from right, passes [f3,f4], stops [f2,f3], passes [f1,f2];
firwin(N, [f1, f2, f3, f4], btype='stop') # multiband filter:
starting from right, stops [f3,f4], passes [f2,f3], stops [f1,f2]
Also works:
firwin(N, f) # low-pass from 0 to f
firwin(N, f, btype='stop') # high-pass from f to 1
Parameters
----------
N -- length of filter (number of coefficients, = filter order + 1)
cutoff -- cutoff frequency of filter (normalized so that 1 corresponds
to
Nyquist or pi radians / sample) OR a list (or array) of cutoff
frequencies (that is, band edges). cutoff should be positive
and
monotonically increasing.
width -- 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.
window -- desired window to use. See get_window function for allowed
values.
btype -- the band type of filter, either 'pass' or 'stop'; specifies
whether the last band in the cutoff list is a passband or a
stopband
scale -- boolean, set to True to scale the coefficients so that the
frequency response is exactly unity at a certain
cutoff-dependent
frequency. The frequency is either:
0 (DC) if the first passband starts at 0
1 (Nyquist) if the first passband ends at 1
center of first passband otherwise
Returns
-------
h -- coefficients of length N fir filter.
"""
assert btype == 'pass' or btype == 'stop'
cutoff = numpy.atleast_1d(cutoff).tolist()
# build up list of pass bands starting from the end (near Nyquist)
bands = []
if btype == 'stop':
cutoff.append(1.0)
while cutoff:
right = cutoff.pop()
if cutoff:
left = cutoff.pop()
else:
left = 0.0
if left != right:
bands.insert(0, PassBand(left, right))
# build up the coefficients
alpha = N//2
m = numpy.arange(0,N) - alpha # time indices of taps
h = 0
for band in bands:
h += band.right*special.sinc(band.right*m)
h -= band.left*special.sinc(band.left*m)
# get and apply 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)
from signaltools import get_window
win = get_window(window,N,fftbins=1)
h *= win
# Now handle scaling if desired
if scale:
firstBand = bands[0]
if firstBand.left == 0:
scale_frequency = 0.0
elif firstBand.right == 1:
scale_frequency = 1.0
else:
scale_frequency = .5*(firstBand.left + firstBand.right)
h /= numpy.sum(h*exp(-1.j*pi*m*scale_frequency))
return h
import warnings
import numpy as np
from numpy.testing import TestCase, assert_array_almost_equal
from scipy.signal import tf2zpk, bessel, BadCoefficients, firwin, freqz
class TestTf2zpk(TestCase):
def test_simple(self):
z_r = np.array([0.5, -0.5])
p_r = np.array([1.j / np.sqrt(2), -1.j / np.sqrt(2)])
# Sort the zeros/poles so that we don't fail the test if the order
# changes
z_r.sort()
p_r.sort()
b = np.poly(z_r)
a = np.poly(p_r)
z, p, k = tf2zpk(b, a)
z.sort()
p.sort()
assert_array_almost_equal(z, z_r)
assert_array_almost_equal(p, p_r)
def test_bad_filter(self):
"""Regression test for #651: better handling of badly conditionned
filter coefficients."""
b, a = bessel(20, 0.1)
warnings.simplefilter("error", BadCoefficients)
try:
try:
z, p, k = tf2zpk(b, a)
raise AssertionError("tf2zpk did not warn about bad "\
"coefficients")
except BadCoefficients:
pass
finally:
warnings.simplefilter("always", BadCoefficients)
class TestFirwin(TestCase):
def check_response(self, h, expected_response, tol=.05):
N = len(h)
alpha = N//2
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, [0, f], btype='pass') # low-pass from 0 to f
self.check_response(h, [(.25,1), (.75,0)])
h = firwin(N+1, [0, f], btype='pass', window='nuttall') # specific
window
self.check_response(h, [(.25,1), (.75,0)])
h = firwin(N+2, [0, f], btype='stop') # 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], btype='pass') # band-pass filter
self.check_response(h, [(.1,0), (.3,1), (.5,0)])
h = firwin(N+4, [f1, f2], btype='stop') # band-stop filter
self.check_response(h, [(.1,1), (.3,0), (.5,1)])
h = firwin(N+5, [f1, f2, f3, f4], btype='pass', scale=False)
self.check_response(h, [(.1,0), (.3,1), (.5,0), (.7,1), (.9,0)])
h = firwin(N+6, [f1, f2, f3, f4], btype='stop') # 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, btype='stop') # 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 = [
([0, .5], (0, 1)),
([0.2, .6], (.4, 1)),
([.5, 1], (1, 1))
]
for cutoff, expected_response in cases:
h = firwin(N, cutoff, scale=False, window='ones')
hs = firwin(N, cutoff, scale=True, window='ones')
self.assertTrue(self.mse(h, [cutoff]) < self.mse(hs, [cutoff]),
'least squares violation')
self.check_response(hs, [expected_response], 1e-12)
--
View this message in context: http://www.nabble.com/firwin-upgrades-tp23246480p23418739.html
Sent from the Scipy-User mailing list archive at Nabble.com.
More information about the SciPy-User
mailing list