diff --git a/src/Filters/design.jl b/src/Filters/design.jl index df485b3f..b733edab 100644 --- a/src/Filters/design.jl +++ b/src/Filters/design.jl @@ -500,7 +500,7 @@ prewarp(ftype::Highpass, fs::Real) = Highpass(prewarp(normalize_freq(ftype.w, fs prewarp(ftype::Bandpass, fs::Real) = Bandpass(prewarp(normalize_freq(ftype.w1, fs)), prewarp(normalize_freq(ftype.w2, fs))) prewarp(ftype::Bandstop, fs::Real) = Bandstop(prewarp(normalize_freq(ftype.w1, fs)), prewarp(normalize_freq(ftype.w2, fs))) # freq in half-samples per cycle -prewarp(f::Real) = 4*tan(pi*f/2) +prewarp(f::Real) = 4 * tanpi(f / 2) # Digital filter design """ @@ -531,7 +531,7 @@ function iirnotch(w::Real, bandwidth::Real; fs=2) bandwidth = normalize_freq(bandwidth, fs) # Eq. 8.2.23 - b = 1/(1+tan(pi*bandwidth/2)) + b = 1 / (1 + tanpi(bandwidth / 2)) # Eq. 8.2.22 cosw0 = cospi(w) Biquad(b, -2b*cosw0, b, -2b*cosw0, 2b-1) diff --git a/src/Filters/filt_order.jl b/src/Filters/filt_order.jl index 5f9e713a..8d22d6ef 100644 --- a/src/Filters/filt_order.jl +++ b/src/Filters/filt_order.jl @@ -196,35 +196,36 @@ end # # Bandstop cost functions and passband minimization # -for filt in (:butterworth, :elliptic, :chebyshev) - @eval begin - function $(Symbol(filt, :_bsfcost))(Wx::Real, uselowband::Bool, Wp::Tuple{Real,Real}, Ws::Tuple{Real,Real}, Rp::Real, Rs::Real) - # override one of the passband edges with the test frequency. - Wpc = uselowband ? (Wx, Wp[2]) : (Wp[1], Wx) +function bsfcost(est_func::F, Wx::Real, uselowband::Bool, Wp::Tuple{Real,Real}, Ws::Tuple{Real,Real}, Rp::Real, Rs::Real) where {F} + # override one of the passband edges with the test frequency. + Wpc = uselowband ? (Wx, Wp[2]) : (Wp[1], Wx) - # get the new warp frequency. - warp = min(abs.((Ws .* (Wpc[1] - Wpc[2])) ./ muladd.(-Wpc[1], Wpc[2], Ws .^ 2))...) + # get the new warp frequency. + warp = min(abs.((Ws .* (Wpc[1] - Wpc[2])) ./ muladd.(-Wpc[1], Wpc[2], Ws .^ 2))...) - # use the new frequency to determine the filter order. - $(Symbol(filt, :_order_estimate))(Rp, Rs, warp) - end + # use the new frequency to determine the filter order. + est_func(Rp, Rs, warp) +end - function $(Symbol(filt, :_bsfmin))(Wp::Tuple{Real,Real}, Ws::Tuple{Real,Real}, Rp::Real, Rs::Real) - # NOTE: the optimization function will adjust the corner frequencies in Wp, returning a new passband tuple. - Δ = eps(typeof(Wp[1]))^(2 / 3) - C₁(w) = $(Symbol(filt, :_bsfcost))(w, true, Wp, Ws, Rp, Rs) - p1 = brent(C₁, Wp[1], Ws[1] - Δ) +function bsfmin(est_func::F, Wp::Tuple{Real,Real}, Ws::Tuple{Real,Real}, Rp::Real, Rs::Real) where {F} + @inline + # NOTE: the optimization function will adjust the corner frequencies in Wp, returning a new passband tuple. + Δ = eps(typeof(Wp[1]))^(2 / 3) + C₁(w) = bsfcost(est_func, w, true, Wp, Ws, Rp, Rs) + p1 = brent(C₁, Wp[1], Ws[1] - Δ) - C₂(w) = $(Symbol(filt, :_bsfcost))(w, false, (p1, Wp[2]), Ws, Rp, Rs) - p2 = brent(C₂, Ws[2] + Δ, Wp[2]) - Wadj = (p1, p2) + C₂(w) = bsfcost(est_func, w, false, (p1, Wp[2]), Ws, Rp, Rs) + p2 = brent(C₂, Ws[2] + Δ, Wp[2]) + Wadj = (p1, p2) - Wa = (Ws .* (p1 - p2)) ./ muladd.(-p1, p2, Ws .^ 2) - min(abs.(Wa)...), Wadj - end - end + Wa = (Ws .* (p1 - p2)) ./ muladd.(-p1, p2, Ws .^ 2) + min(abs.(Wa)...), Wadj end +butterworth_bsfmin(args...) = bsfmin(butterworth_order_estimate, args...) +elliptic_bsfmin(args...) = bsfmin(elliptic_order_estimate, args...) +chebyshev_bsfmin(args...) = bsfmin(chebyshev_order_estimate, args...) + """ (N, ωn) = buttord(Wp::Tuple{Real,Real}, Ws::Tuple{Real,Real}, Rp::Real, Rs::Real; domain=:z) @@ -250,8 +251,8 @@ function buttord(Wp::Tuple{Real,Real}, Ws::Tuple{Real,Real}, Rp::Real, Rs::Real; # pre-warp both components, (if Z-domain specified) if (domain == :z) - Ωp = tan.(π / 2 .* Wps) - Ωs = tan.(π / 2 .* Wss) + Ωp = tanpi.(Wps ./ 2) + Ωs = tanpi.(Wss ./ 2) else Ωp = Wps Ωs = Wss @@ -298,8 +299,8 @@ function buttord(Wp::Real, Ws::Real, Rp::Real, Rs::Real; domain::Symbol=:z) if (domain == :z) # need to pre-warp since we want to use formulae for analog case. - Ωp = tan(π / 2 * Wp) - Ωs = tan(π / 2 * Ws) + Ωp = tanpi(Wp / 2) + Ωs = tanpi(Ws / 2) else # already analog Ωp = Wp @@ -327,76 +328,76 @@ end # # Elliptic/Chebyshev1 Estimation # -for (fcn, est, filt) in ((:ellipord, :elliptic, "Elliptic (Cauer)"), - (:cheb1ord, :chebyshev, "Chebyshev Type I")) +function ordfreq_est(order_estimate, domain::Symbol, Wp::Real, Ws::Real, Rp::Real, Rs::Real) + @inline + ftype = (Wp < Ws) ? Lowpass : Highpass + if (domain == :z) + Ωp = tanpi(Wp / 2) + Ωs = tanpi(Ws / 2) + else + Ωp = Wp + Ωs = Ws + end + # Lowpass/Highpass prototype xform is same as Butterworth. + wa = toprototype(Ωp, Ωs, ftype) + N = ceil(Int, order_estimate(Rp, Rs, wa)) + if (domain == :z) + ωn = (2 / π) * atan(Ωp) + else + ωn = Ωp + end + N, ωn +end + +function ordfreq_est(order_estimate::F, domain::Symbol, + Wp::Tuple{Real,Real}, Ws::Tuple{Real,Real}, Rp::Real, Rs::Real) where {F} + @inline + Wps = sort_W(Wp) + Wss = sort_W(Ws) + (Wps[1] < Wss[1]) != (Wps[2] > Wss[2]) && throw(ArgumentError("Pass and stopband edges must be ordered for Bandpass/Bandstop filters.")) + ftype = (Wps[1] < Wss[1]) ? Bandstop : Bandpass + # pre-warp to analog if z-domain. + (Ωp, Ωs) = (domain == :z) ? (tanpi.(Wps ./ 2), tanpi.(Wss ./ 2)) : (Wps, Wss) + if (ftype == Bandpass) + Wa = muladd.(-Ωp[1], Ωp[2], Ωs .^ 2) ./ (Ωs .* (Ωp[1] - Ωp[2])) + Ωpadj = Ωp + else + (Wa, Ωpadj) = bsfmin(order_estimate, Ωp, Ωs, Rp, Rs) # check scipy. + end + N = ceil(Int, order_estimate(Rp, Rs, min(abs.(Wa)...))) + ωn = (domain == :z) ? Wps : Ωpadj + N, ωn +end + +ellipord(args...; domain::Symbol=:z) = ordfreq_est(elliptic_order_estimate, domain, args...) +cheb1ord(args...; domain::Symbol=:z) = ordfreq_est(chebyshev_order_estimate, domain, args...) + +for (fcn, filt) in ((:ellipord, "Elliptic (Cauer)"), (:cheb1ord, "Chebyshev Type I")) @eval begin """ (N, ωn) = $($fcn)(Wp::Real, Ws::Real, Rp::Real, Rs::Real; domain::Symbol=:z) + (N, ωn) = $($fcn)(Wp::Tuple{Real, Real}, Ws::Tuple{Real, Real}, Rp::Real, Rs::Real; domain::Symbol=:z) Integer and natural frequency order estimation for $($filt) Filters. `Wp` and `Ws` indicate the passband and stopband frequency edges, and `Rp` and `Rs` indicate the - maximum loss in the passband and the minimum attenuation in the stopband, (in dB.)\n + maximum loss in the passband and the minimum attenuation in the stopband, (in dB.) + Based on the ordering of passband and stopband edges, the Lowpass or Highpass filter type is inferred. `N` indicates the smallest integer filter order that achieves the desired specifications, and `ωn` contains the natural frequency of the filter, (in this case, simply the passband edge.)\n If a domain of `:s` is specified, the passband and stopband edges are interpreted - as radians/second, giving an order and natural frequency result for an analog filter. + as radians/second, giving the order and natural frequencies for an analog filter. The default domain is `:z`, interpreting the input frequencies as normalized from 0 to 1, where 1 corresponds to π radians/sample. - """ - function $fcn(Wp::Real, Ws::Real, Rp::Real, Rs::Real; domain::Symbol=:z) - ftype = (Wp < Ws) ? Lowpass : Highpass - if (domain == :z) - Ωp = tan(π / 2 * Wp) - Ωs = tan(π / 2 * Ws) - else - Ωp = Wp - Ωs = Ws - end - # Lowpass/Highpass prototype xform is same as Butterworth. - wa = toprototype(Ωp, Ωs, ftype) - N = ceil(Int, $(Symbol(est, :_order_estimate))(Rp, Rs, wa)) - if (domain == :z) - ωn = (2 / π) * atan(Ωp) - else - ωn = Ωp - end - N, ωn - end + `Wp` and `Ws` can also be 2-element frequency edges for Bandpass/Bandstop cases.\n + `N` is an integer indicating the lowest estimated filter order, with + `ωn` specifying the cutoff or "-3 dB" frequencies. """ - (N, ωn) = $($fcn)(Wp::Tuple{Real, Real}, Ws::Tuple{Real, Real}, Rp::Real, Rs::Real; domain::Symbol=:z) - - Integer and natural frequency order estimation for $($filt) Filters. `Wp` and `Ws` are 2-element - frequency edges for Bandpass/Bandstop cases, with `Rp` and `Rs` representing the ripple maximum loss - in the passband and minimum ripple attenuation in the stopband in dB.\nBased on the ordering of passband - and bandstop edges, the Bandstop or Bandpass filter type is inferred. `N` is an integer indicating the - lowest estimated filter order, with `ωn` specifying the cutoff or "-3 dB" frequencies.\nIf a domain of - `:s` is specified, the passband and stopband frequencies are interpreted as radians/second, giving the - order and natural frequencies for an analog filter. The default domain is `:z`, interpreting the input - frequencies as normalized from 0 to 1, where 1 corresponds to π radians/sample. - """ - function $fcn(Wp::Tuple{Real,Real}, Ws::Tuple{Real,Real}, Rp::Real, Rs::Real; domain::Symbol=:z) - Wps = sort_W(Wp) - Wss = sort_W(Ws) - (Wps[1] < Wss[1]) != (Wps[2] > Wss[2]) && throw(ArgumentError("Pass and stopband edges must be ordered for Bandpass/Bandstop filters.")) - ftype = (Wps[1] < Wss[1]) ? Bandstop : Bandpass - # pre-warp to analog if z-domain. - (Ωp, Ωs) = (domain == :z) ? (tan.(π / 2 .* Wps), tan.(π / 2 .* Wss)) : (Wps, Wss) - if (ftype == Bandpass) - Wa = muladd.(-Ωp[1], Ωp[2], Ωs .^ 2) ./ (Ωs .* (Ωp[1] - Ωp[2])) - Ωpadj = Ωp - else - (Wa, Ωpadj) = $(Symbol(est, :_bsfmin))(Ωp, Ωs, Rp, Rs) # check scipy. - end - N = ceil(Int, $(Symbol(est, :_order_estimate))(Rp, Rs, min(abs.(Wa)...))) - ωn = (domain == :z) ? Wps : Ωpadj - N, ωn - end + $fcn end end - """ (N, ωn) = cheb2ord(Wp::Real, Ws::Real, Rp::Real, Rs::Real; domain::Symbol=:z) @@ -413,7 +414,7 @@ the input frequencies as normalized from 0 to 1, where 1 corresponds to π radia """ function cheb2ord(Wp::Real, Ws::Real, Rp::Real, Rs::Real; domain::Symbol=:z) ftype = (Wp < Ws) ? Lowpass : Highpass - (Ωp, Ωs) = (domain == :z) ? (tan(π / 2 * Wp), tan(π / 2 * Ws)) : (Wp, Ws) + (Ωp, Ωs) = (domain == :z) ? (tanpi(Wp / 2), tanpi(Ws / 2)) : (Wp, Ws) wa = toprototype(Ωp, Ωs, ftype) N = ceil(Int, chebyshev_order_estimate(Rp, Rs, wa)) @@ -444,7 +445,7 @@ function cheb2ord(Wp::Tuple{Real,Real}, Ws::Tuple{Real,Real}, Rp::Real, Rs::Real Wss = sort_W(Ws) (Wps[1] < Wss[1]) != (Wps[2] > Wss[2]) && throw(ArgumentError("Pass and stopband edges must be ordered for Bandpass/Bandstop filters.")) ftype = (Wps[1] < Wss[1]) ? Bandstop : Bandpass - (Ωp, Ωs) = (domain == :z) ? (tan.(π / 2 .* Wps), tan.(π / 2 .* Wss)) : (Wps, Wss) + (Ωp, Ωs) = (domain == :z) ? (tanpi.(Wps ./ 2), tanpi.(Wss ./ 2)) : (Wps, Wss) if (ftype == Bandpass) prod = Ωp[1] * Ωp[2] diff = Ωp[1] - Ωp[2] diff --git a/src/windows.jl b/src/windows.jl index 6311067b..4b3b7272 100644 --- a/src/windows.jl +++ b/src/windows.jl @@ -632,33 +632,31 @@ const IntegerOr2 = Union{Tuple{Integer, Integer}, Integer} const RealOr2 = Union{Tuple{Real, Real}, Real} const BoolOr2 = Union{Tuple{Bool, Bool}, Bool} +function matrix_window(func::F, dims::Tuple{Integer,Integer}, arg::Union{RealOr2,Nothing}=nothing; + padding::IntegerOr2=0, zerophase::BoolOr2=false) where {F} + paddings = argdup(padding) + zerophases = argdup(zerophase) + if isnothing(arg) + w1 = func(dims[1]; padding=paddings[1], zerophase=zerophases[1]) + w2 = func(dims[2]; padding=paddings[2], zerophase=zerophases[2]) + else + args = argdup(arg) + w1 = func(dims[1], args[1]; padding=paddings[1], zerophase=zerophases[1]) + w2 = func(dims[2], args[2]; padding=paddings[2], zerophase=zerophases[2]) + end + return w1 * w2' +end + for func in (:rect, :hanning, :hamming, :cosine, :lanczos, :triang, :bartlett, :bartlett_hann, :blackman) - @eval begin - function $func(dims::Tuple; padding::IntegerOr2=0, - zerophase::BoolOr2=false) - length(dims) == 2 || throw(ArgumentError("`dims` must be length 2")) - paddings = argdup(padding) - zerophases = argdup(zerophase) - w1 = $func(dims[1]; padding=paddings[1], zerophase=zerophases[1]) - w2 = $func(dims[2]; padding=paddings[2], zerophase=zerophases[2]) - w1 * w2' - end + @eval function $func(dims::Tuple{Integer,Integer}; padding::IntegerOr2=0, zerophase::BoolOr2=false) + return matrix_window($func, dims; padding, zerophase) end end for func in (:tukey, :gaussian, :kaiser) - @eval begin - function $func(dims::Tuple, arg::RealOr2; - padding::IntegerOr2=0, zerophase::BoolOr2=false) - length(dims) == 2 || throw(ArgumentError("`dims` must be length 2")) - args = argdup(arg) - paddings = argdup(padding) - zerophases = argdup(zerophase) - w1 = $func(dims[1], args[1]; padding=paddings[1], zerophase=zerophases[1]) - w2 = $func(dims[2], args[2]; padding=paddings[2], zerophase=zerophases[2]) - w1 * w2' - end + @eval function $func(dims::Tuple{Integer,Integer}, arg::RealOr2; padding::IntegerOr2=0, zerophase::BoolOr2=false) + return matrix_window($func, dims, arg; padding, zerophase) end end diff --git a/test/filt_order.jl b/test/filt_order.jl index 3ee67039..f4c8faad 100644 --- a/test/filt_order.jl +++ b/test/filt_order.jl @@ -179,7 +179,7 @@ end # Highpass (nh, Wnh) = cheb2ord(0.21, 0.1, Rp, Rs, domain=:z) @test nh == 8 - @test Wnh == 0.10862150541420543 + @test Wnh == 0.10862150541420544 (nh, Wnh) = cheb2ord(0.21, 0.1, Rp, Rs, domain=:s) @test nh == 8 @test Wnh == 0.10568035411923006 @@ -219,8 +219,8 @@ end # # Using the test-cases highlighted in [^Parks] Figures 8 and 15. # - # [^Parks]: Rabiner, L. R., McClellan, J. H., & Parks, T. W. (1975). - # FIR digital filter design techniques using weighted Chebyshev + # [^Parks]: Rabiner, L. R., McClellan, J. H., & Parks, T. W. (1975). + # FIR digital filter design techniques using weighted Chebyshev # approximation. Proceedings of the IEEE, 63(4), 595-610. # @test remezord(0.41665, 0.49417, 0.0116, 0.0001) == 39 diff --git a/test/windows.jl b/test/windows.jl index 34db2ff6..2473e018 100644 --- a/test/windows.jl +++ b/test/windows.jl @@ -90,7 +90,7 @@ end tukey_jl = tukey(128, 0.4) tukey_ml = readdlm(joinpath(dirname(@__FILE__), "data", "tukey128,0.4.txt"), '\t') @test tukey_jl ≈ tukey_ml - + @test tukey(128, 0) == rect(128) lanczos_jl = lanczos(128) @@ -150,53 +150,38 @@ end @testset "tensor product windows" begin # test all combinations of arguments. Each arg and kwarg can be not present, # a single value, or a 2-tuple - for winf in [zeroarg_wins; onearg_wins] - for arg in (nothing, 0.4, (0.4, 0.5)) - # skip invalid combinations - winf in zeroarg_wins && arg !== nothing && continue - winf in onearg_wins && arg === nothing && continue - for padding in (nothing, 4, (4,5)) - for zerophase in (nothing, true, (true,false)) - w1_expr = :($winf(15)) - w2_expr = :($winf(20)) - w3_expr = :($winf((15,20))) - - if arg isa Real - push!(w1_expr.args, arg) - push!(w2_expr.args, arg) - push!(w3_expr.args, arg) - elseif arg isa Tuple - push!(w1_expr.args, arg[1]) - push!(w2_expr.args, arg[2]) - push!(w3_expr.args, arg) - end - - if padding isa Integer - push!(w1_expr.args, Expr(:kw, :padding, padding)) - push!(w2_expr.args, Expr(:kw, :padding, padding)) - push!(w3_expr.args, Expr(:kw, :padding, padding)) - elseif padding isa Tuple - push!(w1_expr.args, Expr(:kw, :padding, padding[1])) - push!(w2_expr.args, Expr(:kw, :padding, padding[2])) - push!(w3_expr.args, Expr(:kw, :padding, padding)) - end - - if zerophase isa Bool - push!(w1_expr.args, Expr(:kw, :zerophase, zerophase)) - push!(w2_expr.args, Expr(:kw, :zerophase, zerophase)) - push!(w3_expr.args, Expr(:kw, :zerophase, zerophase)) - elseif zerophase isa Tuple - push!(w1_expr.args, Expr(:kw, :zerophase, zerophase[1])) - push!(w2_expr.args, Expr(:kw, :zerophase, zerophase[2])) - push!(w3_expr.args, Expr(:kw, :zerophase, zerophase)) - end - - w1 = eval(w1_expr) - w2 = eval(w2_expr) - w3 = eval(w3_expr) - @test w3 ≈ w1 * w2' - end - end + function push_args!((w1e, w2e, w3e), arg, ::Type{T}, f) where {T} + if arg isa T + push!(w1e.args, f(arg)) + push!(w2e.args, f(arg)) + push!(w3e.args, f(arg)) + elseif arg isa Tuple + push!(w1e.args, f(arg[1])) + push!(w2e.args, f(arg[2])) + push!(w3e.args, f(arg)) + end + end + + for winf in [zeroarg_wins; onearg_wins], + arg in (nothing, 0.4, (0.4, 0.5)) + # skip invalid combinations + winf in zeroarg_wins && arg !== nothing && continue + winf in onearg_wins && arg === nothing && continue + for padding in (nothing, 4, (4,5)), + zerophase in (nothing, true, (true,false)) + w1_expr = :($winf(15)) + w2_expr = :($winf(20)) + w3_expr = :($winf((15,20))) + w_all = (w1_expr, w2_expr, w3_expr) + + push_args!(w_all, arg, Real, identity) + push_args!(w_all, padding, Integer, s -> Expr(:kw, :padding, s)) + push_args!(w_all, zerophase, Bool, s -> Expr(:kw, :zerophase, s)) + + w1 = eval(w1_expr) + w2 = eval(w2_expr) + w3 = eval(w3_expr) + @test w3 ≈ w1 * w2' end end end