[SciPy-user] firwin upgrades

Tom K. tpk@kraussfamily....
Wed May 6 20:33:55 CDT 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