[breeze.signal] IIR filter, Biquad filter

[ktakagaki]

quoted from: https://github.com/almoehi/breeze/wiki/Digital-Signal-Processing#filter-design-thoughts

filter design thoughts

move filter design functions into separate sub packages for fir, iir, sos filters and keep the naming of the actual filter designs consistent across the packges:

import breeze.signal.filter.iir
designButterworth[IIRKernel1D](normCutOff: T, order: Int, flavour: FilterFlavour = FilterFlavour.LoPass, opts: Seq[FilterOpts]: Seq.empty)
designElliptic[IIRKernel1D](...)

import breeze.signal.filter.sos
designButterworth[SOSKernel1D](normCutOff: T, order: Int, flavour: FilterFlavour = FilterFlavour.LoPass, opts: Seq[FilterOpts]: Seq.empty)
designElliptic[SOSKernel1D](...)

this way one can easily swap the filter types by just changing the imported package and the naming is consistent as well. all functions return a FilterKernel suitable for the filter() method

Delegate the actual filter design implementation to an implicit CanDesignFilter[KernelType] Maybe distinguish between CanDesignFilter and CanDesignWindowedFilter - or is it possible to determine a common set of args/parameters for fir(win), iir, sos filters ?

I would prefer having explicit functions for each filter design instead of a generic design() function with lots of parameters (better UX with code completition in my eyes :-))

[ktakagaki]

Hi Hanez, (@almoehi)

Thanks for your patience. Here are my questions/discussion...

• I've never used SOS filters, but wouldn't it just be simpler to return an IIRKernel?
• designButterworth[IIRKernel1D]... what is the purpose of specifying the output type (i.e. is it not always IIRKernel1D)? Ditto designElliptic. One issue with specifying the return type in this way is that one can't specify FIRKernel1D[Double] vs FIRKernel1D[Int] easily. Perhaps kernels other than Double are not generally needed, but for me they are useful, since I often deal with DAQ sampled data which is 12-16 bits integer, and it is sometimes faster to keep them as Int. What do you think? There would be a viable argument to ditch integers and make the whole breeze.signal package Double based only.
• flavour => flavor(American style)? I think the majority of breeze contributors are new-worldian, including me (although I've hopped the pond). Should we consider the scipy style btype or MatLab style ftype, or even just call it type instead (I kind of favor this)?

• I think this is a key point you bring up:

  or is it possible to determine a common set of args/parameters for fir(win), iir, sos filters ?

  I would prefer having explicit functions for each filter design instead of a generic design() function with lots of parameters (better UX with code completition in my eyes :-))

I think for the first half of the issue, it might be possible to an extent (see next bullet), but the second part, I see your point and agree. The filterDesign( ... ) syntax I proposed previously not only becomes an option nightmare, as you point out, it also introduces needless complexities in inferring the return type (i.e. FIRKernel1D, IIRKernel1D, FIRKernel2D, etc...). As long as we avoid half-hazard, <8 character, 20th century names like butter and cheby1, I think it would all work well. So something like this?:

designFiltButterworth(order = 10, omega = 0.2, type = Opt.LowPass)
designFiltButterworth(order = 10, omega = DenseVector(0.2, 0.3), type = Opt.BandPass)
designFiltButterworth(order = 10, omega = 10d, type = Opt.LowPass, samplingRate = 100d)

designFiltElliptical(order = 10, omega = ...., ripplePass = ..., rippleStop = ...)

Re the naming... designFiltButterworth? designButterworth? butterworth/butter? (Personally, I don't like these last two, esp. in light of IDE code completion)

• So for your point whether there is a common set of args/params, I think at least order and omega can be assumed, along with samplingRate (default 2Pi, i.e. Nyquist frequency would be omega=1)

Thoughts? Kenta

[ktakagaki]

P.S.

One challenge which I don't have any clear answer to is how to encorporate stuff like: buttord (Scipy, MatLab into the syntax.)

One way would be to do this all through the filter design function, and specify the order argument as order: OptOrder and to alloworder = (OptOrder.)Automatic. In that case, order = 10 would be implicitly converted to something like order = OptOrder.IntValue(10)

The other way would be to come up with a function like buttorder, but hopefully to make the name sound less odd/ad hoc/20th century...

[ktakagaki]

PPS. I think if we do it this way (especially implementing biquad as a subset of IIR), I think we could put all of the design functions into breeze.signal without providing breeze.signal.iir, etc.?

[dlwh]

FWIW, my general preferences for the signal api are, more or less in this order:

1) Internally consistent 2) consistent with scipy signal processing 3) consistent with breeze.linalg (which mostly just means thinking if things can be ufuncs.)

