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

[Nathaniel Smith]

[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?

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

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

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

-- Your doc string still refers to this function as 'fir2' at one point

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

I don't have any authority or anything, but with those issues fixed
I'd vote for committing it. Tests seem reasonably complete.

-- Nathaniel