diff --git a/src/Filters/coefficients.jl b/src/Filters/coefficients.jl index e6fbe004..dfeac420 100644 --- a/src/Filters/coefficients.jl +++ b/src/Filters/coefficients.jl @@ -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) @@ -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 # @@ -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] @@ -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[] @@ -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)) @@ -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 diff --git a/test/filter_conversion.jl b/test/filter_conversion.jl index 3471c6b7..10535a66 100644 --- a/test/filter_conversion.jl +++ b/test/filter_conversion.jl @@ -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...) @@ -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