[SciPy-user] FIR Filters

[Travis Oliphant]

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.