Reduce macro usage (#609)

* reduce eval usage * prevent `UndefKeywordError` * shorten tests * tanpi more tanpi * prevent possible method ambiguities * Force specialization on function arguments Co-authored-by: Martin Holters <[email protected]> * remove unnecessary exception Co-authored-by: Martin Holters <[email protected]> * inline `bsfmin`, `ordfreq_est` * remove ripple --------- Co-authored-by: Martin Holters <[email protected]>
JuliaDSP · Dec 13, 2024 · d8d423d · d8d423d
1 parent 8607be8
commit d8d423d
Show file tree

Hide file tree

Showing 5 changed files with 139 additions and 155 deletions.
diff --git 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
@@ -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
@@ -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
@@ -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