!(dirname(@__FILE__) in LOAD_PATH) && push!(LOAD_PATH, dirname(@__FILE__))
using DSP, Test, FilterTestHelpers

#
# Butterworth filter prototype
#
@testset "Butterworth" begin
    # Poles of 20 pole Butterworth filter prototype from MATLAB 2013b (buttap(20))
    matlab_p = [-0.07845909572784487+0.996917333733128im,
        -0.07845909572784487-0.996917333733128im,
        -0.2334453638559053+0.9723699203976767im,
        -0.2334453638559053-0.9723699203976767im,
        -0.3826834323650897+0.9238795325112867im,
        -0.3826834323650897-0.9238795325112867im,
        -0.5224985647159488+0.8526401643540923im,
        -0.5224985647159488-0.8526401643540923im,
        -0.6494480483301835+0.760405965600031im,
        -0.6494480483301835-0.760405965600031im,
        -0.7604059656000308+0.6494480483301838im,
        -0.7604059656000308-0.6494480483301838im,
        -0.8526401643540922+0.5224985647159489im,
        -0.8526401643540922-0.5224985647159489im,
        -0.9238795325112867+0.3826834323650899im,
        -0.9238795325112867-0.3826834323650899im,
        -0.9723699203976766+0.2334453638559055im,
        -0.9723699203976766-0.2334453638559055im,
        -0.996917333733128+0.07845909572784507im,
        -0.996917333733128-0.07845909572784507im]
    matlab_butter = ZeroPoleGain(Float64[], matlab_p, 1)

    # Test that our answers are close to MATLAB's and at least as accurate
    butter = Butterworth(20)
    zpkfilter_eq(butter, matlab_butter)
    zpkfilter_accuracy(butter, matlab_butter, Butterworth(BigFloat, 20))

    # Poles of 19 pole Butterworth filter prototype from MATLAB 2013b (buttap(19))
    matlab_p = [-0.08257934547233227+0.9965844930066698im,
        -0.08257934547233227-0.9965844930066698im,
        -0.2454854871407991+0.9694002659393304im,
        -0.2454854871407991-0.9694002659393304im,
        -0.4016954246529694+0.9157733266550574im,
        -0.4016954246529694-0.9157733266550574im,
        -0.5469481581224267+0.8371664782625287im,
        -0.5469481581224267-0.8371664782625287im,
        -0.6772815716257409+0.7357239106731317im,
        -0.6772815716257409-0.7357239106731317im,
        -0.7891405093963935+0.6142127126896679im,
        -0.7891405093963935-0.6142127126896679im,
        -0.8794737512064892+0.4759473930370733im,
        -0.8794737512064892-0.4759473930370733im,
        -0.9458172417006346+0.3246994692046836im,
        -0.9458172417006346-0.3246994692046836im,
        -0.9863613034027224+0.1645945902807336im,
        -0.9863613034027224-0.1645945902807336im,
        -1+0im]

    # Test that our answers are close to MATLAB's and at least as accurate
    zpkfilter_eq(Butterworth(19), ZeroPoleGain(Float64[], matlab_p, 1))
    zpkfilter_accuracy(Butterworth(19), ZeroPoleGain(Float64[], matlab_p, 1), Butterworth(BigFloat, 19))
end

