ENH: Add fs= parameter to filter design functions

[endolith]

This is just a more convenient way of specifying normalized frequencies for digital filters. Instead of having to remember that butter uses fc/(fs/2), while freqz uses fc/(fs/(2*pi)), for instance, you can just specify the frequencies in Hz and the sampling rate in Hz.

Old way:

>>> fc = 1000
>>> fs = 3000
>>> b, a = butter(2, fc/(fs/2))
>>> abs(freqz(b, a, [fc/(fs/(2*pi))])[1])
array([ 0.70710678])

New way:

>>> fc = 1000
>>> fs = 3000
>>> b, a = butter(2, fc, fs=fs)
>>> abs(freqz(b, a, [fc], fs=fs)[1])
array([ 0.70710678])

I proposed this a while ago and no one seemed interested, but it's easy enough to implement so we'll see what people think.

Examples of confusion:

• https://dsp.stackexchange.com/q/46915/29
• https://dsp.stackexchange.com/q/24989/29#comment45751_24992
• https://mail.python.org/pipermail/scipy-user/2008-March/015944.html

[WarrenWeckesser]

Thanks for this, it is a much needed enhancement. I'll try to do a review sometime this weekend (nag me if I don't!).

One thing I'll note now is that unit tests are needed for the new option.

[larsoner]

Agreed, even if it's redundant the new code is much more readable. +1 from me

[larsoner]

I've optimistically marked this for 1.1 (mid-March)

[endolith]

@WarrenWeckesser

One thing I'll note now is that unit tests are needed for the new option.

Yeah I wanted to see if people liked it before doing the rest of the work. Sounds like people do, though.

[endolith]

Instead of default of fs=None, I could make it fs=2 or fs=2*pi or whatever is appropriate for each function? I don't really have an opinion on which is better.

Also can be added to

• iirnotch, iirpeak
• buttord, cheb1ord, etc.
• group_delay
• iirdesign
• band_stop_obj ?

Already had fs:

• firwin, firwin2, firls, remez (all of which use fs=None)

[larsoner]

I prefer explicit fs=2. and similar but IIRC @WarrenWeckesser prefers fs=None as an alias because it can be easier to wrap.

Maybe we can both be happy by using fs=2. as the default but support fs=None as "use the default" (in the code this is the one-liner fs = 2. if fs is None else fs)? That way introspection of the call immediately shows the default (rather than having to look in the docstring; what I like) and functions that wrap these can still pass None to get default behavior (what I think @WarrenWeckesser wants)?

As for what needs to be changed, I'd recommend taking a brief scan of the test_* functions in scipy.signal and seeing where else fs= would make things simpler. But I also don't mind multiple PRs if we find other cases in the future, so I don't think this needs to be exhaustive.

[endolith]

Should freqz, ellipord, etc. return w in fs units if fs specified?

Like if you do

w, h = freqz(b, a, [100, 200, 300], fs=1000)

should w == [100, 200, 300]? What about

w, h = signal.freqz(b, a, 5, fs=8000)

Should w == [0, 800, 1600, 2400, 3200] in Hz? What about

w, h = signal.freqz(b, a, fs=8000)

?

I think it makes sense to be in Hz for all of them. Likewise for buttord, etc. if your inputs are in Hz then it makes sense for output cutoff frequency to also be in Hz, which you can then feed to butter, also using fs parameter.

[endolith]

Added examples that use fs and sos to filter a real signal:

https://5541-1460385-gh.circle-artifacts.com/0/home/ubuntu/scipy/doc/build/html-scipyorg/generated/scipy.signal.butter.html

Also added test for some functions but not all yet

[endolith]

So doctest failed because of

Failed example:
    ax1.axis([0, 1, -2, 2])
Expected nothing
Got:
    [0, 1, -2, 2]

But there are other examples with

>>> ax.axis((10, 1000, -100, 10))

that don't cause failures?

[larsoner]

Looks like ax.axis is explicitly ignored:

https://github.com/scipy/scipy/blob/master/tools/refguide_check.py#L518

This could be improved/extended with e.g. regex that permits names like ax1 and ax2, or just use # doctest:+SKIP or similar with the offending lines.

[larsoner]

