[SciPy-user] FIR Filters

[Travis Oliphant]

> I am trying to get a windowed FIR filter to use in fast convolution 
> filtering.
>  
> I have tried using :
> signal.remez(2048,[300,3000],[1],Hz=44100)
> to generate a 2048 tap filter, which is what I have used previously in 
> Intel IPP. Firstly, this gives me an array full of QNAN errors. Infact 
> the only time I don't get errors is by using a very small number of 
> taps, say 10. Various other values of taps fails to converge or gives 
> QNAN errors. With the 2048 taps it is also very, very slow (about 75 
> seconds on a 1.2GHz). The IPP implementation is sub-second. I don't 
> know if it's throwing an exception on every QNAN or quite what is 
> hapenning.
>  
> I am also confused by not being able to specify a window type and I 
> really need lowpass not bandpass. I can see there is a function 
> get_window() but I don't know how this relates to an FIR filter when 
> the only other parameter is fftbins. Do I do something like tell it 
> the number of fftbins I want the filter to cover (I have 2048 bins of 
> around 10Hz each) then shift it to the correct part of the spectrum?


Bob,

Here is a simple windowed ideal lowpass filter algorithm.  The two 
important parameters are the number of taps and the cutoff frequency 
normalized
  so that 1 corresponds to pi radians/sample (i.e. the Nyquist 
frequency).  If you specify width it will choose a nice Kaiser window 
for you to try and reach
  that width of transistion region (from  passband to stopband).  
Otherwise, you can specify the window type (see signal.get_window for types)

This is now in SciPy CVS.



def firwin(N, cutoff, width=None, window='hamming'):
    """FIR Filter Design using windowed ideal filter method.

    Inputs:
   
      N      -- order of filter (number of taps)
      cutoff -- cutoff frequency of filter (normalized so that 1 
corresponds to
                  Nyquist or pi radians / sample)

      width  -- 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 -- desired window to use.

    Outputs:

      h      -- coefficients of length N fir filter.
    """

    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 = Num.arange(0,N)
    return win*special.sinc(cutoff*(m-alpha))