From 6a904b3b3efed5557f4c4b6391c931a23c807514 Mon Sep 17 00:00:00 2001
From: Martin Holters <martin.holters@hsu-hh.de>
Date: Mon, 21 Jun 2021 17:45:27 +0200
Subject: [PATCH] Make `SecondOrderSections(::ZeroPoleGain)`more robust (#435)

* Make `SecondOrderSections(::ZeroPoleGain)`more robust

The helper function `groupzp` assumes that both zeros and poles are
ordered such that complex entries appear in conjugate pairs
consecutively and further, that these pairs are orderer consistently:
Either always the positive imaginary part first, or always the negative
imaginary part first. The output of `split_real_complex` fulfills this.
However, for the poles, there is an intermediate sorting step which
apparently can interfere with this ordering under certain circumstances.
This can lead to a `BoundsError` (#432), but also to zeros which are not
conjugate pairs being grouped together, thereby silently producing wrong
results. While #433 was a band-aid preventing the `BoundsError`, this
should fix the underlying issue by integrating the sorting into
`split_real_complex` in such a way that the required ordering is
ensured.

* Bump version to 0.7.2
---
 Project.toml                |  2 +-
 src/Filters/coefficients.jl | 17 +++++++----------
 test/filter_conversion.jl   |  8 +++++---
 3 files changed, 13 insertions(+), 14 deletions(-)

diff --git a/Project.toml b/Project.toml
index 0917a287e..7d1c4fe96 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,6 +1,6 @@
 name = "DSP"
 uuid = "717857b8-e6f2-59f4-9121-6e50c889abd2"
-version = "0.7.1"
+version = "0.7.2"
 
 [deps]
 Compat = "34da2185-b29b-5c13-b0c7-acf172513d20"
diff --git a/src/Filters/coefficients.jl b/src/Filters/coefficients.jl
index 016f53083..7aef73bf0 100644
--- a/src/Filters/coefficients.jl
+++ b/src/Filters/coefficients.jl
@@ -322,9 +322,6 @@ function groupzp(z, p)
         groupedz[i] = splice!(z, closest_zero_idx)
         if !isreal(groupedz[i])
             i += 1
-            if closest_zero_idx > length(z)
-                closest_zero_idx = length(z)
-            end
             groupedz[i] = splice!(z, closest_zero_idx)
         end
         i += 1
@@ -337,7 +334,7 @@ end
 # Sort zeros or poles lexicographically (so that poles are adjacent to
 # their conjugates). Handle repeated values. Split real and complex
 # values into separate vectors. Ensure that each value has a conjugate.
-function split_real_complex(x::Vector{T}) where T
+function split_real_complex(x::Vector{T}; sortby=nothing) where T
     # Get counts and store in a Dict
     d = Dict{T,Int}()
     for v in x
@@ -349,7 +346,11 @@ function split_real_complex(x::Vector{T}) where T
 
     c = T[]
     r = typeof(real(zero(T)))[]
-    for k in keys(d)
+    ks = collect(keys(d))
+    if sortby !== nothing
+        sort!(ks, by=sortby)
+    end
+    for k in ks
         if imag(k) != 0
             if !haskey(d, conj(k)) || d[k] != d[conj(k)]
                 # No match for conjugate
@@ -381,13 +382,9 @@ function SecondOrderSections{D,T,G}(f::ZeroPoleGain{D,Z,P}) where {D,T,G,Z,P}
     # Split real and complex poles
     (complexz, realz, matched) = split_real_complex(z)
     matched || throw(ArgumentError("complex zeros could not be matched to their conjugates"))
-    (complexp, realp, matched) = split_real_complex(p)
+    (complexp, realp, matched) = split_real_complex(p; sortby=x->abs(abs(x) - 1))
     matched || throw(ArgumentError("complex poles could not be matched to their conjugates"))
 
-    # Sort poles according to distance to unit circle (nearest first)
-    sort!(complexp, by=x->abs(abs(x) - 1))
-    sort!(realp, by=x->abs(abs(x) - 1))
-
     # Group complex poles with closest complex zeros
     z1, p1 = groupzp(complexz, complexp)
     # Group real poles with remaining complex zeros
diff --git a/test/filter_conversion.jl b/test/filter_conversion.jl
index af56994cc..5b7dc370a 100644
--- a/test/filter_conversion.jl
+++ b/test/filter_conversion.jl
@@ -351,7 +351,9 @@ end
     Fsamp = 180
     responsetype = Bandpass(bp1, bp2; fs = Fsamp)
     designmethod = Elliptic(11, 0.25, 40)
-    bpass = digitalfilter(responsetype, designmethod)
-    H = SecondOrderSections(bpass) # this shouldn't throw
-    @test H isa SecondOrderSections
+    H = digitalfilter(responsetype, designmethod)
+    H′ = ZeroPoleGain(SecondOrderSections(H))
+    @test sort(H.p, by=z->(real(z), imag(z))) ≈ sort(H′.p, by=z->(real(z), imag(z)))
+    @test sort(H.z, by=z->(real(z), imag(z))) ≈ sort(H′.z, by=z->(real(z), imag(z)))
+    @test H.k ≈ H′.k
 end