Accuracy of ZeroPoleGain to PolynomialRatio conversion

[martinholters]

Consider an FIR filter converted to ZPK and back:

julia> Hpf = PolynomialRatio(remez(151, [0, 0.2, 0.25, 0.5], [1, 0]; weight=[1, 2000]), 1);

julia> Hzpk = ZeroPoleGain(Hpf);

julia> Hpf2 = PolynomialRatio(Hzpk);

One may wonder whether that is numerically sound, and it turns out the first conversion is, the second is not:

julia> freqresp(Hpf)[1] ≈ freqresp(Hzpk)[1]
true

julia> freqresp(Hpf2)[1] ≈ freqresp(Hzpk)[1]
false

I.e. the frequency response of the ZPK form is ok, but going back to the polynomial is problematic. And spectacularly so:

julia> extrema(coefb(Hpf))
(-0.08912595567110218, 0.44057088012179)

julia> extrema(coefb(Hpf2))
(-9.006618714365288e13, 8.865389303639448e13)

However, as the frequency response is ok, we can also transform that back to time domain to get the impulse response which at the same time gives the polynomial coefficients. And indeed, that works:

julia> irfft(freqresp(Hzpk, range(0, 2pi, length=length(Hzpk.z)+2)[1:(end+1)÷2]), length(Hzpk.z)+1) ≈ coefb(Hpf)
true

For the IIR case, numerator and denominator could probably be transformed individually in a similar fashion. Does anyone have an insight to share whether that is more generally numerically better or have I just hit an example where that is the case by accident? In any case, it would be nice if PolynomialRatio(::ZeroPoleGain) worker properly in cases like this, where obviously, it could.

[martinholters]

Note that PolynomialRatio(::ZeroPoleGain) currently invokes Polynomials.fromroots, which basically does repeated convolution with [1, -r] for all given roots r. If there is a better strategy, Polynomials.jl might be the better place to implement it than DSP.jl. CC @jverzani although I'm not sure this is actionable outside the (in that respect) somewhat limited scope of DSP.jl.

[jverzani]

If a better algorithm exists, I'm open to including it in Polynomials.jl.

[martinholters]

At least in this case, the order of the roots as returned by roots seems pathologically bad for fromroots:

julia> p = Polynomial(remez(151, [0, 0.2, 0.25, 0.5], [1, 0]; weight=[1, 2000]));

julia> fromroots(roots(p)) * p[end] ≈ p
false

julia> fromroots(shuffle(roots(p))) * p[end] ≈ p
true

julia> fromroots(shuffle(roots(p))) * p[end] ≈ p
true

julia> fromroots(shuffle(roots(p))) * p[end] ≈ p
false

julia> fromroots(shuffle(roots(p))) * p[end] ≈ p
true

julia> fromroots(shuffle(roots(p))) * p[end] ≈ p
true

Some basic monte-carlo testing indicates that fromroots(shuffle(roots(p))) * p[end] ≈ p holds with about 95% probability for this example.

This makes me think that there might be a better reduction strategy (in the sense that fromroot does a reduce(conv, ...) like operation) than the foldl that is currently performed (although not spelled like that). Similar to how pairwise summation is superior to foldl(+, x). But I'm way out of my comfort zone here.

[martinholters]

Some googling later, it appears that Leja Ordering is the answer. Best reference seems to be https://link.springer.com/article/10.1023/A:1025555803588 - unfortunately, behind a paywall, but I have access through my institution. If I can make anything out of it, I'll open a PR a Polynomials.jl.

[jverzani]

I can have a look. Our institution has access. Thanks! I'll let you know if I get anywhere.

[martinholters]

FWIW, I've tried a quick-and-dirty implementation and it did solve the problems reported above. I also have an idea how the mentioned O(n^2) algorithm might work. I'll keep you updated.

[jverzani]

Great. I'd love a PR if you get to that point.

[martinholters]

JuliaMath/Polynomials.jl#586 should fix this, so once there is a new release of Polynomials.jl, this should be resolved.

[jverzani]

I just initiated a tag for Polynomials with Martin's fix.

[martinholters]

Fixed with Polynomials.jl v4.0.12.