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

Warren Weckesser warren.weckesser@enthought....
Wed Nov 24 22:33:36 CST 2010

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

> On Sun, Nov 21, 2010 at 5:28 PM, Warren Weckesser
> <warren.weckesser@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.

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

More information about the SciPy-Dev mailing list