Fix _polyprep, PolynomialRatio, for Number args (#571)

Author: wheeheee

docs/src/internals.md

@@ -11,6 +11,7 @@ DSP.Windows.padplot
 DSP._zeropad
 DSP._zeropad!
 DSP._zeropad_keep_offset
+DSP.Filters._polyprep
 DSP.Filters.freq_eval
 DSP.Filters.build_grid
 DSP.Filters.lagrange_interp

src/Filters/coefficients.jl

@@ -34,13 +34,13 @@ ZeroPoleGain{D}(z::Vector{Z}, p::Vector{P}, k::K) where {D,Z<:Number,P<:Number,K
 Base.promote_rule(::Type{ZeroPoleGain{D,Z1,P1,K1}}, ::Type{ZeroPoleGain{D,Z2,P2,K2}}) where {D,Z1,P1,K1,Z2,P2,K2} =
     ZeroPoleGain{D,promote_type(Z1,Z2),promote_type(P1,P2),promote_type(K1,K2)}
 
-*(f::ZeroPoleGain{D}, g::Number) where {D} = ZeroPoleGain{D}(f.z, f.p, f.k*g)
-*(g::Number, f::ZeroPoleGain{D}) where {D} = ZeroPoleGain{D}(f.z, f.p, f.k*g)
+*(f::ZeroPoleGain{D}, g::Number) where {D} = ZeroPoleGain{D}(f.z, f.p, f.k * g)
+*(g::Number, f::ZeroPoleGain{D}) where {D} = ZeroPoleGain{D}(f.z, f.p, f.k * g)
 *(f1::ZeroPoleGain{D}, f2::ZeroPoleGain{D}) where {D} =
-    ZeroPoleGain{D}([f1.z; f2.z], [f1.p; f2.p], f1.k*f2.k)
+    ZeroPoleGain{D}([f1.z; f2.z], [f1.p; f2.p], f1.k * f2.k)
 *(f1::ZeroPoleGain{D}, fs::ZeroPoleGain{D}...) where {D} =
     ZeroPoleGain{D}(vcat(f1.z, map(f -> f.z, fs)...), vcat(f1.p, map(f -> f.p, fs)...),
-        f1.k*prod(f.k for f in fs))
+        f1.k * prod(f.k for f in fs))
 
 Base.inv(f::ZeroPoleGain{D}) where {D} = ZeroPoleGain{D}(f.p, f.z, inv(f.k))
 
@@ -56,9 +56,9 @@ end
 # Transfer function form
 #
 
-function shiftpoly(p::LaurentPolynomial{T}, i::Integer) where {T<:Number}
+function shiftpoly(p::LaurentPolynomial{T,D}, i::Integer) where {T<:Number,D}
     if !iszero(i)
-        return p * LaurentPolynomial([one(T)], i, indeterminate(p))
+        return p * LaurentPolynomial{T,D}([one(T)], i)
     end
     return p
 end
@@ -122,20 +122,35 @@ PolynomialRatio{:s, Int64}(LaurentPolynomial(3 + 2*s + s²), LaurentPolynomial(4
 
 """
 PolynomialRatio(b, a) = PolynomialRatio{:z}(b, a)
-function PolynomialRatio{:z}(b::LaurentPolynomial{T1}, a::LaurentPolynomial{T2}) where {T1,T2}
-    return PolynomialRatio{:z,typeof(one(T1)/one(T2))}(b, a)
-end
-function PolynomialRatio{:s}(b::LaurentPolynomial{T1}, a::LaurentPolynomial{T2}) where {T1,T2}
-    return PolynomialRatio{:s,promote_type(T1, T2)}(b, a)
-end
+const PolynomialRatioArgTs{T} = Union{T,Vector{T},LaurentPolynomial{T}} where {T<:Number}
+PolynomialRatio{:z}(b::PolynomialRatioArgTs{T1}, a::PolynomialRatioArgTs{T2}) where {T1<:Number,T2<:Number} =
+    PolynomialRatio{:z,typeof(one(T1) / one(T2))}(b, a)
+PolynomialRatio{:s}(b::PolynomialRatioArgTs{T1}, a::PolynomialRatioArgTs{T2}) where {T1<:Number,T2<:Number} =
+    PolynomialRatio{:s,promote_type(T1, T2)}(b, a)
+
+"""
+    _polyprep(D::Symbol, x::Union{T,Vector{T}}, ::Type) where {T<:Number}
+
+Converts `x` to polynomial form. If x is a `Number`, it has to be converted into
+a `Vector`, otherwise `LaurentPolynomial` dispatch goes into stack overflow
+trying to collect a 0-d array into a `Vector`.
 
-# The DSP convention for Laplace domain is highest power first. The Polynomials.jl
-# convention is lowest power first.
-_polyprep(D::Symbol, x, T...) = LaurentPolynomial{T...}(reverse(x), D === :z ? -length(x)+1 : 0, D)
-PolynomialRatio{D,T}(b::Union{Number,Vector{<:Number}}, a::Union{Number,Vector{<:Number}}) where {D,T} =
-    PolynomialRatio{D,T}(_polyprep(D, b, T), _polyprep(D, a, T))
-PolynomialRatio{D}(b::Union{Number,Vector{<:Number}}, a::Union{Number,Vector{<:Number}}) where {D} =
-    PolynomialRatio{D}(_polyprep(D, b), _polyprep(D, a))
+!!! warning
+
+    The DSP convention for Laplace domain is highest power first.\n
+    The Polynomials.jl convention is lowest power first.
+"""
+@inline _polyprep(D::Symbol, x::Union{T,Vector{T}}, ::Type{V}) where {T<:Number,V} =
+    LaurentPolynomial{V,D}(x isa Vector ? reverse(x) : [x], D === :z ? -length(x) + 1 : 0)
+
+function PolynomialRatio{:z,T}(b::Union{Number,Vector{<:Number}}, a::Union{Number,Vector{<:Number}}) where {T<:Number}
+    if isempty(a) || iszero(a[1])
+        throw(ArgumentError("filter must have non-zero leading denominator coefficient"))
+    end
+    return PolynomialRatio{:z,T}(_polyprep(:z, b / a[1], T), _polyprep(:z, a / a[1], T))
+end
+PolynomialRatio{:s,T}(b::Union{Number,Vector{<:Number}}, a::Union{Number,Vector{<:Number}}) where {T<:Number} =
+    PolynomialRatio{:s,T}(_polyprep(:s, b, T), _polyprep(:s, a, T))
 
 PolynomialRatio{D,T}(f::PolynomialRatio{D}) where {D,T} = PolynomialRatio{D,T}(f.b, f.a)
 PolynomialRatio{D}(f::PolynomialRatio{D,T}) where {D,T} = PolynomialRatio{D,T}(f)
@@ -159,16 +174,14 @@ function ZeroPoleGain{D}(f::PolynomialRatio{D,T}) where {D,T}
     return ZeroPoleGain{D}(z, p, f.b[end]/f.a[end])
 end
 
-*(f::PolynomialRatio{D}, g::Number) where {D} = PolynomialRatio{D}(g*f.b, f.a)
-*(g::Number, f::PolynomialRatio{D}) where {D} = PolynomialRatio{D}(g*f.b, f.a)
+*(f::PolynomialRatio{D}, g::Number) where {D} = PolynomialRatio{D}(g * f.b, f.a)
+*(g::Number, f::PolynomialRatio{D}) where {D} = PolynomialRatio{D}(g * f.b, f.a)
 *(f1::PolynomialRatio{D}, f2::PolynomialRatio{D}) where {D} =
-    PolynomialRatio{D}(f1.b*f2.b, f1.a*f2.a)
+    PolynomialRatio{D}(f1.b * f2.b, f1.a * f2.a)
 *(f1::PolynomialRatio{D}, fs::PolynomialRatio{D}...) where {D} =
-    PolynomialRatio{D}(f1.b*prod(f.b for f in fs), f1.a*prod(f.a for f in fs))
+    PolynomialRatio{D}(f1.b * prod(f.b for f in fs), f1.a * prod(f.a for f in fs))
 
-Base.inv(f::PolynomialRatio{D}) where {D} = begin
-    PolynomialRatio{D}(f.a, f.b)
-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}
     if e < 0
@@ -312,7 +325,7 @@ 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)
+    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)
 
@@ -418,7 +431,7 @@ function SecondOrderSections{D,T,G}(f::ZeroPoleGain{D,Z,P}) where {D,T,G,Z,P}
     npairs = length(groupedp) >> 1
     odd = isodd(n)
     for i = 1:npairs
-        pairidx = 2*(npairs-i)
+        pairidx = 2 * (npairs - i)
         biquads[odd+i] = convert(Biquad, ZeroPoleGain{D}(groupedz[pairidx+1:min(pairidx+2, length(groupedz))],
                                                          groupedp[pairidx+1:pairidx+2], one(T)))
     end
@@ -439,12 +452,12 @@ SecondOrderSections{D,T,G}(f::Biquad{D}) where {D,T,G} = SecondOrderSections{D,T
 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)
+*(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)
+    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))
+    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)

src/Filters/remez_fir.jl

@@ -333,8 +333,8 @@ Here we compute the frequency responses and plot them in dB.
 
 ```julia-repl
 julia> using PyPlot
-julia> b = DSP.Filters.PolynomialRatio(bpass, [1.0])
-julia> b2 = DSP.Filters.PolynomialRatio(bpass2, [1.0])
+julia> b = PolynomialRatio(bpass, [1.0])
+julia> b2 = PolynomialRatio(bpass2, [1.0])
 julia> f = range(0, stop=0.5, length=1000)
 julia> plot(f, 20*log10.(abs.(freqresp(b,f,1.0))))
 julia> plot(f, 20*log10.(abs.(freqresp(b2,f,1.0))))

test/filt.jl

@@ -138,7 +138,7 @@ end
     b = [0.4, 1]
     z = [0.4750]
     x  = readdlm(joinpath(dirname(@__FILE__), "data", "spectrogram_x.txt"),'\t')
-    DSP.Filters.filt!(vec(x), b, a, vec(x), z)
+    filt!(vec(x), b, a, vec(x), z)
 
     @test matlab_filt ≈ x
 end
@@ -272,7 +272,7 @@ end
 @testset "fir_filtfilt" begin
     b = randn(10)
     for x in (randn(100), randn(100, 2))
-        @test DSP.Filters.filtfilt(b, x) ≈ DSP.Filters.iir_filtfilt(b, [1.0], x)
-        @test DSP.Filters.filtfilt(b, [2.0], x) ≈ DSP.Filters.iir_filtfilt(b, [2.0], x)
+        @test filtfilt(b, x) ≈ DSP.Filters.iir_filtfilt(b, [1.0], x)
+        @test filtfilt(b, [2.0], x) ≈ DSP.Filters.iir_filtfilt(b, [2.0], x)
     end
 end

test/filter_conversion.jl

@@ -249,11 +249,15 @@ end
     @test_throws InexactError PolynomialRatio{:z,Int}([1, 2], [3, 4])
     # throws because normalization is impossible for a0 = 0
     @test_throws ArgumentError PolynomialRatio{:z,Float64}([1.0, 2.0], [0.0, 4.0])
+    # throws because a0 = 0, in inner constructor
+    @test_throws ArgumentError inv(PolynomialRatio{:z,Float64}(0, 1))
     # does not normalize, uses type from input
     @test @inferred(PolynomialRatio{:s}([1, 2], [3, 4])) isa PolynomialRatio{:s,Int}
     # throws because denominator must not be zero
     @test_throws ArgumentError PolynomialRatio{:s}([1.0, 2.0], [0.0])
     @test_throws ArgumentError PolynomialRatio{:s}([1.0, 2.0], Float64[])
+    # test PolynomialRatio{:s} constructor for LaurentPolynomial input
+    @test inv(PolynomialRatio{:s}(1.0, 2.3)).a == PolynomialRatio{:s}(1.0, 2.3).b
 end
 
 @testset "misc" begin
@@ -286,10 +290,10 @@ end
     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 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 convert(Biquad{:z,Float64}, f1) == convert(Biquad, f1)
-    f = PolynomialRatio(Float64[1.0], Float64[1.0])
+    f = PolynomialRatio(1.0, 1.0)       # doubles as test for Number arguments (PR #571)
 
     @test_throws ArgumentError convert(SecondOrderSections, ZeroPoleGain([0.5 + 0.5im, 0.5 + 0.5im], [0.5 + 0.5im, 0.5 - 0.5im], 1))
     @test_throws ArgumentError convert(SecondOrderSections, ZeroPoleGain([0.5 + 0.5im, 0.5 - 0.5im], [0.5 + 0.5im, 0.5 + 0.5im], 1))

test/remez_fir.jl

@@ -163,19 +163,19 @@ end
 #
 @testset "inverse_sinc_response_function" begin
     L = 64
-    
+
     Fs = 4800*L
     f = range(0, stop=0.5, length=10000)
 
     P = (π*f*Fs/4800) ./ sin.(π*f*Fs/4800)
     Pdb = 20*log10.(abs.(P))
     Pdb[1] = 0.0
-    
+
     g_vec = remez(201, [
                 (    0.0, 2880.0) => (f -> (f==0) ? 1.0 : abs.((π*f/4800) ./ sin.(π*f/4800)), 1.0),
                 (10000.0,  Fs/2) => (0.0, 100.0)
         ]; Hz=Fs)
-    g = DSP.Filters.PolynomialRatio(g_vec, [1.0])
+    g = PolynomialRatio(g_vec, [1.0])
     Gdb = 20*log10.(abs.(freqresp(g, 2π*f)))
 
     passband_indices = (f*Fs) .< 2880.0