Inconsistencies for filters with different number of poles and zeros

[martinholters]

I think for discrete filters, it is unclear whether powers of z or powers of z⁻¹ are being used. If numerator and denominator have the same order, that doesn't matter, but otherwise this gives rise to some inconsistencies:

julia> h = Biquad(0., 1., 0., 0., 0.); # just a unit delay z⁻¹

julia> impz(h, 3) # OK
3-element Array{Float64,1}:
 0.0
 1.0
 0.0

julia> h2 = PolynomialRatio(h)
DSP.Filters.PolynomialRatio{Float64}(Poly(1.0), Poly(1.0*x))

julia> impz(h2, 3) # Oh oh
3-element Array{Float64,1}:
 1.0
 0.0
 0.0

julia> freqz(h, linspace(0, pi)) ≈ freqz(h2, linspace(0, pi)) ≈ exp.(linspace(0, pi)*im) # here it's treated like z instead of z⁻¹
true

Looks like the last one to make any major changes to the involved functions was @simonster three years ago. Any opinions here?

[simonster]

I think that freqz should treat the powers as powers of z⁻¹, and there should be an argument that controls this behavior for PolynomialRatio. I'm open to other ideas, though.

[martinholters]

I think that freqz should treat the powers as powers of z⁻¹

Yes, I agree that's the most reasonable approach, and it looks like that's the case:

julia> h3=PolynomialRatio([1.0, 0.0], [1.0])
DSP.Filters.PolynomialRatio{Float64}(Poly(1.0*x), Poly(1.0))

julia> freqz(h3, linspace(0, pi)) ≈ exp.(-linspace(0, pi)*im)
true

As noted in #135, the order of the coefficients has to be given in reversed order, but then, yes, if x corresponds to z⁻¹, h3 indeed describes a unit delay, and freqz agrees. But:

julia> impz(h3, 3)
3-element Array{Float64,1}:
 1.0
 0.0
 0.0

So as impz uses filt internally, filtering seems broken - there is no reasonable interpretation under which h2 from above and h3 should produce the same impulse response. In fact, I haven't managed to find a PolynomialRatio that makes impz produce a delayed-by-one impulse.

Anyway, before fixing specific issues, it makes sense to straighten out what different representations should mean. I'll give it a try:

constructor	z	Laplace
PolynomialRatio([b0, ..., bN], [a0, ..., aM])	(b0 + ... + bN z⁻ᴺ) / (a0 + ... + aM z⁻ᴹ)	(b0 sᴺ + ... + bN) / (a0 sᴹ + ... + aM)
Biquad(b0, b1, b2, a0, a1, a2)	(b0 + b1z⁻¹ + bN z⁻²) / (a0 + a1z⁻¹ + aM z⁻²)	(b0 s² + b1 s + bN) / (a0 s² + a1 s + aM)
ZeroPoleGain([z1, ..., zN], [p1, ..., pN], g)	g (z-z1)...(z-zN) / ((z-p1)...(z-pN))	g (s-z1)...(s-zN) / ((s-p1)...(s-pN))

Assuming we agree on this, or at least something relatively similar, it is apparent that we cannot unify the z and Laplace columns by introducing an x which is an expression in z or in s, depending on use case. So it looks like the user will have to specify whether he is giving z or Laplace transform parameters. Of course another parameter to the constructors would do, but my current thinking is that a type parameter might be more useful. Especially the conversion between different representations would be more straight forward then.

I'm happy to implement this scheme, but would like to get some feedback first whether it would be accepted. Any others we should CC to make them aware of this discussion?

[martinholters]

And whatever we do, it would probably help a lot if we add customized show methods that leave no doubts. Something like

julia> PolynomialRatio{:z}([3.0], [1, 0.5])
DSP.Filters.PolynomialRatio{:z,Float64}:
     3.0
–––––––––––––
1.0 + 0.5z^-1

julia> PolynomialRatio{:s}([3.0], [1, 0.5])
DSP.Filters.PolynomialRatio{:s,Float64}:
   3.0
––––––––––
1.0s + 0.5

[hamzaelsaawy]

Is there any consensus on this? (Asking in relation to #255)

[martinholters]

If we want to handle z and s domains differently (and I think we do), I still believe we should distinguish this in the type of PolynomialRatio. Obviously, that would be a breaking change, but I don't see any other good options. I had started on an implementation, but lacked motivation (in part due to this issue seemingly not bothering anyone).

[hamzaelsaawy]

Is there a rationale as to changing just PolynomialRatio vs FilterCoefficients? I feel if were going to break the interface, we should go big with the changes, especially if augmenting all filters with domain information standardizes the interface, and allow printing filter coefficients with either z or s... For example, all the IIR design methods return s domain, we could automatically have a bilinear() call if we filter with it since this package only does z.

We could also include sampling rate if that is of value, but I have a feeling it will only complicate things.

I was thinking roughly:

abstract type DomainType end
struct DiscreteDomain{Tsample} <: DomainType end
struct AnalogDomain <: DomainType end

abstract type FilterCoefficients{D<:DomainType} end

struct PolynomialRatio{D<:DomainType, T<:Number} <: FilterCoefficients{D} 
    #...
end

struct ZeroPoleGain{D<:DomainType, Z<:Number, P<:Number, K<:Number} <: FilterCoefficients{D}
    # ...
end

[martinholters]

No rationale. Ob the contrary, changing FilterCoefficients, too, was exactly what I had started doing.