Skip to content

Commit

Permalink
Use Polynomials.jl instead of Polynomial.jl
Browse files Browse the repository at this point in the history
  • Loading branch information
simonster committed Jun 23, 2014
1 parent 4c449be commit fdef6ae
Show file tree
Hide file tree
Showing 4 changed files with 24 additions and 14 deletions.
2 changes: 1 addition & 1 deletion REQUIRE
Original file line number Diff line number Diff line change
@@ -1,2 +1,2 @@
Polynomial
Polynomials
Reexport
23 changes: 14 additions & 9 deletions src/filter_design.jl
Original file line number Diff line number Diff line change
Expand Up @@ -39,7 +39,7 @@
# DAMAGE.

module FilterDesign
using Polynomial
using Polynomials

export ZPKFilter, TFFilter, BiquadFilter, SOSFilter, Butterworth, Lowpass, Highpass, Bandpass,
Bandstop, analogfilter, digitalfilter, Filter
Expand All @@ -66,7 +66,7 @@ end
similarzeros{T}(v::Array{T}, args...) = zeros(T, args...)

# Get coefficients of a polynomial
coeffs{T}(p::Poly{T}) = p.a[1+p.nzfirst:end]
coeffs{T}(p::Poly{T}) = reverse(p.a)

#
# Filter types
Expand All @@ -87,11 +87,16 @@ immutable TFFilter{T} <: Filter
a::Poly{T}

function TFFilter(b::Poly{T}, a::Poly{T})
new(b/a[1], a/a[1])
new(b/a[end], a/a[end])
end
end
TFFilter{T}(b::Poly{T}, a::Poly{T}) = TFFilter{T}(b, a)
TFFilter{T}(b::Vector{T}, a::Vector{T}) = TFFilter{T}(Poly(b), Poly(a))

# The DSP convention is lowest power first. The Polynomials.jl
# convention is highest power first.
TFFilter{T}(b::Vector{T}, a::Vector{T}) =
TFFilter{T}(Poly(b[end:-1:findfirst(b)]), Poly(a[end:-1:findfirst(a)]))

function TFFilter{T,S}(b::Vector{T}, a::Vector{S})
V = promote_type(T, S)
TFFilter(convert(Vector{V}, b), convert(Vector{V}, a))
Expand All @@ -104,7 +109,7 @@ function convert(::Type{TFFilter}, f::ZPKFilter)
end

function convert{T}(::Type{ZPKFilter}, f::TFFilter{T})
k = real(f.b[1])
k = real(f.b[end])
b = f.b / k
z = convert(Vector{Complex{T}}, roots(b))
p = convert(Vector{Complex{T}}, roots(f.a))
Expand Down Expand Up @@ -336,8 +341,8 @@ function transform_prototype(ftype::Bandpass, proto::TFFilter)
bw = ftype.w2 - ftype.w1
wo = sqrt(ftype.w1 * ftype.w2)
tf = convert(TFFilter, proto)
b = tf.b.a
a = tf.a.a
b = coeffs(tf.b)
a = coeffs(tf.a)
D = length(a) - 1
N = length(b) - 1
M = max(N, D)
Expand Down Expand Up @@ -378,8 +383,8 @@ function transform_prototype(ftype::Bandstop, proto::TFFilter)
bw = ftype.w2 - ftype.w1
wo = sqrt(ftype.w1 * ftype.w2)
tf = convert(TFFilter, proto)
b = tf.b.a
a = tf.a.a
b = coeffs(tf.b)
a = coeffs(tf.a)
D = length(a) - 1
N = length(b) - 1
M = max(N, D)
Expand Down
2 changes: 1 addition & 1 deletion src/filter_response.jl
Original file line number Diff line number Diff line change
Expand Up @@ -5,7 +5,7 @@ module FilterResponse

export freqz, freqs

using Polynomial
using Polynomials
using ..FilterDesign


Expand Down
11 changes: 8 additions & 3 deletions test/filter_design.jl
Original file line number Diff line number Diff line change
@@ -1,4 +1,4 @@
using DSP, Base.Test, Polynomial
using DSP, Base.Test, Polynomials
import Base.Sort.Lexicographic
import DSP.FilterDesign: coeffs

Expand Down Expand Up @@ -81,7 +81,7 @@ accurate_a = coeffs(convert(TFFilter, accurate_f).a)

f = convert(ZPKFilter, convert(TFFilter, Butterworth(20)))
@test isempty(f.z)
@test_approx_eq_eps sort(f.p, lt=lt) prototype 1e-7
@test_approx_eq_eps sort(f.p, lt=lt) prototype 1e-6

# Test that our answers are more accurate than MATLAB's
# Output of [z, p, k] = buttap(20); [b, a] = zp2tf(z, p, k); tf2zpk(b, a)
Expand Down Expand Up @@ -221,7 +221,12 @@ for f in (digitalfilter(Lowpass(0.5), Butterworth(1)), digitalfilter(Lowpass(0.5
f2 = convert(ftype1, f)
for ftype2 in (ZPKFilter, TFFilter, BiquadFilter, SOSFilter)
f3 = convert(ftype2, f)
zpkfilter_eq(f, convert(ZPKFilter, f3), eps())
try
zpkfilter_eq(f, convert(ZPKFilter, f3), sqrt(eps()))
catch e
println("Conversion from $ftype1 to $ftype2 failed:")
rethrow(e)
end
end
end
end
Expand Down

0 comments on commit fdef6ae

Please sign in to comment.