[SciPy-Dev] Fwd: [SciPy-User] FIR filter with arbitrary frequency response

[Warren Weckesser]

Forwarded to the mailing list (didn't "reply to all" the first time).


---------- Forwarded message ----------
From: Warren Weckesser <warren.weckesser at enthought.com>
Date: Sun, Nov 21, 2010 at 5:00 PM
Subject: Re: [SciPy-User] FIR filter with arbitrary frequency response
To: Nathaniel Smith <njs at pobox.com>


Nathanial, thanks for the great feedback.  Comments below...


On Sat, Nov 20, 2010 at 11:25 PM, Nathaniel Smith <njs at pobox.com> wrote:

> [Moving this to scipy-dev]
>
> On Sat, Nov 20, 2010 at 8:09 PM, Warren Weckesser
> <warren.weckesser at enthought.com> wrote:
> > There is one implemented as the function firwin2 currently under review
> in
> > this ticket:
> >     http://projects.scipy.org/scipy/ticket/457
> > It will eventually be added to scipy.signal.  Included in the ticket is a
> > patch file that can be applied to the latest version of the scipy source,
> > and also a copy of just the updated file fir_filter_design.py.  You could
> > grab that and try it stand-alone (but you would have to comment out the
> > local import of sigtools, which is used by the remez function in
> > fir_filter_design.py).
> >
> > Feedback would be appreciated, so if you try it, be sure to write back
> with
> > comments or questions.
>
> Ha, what luck! Thanks!
>
> Some quick thoughts just based on reading the source code:
>
> -- Your handling of discontinuities in the desired frequency response
> looks questionable to me (not that I'm an expert). Octave's version of
> this function has a 'ramp_n' parameter that controls how much they
> smooth out discontinuities, documented as "transition width for jumps
> in filter response. Defaults to grid_n/20; a wider ramp gives wider
> transitions but has better stopband characteristics." Your function,
> OTOH, always uses the narrowest possible ramp.
>
> Octave's choice to make ramp_n defined in number-of-grid-point units
> seems bizarre to me; I wouldn't copy that detail. But it probably
> would be good to have some way to specify smarter handling of
> discontinuities?
>
>

As you noticed, currently it does not attempt to smooth the requested
frequency response.  If the user wants a ramp, they can just define their
function with a ramp instead of a discontinuity. :)

To handle transition width issues, I've used the Kaiser design method, which
provides formulas for the number of taps and the Kaiser parameter beta in
terms of the desired attenuation and transition width.  This works
reasonably well.

I have nothing against adding some sort of "auto-smoothing" option.  I'd
prefer something with a sound mathematical basis, rather than an ad hoc
tweak.  But if there is a widely used heuristic that works, it should
probably be included in the function.

In the Octave example, what happens if the user includes a ramp in the
response that isn't as wide as  the transition width defined by grid_n/20?



> -- I suspect the first part of the docstring:
>
> "    From the given set of frequencies and gains, the desired response is
>     constructed in the frequency domain.  The inverse FFT is applied to the
>     desired response to create the associated convolution kernel, and the
>     first `numtaps` coefficients of this kernel, scaled by `window`, are
>     returned."
>
> will not make sense to anyone who doesn't already know how this works.
>
> If it were me, I'd start by saying:
> "Given any desired combination of frequencies (`freq`) and gains
> (`gain`), this function constructs an FIR filter with linear phase and
> (approximately) the given frequency response."
> And I'd put the terse description of how it accomplishes this trick
> lower down in the docstring.
>
>
Thanks--after rereading, I see that it does need to be improved.



> - In info.py, I'd expand the description of firwin to "Windowed FIR
> filter design, with standard high/low/band-pass frequency response",
> to make the contrast between it and firwin2 clear.
>
>
Good idea.



> - I'd also add a reference to http://www.dspguide.com/ch17/1.htm;
> obviously any good DSP book will explain this idea, and the reference
> you already have is fine, but it's much easier to click a link than to
> go to the library... and the explanation is much less terse than the
> one in the docstring :-).
>
>
I'm OK with web links, but the numpy/scipy documentation guidelines state
"Referencing sources of a temporary nature, like web pages, is discouraged."



> -- Your doc string still refers to this function as 'fir2' at one point
>
>
Thanks, I'll fix that.



> -- I have a mild preference for specifying the sampling frequency
> instead of the nyquist frequency, but that's not a big deal. What's
> really bad is that we have no standard for this (it looks like the
> only other filter design function that lets you work in sampling units
> instead of normalized units is remez, and it uses the terribly named
> 'Hz'). Maybe I should start another thread to discuss that...
>
>
Personally, I'm fine with either sample rate or Nyquist rate, but there are
other functions where the default is to express frequencies relative to the
Nyquist rate (e.g. all the IIR functions).   I also noticed that remez has
the unfortunate keyword 'Hz'.   For the sake of consistency, I would like to
deprecate that keyword, and add 'nyq' (or perhaps 'nyquist', or...) as the
preferred keyword.  This option should also be added to firwin(), so all
these FIR filter design functions have a consistent interface.



> I don't have any authority or anything, but with those issues fixed
> I'd vote for committing it. Tests seem reasonably complete.
>
>
Great--thanks again for the comments.

Warren
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.python.org/pipermail/scipy-dev/attachments/20101121/c6650529/attachment.html>