Merge pull request #108 from JuliaDSP/sjk/filter-compose

Allow scaling and composing filter coefficients using *

Author: simonster

src/Filters/coefficients.jl

@@ -2,6 +2,11 @@
 
 abstract FilterCoefficients
 
+realtype(x::DataType) = x
+realtype{T}(::Type{Complex{T}}) = T
+complextype(T::DataType) = Complex{T}
+complextype{T}(::Type{Complex{T}}) = Complex{T}
+
 #
 # Zero-pole gain form
 #
@@ -12,6 +17,19 @@ immutable ZeroPoleGain{Z<:Number,P<:Number,K<:Number} <: FilterCoefficients
     k::K
 end
 
+Base.promote_rule{Z1,P1,K1,Z2,P2,K2}(::Type{ZeroPoleGain{Z1,P1,K1}}, ::Type{ZeroPoleGain{Z2,P2,K2}}) =
+    ZeroPoleGain{promote_type(Z1,Z2),promote_type(P1,P2),promote_type(K1,K2)}
+Base.convert{Z,P,K}(::Type{ZeroPoleGain{Z,P,K}}, f::ZeroPoleGain{Z,P,K}) = f
+Base.convert{Z,P,K}(::Type{ZeroPoleGain{Z,P,K}}, f::ZeroPoleGain) =
+    ZeroPoleGain{Z,P,K}(f.z, f.p, f.k)
+
+*(f::ZeroPoleGain, g::Number) = ZeroPoleGain(f.z, f.p, f.k*g)
+*(g::Number, f::ZeroPoleGain) = ZeroPoleGain(f.z, f.p, f.k*g)
+*(f1::ZeroPoleGain, f2::ZeroPoleGain) =
+    ZeroPoleGain([f1.z; f2.z], [f1.p; f2.p], f1.k*f2.k)
+*(f1::ZeroPoleGain, fs::ZeroPoleGain...) =
+    ZeroPoleGain(vcat(f1.z, [f.z for f in fs]...), vcat(f1.p, [f.p for f in fs]...), f1.k*prod([f.k for f in fs]))
+
 #
 # Transfer function form
 #
