Here I will present some thoughts I have collected about the status of filtering and LTI system analysis in scipy.signal. This issue will hopefully lead to a discussion that can address factors that have frequently arisen in multiple PRs, but are generally of broader scope than can be settled in a single PR. While I have almost certainly overlooked relevant considerations, much of what follows was informed by reviewing the following discussions:

gh-3259, gh-3717, gh-4576, gh-4881

We should keep the following use cases in mind:

Control systems: Design a continuous-time LTI controller, analyze stability, simulate continuous time responses
DSP: Design discrete time FIR/IIR filters, simulate discrete time responses, apply filters, resample
Non-expert use: Just want to lowpass some data in a simple way

I welcome any and all feedback and opinions!

Problems:

• The number of filter/system representations has led to some disarray:
  • Current call signatures use a mishmash of system-tuples, lti objects, and individual coefficient arguments
  • Conversions between representations are not always optimal for a given use case
  • We mostly default to (b, a) form, which is numerically unstable for anything beyond low filter orders
• Inconsistent specification of sampling time/freq in many functions
• Inconsistent specification of angular or regular frequency units in many functions.

Possibilities

Object focused

We have added continuous time LTI objects, and a PR for discrete time objects is well underway. One option is to standardize all of the filter/LTI parts of the module on these LTI objects. I think that the main drawback of this approach would be unfamiliarity to MATLAB/Octave converts, who expect the lfilter(b, a, data) style syntax.

There are a few ways this could go.

One vision, which @endolith has advocated, is a system where the filter representation is mostly hidden from the user. Behind the scenes, the LTI object's methods would automatically use the right representation for the job, and convert as necessary. This allows things like

    fs = 48000
    signal = randn(fs)
    fc = 2000
    lowpass = butter(4, fc, fs=fs, analog=False, output='lti')
    w, h = lowpass.freqresp()  # Uses ZPK internally
    t, y = lowpass.dstep()  # Uses SS internally
    num = lowpass.num  # Automatically calculates TF form internally

Something similar, as @befelix has suggested, would not include automatic conversion, the only instantiable classes are LTI subclasses that define their representation. This would also allow the LTI subclasses to have methods unique to them. The methods would largely be the same, and methods could still internally convert. For example:

    lowpass_zpk = butter(4, fc, fs=fs, analog=False, output='zpk')
    t, y = dstep(lowpass_zpk)  # Converts to SS inside of dstep
    zeros = lowpass_zpk.zeros
    num = lowpass_zpk.as_tf().num

In each of these scenarios, there are multiple possibilities for how a given filter may be applied. For example, there could be one unified filter function that allows different methods (which may be a method of the LTI object and/or accept an LTI object as an argument)

    filter(system, data, method='lfilter')
    filter(system, data, method='filtfilt')
    filter(system, data, method='sosfilt')

Or, the filtering functions could take LTI objects as their arguments:

    lfilter(system, data)
    filtfilt(system, data)
    sosfilt(system, data)

Coefficient focused

Alternatively, maybe people don't want to move towards an object-heavy model. In this paradigm, call signatures would be limited to arguments for specific kinds of coefficients, or use the system-tuple approach.

    z, p, k = butter(4, fc, fs=fs, analog=False, output='zpk')
    b, a = zpk2tf(z, p, k)
    w, mag, ph = bode((z, p, k))
    w, h = freqz(b, a)
    t, out = dstep((A,B,C,D,dt))
    filtered = lfilter(b, a, x)

This is more involved than the object-focused approach, but more familiar to users coming from other languages.

Consistency

Regardless of what API style is agreed up, consistency is of utmost importance. Once an agreement it achieved, all of the filter design functions should be standardized in terms of their arguments.

For example:

• Are discrete-time system frequencies normalized by the nyquist or sampling frequency?
• Should the user specify sampling frequency via fs or sample time via dt? It has been suggested that fs is more likely to be an exact integer, so we can avoid testing equality of floats when comparing rates.
• Do functions like bode, freqresp return a frequency array in angular or normal frequency units?

Open Questions

• What API seems like the best option?
• It'll be a long slog to get to whatever vision is agreed upon with the usual deprecation procedure. What prospect exists for making a clean break?
• What pros/cons do people see in the suggestions above?
• What other forms could the API take?