From 4f6390bd219bf9e1188641da8ccaa6f43481f1f4 Mon Sep 17 00:00:00 2001 From: Jordan Smith <46058723+jsmithnh09@users.noreply.github.com> Date: Fri, 1 Jul 2022 05:43:07 -0400 Subject: [PATCH] Butterworth, Elliptic, and Chebyshev I/II Filter Order Estimation (#453) * Adding filter order routines. * Adding test cases for `cheb1ord` and `ellipord` * Adding `cheb2ord` test cases. * Adding `cheb2ord` docstrings * Adding order estimation routines to the docs. * Adding argument error and `brent` test cases. --- docs/src/filters.md | 21 +- src/Filters/Filters.jl | 8 + src/Filters/filt_order.jl | 460 ++++++++++++++++++++++++++++++++++++++ test/filt_order.jl | 215 ++++++++++++++++++ test/runtests.jl | 2 +- 5 files changed, 704 insertions(+), 2 deletions(-) create mode 100644 src/Filters/filt_order.jl create mode 100644 test/filt_order.jl diff --git a/docs/src/filters.md b/docs/src/filters.md index bbe847bd9..40d97de6b 100644 --- a/docs/src/filters.md +++ b/docs/src/filters.md @@ -84,7 +84,9 @@ or [`Bandstop`](@ref) and includes the edges of the bands. The design method is [`Butterworth`](@ref), [`Chebyshev1`](@ref), [`Chebyshev2`](@ref), [`Elliptic`](@ref), or [`FIRWindow`](@ref), and includes any necessary parameters for the method that affect the shape of the response, -such as filter order, ripple, and attenuation. +such as filter order, ripple, and attenuation. [Filter order estimation methods](@ref order-est-methods) +are available in [`buttord`](@ref), [`cheb1ord`](@ref), [`cheb2ord`](@ref), +and [`ellipord`](@ref) if the corner frequencies for different IIR filter types are known. ```@docs analogfilter @@ -119,6 +121,15 @@ Chebyshev2 Elliptic ``` +### [IIR filter order estimation methods](@id order-est-methods) + +```@docs +buttord +cheb1ord +cheb2ord +ellipord +``` + #### FIR filter design methods ```@docs @@ -180,3 +191,11 @@ responsetype = Lowpass(5; fs=50) designmethod = FIRWindow(hanning(64)) filt(digitalfilter(responsetype, designmethod), x) ``` + +Estimate a Lowpass Elliptic filter order with a normalized +passband cutoff frequency of 0.2, a stopband cutoff frequency of 0.4, +3 dB of passband ripple, and 40 dB attenuation in the stopband: + +```julia +(N, ωn) = ellipord(0.2, 0.4, 3, 40) +``` \ No newline at end of file diff --git a/src/Filters/Filters.jl b/src/Filters/Filters.jl index 541923b09..ed5826db1 100644 --- a/src/Filters/Filters.jl +++ b/src/Filters/Filters.jl @@ -6,6 +6,7 @@ using Polynomials: LaurentPolynomial, Polynomial, coeffs, derivative, fromroots, import Base: * using LinearAlgebra: I, mul!, rmul! using Statistics: middle +using SpecialFunctions: ellipk import ..DSP: filt, filt!, optimalfftfiltlength, os_fft_complexity, SMALL_FILT_CUTOFF import Compat using FFTW @@ -44,6 +45,13 @@ export FilterType, FIRWindow, resample_filter +include("filt_order.jl") +export buttord, + ellipord, + cheb1ord, + cheb2ord + + include("response.jl") export freqresp, phaseresp, diff --git a/src/Filters/filt_order.jl b/src/Filters/filt_order.jl new file mode 100644 index 000000000..1de636415 --- /dev/null +++ b/src/Filters/filt_order.jl @@ -0,0 +1,460 @@ +#==================================================================== + +Filter order estimation routines. +Jordan R. Smith, 2021. + +Filter transformations translated from Scipy to Julia. + +SCIPY license: +Copyright (c) 2001, 2002 Enthought, Inc. +All rights reserved. + +Copyright (c) 2003-2017 SciPy Developers. +All rights reserved. + +Redistribution and use in source and binary forms, with or without +modification, are permitted provided that the following conditions are met: + + a. Redistributions of source code must retain the above copyright notice, + this list of conditions and the following disclaimer. + b. Redistributions in binary form must reproduce the above copyright + notice, this list of conditions and the following disclaimer in the + documentation and/or other materials provided with the distribution. + c. Neither the name of Enthought nor the names of the SciPy Developers + may be used to endorse or promote products derived from this software + without specific prior written permission. + + +THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" +AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE +ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDERS OR CONTRIBUTORS +BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, +OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF +THE POSSIBILITY OF SUCH DAMAGE. + +--- end of scipy license + +Design equations based on [1], Chapter 10. Brent's method with Parabolic interpolation +uses code from [3], modifications from Optim.jl's Brent method in [4]. + +[1] Proakis, J. G., & Manolakis, D. G. (1996). Digital Signal Processing, Fourth Edition. +Prentice Hall, New Jersey. + +[2] Virtanen, P., Gommers, R., Oliphant, T. E., Haberland, M., Reddy, T., +Cournapeau, D., ... & Van Mulbregt, P. (2020). SciPy 1.0: fundamental algorithms for scientific +computing in Python. Nature methods, 17(3), 261-272. + +[3] W. H. Press, et al. Numerical recipes. Third Edition. Cambridge university press, 2007. Chapter 10.3 + +[4] P. K. Mogensen and A. N. Riseth, Optim: A mathematical optimization package for Julia. Journal of Open +Source Software, 2018. https://github.com/JuliaNLSolvers/Optim.jl/blob/master/src/univariate/solvers/brent.jl + +====================================================================# + + + +toprototype(Wp::Real, Ws::Real, ftype::Type{Lowpass}) = Ws / Wp +toprototype(Wp::Real, Ws::Real, ftype::Type{Highpass}) = Wp / Ws +function toprototype(Wp::Tuple{Real,Real}, Ws::Tuple{Real,Real}, ftype::Type{Bandpass}) + # bandpass filter must have two corner frequencies we're computing with + Wa = (Ws.^2 .- Wp[1]*Wp[2]) ./ (Ws .* (Wp[1]-Wp[2])) + minimum(abs.(Wa)) +end +toprototype(Wp::Tuple{Real, Real}, Ws::Tuple{Real, Real}, Rp::Real, Rs::Real, ftype::Type{Bandstop}) = butterworth_bsfmin(Wp, Ws, Rp, Rs) +fromprototype(Wp::Real, Wscale::Real, ftype::Type{Lowpass}) = Wp * Wscale +fromprototype(Wp::Real, Wscale::Real, ftype::Type{Highpass}) = Wp / Wscale + +function fromprototype(Wp::Tuple{Real,Real}, Wscale::Real, ftype::Type{Bandstop}) + Wa = zeros(2) + diff = Wp[2]-Wp[1] + prod = Wp[2]*Wp[1] + Wa[1] = (diff + sqrt(diff^2 + 4*(Wscale^2)*prod)) / (2*Wscale) + Wa[2] = (diff - sqrt(diff^2 + 4*(Wscale^2)*prod)) / (2*Wscale) + sort!(abs.(Wa)) +end + +function fromprototype(Wp::Tuple{Real,Real}, Wscale::Real, ftype::Type{Bandpass}) + Wsc = [-Wscale, Wscale] + Wa = -Wsc .* (Wp[2]-Wp[1])./2 .+ sqrt.( Wsc.^2/4*(Wp[2]-Wp[1]).^2 .+ Wp[1]*Wp[2]) + sort!(abs.(Wa)) +end + +butterworth_order_estimate(Rp::Real, Rs::Real, warp::Real) = (log(db2pow(Rs) - 1) - log(db2pow(Rp) - 1)) / (2*log(warp)) +butterworth_natfreq_estimate(warp::Real, Rs::Real, order::Integer) = warp / (db2pow(Rs) - 1)^(1/(2*order)) + +function elliptic_order_estimate(Rp::Real, Rs::Real, Wa::Real) + # Elliptic integer order estmate. Requires Complete Elliptic integral of first kind. + ϵ = √(db2pow(Rp) - 1) # define selectivity/discrimination parameters. + k₁ = ϵ / √(db2pow(Rs) - 1) + k = Wa^-1 + (k^2 < 1) || throw(DomainError(k^2, "Selectivity parameter specifies too narrow of a transition width.")) # check if in-bounds (k₁ << k <~ 1) + (1-k₁^2 < 1) || throw(DomainError(1-k₁^2, "Discrimination parameter specifies too deep of a stopband.")) + K = tuple(ellipk(k^2), ellipk(1-k^2)) # define the complementary moduli + K₁ = tuple(ellipk(k₁^2), ellipk(1-k₁^2)) + + # other order approach would be using (k₁'/k₁) / (k′/k). + (K[1]*K₁[2])/(K[2]*K₁[1]) +end + +function chebyshev_order_estimate(Rp::Real, Rs::Real, Wa::Real) + # Chebyshev1/2 order estimate. Requires inverse hyperbolic cosine. + ϵs, ϵp = db2pow(Rs) - 1, db2pow(Rp) - 1 + acosh(√(ϵs / ϵp)) / acosh(Wa) # Eq. (10.3.63) in [1] +end + +function brent(f, x1::T, x2::T) where T <: AbstractFloat + a, b = x1, x2 + fa, fb = f(a), f(b) + ϵ, rtol = eps(T), √(eps(T)) + g = 1/2 * (3 - √(5)) # golden ratio + k = zero(T) + k1 = zero(T) + + # first interval calculation. + m = a + g*(b - a) + m1, m2 = m, m + fm = f(m) + fm1, fm2 = fm, fm + iter = 0 + + while(iter < 1_000) + p, q = zero(T), zero(T) + xt = rtol * abs(m) + ϵ + c = (b + a) / 2 + if (abs(m - c) ≤ 2*xt - (b - a)/2) + break + end + iter += 1 + if (abs(k1) > xt) # trial parabolic fit + r = (m - m1) * (fm - fm2) + q = (m - m2) * (fm - fm1) + p = (m - m2) * q - (m - m1) * r + q = 2*(q - r) + if (q > 0) + p = -p + else + q = -q + end + end + if abs(p) < abs(q*k1/2) && (p < q*(b - m)) && (p < q*(m - a)) # parabolic step + k1 = k + k = p/q + if ((m + k) < 2*xt) || ((b - (m + k)) < 2*xt) + k = (m < c) ? xt : -xt + end + else + k1 = (m < c) ? b - m : a - m + k = g * k1 + end + + if (abs(k) > xt) + xn = m + k + else + xn = m + ((k > 0) ? xt : -xt) + end + # re-assign variables for next iteration, ("Housekeeping" in [3]) + fn = f(xn) + if (fn < fm) + if (xn < m) + b = m + else + a = m + end + m2 = m1 + fm2 = fm1 + m1 = m + fm1 = fm + m = xn + fm = fn + else + if (xn < m) + a = xn + else + b = xn + end + if (fn ≤ fm1) || (m1 == m) + m2 = m1 + fm2 = fm1 + m1 = xn + fm1 = fn + elseif (fn ≤ fm2 || m2 == m || m2 == m1) + m2 = xn + fm2 = fn + end + end + end + m +end + +# +# Bandstop cost functions and passband minimization +# +for filt in (:butterworth, :elliptic, :chebyshev) + @eval begin + function $(Symbol(string(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 ? tuple(Wx, Wp[2]) : tuple(Wp[1], Wx) + + # get the new warp frequency. + warp = minimum(abs.((Ws .* (Wpc[1]-Wpc[2])) ./ (Ws.^2 .- (Wpc[1]*Wpc[2])))) + + # use the new frequency to determine the filter order. + $(Symbol(string(filt, "_order_estimate")))(Rp, Rs, warp) + end + + function $(Symbol(string(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(string(filt, "_bsfcost")))(w, true, Wp, Ws, Rp, Rs) + p1 = brent(C₁, Wp[1], Ws[1]-Δ) + + C₂(w) = $(Symbol(string(filt, "_bsfcost")))(w, false, tuple(p1, Wp[2]), Ws, Rp, Rs) + p2 = brent(C₂, Ws[2]+Δ, Wp[2]) + Wadj = tuple(p1, p2) + + Wa = (Ws .* (Wadj[1]-Wadj[2])) ./ (Ws.^2 .- (Wadj[1]*Wadj[2])) + minimum(abs.(Wa)), Wadj + end + end +end + +""" + (N, ωn) = buttord(Wp::Tuple{Real,Real}, Ws::Tuple{Real,Real}, Rp::Real, Rs::Real; domain=:z) + +Butterworth order estimate for bandpass and bandstop filter types. `Wp` and `Ws` are 2-element pass +and stopband frequency edges, with no more than `Rp` dB passband ripple and at least `Rs` dB stopband +attenuation. Based 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. If a domain of `:s` is specified, the passband and stopband frequencies +are interpretted as radians/second, giving an order and natural frequencies for an analog filter. The default +domain is `:z`, interpretting the input frequencies as normalized from 0 to 1, where 1 corresponds to π radians/sample. +""" +function buttord(Wp::Tuple{Real,Real}, Ws::Tuple{Real,Real}, Rp::Real, Rs::Real; domain::Symbol=:z) + + # make sure the band edges are in increasing order. + Wps = (Wp[1] > Wp[2]) ? tuple(Wp[2], Wp[1]) : tuple(Wp[1], Wp[2]) + Wss = (Ws[1] > Ws[2]) ? tuple(Ws[2], Ws[1]) : tuple(Ws[1], Ws[2]) + + # infer filter type based on ordering of edges. + (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 both components, (if Z-domain specified) + if (domain == :z) + Ωp = tan.(π/2 .* Wps) + Ωs = tan.(π/2 .* Wss) + else + Ωp = Wps + Ωs = Wss + end + + if (ftype == Bandstop) + # optimizer modifies passband frequencies. + (wa, wpadj) = toprototype(Ωp, Ωs, Rp, Rs, ftype) + else + # no optimization in BPF case, use original warped passband edges. + wa = toprototype(Ωp, Ωs, ftype) + wpadj = Ωp + end + + # get the integer order estimate. + N = ceil(Int, butterworth_order_estimate(Rp, Rs, wa)) + + wscale = butterworth_natfreq_estimate(wa, Rs, N) + if (domain == :z) + ωn = (2/π).*atan.(fromprototype(wpadj, wscale, ftype)) + else + ωn = fromprototype(wpadj, wscale, ftype) + end + N, ωn +end + +""" + (N, ωn) = buttord(Wp::Real, Ws::Real, Rp::Real, Rs::Real; domain=:z) + +LPF/HPF Butterworth filter order and -3 dB frequency approximation. `Wp` and `Ws` are +the passband and stopband frequencies, whereas Rp and Rs are the passband and stopband +ripple attenuations in dB. If the passband is greater than stopband, the filter type +is inferred to be for estimating the order of a highpass filter. `N` specifies the lowest +possible integer filter order, whereas `ωn` is the cutoff or "-3 dB" frequency. If a domain +of `:s` is specified, the passband and stopband edges are interpretted as radians/second, +giving an order and natural frequency result for an analog filter. The default domain +is `:z`, interpretting the input frequencies as normalized from 0 to 1, where 1 corresponds +to π radians/sample. +""" +function buttord(Wp::Real, Ws::Real, Rp::Real, Rs::Real; domain::Symbol=:z) + # infer which filter type based on the frequency ordering. + ftype = (Wp < Ws) ? Lowpass : Highpass + + if (domain == :z) + # need to pre-warp since we want to use formulae for analog case. + Ωp = tan(π/2 * Wp) + Ωs = tan(π/2 * Ws) + else + # already analog + Ωp = Wp + Ωs = Ws + end + wa = toprototype(Ωp, Ωs, ftype) + + # rounding up fractional order. Using differences of logs instead of division. + N = ceil(Int, butterworth_order_estimate(Rp, Rs, wa)) + + # specifications for the stopband ripple are met precisely. + wscale = butterworth_natfreq_estimate(wa, Rs, N) + + # convert back to the original analog filter + if (domain == :z) + # bilinear xform + ωn = (2/π)*atan(fromprototype(Ωp, wscale, ftype)) + else + # s-domain, no atan call. + ωn = fromprototype(Ωp, wscale, ftype) + end + N, ωn +end + +# +# Elliptic/Chebyshev1 Estimation +# +for (fcn, est, filt) in ((:ellipord, :elliptic, "Elliptic (Cauer)"), + (:cheb1ord, :chebyshev, "Chebyshev Type I")) + @eval begin + """ + (N, ωn) = $($fcn)(Wp::Real, Ws::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.) + 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.) + If a domain of `:s` is specified, the passband and stopband edges are interpretted + as radians/second, giving an order and natural frequency result for an analog filter. + The default domain is `:z`, interpretting 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(string(est, "_order_estimate")))(Rp, Rs, wa)) + if (domain == :z) + ωn = (2/π)*atan(Ωp) + else + ωn = Ωp + end + N, ωn + end + + """ + (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. Based 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. If a domain of + `:s` is specified, the passband and stopband frequencies are interpretted as radians/second, giving an + order and natural frequencies for an analog filter. The default domain is `:z`, interpretting 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 = (Wp[1] > Wp[2]) ? tuple(Wp[2], Wp[1]) : tuple(Wp[1], Wp[2]) + Wss = (Ws[1] > Ws[2]) ? tuple(Ws[2], Ws[1]) : tuple(Ws[1], Ws[2]) + (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 = (Ωs.^2 .- (Ωp[1]*Ωp[2])) ./ (Ωs .* (Ωp[1]-Ωp[2])) + Ωpadj = Ωp + else + (Wa, Ωpadj) = $(Symbol(string(est, "_bsfmin")))(Ωp, Ωs, Rp, Rs) # check scipy. + end + N = ceil(Int, $(Symbol(string(est, "_order_estimate")))(Rp, Rs, minimum(abs.(Wa)))) + ωn = (domain == :z) ? Wps : Ωpadj + N, ωn + end + end +end + + +""" + (N, ωn) = cheb2ord(Wp::Real, Ws::Real, Rp::Real, Rs::Real; domain::Symbol=:z) + +Integer and natural frequency order estimation for Chebyshev Type II (inverse) Filters. `Wp` and `Ws` +are the 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. Based on the ordering of +the passband and stopband edges, the Lowpass or Highpass filter type is inferred. `N` is the integer +indicating the lowest filter order, with `ωn` specifying the "-3 dB" cutoff frequency. If the domain +is specified as `:s`, the passband and stopband frequencies are interpretted as radians/second, giving +the order and natural frequencies for an analog filter. The default domain is `:z`, interpretting the +input frequencies as normalized from 0 to 1, where 1 corresponds to π radians/sample. +""" +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) + wa = toprototype(Ωp, Ωs, ftype) + N = ceil(Int, chebyshev_order_estimate(Rp, Rs, wa)) + + # new frequency for stopband spec. + wnew = 1 / cosh(1 / N * acosh(√(db2pow(Rs) - 1) / √(db2pow(Rp) - 1))) + wa = (ftype == Lowpass) ? Ωp / wnew : Ωp * wnew + ωn = (domain == :z) ? (2/π)*atan(wa) : wa + N, ωn +end + + +""" + (N, ωn) = cheb2ord(Wp::Tuple{Real, Real}, Ws::Tuple{Real, Real}, Rp::Real, Rs::Real; domain::Symbol=:z) + +Integer and natural frequency order estimation for Chebyshev Type II (inverse) 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. Based 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. If a domain of +`:s` is specified, the passband and stopband frequencies are interpretted as radians/second, giving an +order and natural frequencies for an analog filter. The default domain is `:z`, interpretting the input +frequencies as normalized from 0 to 1, where 1 corresponds to π radians/sample. +""" +function cheb2ord(Wp::Tuple{Real, Real}, Ws::Tuple{Real, Real}, Rp::Real, Rs::Real; domain::Symbol=:z) + Wps = (Wp[1] > Wp[2]) ? tuple(Wp[2], Wp[1]) : tuple(Wp[1], Wp[2]) + Wss = (Ws[1] > Ws[2]) ? tuple(Ws[2], Ws[1]) : tuple(Ws[1], Ws[2]) + (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) + if (ftype == Bandpass) + Wa = (Ωs.^2 .- (Ωp[1]*Ωp[2])) ./ (Ωs .* (Ωp[1]-Ωp[2])) + else + (Wa, Ωpadj) = chebyshev_bsfmin(Ωp, Ωs, Rp, Rs) + end + N = ceil(Int, chebyshev_order_estimate(Rp, Rs, minimum(abs.(Wa)))) + + # new frequency for stopband spec. + wnew = 1 / cosh(1 / N * acosh(√(db2pow(Rs) - 1) / √(db2pow(Rp) - 1))) + + # lowpass prototype to analog filter re-map. + Wna = zeros(2) + if (ftype == Bandpass) + Wna[1] = (Ωp[1]-Ωp[2])/(2*wnew) + √( (Ωp[2]-Ωp[1])^2 / (4*wnew^2) + Ωp[1]*Ωp[2]) + Wna[2] = (Ωp[1]*Ωp[2]) / Wna[1] + else + Wna[1] = ((Ωpadj[1]-Ωpadj[2])*wnew)/2 + √((Ωpadj[1]-Ωpadj[2])^2 * wnew^2 / 4 + (Ωpadj[1]*Ωpadj[2])) + Wna[2] = (Ωpadj[1]*Ωpadj[2]) / Wna[1] + end + ωn = (domain == :z) ? tuple((2/π)*atan(Wna[1]), (2/π)*atan(Wna[2])) : tuple(Wna[1], Wna[2]) + N, ωn +end \ No newline at end of file diff --git a/test/filt_order.jl b/test/filt_order.jl new file mode 100644 index 000000000..3a760691a --- /dev/null +++ b/test/filt_order.jl @@ -0,0 +1,215 @@ +using DSP, Test + + +Δ = 1e-3 +@testset "buttord" begin + + # + # https://www.mathworks.com/help/signal/ref/buttord.html#d123e9937 + # Example 1: Lowpass Butterworth Example + # + # Wp = 40/500; + # Ws = 150/500; + # [n, Wn] = buttord(Wp, Ws, 3, 60) + # + + # z-domain test + (n, Wn) = buttord(40/500, 150/500, 3, 60, domain=:z) + @test n == 5 + @test Wn ≈ 0.081038494957764 + + # s-domain test + (ns, Wns) = buttord(40/500, 150/500, 3, 60, domain=:s) + @test ns == 6 + @test Wns ≈ 0.0948683377107 + + # + # Highpass filter example + # Wp = 1200/2000; Ws=600/2000; + # Rs = 3; Rp = 60; + + # z-domain (default) + (nhpf, Wnhpf) = buttord(1200/2000, 600/2000, 3, 60, domain=:z) + @test nhpf == 7 + @test Wnhpf ≈ 0.597905417809 + + + # s-domain test + (nhpfs, Wnhpfs) = buttord(1200/2000, 600/2000, 3, 60, domain=:s) + @test nhpfs == 10 + @test Wnhpfs ≈ 0.598578664562 + + # + # https://www.mathworks.com/help/signal/ref/buttord.html#d123e9937 + # Example 2: Bandpass Butterworth Example + # + # Wp = [100 200]/500; + # Ws = [50 250]/500; + # Rp = 3; + # Rs = 40; + # [n, Wn] = buttord(Wp, Ws, Rp, Rs) + # + + (nbp, Wnbp) = buttord(tuple(100/500, 200/500), tuple(50/500, 250/500), 3, 40, domain=:z) + @test nbp == 8 + @test Wnbp[1] ≈ 0.195101359239 + @test Wnbp[2] ≈ 0.408043633382 + + # s-domain + (nbps, Wnbps) = buttord(tuple(100/500, 200/500), tuple(50/500, 250/500), 3, 40, domain=:s) + @test nbps == 9 + @test Wnbps[1] ≈ 0.198730150231 + @test Wnbps[2] ≈ 0.402555927759 + + # + # Bandstop Example, (44.1 kHz Nyquist) + # Wp = [3200 7800]/22050 + # Ws = [4800 5600]/22050 + # Rp = 2 + # Rs = 60 + + # this test may be more sensitive...MATLAB's implementation of bounded minimization + # will yield different results in comparison to Optim.jl. + (nbs, Wnbs) = buttord(tuple(3200/22050, 7800/22050), tuple(4800/22050, 5600/22050), 2, 60, domain=:z) + @test nbs == 5 + @test ≈(Wnbs[1], 0.172660908966, rtol=Δ) + @test ≈(Wnbs[2], 0.314956388749, rtol=Δ) + + # s-domain + (nbss, Wnbss) = buttord(tuple(3200/22050, 7800/22050), tuple(4800/22050, 5600/22050), 2, 60, domain=:s) + @test ≈(Wnbss[1], 0.173677826752, rtol=Δ) + @test ≈(Wnbss[2], 0.318267164272, rtol=Δ) + +end + +# +# Elliptic/Chebyshev I/II Filter cases tested against SciPy signal results. +# +@testset "ellipord" begin + + Rp, Rs = 3, 40 + Wp = tuple(0.2, 0.7) + Ws = tuple(0.1, 0.8) + + # Lowpass + (nl, Wnl) = ellipord(0.1, 0.2, Rp, Rs, domain=:z) + @test nl == 3 + @test Wnl == 0.1 + + # Highpass + (nh, Wnh) = ellipord(0.3, 0.1, Rp, Rs, domain=:z) + @test nh == 3 + @test Wnh == 0.3 + + # Bandpass (z-domain) + (nbp, Wnbp) = ellipord(Wp, Ws, Rp, Rs, domain=:z) + @test nbp == 4 + @test Wnbp == Wp + + # Bandpass (s-domain) + (nbp, Wnbp) = ellipord(Wp, Ws, Rp, Rs, domain=:s) + @test nbp == 5 + @test Wnbp == Wp + + # Bandstop (z-domain) + (nbs, Wnbs) = ellipord(Ws, Wp, Rp, Rs, domain=:z) + @test nbs == 4 + @test Wnbs == Ws + + # Bandstop (s-domain) + # n, Wn = scipy.signal.ellipord([0.1, 0.8], [0.2, 0.7], 3, 40, True) + (nbs, Wnbs) = ellipord(Ws, Wp, Rp, Rs, domain=:s) + @test nbs == 5 + @test ≈(Wnbs[1], 0.17500000332788998, rtol=Δ) + @test ≈(Wnbs[2], 0.799993389303865, rtol=Δ) + +end + +@testset "cheb1ord" begin + Rp, Rs = 2, 70 + Wp = tuple(0.2, 0.5) + Ws = tuple(0.1, 0.6) + + # Lowpass + (nl, Wnl) = cheb1ord(0.1, 0.21, Rp, Rs, domain=:z) + @test nl == 7 + @test Wnl == 0.1 + + # Highpass + (nh, Wnh) = cheb1ord(0.1, 0.04, Rp, Rs, domain=:z) + @test nh == 6 + @test Wnh == 0.1 + + # Bandpass (z-domain) + (nbp, Wnbp) = cheb1ord(Wp, Ws, Rp, Rs, domain=:z) + @test nbp == 9 + @test Wnbp == (0.2, 0.5) + + + # Bandpass (s-domain) + (nbp, Wnbp) = cheb1ord(Wp, Ws, Rp, Rs, domain=:s) + @test nbp == 10 + @test Wnbp == (0.2, 0.5) + + # Bandstop (z-domain) + (nbs, Wnbs) = cheb1ord(Ws, Wp, Rp, Rs, domain=:z) + @test nbs == 9 + @test Wnbs == Ws + + # Bandstop (s-domain) + (nbs, Wnbs) = cheb1ord(Ws, Wp, Rp, Rs, domain=:s) + @test nbs == 10 + @test ≈(Wnbs[1], 0.166666612185443, rtol=Δ) + @test ≈(Wnbs[2], 0.5999933893038649, rtol=Δ) +end + +@testset "cheb2ord" begin + Rp, Rs = 1.2, 80 + Wp = tuple(0.22, 0.51) + Ws = tuple(0.14, 0.63) + + # Lowpass + (nl, Wnl) = cheb2ord(0.1, 0.21, Rp, Rs, domain=:z) + @test nl == 8 + @test Wnl == 0.19411478246577737 + (nl, Wnl) = cheb2ord(0.1, 0.21, Rp, Rs, domain=:s) + @test nl == 8 + @test Wnl == 0.1987124302811051 + + # Highpass + (nh, Wnh) = cheb2ord(0.21, 0.1, Rp, Rs, domain=:z) + @test nh == 8 + @test Wnh == 0.10862150541420543 + (nh, Wnh) = cheb2ord(0.21, 0.1, Rp, Rs, domain=:s) + @test nh == 8 + @test Wnh == 0.10568035411923006 + + # Bandpass + (nbp, Wnbp) = cheb2ord(Wp, Ws, Rp, Rs, domain=:z) + @test nbp == 9 + @test ≈(Wnbp[1], 0.1608459041132262) + @test ≈(Wnbp[2], 0.6133747025904719) + (nbp, Wnbp) = cheb2ord(Wp, Ws, Rp, Rs, domain=:s) + @test nbp == 11 + @test ≈(Wnbp[1], 0.18262279523940905) + @test ≈(Wnbp[2], 0.6143811338169016) + + # Bandstop + (nbs, Wnbs) = cheb2ord(Ws, Wp, Rp, Rs, domain=:z) + @test nbs == 9 + @test ≈(Wnbs[1], 0.21211425852327126, rtol=Δ) + @test ≈(Wnbs[2], 0.5225427194473862, rtol=Δ) + (nbs, Wnbs) = cheb2ord(Ws, Wp, Rp, Rs, domain=:s) + @test nbs == 11 + @test ≈(Wnbs[1], 0.2159740591083134, rtol=Δ) + @test ≈(Wnbs[2], 0.5195028184932494, rtol=Δ) + +end + +# using some simple examples for testing Brent's method. +@testset "brent" begin + f1(x) = (x+3) * ((x-1)^2) # x³ + x² - 5x + 3 + @test ≈(f1(DSP.Filters.brent(f1, -4.0, 4.0)), 0.0, atol=1e-8) + @test ≈(sin(DSP.Filters.brent(sin, 0.0, 2pi)), -1.0, atol=1e-8) + @test ≈(cos(DSP.Filters.brent(cos, 0.0, 2pi)), -1.0, atol=1e-8) +end diff --git a/test/runtests.jl b/test/runtests.jl index bf7e4a398..8b245c233 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -5,7 +5,7 @@ using Random: seed! testfiles = ["dsp.jl", "util.jl", "windows.jl", "filter_conversion.jl", "diric.jl", - "filter_design.jl", "filter_response.jl", "filt.jl", "filt_stream.jl", + "filter_design.jl", "filter_response.jl", "filt.jl", "filt_order.jl", "filt_stream.jl", "periodograms.jl", "multitaper.jl", "resample.jl", "lpc.jl", "estimation.jl", "unwrap.jl", "remez_fir.jl" ]