@@ -34,77 +52,100 @@ function PolynomialRatio{T<:Number,S<:Number}(b::Union(T,Vector{T}), a::Union(S,
     PolynomialRatio{promote_type(T,S)}(Poly(b[end:-1:findfirst(b)]), Poly(a[end:-1:findfirst(a)]))
 end
 
-function Base.convert(::Type{PolynomialRatio}, f::ZeroPoleGain)
+Base.promote_rule{T,S}(::Type{PolynomialRatio{T}}, ::Type{PolynomialRatio{S}}) = PolynomialRatio{promote_type(T,S)}
+Base.convert{T}(::Type{PolynomialRatio{T}}, f::PolynomialRatio{T}) = f
+Base.convert{T}(::Type{PolynomialRatio{T}}, f::PolynomialRatio) = PolynomialRatio{T}(f.b, f.a)
+
+function Base.convert{T<:Real}(::Type{PolynomialRatio{T}}, f::ZeroPoleGain)
     b = f.k*poly(f.z)
     a = poly(f.p)
-    PolynomialRatio(Poly(real(b.a)), Poly(real(a.a)))
+    PolynomialRatio{T}(Poly(real(b.a)), Poly(real(a.a)))
 end
+Base.convert{Z,P,K}(::Type{PolynomialRatio}, f::ZeroPoleGain{Z,P,K}) =
+    convert(PolynomialRatio{promote_type(realtype(Z),realtype(P),K)}, f)
 
-function Base.convert{T}(::Type{ZeroPoleGain}, f::PolynomialRatio{T})
+function Base.convert{Z,P,K}(::Type{ZeroPoleGain{Z,P,K}}, f::PolynomialRatio)
     k = real(f.b[end])
-    b = f.b / k
-    z = convert(Vector{Complex{T}}, roots(b))
-    p = convert(Vector{Complex{T}}, roots(f.a))
-    ZeroPoleGain(z, p, k)
+    ZeroPoleGain{Z,P,K}(roots(f.b / k), roots(f.a), k)
 end
+Base.convert{T}(::Type{ZeroPoleGain}, f::PolynomialRatio{T}) =
+    convert(ZeroPoleGain{complextype(T),complextype(T),T}, f)
+
+*(f::PolynomialRatio, g::Number) = PolynomialRatio(g*f.b, f.a)
+*(g::Number, f::PolynomialRatio) = PolynomialRatio(g*f.b, f.a)
+*(f1::PolynomialRatio, f2::PolynomialRatio) =
+    PolynomialRatio(f1.b*f2.b, f1.a*f2.a)
+*(f1::PolynomialRatio, fs::PolynomialRatio...) =
+    PolynomialRatio(f1.b*prod([f.b for f in fs]), f1.a*prod([f.a for f in fs]))
 
 coefb(f::PolynomialRatio) = reverse(f.b.a)
 coefa(f::PolynomialRatio) = reverse(f.a.a)
 
 #
 # Biquad filter in transfer function form
-# A separate immutable to improve efficiency of filtering using SecondOrderSectionss
+# A separate immutable to improve efficiency of filtering using SecondOrderSections
 #
 
-immutable Biquad{T} <: FilterCoefficients
+immutable Biquad{T<:Number} <: FilterCoefficients
     b0::T
     b1::T
     b2::T
     a1::T
     a2::T
 end
-Biquad{T}(b0::T, b1::T, b2::T, a0::T, a1::T, a2::T, g::Real=1) =
+Biquad{T}(b0::T, b1::T, b2::T, a0::T, a1::T, a2::T, g::Number=1) =
     Biquad(g*b0/a0, g*b1/a0, g*b2/a0, a1/a0, a2/a0)
 
-Base.convert(::Type{ZeroPoleGain}, f::Biquad) = convert(ZeroPoleGain, convert(PolynomialRatio, f))
+Base.promote_rule{T,S}(::Type{Biquad{T}}, ::Type{Biquad{S}}) = Biquad{promote_type(T,S)}
+Base.convert{T}(::Type{Biquad{T}}, f::Biquad{T}) = f
+Base.convert{T}(::Type{Biquad{T}}, f::Biquad) = Biquad{T}(f.b0, f.b1, f.b2, f.a1, f.a2)
 
-function Base.convert{T}(::Type{PolynomialRatio}, f::Biquad{T})
+Base.convert{Z,P,K}(::Type{ZeroPoleGain{Z,P,K}}, f::Biquad) =
+    convert(ZeroPoleGain{Z,P,K}, convert(PolynomialRatio, f))
+Base.convert(::Type{ZeroPoleGain}, f::Biquad) =
+    convert(ZeroPoleGain, convert(PolynomialRatio, f))
+
+function Base.convert{T}(::Type{PolynomialRatio{T}}, f::Biquad)
     if f.b2 == zero(T) && f.a2 == zero(T)
         if f.b1 == zero(T) && f.a1 == zero(T)
-            b = [f.b0]
-            a = [one(T)]
+            b = T[f.b0]
+            a = T[one(T)]
         else
-            b = [f.b0, f.b1]
-            a = [one(T), f.a1]
+            b = T[f.b0, f.b1]
+            a = T[one(T), f.a1]
         end
     else
-        b = [f.b0, f.b1, f.b2]
-        a = [one(T), f.a1, f.a2]
+        b = T[f.b0, f.b1, f.b2]
+        a = T[one(T), f.a1, f.a2]
     end
 
     PolynomialRatio(b, a)
 end
+Base.convert{T}(::Type{PolynomialRatio}, f::Biquad{T}) = convert(PolynomialRatio{T}, f)
 
-Base.convert(::Type{Biquad}, f::ZeroPoleGain) = convert(Biquad, convert(PolynomialRatio, f))
-
-function Base.convert{T}(::Type{Biquad}, f::PolynomialRatio{T})
+function Base.convert{T}(::Type{Biquad{T}}, f::PolynomialRatio)
     a, b = f.a, f.b
     xs = max(length(b), length(a))
 
     if xs == 3
-        Biquad(b[2], b[1], b[0], a[1], a[0])
+        Biquad{T}(b[2], b[1], b[0], a[1], a[0])
     elseif xs == 2
-        Biquad(b[1], b[0], zero(T), a[0], zero(T))
+        Biquad{T}(b[1], b[0], zero(T), a[0], zero(T))
     elseif xs == 1
-        Biquad(b[0], zero(T), zero(T), zero(T), zero(T))
+        Biquad{T}(b[0], zero(T), zero(T), zero(T), zero(T))
     elseif xs == 0
         throw(ArgumentError("cannot convert an empty PolynomialRatio to Biquad"))
     else
         throw(ArgumentError("cannot convert a filter of length > 3 to Biquad"))
     end
 end
+Base.convert{T}(::Type{Biquad}, f::PolynomialRatio{T}) = convert(Biquad{T}, f)
+
+Base.convert{T}(::Type{Biquad{T}}, f::ZeroPoleGain) = convert(Biquad{T}, convert(PolynomialRatio, f))
+Base.convert(::Type{Biquad}, f::ZeroPoleGain) = convert(Biquad, convert(PolynomialRatio, f))
 
 *(f::Biquad, g::Number) = Biquad(f.b0*g, f.b1*g, f.b2*g, f.a1, f.a2)
+*(g::Number, f::Biquad) = Biquad(f.b0*g, f.b1*g, f.b2*g, f.a1, f.a2)
 
 #
 # Second-order sections (array of biquads)
@@ -115,27 +156,40 @@ immutable SecondOrderSections{T,G} <: FilterCoefficients
     g::G
 end
 
-realtype(x::DataType) = x
-realtype{T}(::Type{Complex{T}}) = T
-complextype(T::DataType) = Complex{T}
-complextype{T}(::Type{Complex{T}}) = Complex{T}
+Base.promote_rule{T1,G1,T2,G2}(::Type{SecondOrderSections{T1,G1}}, ::Type{SecondOrderSections{T2,G2}}) =
+    SecondOrderSections{promote_type(T1,T2),promote_type(G1,G2)}
+Base.convert{T,G}(::Type{SecondOrderSections{T,G}}, f::SecondOrderSections{T,G}) = f
+Base.convert{T,G}(::Type{SecondOrderSections{T,G}}, f::SecondOrderSections) =
+    SecondOrderSections{T,G}(f.biquads, f.g)
 
-function Base.convert{T}(::Type{ZeroPoleGain}, f::SecondOrderSections{T})
-    t = complextype(T)
-    z = t[]
-    p = t[]
+function Base.convert{Z,P,K}(::Type{ZeroPoleGain{Z,P,K}}, f::SecondOrderSections)
+    z = Z[]
+    p = P[]
     k = f.g
     for biquad in f.biquads
         biquadzpk = convert(ZeroPoleGain, biquad)
         append!(z, biquadzpk.z)
         append!(p, biquadzpk.p)
         k *= biquadzpk.k
     end
-    ZeroPoleGain(z, p, k)
+    ZeroPoleGain{Z,P,K}(z, p, k)
 end
+Base.convert{T,G}(::Type{ZeroPoleGain}, f::SecondOrderSections{T,G}) =
+    convert(ZeroPoleGain{complextype(T),complextype(T),G}, f)
 
-Base.convert(to::Union(Type{PolynomialRatio}, Type{Biquad}), f::SecondOrderSections) =
-    convert(to, convert(ZeroPoleGain, f))
+function Base.convert{T}(::Type{Biquad{T}}, f::SecondOrderSections)
+    if length(f.biquads) != 1
+        throw(ArgumentError("only a single second order section may be converted to a biquad"))
+    end
+    convert(Biquad{T}, f.biquads[1]*f.g)
+end
+Base.convert{T,G}(::Type{Biquad}, f::SecondOrderSections{T,G}) =
+    convert(Biquad{promote_type(T,G)}, f)
+
+Base.convert{T}(::Type{PolynomialRatio{T}}, f::SecondOrderSections) =
+    convert(PolynomialRatio{T}, convert(ZeroPoleGain, f))
+Base.convert(::Type{PolynomialRatio}, f::SecondOrderSections) =
+    convert(PolynomialRatio, convert(ZeroPoleGain, f))
 
 # Group each pole in p with its closest zero in z
 # Remove paired poles from p and z
@@ -252,4 +306,18 @@ function Base.convert{Z,P}(::Type{SecondOrderSections}, f::ZeroPoleGain{Z,P})
     SecondOrderSections(biquads, f.k)
 end
 
+Base.convert{T,G}(::Type{SecondOrderSections{T,G}}, f::Biquad) = SecondOrderSections{T,G}([f], one(G))
+Base.convert{T}(::Type{SecondOrderSections}, f::Biquad{T}) = convert(SecondOrderSections{T,Int}, f)
 Base.convert(::Type{SecondOrderSections}, f::FilterCoefficients) = convert(SecondOrderSections, convert(ZeroPoleGain, f))
+
+*(f::SecondOrderSections, g::Number) = SecondOrderSections(f.biquads, f.g*g)
+*(g::Number, f::SecondOrderSections) = SecondOrderSections(f.biquads, f.g*g)
+*(f1::SecondOrderSections, f2::SecondOrderSections) =
+    SecondOrderSections([f1.biquads; f2.biquads], f1.g*f2.g)
+*(f1::SecondOrderSections, fs::SecondOrderSections...) =
+    SecondOrderSections(vcat(f1.biquads, [f.biquads for f in fs]...), f1.g*prod([f.g for f in fs]))
+
+*(f1::Biquad, f2::Biquad) = SecondOrderSections([f1, f2], 1)
+*(f1::Biquad, fs::Biquad...) = SecondOrderSections([f1, fs...], 1)
+*(f1::SecondOrderSections, f2::Biquad) = SecondOrderSections([f1.biquads; f2], f1.g)
+*(f1::Biquad, f2::SecondOrderSections) = SecondOrderSections([f1; f2.biquads], f2.g)

test/filter_conversion.jl

@@ -1,5 +1,5 @@
 require(joinpath(dirname(@__FILE__), "FilterTestHelpers.jl"))
-using DSP, Base.Test, FilterTestHelpers
+using DSP, Base.Test, FilterTestHelpers, Polynomials
 
 # Test conversion to SOS against MATLAB
 # Test poles/zeros generated with:
@@ -131,23 +131,72 @@ for proto in (Butterworth(3), Chebyshev1(3, 1), Chebyshev2(3, 1))
     end
 end
 
+# Test setting filter gain
+x = randn(100)
+f1 = digitalfilter(Lowpass(0.3), Butterworth(2))
+y = filt(f1, x)
+for ty in (ZPKFilter, PolynomialRatio, Biquad, SecondOrderSections)
+    @test_approx_eq filt(3*convert(ty, f1), x) 3*y
+    @test_approx_eq filt(convert(ty, f1)*3, x) 3*y
+end
+
+# Test composing filters
+f2 = digitalfilter(Highpass(0.5), Butterworth(1))
+y = filt(f2, y)
+for ty in (ZPKFilter, PolynomialRatio, Biquad, SecondOrderSections)
+    @test_approx_eq filt(convert(ty, f1)*convert(ty, f2), x) y
+end
+
+f3 = digitalfilter(Bandstop(0.35, 0.4), Butterworth(1))
+y = filt(f3, y)
+for ty in (ZPKFilter, PolynomialRatio, Biquad, SecondOrderSections)
+    @test_approx_eq filt(convert(ty, f1)*convert(ty, f2)*convert(ty, f3), x) y
+end
+@test_approx_eq filt(convert(Biquad, f1)*(convert(Biquad, f2)*convert(Biquad, f3)), x) y
+@test_approx_eq filt((convert(Biquad, f1)*convert(Biquad, f2))*convert(Biquad, f3), x) y
+
 # Test some otherwise untested code paths
+@test promote_type(ZeroPoleGain{Complex64,Complex128,Float32}, ZeroPoleGain{Complex128,Complex64,Float64}) == ZeroPoleGain{Complex128,Complex128,Float64}
+@test convert(ZeroPoleGain{Float64,Complex128,Float64}, f1) === f1
+f1f = convert(ZeroPoleGain{Complex64,Complex64,Float32}, f1)
+@test f1f.z == convert(Vector{Complex64}, f1.z)
+@test f1f.p == convert(Vector{Complex64}, f1.p)
+@test f1f.k == convert(Float32, f1.k)
+
+@test_throws ArgumentError PolynomialRatio(Float64[], Float64[])
+@test promote_type(PolynomialRatio{Float32}, PolynomialRatio{Int}) == PolynomialRatio{Float32}
+f1f = convert(PolynomialRatio{Float32}, f1)
+f1p = convert(PolynomialRatio, f1)
+@test f1f.b == convert(Poly{Float32}, f1p.b)
+@test f1f.a == convert(Poly{Float32}, f1p.a)
+
 b = Biquad(0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 2)
 @test b.b0 === 0.0
 @test b.b1 === 2*1/3
 @test b.b2 === 2*2/3
 @test b.a1 === 4/3
 @test b.a2 === 5/3
+@test promote_type(Biquad{Float32}, Biquad{Int}) == Biquad{Float32}
+@test convert(Biquad{Float32}, b) == Biquad{Float32}(0.0, 2*1/3, 2*2/3, 4/3, 5/3)
+zpkfilter_eq(convert(ZeroPoleGain{Complex128,Complex128,Float64}, b), convert(ZeroPoleGain, b), 0.0)
 f = convert(PolynomialRatio, Biquad(2.0, 0.0, 0.0, 0.0, 0.0))
 @test coefb(f) == [2.0]
 @test coefa(f) == [1.0]
 @test convert(Biquad, PolynomialRatio([4.0], [2.0])) == Biquad(2.0, 0.0, 0.0, 0.0, 0.0)
 @test Biquad(2.0, 0.0, 0.0, 0.0, 0.0)*2 == Biquad(4.0, 0.0, 0.0, 0.0, 0.0)
-
-@test_throws ArgumentError PolynomialRatio(Float64[], Float64[])
+@test convert(Biquad{Float64}, f1) == convert(Biquad, f1)
 f = PolynomialRatio(Float64[1.0], Float64[1.0])
 empty!(f.b.a)
 empty!(f.a.a)
 @test_throws ArgumentError convert(Biquad, f)
 @test_throws ArgumentError convert(SOSFilter, ZPKFilter([0.5 + 0.5im, 0.5 + 0.5im], [0.5 + 0.5im, 0.5 - 0.5im], 1))
 @test_throws ArgumentError convert(SOSFilter, ZPKFilter([0.5 + 0.5im, 0.5 - 0.5im], [0.5 + 0.5im, 0.5 + 0.5im], 1))
+
+@test promote_type(SecondOrderSections{Float64,Float32}, SecondOrderSections{Float32,Float64}) == SecondOrderSections{Float64,Float64}
+@test convert(SecondOrderSections{Float32,Float32}, convert(SecondOrderSections, b)).biquads == convert(SecondOrderSections, convert(Biquad{Float32}, b)).biquads
+@test convert(Biquad, convert(SecondOrderSections, b)) == b
+@test_throws ArgumentError convert(Biquad, convert(SecondOrderSections, f1*f2))
+f1p = convert(PolynomialRatio{Float32}, convert(PolynomialRatio, convert(SecondOrderSections, f1)))
+f1f = convert(PolynomialRatio{Float32}, convert(SecondOrderSections, f1))
+@test f1p.b == f1f.b
+@test f1p.a == f1f.a