Changes Unknown when pulling 1f08fd4 on Eric89GXL:sosic into * on scipy:master*.

For compatibility with python 3 use xrange from scipy.lib.six, or just use range?

numpy.pad seems to have been added in 1.7.0 but scipy currently supports 1.5.1+

Changes Unknown when pulling 0288b4d on Eric89GXL:sosic into * on scipy:master*.

The weather in the cloud seems to be causing the reported error https://status.heroku.com/incidents/636. Or that's what I guess, because the TravisCI status page reports Heroku problems, and the builds do not usually exceed the 50 min limit.

Changes Unknown when pulling 01f30bb on Eric89GXL:sosic into * on scipy:master*.

Changes Unknown when pulling e4f4aa0 on Eric89GXL:sosic into * on scipy:master*.

I think having tf2sos is good because even though the coefficients are not optimal, implementing it as sos is still better than implementing it as a single big filter? also for compatibility with octave and matlab

Octave and Matlab both include a scaling/gain vector also, like [z, p, k] = sos2zp(sos, g) and [sos, g] = zp2sos(z, p, k), but it can be added at a later date without breaking anything by def sos2zpk(sos, g=1.0) and making g output optional:

sos = zpk2sos(z, p, k, scale=None)
sos, g = zpk2sos(z, p, k, scale='two')

http://www.mathworks.com/help/signal/ref/tf2sos.html

Matlab also lets you sort the sections by distance from the unit circle, but that could also be added later.

If the filter is already in tf form, is it still better to do SOS? I had assumed that it wouldn't be since (another) conversion seems like it would introduce more errors, but I haven't checked...

I was avoiding the scale-factor stuff because I wasn't sure if it was necessary here. I thought it was mostly for cases where you needed to stay within a particular range (fixed-point), but wasn't too useful for floating-point arithmetic. It simplifies a the code in a few places not to have to deal with it. That being said, if it seems more future-compatible to put it back in, I can.

Thinking about it more, I think you're probably right about the sos, g output at least being optional. I'll take a stab at implementing that on Monday.

From a quick browse, this looks promising. I agree with your comments about filtfilt and lfiltic.

If the filter is already in tf form, is it still better to do SOS? I had assumed that it wouldn't be since (another) conversion seems like it would introduce more errors, but I haven't checked...

Actually I guess you're right. Blue is sosfilt(tf2sos(ba)), green is lfilter(ba):

I was avoiding the scale-factor stuff because I wasn't sure if it was necessary here. I thought it was mostly for cases where you needed to stay within a particular range (fixed-point), but wasn't too useful for floating-point arithmetic.

Yes, I think that's the main purpose, which is why it's not used in matlab's sosfilt, but you can still use these functions to design fixed-point filters for use on other platforms. (lfilter can sort of simulate fixed-point filtering by using an object array of ints, but it can't operate on an actual int array?) But as I said it can always be added later. The g output from zpk2sos should be optional regardless, so the interface will not change.

The docstring of this function would benefit from an example, especially one that shows a problem where you need sos format to get good numerical precision.

freqz doesn't support sos yet, though, so how to show this?

z, p, k = signal.cheby2(9, 60, [0.25, 0.28], 'band', output='zpk')
sos = zpk2sos(z, p, k)

b, a = signal.cheby2(9, 60, [0.25, 0.28], 'band', output='ba')

nsamps = 2**14 # Number of impulse-response samples

x = np.concatenate(([1], np.zeros(nsamps-1))) # impulse

y1 = sosfilt(sos, x)
y2 = lfilter(b, a, x)

# Plot amplitude of frequency response 
Y1 = rfft(y1) # sampled frequency response
Y2 = rfft(y2) # sampled frequency response
f = np.linspace(0, 1, nsamps/2+1, endpoint=True)

plt.plot(f, 20 * np.log10(Y1))
plt.plot(f, 20 * np.log10(Y2))
plt.grid(True)
plt.xlim(0.2, 0.33)
plt.ylim(-100, 10)

or could show that the impulse response blows up for the ba form?