Maybe the thing to do is to change it to .axis(from ax.axis. It should be general enough. Same with plt.plot(, it seems like this should be .plot(.

This should probably be a separate PR, though, so that it gets proper discussion as it modifies the test suite.

[ev-br]

Please no doctest:+SKIP, just fix the refguide-checker.

[WarrenWeckesser]

I can live with fs=2.0 (or whatever the appropriate constant is for the function being changed) for the default keyword argument.

[WarrenWeckesser]

@endolith asked

Should freqz, ellipord, etc. return w in fs units if fs specified?

My expectation was that fs would set the units for everything--input and output--so the frequencies returned would be in the same units as fs. I'm trying to think of a reason for it to be otherwise, but so far I'm not succeeding. Can anyone make a compelling argument for the returned frequencies to not match fs?

[larsoner]

Can anyone make a compelling argument for the returned frequencies to not match fs?

As long as the new defaults give the same old results, we should be good to go on this. I have the same expectation as you, that all outputs would be modified appropriately.

[WarrenWeckesser]

Another observation: the functions periodogram, welch, csd, coherence, spectrogram and stft all accept the argument fs already, and the frequencies that they return are in the same units as fs. Given that existing behavior, for API consistency the functions updated in this pull request should also return frequencies in the same units as fs.

[endolith]

I'm not sure what to do about the failures.

operator.index(array([123])) is raising DeprecationWarning: converting an array with ndim > 0 to an index will result in an error in the future

but it's supposed to produce an error, and does in Python 3.6.3, but not in Python 3.4.6 the test is running in?

[larsoner]

but it's supposed to produce an error, and does in Python 3.6.3, but not in Python 3.4.6 the test is running in?

Different Travis runs use different NumPy versions, this is probably to blame. The 3.4.6 appears to use NumPy 1.8.2. The code for freqz probably needs to be improved to deal with this case.

[endolith]

So this goes back to #7351 and how to verify int parameters.

I made a chart of how different types are handled here: https://stackoverflow.com/a/48940855/125507

I did this using .index() since that seems to be what people want.

[larsoner]

operator.index is indeed what we settled on for int checks.

I made some small comments on the code, let's see what @WarrenWeckesser says about the fs default, and once that's settled I'll take a look at the examples and tests.

[larsoner]

You could save on branching/lines a little bit with else instead of elif here:

else:
    try:
        N = operator.index(worN)
    except TypeError:
        w = worN
    else:
        w = findfreqs(b, a, N)

[endolith]

It needs to do the ndim check to interpret array([8]) as "measure at point 8", rather than "measure at 8 points" since apparently older numpy versions returned .index() on 1-element N-D arrays.

[WarrenWeckesser]

Using, for example, worN=3.0, results in UnboundLocalError: local variable 'w' referenced before assignment

[endolith]

I modified it to fix the unbound error. Now it calls index twice, but that's not a big deal? Otherwise I'd have to duplicate other things and it would be harder to read. Still need to write tests for this condition

[larsoner]

Calling index twice seems fine if it leads to the simplest code.

[larsoner]

Need an extra newline before the .. directive otherwise it doesn't render properly (not seen as separate paragraph)

https://5706-1460385-gh.circle-artifacts.com/0/home/ubuntu/scipy/doc/build/html-scipyorg/generated/scipy.signal.freqz.html#scipy.signal.freqz

[larsoner]

This is a case where setting fs=2*pi makes things simpler. Here you can just say they are in units of the sample rate, and you can get rid of the conditionals below. We can make fs=None behave like 2 * pi if this would make achieving "use-the-default" behavior easier for people. @WarrenWeckesser WDYT?

[endolith]

The worN=None arguments follow the same pattern, right? They actually default to numbers, too.

[larsoner]

Yes, I am arguing against that design choice as it hurts argument introspection. I think we should ideally change this everywhere, so at least avoid propagating the problem here.

[larsoner]

I just remembered this was discussed in #7788. I'll try to move that forward so we can proceed here.

[endolith]

Maybe that should be done in a separate PR then? (And changing worN=512 as well?) Changing these is 100% backwards compatible, right?

[WarrenWeckesser]

Sorry for disappearing from this conversation. If no one else objects to fs=2*np.pi, then go ahead. (Of course, introspection of the signature will then show fs=6.283185307179586.)

[endolith]

@WarrenWeckesser We could make an object that has repr of 2*pi but otherwise behaves as a float.

class Foo(float):
    def __new__(self, value):
        return float.__new__(self, value)

    def __init__(self, value):
        float.__init__(value)

    def __repr__(self):
        return '2*pi'


foo = Foo(1)

def freqz(fs=foo):
    print(fs)

[WarrenWeckesser]

That's a neat idea, but I don't think the extra code is worth it just to fix a cosmetic issue.

[larsoner]

Back to this line of code -- we should just say that the frequencies are in units of fs (and if you want, add that the default fs=2*pi implies the units are radians/sample)

[endolith]

Ok updated the version number

[ilayn]

I guess everybody's happy now. In it goes then. Thanks a lot @endolith !!

[rgommers]

Nice:) @endolith would you mind adding a few sentences to the release notes at https://github.com/scipy/scipy/wiki/Release-note-entries-for-SciPy-1.2.0?