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

[Warren Weckesser]

On Sun, Nov 21, 2010 at 8:18 PM, Nathaniel Smith <njs at pobox.com> wrote:

> On Sun, Nov 21, 2010 at 5:28 PM, Warren Weckesser
> <warren.weckesser at enthought.com> wrote:
> >> -- 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. :)
>
> Certainly :-) But some convenience handling might be useful...
>
> > 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.
>
> You mean, you've used this in the past, not, you've used that in
> coding up this function, right?
>


Yes, that's right.



>
> > 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.
>
> Matlab appears to use a similar tweak (see
> http://www.mathworks.com/help/toolbox/signal/fir2.html#f7-957922), and
> I guess you don't get more widely used than that. Not that this is a
> guarantee of quality... (In other news, this strange choice of
> specifying the ramp in terms of grid-units is probably their fault
> too.)
>
> > 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?
>
> If I'm reading the code right, Octave's algorithm is:
>  -- find all the jump discontinuities
>  -- for each one, subtract ramp_width/2 from the left edge, and add
> ramp_width/2 to the right edge (where ramp_width is ramp_n/grid_n)
>  -- use interpolation to recalculate the corresponding *gains* by
> calculating what the gain *would* have been at the *new* frequency
> points if we hadn't done this ramping stuff.
>  -- then do the linear interpolation step as usual
>
> So, AFAICT the code says that this only happens for jump
> discontinuities. Of course, the comments say it happens for "any
> transition less than ramp_n units". Who you gonna trust...
>
> >> - 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."
>
> I think the intention of that rule is more like "don't make your
> primary reference some random schmuck's blog entry" -- not to force
> people to go trek to the library just because. That's a real published
> book
>  http://www.dspguide.com/editions.htm
> that happens to also be online, and the link's been stable since 2006
>
> http://web.archive.org/web/20060612222046/http://www.dspguide.com/ch17/1.htm
> So I'd add it :-).
>
> >> -- 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.
>
> There are two orthogonal issues here -- whether people specify
> sampling or nyquist rate, and compatibility. If we let people specify
> the sampling rate instead of the Nyquist frequency, and we make the
> default sampling rate be "2", then compatibility will be preserved.
> But really *anything* would be better than what we have now (if
> nothing else, this would force it to be clear what the IIR functions
> expect -- right now it's just undocumented!).
>


Some of the IIR functions (e.g. ellipord, buttord) say that the frequencies
are "normalized from 0 to 1 (1 corresponds to pi radians / sample)".
Whether one immediately sees that that means they are normalized by the
Nyquist rate depends on one's background.  The docstrings of other IIR
functions (e.g. ellip, butter) definitely need work.


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