The one problem I'm hitting with the k=1. argument is in sosfilt. The current signature is:

sosfilt(sos, x, axis=-1, zi=None)

Adding k, we have:

sosfilt(sos, k, x, axis=-1, zi=None)

This isn't great, since then users must specify k. The other option is to keep the current signature, but allow for something like:

sosfilt((sos, k), x, axis=-1, zi=None)

The disadvantage is it's inconsistent with the other filtering functions, where you'd do lfilter(*zpk2tf(*zpk), x), here you'd do sosfilt(zpk2sos(*zpk), x), as it loses the * preceding the first argument. This will probably become moot with #3259, but @endolith did you have an idea for this?

Well, does it need to be in sosfilt if that's a floating-point implementation? Or should sosfilt accept object arrays, etc, too? (Oh it does already.) Matlab and octave's sosfilt don't have the k argument. Maybe make it the last argument?

No, I suppose not. We can add a note to zpk2sos, etc. that optionally return k making it clear that it shouldn't be separate if sosfilt (or eventually filtfilt) is to be used.

Come to think of it, how about we just make it so that there's a function to extract the overall gain from or apply an overall gain to a given sos system? This way, everything can use just sos directly (which should cover 95%+ of use cases, probably), and people who need the gain separate (e.g., for designing fixed-point systems) can still get it?

I think would be good to follow matlab's options for zpk2sos, except that scale=None wouldn't produce a g vector at all, and is the default:

def zpk2sos(z, p, k, order='up', scale=None, zeroflag=False):

so for example:

sos = zpk2sos(z, p, k)
sos = zpk2sos(z, p, k, order='down')
sos, g = zpk2sos(z, p, k, order='up', scale='inf')
sos, g = zpk2sos(z, p, k, 'down', 'two', True)

and yes, leave out the k/g term from sosfilt and the docstring for zpk2sos could say for use with sosfilt, use scale=None and if people want to model fixed-point processing they can just make their own loop of lfilters?

Or add it to the end of sosfilt for the unusual cases that need it:

def sosfilt(sos, x, axis=-1, zi=None, g=1):

?

Those options all make sense. However, it just occurred to me that adding parameters for scale (including having a second output for zpk2sos), order, and zeroflag can be added to the appropriate functions in a future PR, so the current iteration is future-compatible in that sense. Given this idea, I'd prefer to wait to add these features, and instead focus on ensuring the majority-of-potential-use-cases functionality (zpk2sos for use in sosfilter in scipy.signal) implemented here works properly first. It should reduce the probability of errors in implementation (esp. WRT corner cases) and reduce reviewing load. Does that sound reasonable?

Yes, that's what I meant by "can always be added later" :) The only exception is g parameter in sosfilt, but if it's placed as the last parameter, then that could be left for another PR too, if it's even needed.

Yeah, I can't think of a better place for it to go than last. I suspect it
would be the least likely of those parameters to be used.

Comments addressed, added the following example showing sosfilt stability compared to tf (the simplest example I could come up with off the top of my head):

    >>> [b, a] = butter(6, [0.10, 0.105], 'band', output='ba')
    >>> sos = butter(6, [0.10, 0.105], 'band', output='sos')
    >>> x = zeros(1500)
    >>> x[0] = 1.
    >>> y_tf = lfilter(b, a, x)
    >>> y_sos = sosfilt(sos, x)
    >>> plot(y_tf, 'r')
    >>> plot(y_sos, 'k')

Changes Unknown when pulling d30f250 on Eric89GXL:sosic into * on scipy:master*.

This is probably not right:

sos2tf(zpk2sos(z, p, k))
Out[53]: 
(array([ 0.1667,  0.5   ,  0.5   ,  0.1667,  0.    ]),
 array([ 1.    ,  0.    ,  0.3333,  0.    ,  0.    ]))

zpk2tf(z, p, k)
Out[54]: 
(array([ 0.1667,  0.5   ,  0.5   ,  0.1667]),
 array([ 1.    ,  0.    ,  0.3333,  0.    ]))