(and yes, american spelling, please.)

On Tue, May 27, 2014 at 10:11 PM, Kenta Takagaki notifications@github.comwrote:

PPS. I think if we do it this way (especially implementing biquad as a subset of IIR), I think we could put all of the design functions into breeze.signal without providing breeze.signal.iir, etc.?

— Reply to this email directly or view it on GitHubhttps://github.com/scalanlp/breeze/issues/255#issuecomment-44342234 .

[ktakagaki]

Hi David,

Re 1) and 2), I know we've discussed this a bit in the past... https://github.com/scalanlp/breeze/issues/145

Do you prefer function names like butter, cheby1, ellip for designing filters? If you really insist on this, I would be OK with it personally, I mean these basic names at least are pretty standard in MatLab/SciPy-like packages. (((I would just point out that the abbreviation scheme is pretty incoherent, and I kind of like the idea of starting to type designFilt... or something like that, and the IDE giving me all the completion options.)))

Re 3), I think UFuncs are a tough call if we want to use option case objects/classes and implicit conversions to realize something like the polymorphic argument types in MatLab/Scipy. My understanding is that double implicits will never happen in Scala (i.e. implicit conversion of arguments sent to the implicit internal UFunc implementations), so that would leave us with either (1) the old CanXXX syntax or (2) perhaps making implementation UFuncs like breeze.signal.support.canDesignButterworth extends UFunc which would be called by the non-UFunc (e.g. designButterworth), with all of its option objects, etc. I'm thinking (2) would serve no purpose, and would be overkill...???

Cheers, Kenta

[almoehi]

Hi Kenta @ktakagaki

re SOS: focus would be on IIR in first place for sure, SOS is just a personal curiosity I would like to contribute by some time

re 2) the type param was just for illustration purpose to indicate what these functions should return - should of course be generic for Doubles, Ints ...

re 3) british naming wasn't intentional - I totaly agree with american style :-)

I would also prefer specific naming instead of generic design() + option hell. I think we should put the common parameters as you outlined to be the first parameters and alll optional parameters after them - so we also have a consistent parameter ordering

re buttord stuff: just want to throw in this idea - we could use companion objects which can have additional methods like "buttord" for butterworth filters I'm no sure right now if this would work, but we could implement the designFiltButterworth as (case) object where the apply method(s) serve as design function and addition (filter specific) functions can be attached to the object - s.th. like this:

object designFiltButterworth {
  def apply(order = 10, omega = 0.2, type = Opt.LowPass)
  def buttord() ....
}

downside: this would only be possible when not using UFuncs I guess ....

but how about the filter() function, I guess this could be implemented as UFunc ?

[ktakagaki]

Hi Hanez,

call me stupid, but isn't an SOS/biquad filter just a short 2-order IIR filter? Sorry, I don't really know much about them.

Re filter, you are probably right, it can probably be made into a UFunc...

  def filter[Input, Kernel, Output](data: Input, kernel: Kernel,
                            overhang: OptOverhang = OptOverhang.PreserveLength,
                            padding: OptPadding = OptPadding.Zero)
        (implicit canFilter: CanFilter[Input, Kernel, Output]): Output =
    canFilter(data, kernel, overhang, padding)

I am kind of wedded to the overhang and padding options: there would be 3 ways to deal with them in a UFunc:

