From 08be44b309b7f5d88b8c3bb71c502b935fa536e4 Mon Sep 17 00:00:00 2001 From: wheeheee <104880306+wheeheee@users.noreply.github.com> Date: Sun, 6 Oct 2024 13:38:19 +0800 Subject: [PATCH 1/9] reduce eval usage --- src/Filters/filt_order.jl | 146 ++++++++++++++++++-------------------- src/windows.jl | 41 ++++++----- 2 files changed, 91 insertions(+), 96 deletions(-) diff --git a/src/Filters/filt_order.jl b/src/Filters/filt_order.jl index 5f9e713a7..58a0da113 100644 --- a/src/Filters/filt_order.jl +++ b/src/Filters/filt_order.jl @@ -196,35 +196,35 @@ 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, 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) - # 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::T, Wp::Tuple{Real,Real}, Ws::Tuple{Real,Real}, Rp::Real, Rs::Real) where T + # 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) @@ -327,76 +327,72 @@ end # # Elliptic/Chebyshev1 Estimation # -for (fcn, est, filt) in ((:ellipord, :elliptic, "Elliptic (Cauer)"), - (:cheb1ord, :chebyshev, "Chebyshev Type I")) +function ordfreq_est(order_estimate, 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, order_estimate(Rp, Rs, wa)) + if (domain == :z) + ωn = (2 / π) * atan(Ωp) + else + ωn = Ωp + end + N, ωn +end + +function ordfreq_est(order_estimate, 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) = 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) = ordfreq_est(elliptic_order_estimate, args...; domain) +cheb1ord(args...; domain) = ordfreq_est(chebyshev_order_estimate, args...; domain) + +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 ripple loss in the passband and the minimum ripple attenuation in the stopband, (in dB.)\n 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) diff --git a/src/windows.jl b/src/windows.jl index 6311067b7..fcf5873af 100644 --- a/src/windows.jl +++ b/src/windows.jl @@ -632,33 +632,32 @@ 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, dims::Tuple, arg::Union{RealOr2,Nothing}=nothing; + padding::IntegerOr2=0, zerophase::BoolOr2=false) + length(dims) == 2 || throw(ArgumentError("`dims` must be length 2")) + 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; 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, arg::RealOr2; padding::IntegerOr2=0, zerophase::BoolOr2=false) + return matrix_window($func, dims, arg; padding, zerophase) end end From a94a2243d54699e6f0d8ad342b0250870a28be2e Mon Sep 17 00:00:00 2001 From: wheeheee <104880306+wheeheee@users.noreply.github.com> Date: Sun, 6 Oct 2024 21:33:24 +0800 Subject: [PATCH 2/9] prevent `UndefKeywordError` --- src/Filters/filt_order.jl | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/src/Filters/filt_order.jl b/src/Filters/filt_order.jl index 58a0da113..21a1bd07b 100644 --- a/src/Filters/filt_order.jl +++ b/src/Filters/filt_order.jl @@ -327,7 +327,7 @@ end # # Elliptic/Chebyshev1 Estimation # -function ordfreq_est(order_estimate, Wp::Real, Ws::Real, Rp::Real, Rs::Real; domain::Symbol=:z) +function ordfreq_est(order_estimate, domain::Symbol, Wp::Real, Ws::Real, Rp::Real, Rs::Real) ftype = (Wp < Ws) ? Lowpass : Highpass if (domain == :z) Ωp = tan(π / 2 * Wp) @@ -347,7 +347,8 @@ function ordfreq_est(order_estimate, Wp::Real, Ws::Real, Rp::Real, Rs::Real; dom N, ωn end -function ordfreq_est(order_estimate, Wp::Tuple{Real,Real}, Ws::Tuple{Real,Real}, Rp::Real, Rs::Real; domain::Symbol=:z) +function ordfreq_est(order_estimate, domain::Symbol, + Wp::Tuple{Real,Real}, Ws::Tuple{Real,Real}, Rp::Real, Rs::Real) 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.")) @@ -365,8 +366,8 @@ function ordfreq_est(order_estimate, Wp::Tuple{Real,Real}, Ws::Tuple{Real,Real}, N, ωn end -ellipord(args...; domain) = ordfreq_est(elliptic_order_estimate, args...; domain) -cheb1ord(args...; domain) = ordfreq_est(chebyshev_order_estimate, args...; domain) +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 From 53fdb96c863d66691dc4db4871d8a848af2535da Mon Sep 17 00:00:00 2001 From: wheeheee <104880306+wheeheee@users.noreply.github.com> Date: Mon, 7 Oct 2024 17:59:53 +0800 Subject: [PATCH 3/9] shorten tests --- test/windows.jl | 81 ++++++++++++++++++++----------------------------- 1 file changed, 33 insertions(+), 48 deletions(-) diff --git a/test/windows.jl b/test/windows.jl index 34db2ff66..2473e018c 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 From 98c4ba1861e95fa59a90b993decc51708c026162 Mon Sep 17 00:00:00 2001 From: wheeheee <104880306+wheeheee@users.noreply.github.com> Date: Mon, 18 Nov 2024 22:31:31 +0800 Subject: [PATCH 4/9] tanpi more tanpi --- src/Filters/design.jl | 4 ++-- src/Filters/filt_order.jl | 18 +++++++++--------- test/filt_order.jl | 6 +++--- 3 files changed, 14 insertions(+), 14 deletions(-) diff --git a/src/Filters/design.jl b/src/Filters/design.jl index df485b3fe..b733edabb 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 21a1bd07b..a01b551d7 100644 --- a/src/Filters/filt_order.jl +++ b/src/Filters/filt_order.jl @@ -250,8 +250,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 +298,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 @@ -330,8 +330,8 @@ end function ordfreq_est(order_estimate, domain::Symbol, Wp::Real, Ws::Real, Rp::Real, Rs::Real) ftype = (Wp < Ws) ? Lowpass : Highpass if (domain == :z) - Ωp = tan(π / 2 * Wp) - Ωs = tan(π / 2 * Ws) + Ωp = tanpi(Wp / 2) + Ωs = tanpi(Ws / 2) else Ωp = Wp Ωs = Ws @@ -354,7 +354,7 @@ function ordfreq_est(order_estimate, domain::Symbol, (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) + (Ω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 @@ -410,7 +410,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)) @@ -441,7 +441,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/test/filt_order.jl b/test/filt_order.jl index 3ee67039a..f4c8faad7 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 From eba2957adbb65e5a42149ac99ea3d1f71acf33e6 Mon Sep 17 00:00:00 2001 From: wheeheee <104880306+wheeheee@users.noreply.github.com> Date: Sat, 30 Nov 2024 11:19:16 +0800 Subject: [PATCH 5/9] prevent possible method ambiguities --- src/windows.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/windows.jl b/src/windows.jl index fcf5873af..da67e04f4 100644 --- a/src/windows.jl +++ b/src/windows.jl @@ -632,7 +632,7 @@ 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, dims::Tuple, arg::Union{RealOr2,Nothing}=nothing; +function matrix_window(func, dims::Tuple{Integer,Integer}, arg::Union{RealOr2,Nothing}=nothing; padding::IntegerOr2=0, zerophase::BoolOr2=false) length(dims) == 2 || throw(ArgumentError("`dims` must be length 2")) paddings = argdup(padding) @@ -650,13 +650,13 @@ end for func in (:rect, :hanning, :hamming, :cosine, :lanczos, :triang, :bartlett, :bartlett_hann, :blackman) - @eval function $func(dims; padding::IntegerOr2=0, zerophase::BoolOr2=false) + @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 function $func(dims, arg::RealOr2; padding::IntegerOr2=0, zerophase::BoolOr2=false) + @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 From 8f20f0acd4f6c54cdfce1a0b89fb5c8011b642e9 Mon Sep 17 00:00:00 2001 From: wheeheee <104880306+wheeheee@users.noreply.github.com> Date: Fri, 13 Dec 2024 17:10:35 +0800 Subject: [PATCH 6/9] Force specialization on function arguments Co-authored-by: Martin Holters --- src/Filters/filt_order.jl | 8 ++++---- src/windows.jl | 4 ++-- 2 files changed, 6 insertions(+), 6 deletions(-) diff --git a/src/Filters/filt_order.jl b/src/Filters/filt_order.jl index a01b551d7..19e8a41e2 100644 --- a/src/Filters/filt_order.jl +++ b/src/Filters/filt_order.jl @@ -196,7 +196,7 @@ end # # Bandstop cost functions and passband minimization # -function bsfcost(est_func, Wx::Real, uselowband::Bool, Wp::Tuple{Real,Real}, Ws::Tuple{Real,Real}, Rp::Real, Rs::Real) +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) @@ -207,7 +207,7 @@ function bsfcost(est_func, Wx::Real, uselowband::Bool, Wp::Tuple{Real,Real}, Ws: est_func(Rp, Rs, warp) end -function bsfmin(est_func::T, Wp::Tuple{Real,Real}, Ws::Tuple{Real,Real}, Rp::Real, Rs::Real) where T +function bsfmin(est_func::F, Wp::Tuple{Real,Real}, Ws::Tuple{Real,Real}, Rp::Real, Rs::Real) where {F} # 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) @@ -347,8 +347,8 @@ function ordfreq_est(order_estimate, domain::Symbol, Wp::Real, Ws::Real, Rp::Rea N, ωn end -function ordfreq_est(order_estimate, domain::Symbol, - Wp::Tuple{Real,Real}, Ws::Tuple{Real,Real}, Rp::Real, Rs::Real) +function ordfreq_est(order_estimate::F, domain::Symbol, + Wp::Tuple{Real,Real}, Ws::Tuple{Real,Real}, Rp::Real, Rs::Real) where {F} 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.")) diff --git a/src/windows.jl b/src/windows.jl index da67e04f4..3315e2d23 100644 --- a/src/windows.jl +++ b/src/windows.jl @@ -632,8 +632,8 @@ 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, dims::Tuple{Integer,Integer}, arg::Union{RealOr2,Nothing}=nothing; - padding::IntegerOr2=0, zerophase::BoolOr2=false) +function matrix_window(func::F, dims::Tuple{Integer,Integer}, arg::Union{RealOr2,Nothing}=nothing; + padding::IntegerOr2=0, zerophase::BoolOr2=false) where {F} length(dims) == 2 || throw(ArgumentError("`dims` must be length 2")) paddings = argdup(padding) zerophases = argdup(zerophase) From 26f4e4aedbefcde9b36c152e0c3cd455a01fad0f Mon Sep 17 00:00:00 2001 From: wheeheee <104880306+wheeheee@users.noreply.github.com> Date: Fri, 13 Dec 2024 17:13:53 +0800 Subject: [PATCH 7/9] remove unnecessary exception Co-authored-by: Martin Holters --- src/windows.jl | 1 - 1 file changed, 1 deletion(-) diff --git a/src/windows.jl b/src/windows.jl index 3315e2d23..4b3b72723 100644 --- a/src/windows.jl +++ b/src/windows.jl @@ -634,7 +634,6 @@ 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} - length(dims) == 2 || throw(ArgumentError("`dims` must be length 2")) paddings = argdup(padding) zerophases = argdup(zerophase) if isnothing(arg) From f61131859a4de81d0e73b1c98c9ca55167bad110 Mon Sep 17 00:00:00 2001 From: wheeheee <104880306+wheeheee@users.noreply.github.com> Date: Fri, 13 Dec 2024 17:37:51 +0800 Subject: [PATCH 8/9] inline `bsfmin`, `ordfreq_est` --- src/Filters/filt_order.jl | 3 +++ 1 file changed, 3 insertions(+) diff --git a/src/Filters/filt_order.jl b/src/Filters/filt_order.jl index 19e8a41e2..a08c571fd 100644 --- a/src/Filters/filt_order.jl +++ b/src/Filters/filt_order.jl @@ -208,6 +208,7 @@ function bsfcost(est_func::F, Wx::Real, uselowband::Bool, Wp::Tuple{Real,Real}, end 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) @@ -328,6 +329,7 @@ end # Elliptic/Chebyshev1 Estimation # 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) @@ -349,6 +351,7 @@ 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.")) From 3197b9825a52bd663014793c3e3fb8a7df2952c8 Mon Sep 17 00:00:00 2001 From: wheeheee <104880306+wheeheee@users.noreply.github.com> Date: Fri, 13 Dec 2024 17:48:36 +0800 Subject: [PATCH 9/9] remove ripple --- src/Filters/filt_order.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/Filters/filt_order.jl b/src/Filters/filt_order.jl index a08c571fd..8d22d6ef0 100644 --- a/src/Filters/filt_order.jl +++ b/src/Filters/filt_order.jl @@ -380,7 +380,8 @@ for (fcn, filt) in ((:ellipord, "Elliptic (Cauer)"), (:cheb1ord, "Chebyshev Type 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 ripple loss in the passband and the minimum ripple 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