[SciPy-user] FIR Filters

Bob.Cowdery at CGI-Europe.com Bob.Cowdery at CGI-Europe.com
Thu Feb 26 13:27:40 CST 2004


Travis

Thanks very much for the responses. Don't confuse me with someone who knows
what they are doing (but I am working on it!). To give you some background I
am converting a .NET application that uses Intel IPP to Python. The
application is for a software defined radio. At the moment I have audio
capture which contains my quadrature samples in blocks of 4096 (2048 each
for I&Q) with 44100 sampling rate (46ms per block). I have a simple dsp
process to test each signal processing call. At the moment this should
output exactly what it receives, i.e just the audio input as the demodulator
needs the FIR filter and some magic applied at the point indicated in the
code. When I have written the demodulator I can test firwin() properly. At
the moment I can see it outputs valid data, and it's fast which is good. The
code below is the very simple path. I have had to write two functions for
polar/cartesian conversions as I didn't find any in the library. I don't
know if these work right as I can hear clicking sounds on the audio not
there when they are commented out. However, it could be audio overrun as I
am up at 75% processor. The functions are also quite slow (adding 30% to my
processor loading). Is there a more efficient way to do this?

    def sigproc( self, samples ) :        
        
        # convert to a numarray array of type double
        x = asarray(samples, typecode=Float64)
        
        # first to do is sort into real and image arrays from interleaved
array
        realin = x[1::2]  # every other element starting at index 1
        imagin = x[::2]   # every other element starting at index 0
                
        # convert to frequency domain
        realfd = fft(realin)
        imagfd = fft(imagin)
        
        # cartesian to polar in frequency domain
        magn,phase = self.cartesian_to_polar(realfd, imagfd)
        
        # apply filter and demodulate
                
        # polar to cartesian in frequency domain
        realfd, imagfd = self.polar_to_cartesian(magn, phase)
        
        # convert to time domain
        realtd = ifft(realfd)
        imagtd = ifft(imagfd)
                
        # convert back to interleaved array
        x[::2] = realtd.astype(Float64)
        x[1::2] = imagtd.astype(Float64)
        
        # convert back to an integer array
        x = x.astype('l')
        
        # return python list of processed samples
        return x.tolist()
    
    # utility functions not in the library
    def polar_to_cartesian(self, magn, phase):
        return magn * cos(phase), magn * sin(phase)

    def cartesian_to_polar(self, real, imag):
        # convert the arrays real, imag into the polar system
        # magnitude array
        magn = sqrt(real**2 + imag**2)
        # phase array
        phase =
choose(equal(real,0),(arctan(imag/real),(choose(less(imag,0),((pi/2),-(pi/2)
)))))
        return magn, phase


-----Original Message-----
From: Travis Oliphant [mailto:oliphant at ee.byu.edu] 
Sent: 26 February 2004 07:10
To: SciPy Users List
Subject: Re: [SciPy-user] FIR Filters



> 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))
   
   

_______________________________________________
SciPy-user mailing list
SciPy-user at scipy.net
http://www.scipy.net/mailman/listinfo/scipy-user

*** Confidentiality Notice *** Proprietary/Confidential
Information belonging to CGI Group Inc. and its affiliates
may be contained in this message. If you are not a recipient
indicated or intended in this message (or responsible for
delivery of this message to such person), or you think for
any reason that this message may have been addressed to you
in error, you may not use or copy or deliver this message
to anyone else.  In such case, you should destroy this
message and are asked to notify the sender by reply email.
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://www.scipy.net/pipermail/scipy-user/attachments/20040226/c79bdfc0/attachment-0001.htm


More information about the SciPy-user mailing list