1. create separate implementations for each permutation, to allow for each default option:: Impl4[Input, Kernel, OptOverhang, OptPadding, Input], Impl3[Input, Kernel, OptOverhang, Input], Impl2[Input, Kernel, Input], ... (the things lost here would be the ability to specify options by name overhang = Full and the convenience of implicits padding = 0 => `padding = OptPadding.IntValue(0)
2. just get rid of these overhang and padding options (which I think would be unfortunate)
3. incorporate the padding options into the Kernel specification, so that the filter UFunc would just read the options from the kernel object.

I find your idea of dealing with butterord-parity really smart in an OOB-Scala sort of way. I would just change the name to order:

object designFiltButterworth {
  def apply(order = 10, omega = 0.2, type = Opt.LowPass): IIRKernel1D = {...}
  def order(omega: Double or (Double, Double),  attenuationPass: Double, attenuationStop: Double, type: OptType): Int = { .... }
}

Lets see what @dlwh says, he's the final arbiter of how breeze should be. Then, I think I'll set up a "Filtering data with breeze.signal" Wiki page, so that we can finalize the filter design syntax stuff before we dive in. (I wanted to expand the FIR filter design too, beyond firwin). The cheat sheets with Matlab/Scipy comparisons are all fine, but the table format isn't really working for breeze.signal.

[dlwh]

It's possible that making things UFuncs won't work. I ran into this when I was making the optimize main function in the optimize package. I like the def filter that @ktakagaki proposed. How bad would it be to have, e.g. butterworth(padding=3, ...) as the kernel argument like you did with omega etc.? Does it make sense to say that the padding and overhang is a property of the kernel, or is it something different?

Remember you can't use type as a name in Scala. You can use type, but... tpe is the common variant (on analogy to using clss.)

-- David

[almoehi]

Hi, @ktakagaki yes, that's probably what SOS is - haven't had a detailed look too yet but came across it in some paper I read a while ago - no priority, could just be a bullet point on the feature list :-)

I guess moving all the padding / overhang options into the kernel is a good idea - that way the filter() method is straight forward: data + kernel => filtered data

If we use the scala OOB way we could provide several apply() methods to deal with several optional parameters - this can be considered like a factory pattern for creating filter kernels

I think this could make the filter subsystem much more cleaner and easier to use + filter() can be designed as a simple UFunc

cheers!

[ktakagaki]

Hi All,

So I've summarized what we discussed into a Wiki specification, please provide input: https://github.com/ktakagaki/breeze/wiki/%5Bbreeze.signal%5D-Designing-Filters

1. I see the merit in making filter a UFunc and pushing the padding and overhang options into Kernel via filter design arguments.
2. Question: maybe we should push the OOB a bit further, and change the names designFiltWindowed, designFiltButterworth, etc... to kernelWindowed, kernelButterworth, etc...???

[almoehi]

Thanks for the summarization ! Looks good to me so far.

reg 2.) I thought about this as well - but if I were a first-time user and want to create a filter, I would probably try to type something like "filter..." or "createFilter..." or "design..." and see what the IDE code completion comes up with - don't know if I would ever try to type "kernel...." to see waht comes up :-)

[dlwh]

Where are we on this?

[almoehi]

sry, I got kind of "interrupted" by my thesis deadline - I think we're just missing mutual agreement on which naming convention we'd like to use. After that we could start with the refactoring as discussed above (filter -> ufunc, moving parameters to kernels ...) I'm currently at siggraph but will have a look at once I'm back to see where I can start - not sure what @ktakagaki workload looks like right now ?

cheers!

[dlwh]

no apologies necessary! I was just checking in on some old tickets. Hope your thesis went well!

On Mon, Aug 11, 2014 at 11:14 PM, almoehi notifications@github.com wrote:

sry, I got kind of "interrupted" by my thesis deadline - I think we're just missing mutual agreement on which naming conection we'd like to use. After that we could start with the refactoring as discussed above (filter -> ufunc, mogin parameters to kernels ...) I'm currently at siggraph but will have a look at once I'm back to see where I can start - not sure what @ktakagaki https://github.com/ktakagaki workload looks like right now.

cheers!

— Reply to this email directly or view it on GitHub https://github.com/scalanlp/breeze/issues/255#issuecomment-51877316.

[ktakagaki]

Hi All,

Before I took a hiatus, I rethought and now believe that the "padding" and "overhang" options should really be specified within the filter function (as is currently), instead of incorporating into the filter design functions. The main reason is that padding/overhang will affect the output of the filter differently, depending upon the filter method. For example, if you set the overhang option to complete and you choose a bidirectional convolution for your filter method, then there will be two full filter windows worth of overhang. Therefore, I think these options really belong conceptually with the filter function, instead of the kernel design functions. This would mean that we can't use UFuncs meaningfully for this.

I have adjusted the summary document accordingly https://github.com/ktakagaki/breeze/wiki/%5Bbreeze.signal%5D-Designing-Filters

The other question is regarding the object-based filter design syntax that @almoehi proposed. I was considering something like

FilterButterworth.design(order = 10, omega = 0.2, tpe = Opt.LowPass): IIRKernel1D

as a slight modification of @almoehi 's initial proposal of

designFilterButterworth(order = 10, omega = 0.2, tpe = Opt.LowPass): IIRKernel1D

Please let me know what you think

Sincerely, Kenta

[dlwh]

Isn't that an argument for putting it in the design? If the behavior of padding etc. changes based on what kind of filter you're using, shouldn't you be specifying it when you create the filter?

Or maybe I misunderstand.

On Thu, Aug 14, 2014 at 5:28 AM, Kenta Takagaki notifications@github.com wrote:

Hi All,

Before I took a hiatus, I rethought and now believe that the "padding" and "overhang" options should be specified within the filter function, as is, instead of incorporating into the filter design functions. The main reason is that padding/overhang will affect the output of the filter differently, depending upon the filter method. For example, if you set the overhang option to complete and you choose a bidirectional convolution for your filter method, then there will be two full filter windows worth of overhang. Therefore, I think these options really belong conceptually with the filter function, instead of the kernel design functions.

I have adjusted the summary document accordingly

https://github.com/ktakagaki/breeze/wiki/%5Bbreeze.signal%5D-Designing-Filters

The other question is regarding the object-based filter design syntax that @almoehi https://github.com/almoehi proposed. I was considering something like

FilterButterworth.design(order = 10, omega = 0.2, tpe = Opt.LowPass): IIRKernel1D

as a slight modification of @almoehi https://github.com/almoehi 's initial proposal of

designFilterButterworth(order = 10, omega = 0.2, tpe = Opt.LowPass): IIRKernel1D

Please let me know what you think

Sincerely, Kenta

— Reply to this email directly or view it on GitHub https://github.com/scalanlp/breeze/issues/255#issuecomment-52176800.

[ktakagaki]

I realize that there are arguments either way, and I'm somewhat fine either way as long as it's not a mess, but I think I didn't convey my argument clearly.

I think that the same FIRKernel1D or IIRKernel1D object, even if it included encapsulated fixed overhang and padding options, would give different output lengths depending upon whether you filtered data with the equivalent of, say filter or filt or filtfilt, or whatever... Put in other words, the same "full overhang" option would mean different things to say a filt equivalent vs a filtfilt equivalent, etc. Given the cyclical nature of fftfilt, for it, the padding would also probably mean something different. Therefore, I thought that these options should probably be included in the verb part, i.e. the filtering function itself, and not the upstream design functions, so that it is clearer at the time of use what the behavior of the filtering would be.

[dlwh]

I see, ok. I'll go with your judgment on this one.

On Sun, Aug 17, 2014 at 6:50 AM, Kenta Takagaki notifications@github.com wrote:

I realize that there are arguments either way, and I'm somewhat fine either way as long as it's not a mess, but I think I didn't convey my argument clearly.

I think that the same FIRKernel1D or IIRKernel1D object, even if it included encapsulated fixed overhang and padding options, would give different output lengths depending upon whether you filtered data with the equivalent of, say filter or filt or filtfilt, or whatever... Put in other words, the same "full overhang" option would mean different things to say a filt equivalent vs a filtfilt equivalent, etc. Given the cyclical nature of fftfilt, for it, the padding would also probably mean something different. Therefore, I thought that these options should probably be included in the verb part, i.e. the filtering function itself, and not the upstream design functions, so that it is clearer at the time of use what the behavior of the filtering would be.

— Reply to this email directly or view it on GitHub https://github.com/scalanlp/breeze/issues/255#issuecomment-52422360.

[almoehi]

I commited a quick draft of the refactoring of the filter design stuff. For now, I moved everything into a breeze.signal.filter subpackage to make it easer to review, except for the "Opt" objects which I added to the existing options.scala please review them here: https://github.com/almoehi/breeze/tree/master/math/src/main/scala/breeze/signal/filter

With the refactoring in place I think we can remove all of the filterLP/HP/BP functions from the breeze.signal package object - right ?

I propse to move all the filter stuff and support into breeze.signal.filter and breeze.signal.filter.support sub packages. I expect the breeze.signal package to grow over time and it will eventually end up in a big package object "blob" or something.

And we are keeping the filter function as is and do not refactor it to a UFunc right ?

[ktakagaki]

Hi Hannes, thanks for your efforts, I've certainly dropped the ball on this. I quickly made the code compile and the tests pass, so that we can start coding. I posted it as an incomplete pull request, so maybe we can switch the discussion thread to there. https://github.com/scalanlp/breeze/pull/317 Kenta