[SciPy-user] FIR Filters

Bob.Cowdery at CGI-Europe.com Bob.Cowdery at CGI-Europe.com
Fri Feb 27 14:10:39 CST 2004


Travis

With your help I am making progress. I started off using numarray, then ran
into problems so I am now using Numeric. I only import scipy so I am using
whatever it imports which I believe is Numeric. I have tried to follow your
advice and move to complex types as I can see it makes sense. I don't think
I have done a complete job yet. The audio output has a pulse type modulation
on it which should not be there, however selecting different passbands is
definitely having the desired effect. I have copied all the code in below
for context. Yes the filtering is in the frequency domain which I am told
gives better stopband attenuation and good shape. The way the filter and
demodulation works is still a bit of a mystery to me as I didn't write the
original code. I would love to work out exactly what it does but I'm not
quite there yet. I can't really ask very targetted questions yet until I
understand more. The code cleary still has problems but I don't know why or
where. The good news is the processor is down to 50% now from 80%. If you do
have a few moments perhaps you could point out any more obvious errors.

Thanks
Bob

import time
import threading
import fastaudio as f
from scipy import *

class dsp( threading.Thread ):

    # signal processing
    def sig( self, samples ) :
        
        # allocate the filter arrays
        global m_filterMagn
        global m_filterPhase
        
        # convert to a numeric array of type double
        x = asarray(samples, typecode=Float64)
        
        # get a complex array ordered from the I,Q source
        complexin = x[1::2] + 1j*x[::2]
        
        # convert to complex frequency domain
        cart_complex_in_fd = fft(complexin)
        
        # apply filter to magnitude and phase
        magn = abs(complexin)*m_filterMagn
        phase = angle(complexin)*m_filterPhase
        
        # polar to cartesian in frequency domain
        demodulated = magn * cos(phase) + 1j*(magn * sin(phase))
        
        # convert to time domain
        complexout = ifft(demodulated)
        
        # convert back to interleaved array
        x[1::2] = real(complexout)
        x[::2] = real(complexout)
        
        # convert back to an integer array
        x = x.astype('l')
        
        # return sequence of processed samples
        return x
    
    def firwin(self, N, cutoff, width=None, window='blackman'):
        """
        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.
        """

        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 = signal.get_window(window,N,fftbins=1)
        alpha = N//2
        m = signal.Num.arange(0,N)
        return win*special.sinc(cutoff*(m-alpha))
    
    # calculate the complex filter
    def filter_calc( self, fLow, fHigh ):
    
        global m_filterMagn
        global m_filterPhase

        m_filterMagn = zeros(2048)
        m_filterPhase = zeros(2048)
        
        # create empty arrays to hold results
        Rh = zeros(2048)
        Ih = zeros(2048)
        
        # calculate the high/low cutoff freq as a fraction of the sample
rate
        fhc = fHigh/44100.0
        flc = fLow/44100.0
        
        # calculate fft bin size
        bin_size = 44100.0/4096.0
        
        # calculate no of bins in filter width
        bins = fHigh/bin_size
        
        # calculate FIR lowpass filter coefficients
        Rh = self.firwin(2048.0,(fhc-flc)/2.0)
        
        # compute USB filter by modulating the lowpass filter up
        filterPhase = 0.0
        filterFreq = (flc+fhc)*pi
        phaseSin = 0
        phaseCos = 0

        # modulate to USB centre frequency
        i = 0
        while i< 2048 :
            phaseSin = sin(filterPhase)
            phaseCos = cos(filterPhase)
            Ih[i] = Rh[i]*phaseSin
            Rh[i] = Rh[i]*phaseCos
            filterPhase = filterPhase+filterFreq
            if filterPhase >= 6.28318530717959 :
                filterPhase = filterPhase - 6.28318530717959
            i = i+1
            
        # for LSB create conjugate of USB filter
        # comment out for LSB
        Ih = -Ih
        
        # compute complex frequency domain response of LSB or USB filter
        complexin = Rh + 1j*Ih
        complexfd = fft(complexin)
        
        # Cartesian to polar gives the magnitude and phase components of the
fft result
        m_filterMagn = abs(complexfd)
        m_filterPhase = angle(complexfd)
                  
        
    # start streaming
    def run( self ):

        global running
        running = 1
        
        #s = f.stream(7500,       # sample rate
        #             2,          # number of channels
        #             'int16',    # sample format, choose from 'int8',
'int16', 'int32'.
        #             4096,       # frames per buffer
        #             16)         # maximum number of input buffers

        # open stream for input/output
        print "Create, open and start stream"
        global s
        s = f.stream(44100,2,'int16',2048,4)

        # open and start streaming
        s.open()
        s.start()

        print "Thread running..."
        
        # calculate the SSB filter
        self.filter_calc(300.0,3000.0)
        
        while( running ):

            # read next block of data from the input buffer
            samples = s.read()

            if( len( samples ) > 0 ):
                demod = self.sig( samples )
                # now, write the data back
                s.write( demod )

            # sleep for 10 ms
            # a stero block of 4096 samples takes 46 ms so 10ms should avoid
dropped blocks
            time.sleep(0.01)
            
        print "Thread exiting..."

def start():
    global DSP
    DSP = dsp()
    DSP.start()
    
def stop():
    # stop the stream
    global running
    running = 0
    global s
    s.stop()
    # and close it
    s.close()
    
def filter(id):
    # set the filter passband
    global DSP
    if id == 0:
        DSP.filter_calc(300.0,3000.0)
    elif id == 1:
        DSP.filter_calc(300.0,2000.0)
    elif id == 2:
        DSP.filter_calc(500.0,700.0)  

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


Bob.Cowdery at CGI-Europe.com wrote:

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

All the best in your education...

First,

Are you using numarray?  if so then that might cause a slower 
application than using Numeric (older but often faster arraymodule)

Comments on your code below:

>     def sigproc( self, samples ) :       
>        
>         # convert to a numarray array of type double
>         x = asarray(samples, typecode=Float64)
>
            is this numarray.asarray or Numeric.asarray ?

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

          Normally, real and imaginary parts of the data are handled 
together when a complex data type is available.    I would do
          in = x[1::2] + 1j*x[::2]
          fd = fft(in)

>        
>         # cartesian to polar in frequency domain
>         magn,phase = self.cartesian_to_polar(realfd, imagfd)
>        
>
          This seems strange to me, why are you filtering the magnitude 
and phase separately?  This is very non-standard.

>
>         # apply filter and demodulate
>
           I doubt you want a 2048 tap filter at this point.  Probably 
something a lot less.
           What kind of modulation are you demodulating?   Are you doing 
frequency-domain filtering then?

>                
>         # 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()
>
     Why are you going back to a python list?  A numeric array is a 
valid sequence.  Why does it need to be a list.  You have a lot of 
memory copying going on which is going to slow down everything.


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

If you have a complex-valued array, then abs(x) is the magnitude and 
angle(x) is the phase.

My colleagues down the hall do software radio and I have a lot of 
background in signal processing.  Let me know if you have more specific 
questions.

Thanks,

-Travis O.


_______________________________________________
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/20040227/d67b11e1/attachment-0001.htm


More information about the SciPy-user mailing list