[SciPy-user] FIR Filters
Travis Oliphant
oliphant at ee.byu.edu
Thu Feb 26 15:38:17 CST 2004
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.
More information about the SciPy-user
mailing list