[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