The first is just the second times z^-1/z^-1, so they're equivalent, but the extra 0s should probably not be there.

I'm not sure the best way to deal with that. On the one hand, you might expect sos2tf to always give a multiple of three for poles and zeros. On the other hand, the extras are unnecessary. I guess it sos2zpk we could remove any zero-valued poles or zeros?

No I don't see a conflict with zeroflag other than thinking about the name and both options being about "pairing" in some way.

Maybe "nearest" instead of "simple"?

I'm okay with nearest and keep_odd -- @WarrenWeckesser are you happy with those too?

Or should it also be a Boolean flag? keepodd = true, false?

I'm okay with the current API -- I think that string names for the algorithms are a good way to go for now. If we want to expand to zeroflag in the future, it could be done as pairing='nearest_zeroflag, pairing='nearest', zeroflag=True, pairing=dict(algorithm='nearest', zeroflag=True) and these are all compatible with the current iteration. I can't think of a better / more future compatible option in any case...

Last amended to change the algorithm name from simple to nearest.

I guess I'm more or less ok with that, too. With my previous comment I meant that if "nearest" is the default algorithm (there will never be a pairing='farthest'), and "zeroflag" or "keepodd" are additional constraints on top of that, it might make more sense for it to be zpk2sos(z, p, k, keep_odd=False, zero_flag=False).

But combining pairing flags into a single parameter is ok, too. It could even be updated to something like zpk2sos(z, p, k, pairing={'keep_odd': True, 'zero_flag': True) or zpk2sos(z, p, k, pairing=('keep_odd', 'zero_flag')) without breaking backwards compatibility, so it doesn't really matter.

I'm not trying to be a pain; it just seems difficult to change an API once it's been created. :)

Making API changes is indeed a pain, and I agree it's important to do our due diligence in these design choices, so I do appreciate the discussion and progress we've made. At the same time, we have missed one release cycle (0.15) by debating these and related issues, so I'm starting to think that it's likely that the feedback we will get through people actually using this tool will out-weigh our ability to accurately predict what API or practical changes will be most useful to end users. I am thus a bit wary of spending more time tweaking the API, assuming that we are agreed that the API is pretty future compatible, or at least that there isn't a better formulated proposal available...

Yes, that's why I wanted to leave pairing/scaling/ordering for other PRs, since the filters mostly work fine without them. But yes, I'm happy with it the way it is. A single parameter called pairing that controls all pairing options is fine, and the options can be changed/reworded in the future in a backwards-compatible way if needed. What else needs to be done to merge this?

Nothing from my end. If you're happy, then we just need to convince @WarrenWeckesser :)

Consider me convinced. I'll have a couple follow-up PRs, but this PR has cooked enough.

@Eric89GXL and @endolith: Thanks for all the hard work, and for your patience. :)

Thanks for so many rounds of review, they were useful. Feel free to ping me for review / comments on follow-up or other signals-related PRs.

Yay! Thanks @Eric89GXL and @WarrenWeckesser

Great job guys! This must be one of the most thoroughly hashed out PRs in this repo

First follow-up: #4447

Adding sosfreqz in #4465

I found a problem with the handling of the zi argument of sosfilt. Fixed here: #4473

We may find that once this code is released, more experienced engineers will want a different ordering for the analog case

Actually, I just realized this isn't just a matter of pairing of roots/ordering of stages. It's also the coefficient position:

b, a = butter(1, 1, analog=True) produces b = [1], a = [1, 1] → H(s) = 1/(s+1)

sos = tf2sos(b, a) produces [1, 0, 0, 1, 1, 0]

this stage is b = [1, 0, 0], a = [1, 1, 0] → H(s) = s^2/(s^2+s) = s/(s+1) which is not the same as above. For digital filters, trailing zeros don't matter, but for analog filters, leading zeros don't matter. The correct result should be [0, 0, 1, 0, 1, 1]

(I wanted to use this to evaluate high-order analog filters accurately, like an sosfreqs, and was confused as to why it was changing lowpass into highpass filters. But that's a possible use for it.)

So docstring should probably say it's not meant for analog at all for now?