#
# Chebyshev type I filter
#
@testset "Chebyshev I" begin
    # Poles of 20 pole Butterworth filter prototype with 1 dB passband ripple from MATLAB 2013b:
    #=
        [z, p, k] = cheb1ap(20, 1)
        sprintf('%.16g%+.16gim,', [real(p) imag(p)]')
    =#
    matlab_p = [-0.005606643513655412+0.9994594480354867im,
        -0.01668187637027696+0.9748494394091176im,
        -0.02734634606668882+0.9262354022446964im,
        -0.03733745796561948+0.854814375477951im,
        -0.04640919773351807+0.7623449818206758im,
        -0.05433818902992066+0.6511041246330369im,
        -0.06092919377429644+0.5238309230062319im,
        -0.06601991955554824+0.3836592655598524im,
        -0.06948501580979141+0.2340406437032877im,
        -0.07123916036725134+0.07865916446238384im,
        -0.07123916036725134-0.07865916446238384im,
        -0.06948501580979141-0.2340406437032877im,
        -0.06601991955554824-0.3836592655598524im,
        -0.06092919377429644-0.5238309230062319im,
        -0.05433818902992066-0.6511041246330369im,
        -0.04640919773351807-0.7623449818206758im,
        -0.03733745796561948-0.854814375477951im,
        -0.02734634606668882-0.9262354022446964im,
        -0.01668187637027696-0.9748494394091176im,
        -0.005606643513655412-0.9994594480354867im]
    matlab_k = 3.748372513504540e-06

    # Test that our answers are close to MATLAB's and at least as accurate
    zpkfilter_eq(Chebyshev1(20, 1), ZeroPoleGain(Float64[], matlab_p, matlab_k))
    zpkfilter_accuracy(Chebyshev1(20, 1), ZeroPoleGain(Float64[], matlab_p, matlab_k), Chebyshev1(BigFloat, 20, 1))

    # Poles of 19 pole Butterworth filter prototype with 1 dB passband ripple from MATLAB 2013b:
    #=
        [z, p, k] = cheb1ap(19, 1)
        sprintf('%.16g%+.16gim,', [real(p) imag(p)]')
        sprintf('%.16g', k)
    =#
    matlab_p = [-0.006212227114275604+0.9994004289494353im,
        -0.0184672279812169+0.9721393904901743im,
        -0.03021849100930104+0.9183609236365976im,
        -0.04114547237637834+0.8395319647744703im,
        -0.05095011251526942+0.7378027624097928im,
        -0.0593649664017747+0.6159482239948341im,
        -0.06616049875775722+0.4772922236864799im,
        -0.07115134517517649+0.3256169357239265im,
        -0.07420136837393095+0.1650596665748249im,
        -0.07522737167197566+0im,
        -0.07420136837393095-0.1650596665748249im,
        -0.07115134517517649-0.3256169357239265im,
        -0.06616049875775722-0.4772922236864799im,
        -0.0593649664017747-0.6159482239948341im,
        -0.05095011251526942-0.7378027624097928im,
        -0.04114547237637834-0.8395319647744703im,
        -0.03021849100930104-0.9183609236365976im,
        -0.0184672279812169-0.9721393904901743im,
        -0.006212227114275604-0.9994004289494353im]
    matlab_k = 7.496745027009062e-06

    # Test that our answers are close to MATLAB's and at least as accurate
    zpkfilter_eq(Chebyshev1(19, 1), ZeroPoleGain(Float64[], matlab_p, matlab_k))
    zpkfilter_accuracy(Chebyshev1(19, 1), ZeroPoleGain(Float64[], matlab_p, matlab_k), Chebyshev1(BigFloat, 19, 1))
end

#
# Chebyshev type II filter
#
@testset "Chebyshev II" begin
    # 20 pole Butterworth filter prototype with 1 dB passband ripple from MATLAB 2013b:
    #=
        [z, p, k] = cheb2ap(20, 1)
        sprintf('%.16g%+.16gim,', [real(z) imag(z)]')
        sprintf('%.16g%+.16gim,', [real(p) imag(p)]')
        sprintf('%.16g', k)
    =#
    matlab_z = [0+1.003092198482826im, 0-1.003092198482826im,
        0+1.028415193665209im, 0-1.028415193665209im,
        0+1.082392200292394im, 0-1.082392200292394im,
        0+1.172827696614009im, 0-1.172827696614009im,
        0+1.315086999890785im, 0-1.315086999890785im,
        0+1.539769043222366im, 0-1.539769043222366im,
        0+1.913880855430943im, 0-1.913880855430943im,
        0+2.613125929752753im, 0-2.613125929752753im,
        0+4.283657569731185im, 0-4.283657569731185im,
        0+12.74549484318238im, 0-12.74549484318238im]
    matlab_p = [-0.001929675700544753-1.002788598382701im,
        -0.006034873757627809-1.028072311003105im,
        -0.010957850142167-1.081957626858823im,
        -0.01756373888305481-1.172213900920771im,
        -0.02744258171352171-1.314120759169777im,
        -0.04403133250736632-1.538048178699878im,
        -0.07621942506983267-1.910267531994771im,
        -0.1536690174908018-2.6032737531287im,
        -0.4316587213861887-4.238414904224276im,
        -3.610083779585366-11.62012078608171im,
        -3.610083779585366+11.62012078608171im,
        -0.4316587213861887+4.238414904224276im,
        -0.1536690174908018+2.6032737531287im,
        -0.07621942506983267+1.910267531994771im,
        -0.04403133250736632+1.538048178699878im,
        -0.02744258171352171+1.314120759169777im,
        -0.01756373888305481+1.172213900920771im,
        -0.010957850142167+1.081957626858823im,
        -0.006034873757627809+1.028072311003105im,
        -0.001929675700544753+1.002788598382701im]
    matlab_k = 0.8912509381337452
    matlab_f = ZeroPoleGain(matlab_z, matlab_p, matlab_k)

    # The gain shows a greater discrepancy with the gain of the enhanced
    # precision filter than MATLAB's, but AFAICT our filter is still more
    # accurate: the poles and zeros are each more accurate, and computing
    # the gain in enhanced precision yields the same value.

    f = Chebyshev2(20, 1)
    zpkfilter_eq(f, matlab_f)
    zpkfilter_accuracy(f, matlab_f, Chebyshev2(BigFloat, 20, 1), relerr=4)

    # 19 pole Butterworth filter prototype with 1 dB passband ripple from MATLAB 2013b:
    #=
        [z, p, k] = cheb2ap(19, 1)
        sprintf('%.16g%+.16gim,', [real(z) imag(z)]')
        sprintf('%.16g%+.16gim,', [real(p) imag(p)]')
        sprintf('%.16g', k)
    =#
    matlab_z = [0+1.003427212662145im, 0-1.003427212662145im,
    0+1.031565634068626im, 0-1.031565634068626im,
    0+1.091973276457601im, 0-1.091973276457601im,
    0+1.194505544554793im, 0-1.194505544554793im,
    0+1.359205519207709im, 0-1.359205519207709im,
    0+1.628100459889458im, 0-1.628100459889458im,
    0+2.101072544213108im, 0-2.101072544213108im,
    0+3.079770972368366im, 0-3.079770972368366im,
    0+6.075533820974263im, 0-6.075533820974263im]
    matlab_p = [-0.002139218568677465-1.003090263923361im,
        -0.006720708432694556-1.031180124433753im,
        -0.01232195186345456-1.091472452757043im,
        -0.02007306444044342-1.193772338103026im,
        -0.03217434756636293-1.357992987527419im,
        -0.05375960613144559-1.625783486413104im,
        -0.09966387580015423-2.09563676642751im,
        -0.2295218019875045-3.061543704346429im,
        -0.9149742173123645-5.932401750359961im,
        -38.84170469908301+0im,
        -0.9149742173123645+5.932401750359961im,
        -0.2295218019875045+3.061543704346429im,
        -0.09966387580015423+2.09563676642751im,
        -0.05375960613144559+1.625783486413104im,
        -0.03217434756636293+1.357992987527419im,
        -0.02007306444044342+1.193772338103026im,
        -0.01232195186345456+1.091472452757043im,
        -0.006720708432694556+1.031180124433753im,
        -0.002139218568677465+1.003090263923361im]
    matlab_k = 37.33930783884512
    matlab_f = ZeroPoleGain(matlab_z, matlab_p, matlab_k)

    f = Chebyshev2(19, 1)
    zpkfilter_eq(f, matlab_f)
    zpkfilter_accuracy(f, matlab_f, Chebyshev2(BigFloat, 19, 1), relerr=2)
end

#
# Elliptic filter
#
@testset "elliptic" begin
    # 20 pole elliptic filter prototype with 0.1 dB passband ripple and 10
    # dB stopband ripple from MATLAB 2013b:
    #=
        [z, p, k] = ellipap(20, 0.1, 10)
        sprintf('%.16g%+.16gim,', [real(z) imag(z)]')
        sprintf('%.16g%+.16gim,', [real(p) imag(p)]')
        sprintf('%.16g', k)
    =#
    matlab_z = [0-1.953252853757711im, 0+1.953252853757711im,
        0-1.069599937693626im, 0+1.069599937693626im,
        0-1.007032209276402im, 0+1.007032209276402im,
        0-1.000730367252677im, 0+1.000730367252677im,
        0-1.000076070678288im, 0+1.000076070678288im,
        0-1.000007925871874im, 0+1.000007925871874im,
        0-1.000000826312751im, 0+1.000000826312751im,
        0-1.000000086632384im, 0+1.000000086632384im,
        0-1.000000009576065im, 0+1.000000009576065im,
        0-1.000000001633984im, 0+1.000000001633984im]
    matlab_p = [-0.6665552331000151-0.9586923119289487im,
        -0.6665552331000151+0.9586923119289487im,
        -0.06598422536270995-1.015686365686964im,
        -0.06598422536270995+1.015686365686964im,
        -0.006773291804375061-1.001821704927991im,
        -0.006773291804375061+1.001821704927991im,
        -0.0007045345890457818-1.000191780070487im,
        -0.0007045345890457818+1.000191780070487im,
        -7.339097059946433e-05-1.000020003037497im,
        -7.339097059946433e-05+1.000020003037497im,
        -7.646273200288023e-06-1.000002084836051im,
        -7.646273200288023e-06+1.000002084836051im,
        -7.966432964216136e-07-1.000000217756124im,
        -7.966432964216136e-07+1.000000217756124im,
        -8.299912436196738e-08-1.000000023227775im,
        -8.299912436196738e-08+1.000000023227775im,
        -8.637784231770464e-09-1.000000002962719im,
        -8.637784231770464e-09+1.000000002962719im,
        -8.070974080610083e-10-1.000000000874029im,
        -8.070974080610083e-10+1.000000000874029im]
    matlab_k = 0.3162277662398871
    matlab_f = ZeroPoleGain(matlab_z, matlab_p, matlab_k)

    f = Elliptic(20, 0.1, 10)
    zpkfilter_eq(f, matlab_f)
    zpkfilter_accuracy(f, matlab_f, Elliptic(BigFloat, 20, 0.1, 10), eps=1e-9)

    # 19 pole elliptic filter prototype with 0.1 dB passband ripple and 10
    # dB stopband ripple from MATLAB 2013b:
    #=
        [z, p, k] = ellipap(19, 0.1, 10)
        sprintf('%.16g%+.16gim,', [real(z) imag(z)]')
        sprintf('%.16g%+.16gim,', [real(p) imag(p)]')
        sprintf('%.16g', k)
    =#
    matlab_z = [0-1.232609672486912im, 0+1.232609672486912im,
        0-1.021948255265862im, 0+1.021948255265862im,
        0-1.002264470627064im, 0+1.002264470627064im,
        0-1.000235691798532im, 0+1.000235691798532im,
        0-1.000024555183434im, 0+1.000024555183434im,
        0-1.000002559984621im, 0+1.000002559984621im,
        0-1.000000268394018im, 0+1.000000268394018im,
        0-1.00000002966741im, 0+1.00000002966741im,
        0-1.000000005062212im, 0+1.000000005062212im]
    matlab_p = [-0.2102804655411482-1.034611529620449im,
        -0.2102804655411482+1.034611529620449im,
        -0.0210635609742992-1.005504639644632im,
        -0.0210635609742992+1.005504639644632im,
        -0.002183581313750937-1.00059265523267im,
        -0.002183581313750937+1.00059265523267im,
        -0.0002273805846473938-1.000061954788962im,
        -0.0002273805846473938+1.000061954788962im,
        -2.368886404726722e-05-1.000006458811389im,
        -2.368886404726722e-05+1.000006458811389im,
        -2.468065209676509e-06-1.00000067462382im,
        -2.468065209676509e-06+1.00000067462382im,
        -2.571378798158274e-07-1.000000071961468im,
        -2.571378798158274e-07+1.000000071961468im,
        -2.676054057433227e-08-1.000000009178736im,
        -2.676054057433227e-08+1.000000009178736im,
        -2.500451772444974e-09-1.00000000270781im,
        -2.500451772444974e-09+1.00000000270781im,
        -1.309071549907812+0im]
    matlab_k = 0.9266824319626478
    matlab_f = ZeroPoleGain(matlab_z, matlab_p, matlab_k)

    f = Elliptic(19, 0.1, 10)
    zpkfilter_eq(f, matlab_f)
    zpkfilter_accuracy(f, matlab_f, Elliptic(BigFloat, 19, 0.1, 10), eps=4e-9)
end

#
# Conversion between zpk and tf
#
@testset "ZPK to TF" begin
    # Output of [z, p, k] = buttap(20); [b, a] = zp2tf(z, p, k)
    m_a = [1, 12.74549484318237, 81.22381939879423, 343.6513712403923,
        1081.352361133001, 2687.409807920676, 5468.931438945091,
        9326.061201886809, 13528.36656744904, 16852.27707949905,
        18122.54155403868, 16852.27707949905, 13528.36656744904,
        9326.061201886809, 5468.931438945092, 2687.409807920676,
        1081.352361133001, 343.6513712403923, 81.22381939879423,
        12.74549484318237, 1]

    f = convert(PolynomialRatio, Butterworth(20))
    @test coefb(f) ≈ [1]
    @test coefa(f) ≈ m_a

    # Test that our answers are more accurate than MATLAB's
    accurate_a = coefa(convert(PolynomialRatio, Butterworth(BigFloat, 20)))
    @test coefa(f) ≈ map(Float64, accurate_a)
    @test sum(abs.(coefa(f) - accurate_a)) <= sum(abs.(m_a - accurate_a))
end

#
# Conversion between tf and zpk
#
@testset "TF to ZPK" begin
    f = convert(ZeroPoleGain, convert(PolynomialRatio, Butterworth(20)))
    zpkfilter_eq(f, Butterworth(20), 1e-6)

    # TODO: Make this more accurate!

    # Output of [z, p, k] = buttap(20); [b, a] = zp2tf(z, p, k); tf2zpk(b, a)
    m_p = [-0.07845909573254482+0.9969173337335029im,
        -0.07845909573254482-0.9969173337335029im,
        -0.2334453637958131+0.9723699203918822im,
        -0.2334453637958131-0.9723699203918822im,
        -0.3826834327796701+0.9238795325396184im,
        -0.3826834327796701-0.9238795325396184im,
        -0.5224985628488221+0.8526401643454914im,
        -0.5224985628488221-0.8526401643454914im,
        -0.6494480541398985+0.7604059651905597im,
        -0.6494480541398985-0.7604059651905597im,
        -0.760405952587916+0.6494480502272874im,
        -0.760405952587916-0.6494480502272874im,
        -0.8526401859847815+0.5224985598169277im,
        -0.8526401859847815-0.5224985598169277im,
        -0.9238795057196649+0.3826834415328767im,
        -0.9238795057196649-0.3826834415328767im,
        -0.9969173244841298+0.07845911266921719im,
        -0.9969173244841298-0.07845911266921719im,
        -0.97236994351794+0.2334453500964366im,
        -0.97236994351794-0.2334453500964366im]
    # println(ComplexF64([sort(f.p, lt=lt) - sort(Butterworth(BigFloat, 20).p, lt=lt) sort(m_p, lt=lt) - sort(Butterworth(BigFloat, 20).p, lt=lt)]))
    zpkfilter_accuracy(f, ZeroPoleGain(Float64[], m_p, 1), Butterworth(BigFloat, 20); eps=1e-6, relerr=9)
end

#
# Frequency scaling
#
@testset "freq. scaling" begin
    # Output of:
    #=
        [z, p, k] = butter(20, 0.5, 's')
        sprintf('%.16g%+.16gim,', [real(z) imag(z)]')
        sprintf('%.16g%+.16gim,', [real(p) imag(p)]')
        sprintf('%.16g', k)
    =#
    m_z = Float64[]
    m_p = [-0.498458666866501+0.03922954786337574im,
        -0.498458666866501-0.03922954786337574im,
        -0.4861849601989526+0.1167226819283556im,
        -0.4861849601989526-0.1167226819283556im,
        -0.4619397662555457+0.191341716182291im,
        -0.4619397662555457-0.191341716182291im,
        -0.4263200821771486+0.2612492823580798im,
        -0.4263200821771486-0.2612492823580798im,
        -0.3802029827999479+0.3247240241651102im,
        -0.3802029827999479-0.3247240241651102im,
        -0.3247240241650893+0.3802029827999696im,
        -0.3247240241650893-0.3802029827999696im,
        -0.2612492823580033+0.4263200821770542im,
        -0.2612492823580033-0.4263200821770542im,
        -0.1913417161825348+0.461939766255657im,
        -0.1913417161825348-0.461939766255657im,
        -0.1167226819279489+0.486184960198832im,
        -0.1167226819279489-0.486184960198832im,
        -0.03922954786392431+0.4984586668665641im,
        -0.03922954786392431-0.4984586668665641im]
    m_k = 9.5367431640625e-07
    m_f = ZeroPoleGain(m_z, m_p, m_k)

    f = analogfilter(Lowpass(0.5), Butterworth(20))
    zpkfilter_eq(f, m_f)

    # Test that our answers are more accurate than MATLAB's
    accurate_ft = analogfilter(Lowpass(0.5), Butterworth(BigFloat, 20))
    zpkfilter_eq(f, accurate_ft)
    zpkfilter_accuracy(f, m_f, accurate_ft)

    # Output of
    #=
        [z, p, k] = ellip(20, 0.5, 20, 0.3, 's')
        sprintf('%.16g%+.16gim,', [real(z) imag(z)]')
        sprintf('%.16g%+.16gim,', [real(p) imag(p)]')
        sprintf('%.16g', k)
    =#
    m_z = [0+0.6269347368248613im, 0-0.6269347368248613im,
        3.244237104018702e-18+0.3275342448365736im,
        3.244237104018702e-18-0.3275342448365736im,
        -6.938893903907228e-18+0.3032929712822612im,
        -6.938893903907228e-18-0.3032929712822612im,
        -1.387778780781446e-17+0.3004076906327037im,
        -1.387778780781446e-17-0.3004076906327037im,
        -4.163336342344337e-17+0.3000506878920809im,
        -4.163336342344337e-17-0.3000506878920809im,
        6.938893903907228e-18+0.3000063060899686im,
        6.938893903907228e-18-0.3000063060899686im,
        -6.938893903907228e-18+0.3000007854095975im,
        -6.938893903907228e-18-0.3000007854095975im,
        -1.387778780781446e-17+0.30000009864114im,
        -1.387778780781446e-17-0.30000009864114im,
        1.387778780781446e-17+0.3000000132267227im,
        1.387778780781446e-17-0.3000000132267227im,
        2.775557561562891e-17+0.3000000027630185im,
        2.775557561562891e-17-0.3000000027630185im]
    m_p = [-0.1398677695014979+0.188138982052175im,
        -0.1398677695014979-0.188138982052175im,
        -0.02345522399724158+0.2890911154818891im,
        -0.02345522399724158-0.2890911154818891im,
        -0.003009381182870904+0.2987256455960404im,
        -0.003009381182870904-0.2987256455960404im,
        -0.0003757590363774529+0.2998428279105927im,
        -0.0003757590363774529-0.2998428279105927im,
        -4.676623056761031e-05+0.299980469734862im,
        -4.676623056761031e-05-0.299980469734862im,
        -5.818096677288886e-06+0.2999975716762768im,
        -5.818096677288886e-06-0.2999975716762768im,
        -7.237817530730828e-07+0.2999996988517393im,
        -7.237817530730828e-07-0.2999996988517393im,
        -9.003660077222531e-08+0.2999999634684986im,
        -9.003660077222531e-08-0.2999999634684986im,
        -1.117941275663448e-08+0.2999999963792452im,
        -1.117941275663448e-08-0.2999999963792452im,
        -1.220070558471953e-09+0.300000000410979im,
        -1.220070558471953e-09-0.300000000410979im]
    m_k = 0.1000000000276181
    m_f = ZeroPoleGain(m_z, m_p, m_k)

    f = convert(ZeroPoleGain, analogfilter(Lowpass(0.3), Elliptic(20, 0.5, 20)))
    zpkfilter_eq(f, m_f)
end

#
# High pass filter construction
#
@testset "highpass" begin
    # Output of
    #=
        [z, p, k] = butter(20, 0.5, 'high', 's')
        sprintf('%.16g%+.16gim,', [real(z) imag(z)]')
        sprintf('%.16g%+.16gim,', [real(p) imag(p)]')
        sprintf('%.16g', k)
    =#
    m_z = zeros(ComplexF64, 20)
    m_p = [-0.4984586668666513+0.03922954786370267im,
        -0.4984586668666513-0.03922954786370267im,
        -0.4861849601986916+0.1167226819280185im,
        -0.4861849601986916-0.1167226819280185im,
        -0.4619397662557133+0.191341716182565im,
        -0.4619397662557133-0.191341716182565im,
        -0.4263200821770351+0.2612492823579622im,
        -0.4263200821770351-0.2612492823579622im,
        -0.3802029828000168+0.3247240241650909im,
        -0.3802029828000168-0.3247240241650909im,
        -0.3247240241650923+0.3802029828000145im,
        -0.3247240241650923-0.3802029828000145im,
        -0.03922954786392298+0.4984586668665641im,
        -0.03922954786392298-0.4984586668665641im,
        -0.2612492823579744+0.4263200821770468im,
        -0.2612492823579744-0.4263200821770468im,
        -0.1913417161825449+0.4619397662556436im,
        -0.1913417161825449-0.4619397662556436im,
        -0.1167226819279526+0.4861849601988384im,
        -0.1167226819279526-0.4861849601988384im]
    m_k = 1
    m_f = ZeroPoleGain(m_z, m_p, m_k)

    f = analogfilter(Highpass(0.5), Butterworth(20))
    zpkfilter_eq(f, m_f)

    # Test that our answers are more accurate than MATLAB's
    accurate_ft = analogfilter(Highpass(0.5), Butterworth(BigFloat, 20))
    zpkfilter_eq(f, accurate_ft)
    zpkfilter_accuracy(f, m_f, accurate_ft)

    # Output of
    #=
        [z, p, k] = ellip(20, 0.5, 20, 0.4, 'high', 's')
        sprintf('%.16g%+.16gim,', [real(z) imag(z)]')
        sprintf('%.16g%+.16gim,', [real(p) imag(p)]')
        sprintf('%.16g', k)
    =#
    m_z = [4.289315073341702e-17+0.1914074830304432im,
        4.289315073341702e-17-0.1914074830304432im,
        -1.734723475976807e-17+0.3663739040779542im,
        -1.734723475976807e-17-0.3663739040779542im,
        -1.387778780781446e-17+0.3956570424057772im,
        -1.387778780781446e-17-0.3956570424057772im,
        -2.775557561562891e-17+0.3994571502056489im,
        -2.775557561562891e-17-0.3994571502056489im,
        2.775557561562891e-17+0.3999324275609072im,
        2.775557561562891e-17-0.3999324275609072im,
        8.673617379884035e-18+0.3999915920567793im,
        8.673617379884035e-18-0.3999915920567793im,
        4.163336342344337e-17+0.399998952789945im,
        4.163336342344337e-17-0.399998952789945im,
        1.387778780781446e-17+0.3999998684785239im,
        1.387778780781446e-17-0.3999998684785239im,
        -2.775557561562891e-17+0.3999999963159753im,
        -2.775557561562891e-17-0.3999999963159753im,
        6.938893903907228e-18+0.3999999823643708im,
        6.938893903907228e-18-0.3999999823643708im]
    m_p = [-0.3053922020601836+0.4107892635097733im,
        -0.3053922020601836-0.4107892635097733im,
        -0.03345814372376468+0.4123794380387859im,
        -0.03345814372376468-0.4123794380387859im,
        -0.004046405082073851+0.4016656239382425im,
        -0.004046405082073851-0.4016656239382425im,
        -0.0005015366410947575+0.4002090441157958im,
        -0.0005015366410947575-0.4002090441157958im,
        -6.236309209528557e-05+0.400026032326622im,
        -6.236309209528557e-05-0.400026032326622im,
        -7.757587819075262e-06+0.4000032376407234im,
        -7.757587819075262e-06-0.4000032376407234im,
        -9.650442748948551e-07+0.4000004015290894im,
        -9.650442748948551e-07-0.4000004015290894im,
        -1.200488304875513e-07+0.4000000487086385im,
        -1.200488304875513e-07-0.4000000487086385im,
        -1.490588337482723e-08+0.4000000048276729im,
        -1.490588337482723e-08-0.4000000048276729im,
        -1.626761100825824e-09+0.3999999994520274im,
        -1.626761100825824e-09-0.3999999994520274im]
    m_k = 0.9440608762859236
    m_f = ZeroPoleGain(m_z, m_p, m_k)

    f = convert(ZeroPoleGain, analogfilter(Highpass(0.4), Elliptic(20, 0.5, 20)))
    zpkfilter_eq(f, m_f)
end

#
# Band pass filter construction
#
@testset "bandpass" begin
    # Output of
    #=
        [z, p, k] = butter(10, [0.25 0.75], 's')
        sprintf('%.16g%+.16gim,', [real(z) imag(z)]')
        sprintf('%.16g%+.16gim,', [real(p) imag(p)]')
        sprintf('%.16g', k)
    =#
    m_z = zeros(ComplexF64, 10)
    m_p = [-0.05852661385905927+0.7442329445709392im,
        -0.05852661385905927-0.7442329445709392im,
        -0.1665549471789192+0.6992510095482473im,
        -0.1665549471789192-0.6992510095482473im,
        -0.2479888914449238+0.615606020752747im,
        -0.2479888914449238-0.615606020752747im,
        -0.2869902139927978+0.5070583156429217im,
        -0.2869902139927978-0.5070583156429217im,
        -0.2738314211701858+0.3979722116506606im,
        -0.2738314211701858-0.3979722116506606im,
        -0.2200127491273796+0.3197549791305416im,
        -0.2200127491273796-0.3197549791305416im,
        -0.1585130481013887+0.280063065773154im,
        -0.1585130481013887-0.280063065773154im,
        -0.1055644991483508+0.2620526301594698im,
        -0.1055644991483508-0.2620526301594698im,
        -0.06044030269085397+0.2537477474540631im,
        -0.06044030269085397-0.2537477474540631im,
        -0.01969061866105566+0.2503887742733703im,
        -0.01969061866105566-0.2503887742733703im]
    m_k = 0.0009765624999999997
    m_f = ZeroPoleGain(m_z, m_p, m_k)

    f = analogfilter(Bandpass(0.25, 0.75), Butterworth(10))
    zpkfilter_eq(f, m_f)

    # Test that our answers are more accurate than MATLAB's
    accurate_ft = analogfilter(Bandpass(0.25, 0.75), Butterworth(BigFloat, 10))
    zpkfilter_eq(f, accurate_ft)
    zpkfilter_accuracy(f, m_f, accurate_ft)

    # Output of
    #=
        [z, p, k] = ellip(10, 0.5, 20, [0.2 0.8], 's')
        sprintf('%.16g%+.16gim,', [real(z) imag(z)]')
        sprintf('%.16g%+.16gim,', [real(p) imag(p)]')
        sprintf('%.16g', k)
    =#
    m_z = [1.387778780781446e-16+1.370743857088904im,
        1.387778780781446e-16-1.370743857088904im,
        -5.204170427930421e-18+0.1167249440313361im,
        -5.204170427930421e-18-0.1167249440313361im,
        0+0.8445880526167882im, 0-0.8445880526167882im,
        -1.734723475976807e-18+0.8053338099984824im,
        -1.734723475976807e-18-0.8053338099984824im,
        -1.387778780781446e-17+0.1894414673571001im,
        -1.387778780781446e-17-0.1894414673571001im,
        0+0.1986753790956593im, 0-0.1986753790956593im,
        5.551115123125783e-17+0.800148411646084im,
        5.551115123125783e-17-0.800148411646084im,
        1.108054620280186e-16+0.8007109707175198im,
        1.108054620280186e-16-0.8007109707175198im,
        8.673617379884035e-18+0.199822415142662im,
        8.673617379884035e-18-0.199822415142662im,
        4.85722573273506e-17+0.1999629039703328im,
        4.85722573273506e-17-0.1999629039703328im]
    m_p = [-0.2019763467382846+0.6120628823384283im,
        -0.2019763467382846-0.6120628823384283im,
        -0.03721004055421873+0.7823115667177283im,
        -0.03721004055421873-0.7823115667177283im,
        -0.004810677763557675+0.7980130970318471im,
        -0.004810677763557675-0.7980130970318471im,
        -0.0006000695021327641+0.7998051371015736im,
        -0.0006000695021327641-0.7998051371015736im,
        -6.552554548894252e-05+0.8000220660203423im,
        -6.552554548894252e-05-0.8000220660203423im,
        -0.07779252812598642+0.235740074311183im,
        -0.07779252812598642-0.235740074311183im,
        -0.009705975806501373+0.2040604370920958im,
        -0.009705975806501373-0.2040604370920958im,
        -0.001208621817121024+0.2004906765377495im,
        -0.001208621817121024-0.2004906765377495im,
        -0.0001500903998086778+0.2000486149853388im,
        -0.0001500903998086778-0.2000486149853388im,
        -1.638048261950673e-05+0.1999944823054316im,
        -1.638048261950673e-05-0.1999944823054316im]
    m_k = 0.1000000000000025
    m_f = ZeroPoleGain(m_z, m_p, m_k)

    f = convert(ZeroPoleGain, analogfilter(Bandpass(0.2, 0.8), Elliptic(10, 0.5, 20)))
    zpkfilter_eq(f, m_f)
end

#
# Band stop filter construction
#
@testset "bandstop" begin
    # Output of
    #=
        [z, p, k] = butter(10, [0.25 0.75], 'stop', 's')
        sprintf('%.16g%+.16gim,', [real(z) imag(z)]')
        sprintf('%.16g%+.16gim,', [real(p) imag(p)]')
        sprintf('%.16g', k)
    =#
    m_z = [0+0.4330127018922193im, -0-0.4330127018922193im,
        0+0.4330127018922193im, -0-0.4330127018922193im,
        0+0.4330127018922193im, -0-0.4330127018922193im,
        0+0.4330127018922193im, -0-0.4330127018922193im,
        0+0.4330127018922193im, -0-0.4330127018922193im,
        0+0.4330127018922193im, -0-0.4330127018922193im,
        0+0.4330127018922193im, -0-0.4330127018922193im,
        0+0.4330127018922193im, -0-0.4330127018922193im,
        0+0.4330127018922193im, -0-0.4330127018922193im,
        0+0.4330127018922193im, -0-0.4330127018922193im]
    m_p = [-0.05852661385905998+0.7442329445709392im,
        -0.05852661385905998-0.7442329445709392im,
        -0.1665549471789191+0.6992510095482469im,
        -0.1665549471789191-0.6992510095482469im,
        -0.2479888914449218+0.6156060207527428im,
        -0.2479888914449218-0.6156060207527428im,
        -0.2869902139927988+0.5070583156429307im,
        -0.2869902139927988-0.5070583156429307im,
        -0.2738314211701937+0.3979722116506589im,
        -0.2738314211701937-0.3979722116506589im,
        -0.2200127491273721+0.3197549791305355im,
        -0.2200127491273721-0.3197549791305355im,
        -0.1585130481013898+0.2800630657731571im,
        -0.1585130481013898-0.2800630657731571im,
        -0.1055644991483515+0.2620526301594707im,
        -0.1055644991483515-0.2620526301594707im,
        -0.06044030269085383+0.2537477474540633im,
        -0.06044030269085383-0.2537477474540633im,
        -0.01969061866105538+0.2503887742733704im,
        -0.01969061866105538-0.2503887742733704im]
    m_k = 1
    m_f = ZeroPoleGain(m_z, m_p, m_k)

    f = analogfilter(Bandstop(0.25, 0.75), Butterworth(10))
    zpkfilter_eq(f, m_f)

    # Test that our answers are more accurate than MATLAB's
    accurate_ft = analogfilter(Bandstop(0.25, 0.75), Butterworth(BigFloat, 10))
    zpkfilter_eq(f, accurate_ft)
    zpkfilter_accuracy(f, m_f, accurate_ft)

    #=
        [z, p, k] = ellip(10, 0.5, 20, [0.2 0.8], 'stop', 's')
        sprintf('%.16g%+.16gim,', [real(z) imag(z)]')
        sprintf('%.16g%+.16gim,', [real(p) imag(p)]')
        sprintf('%.16g', k)
    =#
    m_z = [-1.402517154063671e-17+0.5685129776423937im,
        -1.402517154063671e-17-0.5685129776423937im,
        -6.156235460644255e-18+0.2814359676774932im,
        -6.156235460644255e-18-0.2814359676774932im,
        -2.775557561562891e-17+0.760016963542146im,
        -2.775557561562891e-17-0.760016963542146im,
        -8.326672684688674e-17+0.7947386855861924im,
        -8.326672684688674e-17-0.7947386855861924im,
        5.551115123125783e-17+0.799851645232653im,
        5.551115123125783e-17-0.799851645232653im,
        8.326672684688674e-17+0.7992903327101535im,
        8.326672684688674e-17-0.7992903327101535im,
        0+0.2105216168522108im, 0-0.2105216168522108im,
        -4.423544863740858e-17+0.2013240363176547im,
        -4.423544863740858e-17-0.2013240363176547im,
        6.938893903907228e-18+0.200037095570968im,
        6.938893903907228e-18-0.200037095570968im,
        5.551115123125783e-17+0.2001775743458425im,
        5.551115123125783e-17-0.2001775743458425im]
    m_p = [-0.3778783653696471+0.7819846151027801im,
        -0.3778783653696471-0.7819846151027801im,
        -0.04044889726976308+0.8144610469677399im,
        -0.04044889726976308-0.8144610469677399im,
        -0.004859687206700769+0.8019364482578354im,
        -0.004859687206700769-0.8019364482578354im,
        -0.0006006726587018951+0.8001940294215477im,
        -0.0006006726587018951-0.8001940294215477im,
        -6.551807382872143e-05+0.799977924147539im,
        -6.551807382872143e-05-0.799977924147539im,
        -0.08015534922525629+0.1658741427311564im,
        -0.08015534922525629-0.1658741427311564im,
        -0.009732315617052817+0.1959655887284332im,
        -0.009732315617052817-0.1959655887284332im,
        -0.00120901710644529+0.1995097303565812im,
        -0.00120901710644529-0.1995097303565812im,
        -0.0001500952639848091+0.1999513917362522im,
        -0.0001500952639848091-0.1999513917362522im,
        -1.638042236460496e-05+0.2000055177738602im,
        -1.638042236460496e-05-0.2000055177738602im]
    m_k = 0.9440608762859234
    m_f = ZeroPoleGain(m_z, m_p, m_k)

    f = analogfilter(Bandstop(0.2, 0.8), Elliptic(10, 0.5, 20))
    zpkfilter_eq(f, m_f)
end

#
# IIR digital filter creation
#
@testset "digital IIR" begin
    # Output of [b, a] = butter(20, 0.5)
    m_b = [8.555155699386467e-06, 0.0001711031139877293,
        0.001625479582883429, 0.009752877497300572,
        0.04144972936352743, 0.1326391339632878,
        0.3315978349082195, 0.663195669816439,
        1.077692963451713, 1.436923951268951,
        1.580616346395846, 1.436923951268951,
        1.077692963451713, 0.663195669816439,
        0.3315978349082195, 0.1326391339632878,
        0.04144972936352743, 0.009752877497300572,
        0.001625479582883429, 0.0001711031139877293,
        8.555155699386467e-06]
    m_a = [1, -9.216626485446244e-15, 2.719352210996272,
        -2.276898740746355e-14, 2.954502383565172,
        -2.181204855972978e-14, 1.663404391798304,
        -9.886771142150971e-15, 0.5272563634534938,
        -1.578399203426284e-15, 0.09588816575455775,
        3.201156931050434e-16, 0.009788250615318749,
        1.651392612127069e-16, 0.0005259945756568433,
        3.034673947897289e-17, 1.306795536212839e-05,
        3.205637858266449e-18, 1.137794123801438e-07,
        9.233765069631315e-20, 1.463813780864114e-10]
    m_f = PolynomialRatio(m_b, m_a)

    f = convert(PolynomialRatio, digitalfilter(Lowpass(0.5), Butterworth(20)))
    tffilter_eq(f, m_f)

    # Output of [b, a] = butter(20, 0.3, 'high')
    m_b = [0.001797932208654781, -0.03595864417309561, 0.3416071196444083,
        -2.04964271786645, 8.710981550932413, -27.87514096298372,
        69.6878524074593, -139.3757048149186, 226.4855203242427,
        -301.980693765657, 332.1787631422226, -301.980693765657,
        226.4855203242427, -139.3757048149186, 69.6878524074593,
        -27.87514096298372, 8.710981550932413, -2.04964271786645,
        0.3416071196444083, -0.03595864417309561, 0.001797932208654781]
    m_a = [1, -7.99378842342338, 32.24776489137069, -86.15266953725099,
        169.613151737271, -260.1393113842562, 321.1838937487486,
        -325.8437287835835, 275.1584488419961, -194.8984651152305,
        116.2180984348268, -58.35420576527059, 24.59604119108612,
        -8.645276014889031, 2.506713994242244, -0.5898102957633149,
        0.1098868139124379, -0.01561185396257133, 0.001590100744530764,
        -0.0001034620065564712, 3.232560226924501e-06]
    m_f = PolynomialRatio(m_b, m_a)

    f = convert(PolynomialRatio, digitalfilter(Highpass(0.3), Butterworth(20)))
    tffilter_eq(f, m_f)

    # Output of [b, a] = butter(7, [0.1 0.4])
    m_b = [0.0009628940476886417, 0, -0.006740258333820491, 0,
        0.02022077500146148, 0, -0.03370129166910246,
        0, 0.03370129166910246, 0, -0.02022077500146148,
        0, 0.006740258333820491, 0, -0.0009628940476886417]
    m_a = [1, -7.763446908073567, 29.02150968091795, -69.51575164764138,
        119.5092528789707, -156.1341322158939, 159.8101214681526,
        -130.0847734566226, 84.56503905088692, -43.67462798637182,
        17.64013826588496, -5.405682857208016, 1.18954580728982,
        -0.1686342472129378, 0.01166272804918956]
    m_f = PolynomialRatio(m_b, m_a)

    f = convert(PolynomialRatio, digitalfilter(Bandpass(0.1, 0.4), Butterworth(7)))
    tffilter_eq(f, m_f)

    # Output of [b, a] = butter(4, [0.2 0.25], 'stop')
    m_b = [0.8142545568862771, -4.968628804996451, 14.62659178092118,
        -26.46885325790832, 32.03454308505385, -26.46885325790833,
        14.62659178092118, -4.968628804996451, 0.8142545568862771]
    m_a = [1, -5.789124865012897, 16.16849931613174, -27.77244867328942,
        31.92014561624699, -25.06018377653858, 13.16458034390183,
        -4.253206810966346, 0.6630104843858931]
    m_f = PolynomialRatio(m_b, m_a)

    f = convert(PolynomialRatio, digitalfilter(Bandstop(0.2, 0.25), Butterworth(4)))
    tffilter_eq(f, m_f)

    # Output of [b, a] = ellip(10, 0.7, 13, 0.4)
    m_b = [0.2797942479261318, -0.546743606764335, 1.704005589668604,
        -2.159911306584612, 3.642222926761125, -3.210370872716824,
        3.642222926761128, -2.159911306584614, 1.704005589668607,
        -0.5467436067643362, 0.2797942479261327]
    m_a = [1, -3.029857317356546, 7.982710263075236, -12.64912738840943,
        17.41532979440429, -17.06917827530357, 14.5191214973975,
        -8.798609366466101, 4.541539167299845, -1.419553904214128, 0.3565809079146795]
    m_f = PolynomialRatio(m_b, m_a)

    f = convert(PolynomialRatio, digitalfilter(Lowpass(0.4), Elliptic(10, 0.7, 13)))
    tffilter_eq(f, m_f)

    # Output of [b, a] = ellip(11, 0.7, 17, 0.2, 'high')
    m_b = [0.5737867105773529, -5.288884465999311, 23.07976318423696,
        -62.6831981311594, 117.4527496444787, -159.209133733912,
        159.2091337339121, -117.4527496444788, 62.68319813115945,
        -23.079763184237, 5.288884465999319, -0.5737867105773541]
    m_a = [1, -8.215565996888206, 32.26875224082856, -79.26375768096621,
        134.6328215735552, -165.4634116073934, 149.7532535390064,
        -99.57881504158784, 47.56023313164873, -15.48547977535943,
        3.075141867250689, -0.2777992862430019]
    m_f = PolynomialRatio(m_b, m_a)

    f = convert(PolynomialRatio, digitalfilter(Highpass(0.2), Elliptic(11, 0.7, 17)))
    tffilter_eq(f, m_f)

    # Output of [b, a] = ellip(10, 0.5, 2, [0.2 0.3])
    m_b = [0.7571867287399742, -10.55194327196835, 73.54073467886568,
        -338.6101060102045, 1150.502868412841, -3057.625219554799,
        6581.673203566225, -11732.90330745444, 17574.06211811366,
        -22320.04109078996, 24158.39342879073, -22320.04109078995,
        17574.06211811365, -11732.90330745443, 6581.673203566216,
        -3057.625219554793, 1150.502868412839, -338.6101060102037,
        73.54073467886548, -10.55194327196831, 0.7571867287399715]
    m_a = [1, -13.88122062422027, 96.34606972585827, -441.7050045576331,
        1494.048255114342, -3952.085137071483, 8465.769048060451,
        -15015.81246941687, 22374.64201809294, -28265.00095948501,
        30424.52053126076, -27950.27646663918, 21879.36223420588,
        -14520.36940748522, 8095.785094725445, -3737.670491911818,
        1397.479471832041, -408.6480215434079, 88.17084042962101,
        -12.56723162073412, 0.8957646298422474]
    m_f = PolynomialRatio(m_b, m_a)

    f = convert(PolynomialRatio, digitalfilter(Bandpass(0.2, 0.3), Elliptic(10, 0.5, 2)))
    tffilter_eq(f, m_f)

    # Output of [b, a] = ellip(15, 1.8, 52, [0.6 0.7], 'stop')
    m_b = [0.4754470525338341, 6.428928219976483, 47.47972043083799,
        245.8985324356427, 989.107195159456, 3259.995646867947,
        9099.649115772192, 21987.91200123768, 46709.48715187636,
        88217.96445932264, 149368.0297431405, 228141.9056208539,
        315792.8309895032, 397458.2884412499, 455878.980369108,
        477133.9777299722, 455878.9803691082, 397458.2884412502,
        315792.8309895035, 228141.9056208542, 149368.0297431408,
        88217.96445932281, 46709.48715187647, 21987.91200123774,
        9099.649115772219, 3259.995646867959, 989.1071951594597,
        245.8985324356436, 47.4797204308382, 6.428928219976513,
        0.475447052533836]
    m_a = [1, 12.89634891736776, 90.77856214451822, 447.9765044361528,
        1716.574771803699, 5388.380608104304, 14321.19101752449,
        32940.21064000877, 66586.59389679749, 119619.3933180238,
        192555.4138276511, 279453.5193544072, 367301.6888247615,
        438620.3975638363, 476900.358709065, 472649.3359368506,
        427104.9678341324, 351676.250633706, 263451.4609235193,
        179107.598840823, 110104.024869051, 60897.74794125407,
        30103.22741782092, 13182.01328029648, 5052.946545447418,
        1668.118080642225, 463.5003850459777, 104.731419523004,
        18.20949251383043, 2.194519514615899, 0.1423868039113367]
    m_f = PolynomialRatio(m_b, m_a)

    f = convert(PolynomialRatio, digitalfilter(Bandstop(0.6, 0.7), Elliptic(15, 1.8, 52)))
    tffilter_eq(f, m_f)
end

#
# Special filters
#
@testset "special" begin
    # Output of [num, den] = iirnotch(0.1, 0.1)
    notch = iirnotch(0.1, 0.1)
    @test notch.b0 ≈ 0.863271264002681
    @test notch.b1 ≈ -1.642039521920206
    @test notch.b2 ≈ 0.863271264002681
    @test notch.a1 ≈ -1.642039521920206
    @test notch.a2 ≈ 0.726542528005361

    notch = iirnotch(60, 1, fs=1000)
    @test notch.b0 ≈ 0.996868235770807
    @test notch.b1 ≈ -1.853729290297204
    @test notch.b2 ≈ 0.996868235770807
    @test notch.a1 ≈ -1.853729290297204
    @test notch.a2 ≈ 0.993736471541615
end

#
# IIR digital filter error conditions
#
@testset "errors" begin
    # zero pole
    @test_throws ErrorException Butterworth(0)
    @test_throws ErrorException Chebyshev1(0, 1)
    @test_throws ErrorException Chebyshev2(0, 1)
    @test_throws ErrorException Elliptic(0, 0.1, 0.5)

    # too much passband ripple
    @test_throws ErrorException Elliptic(2, 0.5, 0.1)

    # negative ripple
    @test_throws ErrorException Chebyshev1(2, -1)
    @test_throws ErrorException Chebyshev2(2, -1)
    @test_throws ErrorException Elliptic(2, 0.1, -0.5)
end

#
# Specification of sample rate
#
@testset "sample rate" begin
    @test Lowpass(20, fs=100) == Lowpass(0.4)
    @test Highpass(20, fs=100) == Highpass(0.4)
    @test Bandpass(20, 21, fs=100) == Bandpass(0.4, 0.42)
    @test Bandstop(20, 21, fs=100) == Bandstop(0.4, 0.42)
end

#
# kaiserord
#
@testset "kaiserord" begin
    kaiserord_jl    = kaiserord(0.1, 60)
    kaiserord_scipy = (74, 1.7994885471673767)    # (n,beta)=kaiserord(60,.1); (n,beta/pi)
    @test kaiserord_jl[1] == kaiserord_scipy[1]
    @test kaiserord_jl[2] ≈ kaiserord_scipy[2]

    kaiserord_jl    = kaiserord(0.1, 21.0 + eps(21.0))
    kaiserord_scipy = (20, 3.087704731700782e-07) # (n,beta)=kaiserord(21.000000000000004,.1); (n,beta/pi)
    @test kaiserord_jl[1] == kaiserord_scipy[1]
    @test kaiserord_jl[2] ≈ kaiserord_scipy[2]

    kaiserord_jl    = kaiserord(0.1, 15)
    kaiserord_scipy = (11, 0.0)                   # (n,beta)=kaiserord(15,.1); (n,beta/pi)
    @test kaiserord_jl[1] == kaiserord_scipy[1]
    @test kaiserord_jl[2] ≈ kaiserord_scipy[2]
end

#
# Window FIR filter taps
#
@testset "window FIR" begin
    winfirtaps_jl    = digitalfilter(Lowpass(0.25; fs=1),FIRWindow(hamming(128), scale=false))
    # firwin(128, 0.25, nyq=.5, scale=False)
    winfirtaps_scipy = readdlm(joinpath(dirname(@__FILE__), "data", "digitalfilter_hamming_128_lowpass_fc0.25_fs1.0.txt"),'\t')
    @test winfirtaps_jl ≈ winfirtaps_scipy

    winfirtaps_jl    = digitalfilter(Lowpass(0.25; fs=1),FIRWindow(hamming(129), scale=false))
    # firwin(129, 0.25, nyq=.5, scale=False)
    winfirtaps_scipy = readdlm(joinpath(dirname(@__FILE__), "data", "digitalfilter_hamming_129_lowpass_fc0.25_fs1.0.txt"),'\t')
    @test winfirtaps_jl ≈ winfirtaps_scipy

    @test_throws ArgumentError digitalfilter(Highpass(0.25; fs=1),FIRWindow(hamming(128), scale=false))

    winfirtaps_jl    = digitalfilter(Highpass(0.25; fs=1),FIRWindow(hamming(129), scale=false))
    # firwin(129, 0.25, nyq=.5, scale=False, pass_zero=False)
    winfirtaps_scipy = readdlm(joinpath(dirname(@__FILE__), "data", "digitalfilter_hamming_129_highpass_fc0.25_fs1.0.txt"),'\t')
    @test winfirtaps_jl ≈ winfirtaps_scipy

    winfirtaps_jl    = digitalfilter(Bandpass(0.1, 0.2; fs=1),FIRWindow(hamming(128), scale=false))
    # firwin(128, [0.1, 0.2], nyq=.5, scale=False, pass_zero=False)
    winfirtaps_scipy = readdlm(joinpath(dirname(@__FILE__), "data", "digitalfilter_hamming_128_bandpass_fc0.1_0.2_fs1.0.txt"),'\t')
    @test winfirtaps_jl ≈ winfirtaps_scipy

    winfirtaps_jl    = digitalfilter(Bandpass(0.1, 0.2; fs=1),FIRWindow(hamming(129), scale=false))
    # firwin(129, [0.1, 0.2], nyq=.5, scale=False, pass_zero=False)
    winfirtaps_scipy = readdlm(joinpath(dirname(@__FILE__), "data", "digitalfilter_hamming_129_bandpass_fc0.1_0.2_fs1.0.txt"),'\t')
    @test winfirtaps_jl ≈ winfirtaps_scipy

    @test_throws ArgumentError digitalfilter(Bandstop(0.1, 0.2; fs=1),FIRWindow(hamming(128), scale=false))

    winfirtaps_jl    = digitalfilter(Bandstop(0.1, 0.2; fs=1),FIRWindow(hamming(129), scale=false))
    # firwin(129, [0.1, 0.2], nyq=.5, scale=False)
    winfirtaps_scipy = readdlm(joinpath(dirname(@__FILE__), "data", "digitalfilter_hamming_129_bandstop_fc0.1_0.2_fs1.0.txt"),'\t')
    @test winfirtaps_jl ≈ vec(winfirtaps_scipy)

    winfirtaps_jl    = digitalfilter(Lowpass(0.25; fs=1),FIRWindow(hamming(128), scale=true))
    # firwin(128, 0.25, nyq=.5)
    winfirtaps_scipy = readdlm(joinpath(dirname(@__FILE__), "data", "digitalfilter_hamming_128_lowpass_scaled_fc0.25_fs1.0.txt"),'\t')
    @test winfirtaps_jl ≈ winfirtaps_scipy

    winfirtaps_jl    = digitalfilter(Lowpass(0.25; fs=1),FIRWindow(hamming(129), scale=true))
    # firwin(129, 0.25, nyq=.5)
    winfirtaps_scipy = readdlm(joinpath(dirname(@__FILE__), "data", "digitalfilter_hamming_129_lowpass_scaled_fc0.25_fs1.0.txt"),'\t')
    @test winfirtaps_jl ≈ winfirtaps_scipy

    winfirtaps_jl    = digitalfilter(Highpass(0.25; fs=1),FIRWindow(hamming(129), scale=true))
    # firwin(129, 0.25, nyq=.5, pass_zero=False)
    winfirtaps_scipy = readdlm(joinpath(dirname(@__FILE__), "data", "digitalfilter_hamming_129_highpass_scaled_fc0.25_fs1.0.txt"),'\t')
    @test winfirtaps_jl ≈ winfirtaps_scipy

    winfirtaps_jl    = digitalfilter(Bandpass(0.1, 0.2; fs=1),FIRWindow(hamming(128), scale=true))
    # firwin(128, [0.1, 0.2], nyq=.5, pass_zero=False)
    winfirtaps_scipy = readdlm(joinpath(dirname(@__FILE__), "data", "digitalfilter_hamming_128_bandpass_scaled_fc0.1_0.2_fs1.0.txt"),'\t')
    @test winfirtaps_jl ≈ winfirtaps_scipy

    winfirtaps_jl    = digitalfilter(Bandpass(0.1, 0.2; fs=1),FIRWindow(hamming(129), scale=true))
    # firwin(129, [0.1, 0.2], nyq=.5, scale=False, pass_zero=False)
    winfirtaps_scipy = readdlm(joinpath(dirname(@__FILE__), "data", "digitalfilter_hamming_129_bandpass_scaled_fc0.1_0.2_fs1.0.txt"),'\t')
    @test winfirtaps_jl ≈ winfirtaps_scipy

    winfirtaps_jl    = digitalfilter(Bandstop(0.1, 0.2; fs=1),FIRWindow(hamming(129), scale=true))
    # firwin(129, [0.1, 0.2], nyq=.5, scale=False)
    winfirtaps_scipy = readdlm(joinpath(dirname(@__FILE__), "data", "digitalfilter_hamming_129_bandstop_scaled_fc0.1_0.2_fs1.0.txt"),'\t')
    @test winfirtaps_jl ≈ vec(winfirtaps_scipy)
end

# make sure typing is preserved
# ie, user passes in F32, should get ZPK{CF32, CF32, F32}
@testset "typing" begin
    @testset "prototypes" begin
        @test_throws MethodError Butterworth(ComplexF64, 4)
        @test_throws MethodError Chebyshev1(ComplexF64, 4, 1.01)
        @test_throws MethodError Chebyshev2(ComplexF64, 4, 1.01)
        @test_throws MethodError Elliptic(ComplexF64, 4, 0.5, 1.01)

        @test_throws MethodError Butterworth(Char, 4)
        @test_throws MethodError Chebyshev1(Char, 4, 1.01)
        @test_throws MethodError Chebyshev2(Char, 4, 1.01)
        @test_throws MethodError Elliptic(Char, 4, 0.5, 1.01)

        @testset "$T" for T in [BigFloat, Float16, Float32, Float64]
            @test Butterworth(T, 4) isa ZeroPoleGain{:s, Complex{T}, Complex{T}, T}
            @test Chebyshev1(T, 5, 1.16) isa ZeroPoleGain{:s, Complex{T}, Complex{T}, T}
            @test Chebyshev2(T, 5, 1.16) isa ZeroPoleGain{:s, Complex{T}, Complex{T}, T}
            @test Elliptic(T, 5, 0.31, 1.16) isa ZeroPoleGain{:s, Complex{T}, Complex{T}, T}
        end
    end

    @testset "transformed types" begin
        @testset "$T & $R" for T in [BigFloat, Float16, Float32, Float64],
                R in [Lowpass(0.4), Highpass(0.6), Bandpass(0.3, 0.6), Bandstop(0.7, 0.8)]
            @test digitalfilter(R, Butterworth(T, 5)) isa ZeroPoleGain{:z, Complex{T}, Complex{T}, T}
            @test digitalfilter(R, Chebyshev1(T, 5, 1.16)) isa ZeroPoleGain{:z, Complex{T}, Complex{T}, T}
            @test digitalfilter(R, Chebyshev2(T, 5, 1.16)) isa ZeroPoleGain{:z, Complex{T}, Complex{T}, T}
            @test digitalfilter(R, Elliptic(T, 5, 0.31, 1.16)) isa ZeroPoleGain{:z, Complex{T}, Complex{T}, T}
        end
    end
end
