Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Reduce macro usage #609

Merged
merged 9 commits into from
Dec 13, 2024
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
4 changes: 2 additions & 2 deletions src/Filters/design.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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
"""
Expand Down Expand Up @@ -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)
Expand Down
163 changes: 82 additions & 81 deletions src/Filters/filt_order.jl
Original file line number Diff line number Diff line change
Expand Up @@ -196,35 +196,36 @@
#
# 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...)

Check warning on line 226 in src/Filters/filt_order.jl

View check run for this annotation

Codecov / codecov/patch

src/Filters/filt_order.jl#L226

Added line #L226 was not covered by tests
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)

Expand All @@ -250,8 +251,8 @@

# 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
Expand Down Expand Up @@ -298,8 +299,8 @@

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
Expand Down Expand Up @@ -327,76 +328,76 @@
#
# 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

Check warning on line 339 in src/Filters/filt_order.jl

View check run for this annotation

Codecov / codecov/patch

src/Filters/filt_order.jl#L338-L339

Added lines #L338 - L339 were not covered by tests
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

Check warning on line 347 in src/Filters/filt_order.jl

View check run for this annotation

Codecov / codecov/patch

src/Filters/filt_order.jl#L347

Added line #L347 was not covered by tests
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)

Expand All @@ -413,7 +414,7 @@
"""
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))

Expand Down Expand Up @@ -444,7 +445,7 @@
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]
Expand Down
40 changes: 19 additions & 21 deletions src/windows.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand Down
6 changes: 3 additions & 3 deletions test/filt_order.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down
Loading
Loading