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

Nathaniel Smith njs@pobox....
Sun Nov 21 20:18:00 CST 2010


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?

> 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!).

> Great--thanks again for the comments.

np.

-- Nathaniel


More information about the SciPy-Dev mailing list