Skip to content

Commit

Permalink
Remove recursion from Base.:^ (#618)
Browse files Browse the repository at this point in the history 
Co-authored-by: Martin Holters <[email protected]>
  • Loading branch information
@wheeheee @martinholters
wheeheee and martinholters authored Jan 16, 2025
1 parent 1f70965 commit 1b76b87
Show file tree
Hide file tree
Showing 2 changed files with 71 additions and 52 deletions.
102 changes: 51 additions & 51 deletions src/Filters/coefficients.jl
View file
Original file line number Diff line number Diff line change
@@ -1,5 +1,7 @@
# Filter types and conversions

using Base: uabs

abstract type FilterCoefficients{Domain} end

Base.convert(::Type{T}, f::FilterCoefficients) where {T<:FilterCoefficients} = T(f)
Expand Down Expand Up @@ -45,11 +47,9 @@ Base.promote_rule(::Type{ZeroPoleGain{D,Z1,P1,K1}}, ::Type{ZeroPoleGain{D,Z2,P2,
Base.inv(f::ZeroPoleGain{D}) where {D} = ZeroPoleGain{D}(f.p, f.z, inv(f.k))

function Base.:^(f::ZeroPoleGain{D}, e::Integer) where {D}
if e < 0
return inv(f^-e)
else
return ZeroPoleGain{D}(repeat(f.z, e), repeat(f.p, e), f.k^e)
end
ae = uabs(e)
z, p = repeat(f.z, ae), repeat(f.p, ae)
return e < 0 ? ZeroPoleGain{D}(p, z, inv(f.k)^ae) : ZeroPoleGain{D}(z, p, f.k^ae)
end

#
Expand Down Expand Up @@ -184,11 +184,12 @@ end
Base.inv(f::PolynomialRatio{D}) where {D} = PolynomialRatio{D}(f.a, f.b)

function Base.:^(f::PolynomialRatio{D,T}, e::Integer) where {D,T}
ae = uabs(e)
b, a = f.b^ae, f.a^ae
if e < 0
return PolynomialRatio{D}(f.a^-e, f.b^-e)
else
return PolynomialRatio{D}(f.b^e, f.a^e)
b, a = a, b
end
return PolynomialRatio{D}(b, a)
end

coef_s(p::LaurentPolynomial) = p[end:-1:0]
Expand Down Expand Up @@ -304,7 +305,49 @@ Base.promote_rule(::Type{SecondOrderSections{D,T1,G1}}, ::Type{SecondOrderSectio

SecondOrderSections{D,T,G}(f::SecondOrderSections) where {D,T,G} =
SecondOrderSections{D,T,G}(f.biquads, f.g)
SecondOrderSections{D,T,G}(f::Biquad{D}) where {D,T,G} = SecondOrderSections{D,T,G}([f], one(G))

SecondOrderSections{D}(f::SecondOrderSections{D,T,G}) where {D,T,G} = SecondOrderSections{D,T,G}(f)
SecondOrderSections{D}(f::Biquad{D,T}) where {D,T} = SecondOrderSections{D,T,Int}(f)

*(f::SecondOrderSections{D}, g::Number) where {D} = SecondOrderSections{D}(f.biquads, f.g * g)
*(g::Number, f::SecondOrderSections{D}) where {D} = SecondOrderSections{D}(f.biquads, f.g * g)
*(f1::SecondOrderSections{D}, f2::SecondOrderSections{D}) where {D} =
SecondOrderSections{D}([f1.biquads; f2.biquads], f1.g * f2.g)
*(f1::SecondOrderSections{D}, fs::SecondOrderSections{D}...) where {D} =
SecondOrderSections{D}(vcat(f1.biquads, map(f -> f.biquads, fs)...), f1.g * prod(f.g for f in fs))

*(f1::Biquad{D}, f2::Biquad{D}) where {D} = SecondOrderSections{D}([f1, f2], 1)
*(f1::Biquad{D}, fs::Biquad{D}...) where {D} = SecondOrderSections{D}([f1, fs...], 1)
*(f1::SecondOrderSections{D}, f2::Biquad{D}) where {D} =
SecondOrderSections{D}([f1.biquads; f2], f1.g)
*(f1::Biquad{D}, f2::SecondOrderSections{D}) where {D} =
SecondOrderSections{D}([f1; f2.biquads], f2.g)

Base.inv(f::SecondOrderSections{D}) where {D} = SecondOrderSections{D}(inv.(f.biquads), inv(f.g))

function Base.:^(f::SecondOrderSections{D}, e::Integer) where {D}
ae = uabs(e)
if e < 0
inv_f = inv(f)
return SecondOrderSections{D}(repeat(inv_f.biquads, ae), inv_f.g^ae)
else
return SecondOrderSections{D}(repeat(f.biquads, ae), f.g^ae)
end
end

function Base.:^(f::Biquad{D}, e::Integer) where {D}
ae = uabs(e)
return SecondOrderSections{D}(fill(e < 0 ? inv(f) : f, ae), 1)
end

function Biquad{D,T}(f::SecondOrderSections{D}) where {D,T}
if length(f.biquads) != 1
throw(ArgumentError("only a single second order section may be converted to a biquad"))
end
Biquad{D,T}(f.biquads[1] * f.g)
end
Biquad{D}(f::SecondOrderSections{D,T,G}) where {D,T,G} = Biquad{D,promote_type(T, G)}(f)

function ZeroPoleGain{D,Z,P,K}(f::SecondOrderSections{D}) where {D,Z,P,K}
z = Z[]
Expand All @@ -321,14 +364,6 @@ end
ZeroPoleGain{D}(f::SecondOrderSections{D,T,G}) where {D,T,G} =
ZeroPoleGain{D,complex(T),complex(T),G}(f)

function Biquad{D,T}(f::SecondOrderSections{D}) where {D,T}
if length(f.biquads) != 1
throw(ArgumentError("only a single second order section may be converted to a biquad"))
end
Biquad{D,T}(f.biquads[1] * f.g)
end
Biquad{D}(f::SecondOrderSections{D,T,G}) where {D,T,G} = Biquad{D,promote_type(T,G)}(f)

PolynomialRatio{D,T}(f::SecondOrderSections{D}) where {D,T} = PolynomialRatio{D,T}(ZeroPoleGain(f))
PolynomialRatio{D}(f::SecondOrderSections{D}) where {D} = PolynomialRatio{D}(ZeroPoleGain(f))

Expand Down Expand Up @@ -447,39 +482,4 @@ end
SecondOrderSections{D}(f::ZeroPoleGain{D,Z,P,K}) where {D,Z,P,K} =
SecondOrderSections{D,promote_type(real(Z), real(P)), K}(f)


SecondOrderSections{D,T,G}(f::Biquad{D}) where {D,T,G} = SecondOrderSections{D,T,G}([f], one(G))
SecondOrderSections{D}(f::Biquad{D,T}) where {D,T} = SecondOrderSections{D,T,Int}(f)
SecondOrderSections{D}(f::FilterCoefficients{D}) where {D} = SecondOrderSections{D}(ZeroPoleGain(f))

*(f::SecondOrderSections{D}, g::Number) where {D} = SecondOrderSections{D}(f.biquads, f.g * g)
*(g::Number, f::SecondOrderSections{D}) where {D} = SecondOrderSections{D}(f.biquads, f.g * g)
*(f1::SecondOrderSections{D}, f2::SecondOrderSections{D}) where {D} =
SecondOrderSections{D}([f1.biquads; f2.biquads], f1.g * f2.g)
*(f1::SecondOrderSections{D}, fs::SecondOrderSections{D}...) where {D} =
SecondOrderSections{D}(vcat(f1.biquads, map(f -> f.biquads, fs)...), f1.g * prod(f.g for f in fs))

*(f1::Biquad{D}, f2::Biquad{D}) where {D} = SecondOrderSections{D}([f1, f2], 1)
*(f1::Biquad{D}, fs::Biquad{D}...) where {D} = SecondOrderSections{D}([f1, fs...], 1)
*(f1::SecondOrderSections{D}, f2::Biquad{D}) where {D} =
SecondOrderSections{D}([f1.biquads; f2], f1.g)
*(f1::Biquad{D}, f2::SecondOrderSections{D}) where {D} =
SecondOrderSections{D}([f1; f2.biquads], f2.g)

Base.inv(f::SecondOrderSections{D}) where {D} = SecondOrderSections{D}(inv.(f.biquads), inv(f.g))

function Base.:^(f::SecondOrderSections{D}, e::Integer) where {D}
if e < 0
return inv(f)^-e
else
return SecondOrderSections{D}(repeat(f.biquads, e), f.g^e)
end
end

function Base.:^(f::Biquad{D}, e::Integer) where {D}
if e < 0
return inv(f)^-e
else
return SecondOrderSections{D}(fill(f, e), 1)
end
end
21 changes: 20 additions & 1 deletion test/filter_conversion.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -219,7 +219,7 @@ end
p = rand(ComplexF64, Npc) .- (0.5 + 0.5im);
p = [p; conj(p); rand(Npr) .- 0.5; zeros(max(2Nzc+Nzr-2Npc-Npr, 0))]
H′ = ZeroPoleGain(z, p,
(rand() + 0.5) * rand([-1, 1]), # non-zero gain with random sign
(rand() + 0.5) * rand((-1, 1)), # non-zero gain with random sign
)
maybe_biquad = length(z) ≤ 2 && length(p) ≤ 2 ? (Biquad,) : ()
for T ∈ (PolynomialRatio, ZeroPoleGain, SecondOrderSections, maybe_biquad...)
Expand All @@ -236,6 +236,25 @@ end
@test filt(H^0, x) ≈ x
end
end
# test that ^ doesn't cause stack overflow
let H = PolynomialRatio([1.0], [2.0])^typemin(Int8)
@test coefb(H) == [2.0^128]
@test coefa(H) == [1.0]
end
let zpg = ZeroPoleGain([1], [2], 3)^typemin(Int8)
@test length(zpg.z) == length(zpg.p) == 128
@test all(==(2), zpg.z)
@test all(==(1), zpg.p)
@test zpg.k ≈ inv(3)^128
end
let bq = Biquad(1:5...)
sos1 = bq^typemin(Int8)
sos2 = SecondOrderSections(bq)^typemin(Int8)
@test length(sos1.biquads) == length(sos2.biquads) == 128
@test all(==(inv(bq)), sos1.biquads)
@test all(==(inv(bq)), sos2.biquads)
@test sos1.g == sos2.g == 1
end
end

@testset "types" begin
Expand Down

0 comments on commit 1b76b87

Please sign in to comment.