Expose resample_filter optional arguments in resample, FIRFilter (

#621) * expose resample_filter options in resample * document `FIRFilter`, `resample_filter` _resample_filter * const fields for `FIRKernel`s * misc, spacing remove inbounds, correct typo * fix `filtfilt` pad length * Apply suggestions from code review Co-authored-by: Martin Holters <[email protected]> * add tests * add more tests * reuse `FIRFilter`s for arrays * interpolation * test resample with resample_filter * Specify Int64, use `@test` broken argument Co-authored-by: Martin Holters <[email protected]> --------- Co-authored-by: Martin Holters <[email protected]>
JuliaDSP · Jan 30, 2025 · cd2912c · cd2912c
1 parent 06f8776
commit cd2912c
Show file tree

Hide file tree

Showing 7 changed files with 199 additions and 163 deletions.
diff --git a/docs/src/filters.md b/docs/src/filters.md
@@ -59,6 +59,7 @@ parameters.
 
 ```@docs
 FIRFilter
+resample_filter
 ```
 
 ## Filter application
@@ -196,7 +197,7 @@ coefa
 ## Examples
 
 Construct a 4th order elliptic lowpass filter with normalized cutoff
-frequency 0.2, 0.5 dB of passband ripple, and 30 dB attentuation in
+frequency 0.2, 0.5 dB of passband ripple, and 30 dB attenuation in
 the stopband and extract the coefficients of the numerator and
 denominator of the transfer function:
 

diff --git a/src/Filters/design.jl b/src/Filters/design.jl
@@ -124,7 +124,7 @@ function landen(k::Real)
     kn = Vector{typeof(k)}(undef, niter)
     # Eq. (50)
     for i = 1:niter
-        kn[i] = k = abs2(k/(1+sqrt(1-abs2(k))))
+        kn[i] = k = abs2(k / (1 + sqrt(1 - abs2(k))))
     end
     kn
 end
@@ -223,7 +223,7 @@ end
     Elliptic(n::Integer, rp::Real, rs::Real)
 
 `n` pole elliptic (Cauer) filter with `rp` dB ripple in the
-passband and `rs` dB attentuation in the stopband.
+passband and `rs` dB attenuation in the stopband.
 """
 Elliptic(n::Integer, rp::Real, rs::Real) = Elliptic(Float64, n, rp, rs)
 
@@ -473,25 +473,25 @@ Here, `m` and `n` are respectively the numbers of zeros and poles in the s-domai
 then additional `n-m` zeros are added at `z = -1`.
 """
 function bilinear(f::ZeroPoleGain{:s,Z,P,K}, fs::Real) where {Z,P,K}
-    ztype = typeof(0 + zero(Z)/fs)
+    ztype = typeof(0 + zero(Z) / fs)
     z = fill(convert(ztype, -1), max(length(f.p), length(f.z)))
 
     ptype = typeof(0 + zero(P) / fs)
     p = Vector{ptype}(undef, length(f.p))
 
     num = one(one(fs) - one(Z))
     for i = 1:length(f.z)
-        z[i] = (2 + f.z[i] / fs)/(2 - f.z[i] / fs)
+        z[i] = (2 + f.z[i] / fs) / (2 - f.z[i] / fs)
         num *= (2 * fs - f.z[i])
     end
 
     den = one(one(fs) - one(P))
     for i = 1:length(f.p)
-        p[i] = (2 + f.p[i] / fs)/(2 - f.p[i]/fs)
+        p[i] = (2 + f.p[i] / fs) / (2 - f.p[i] / fs)
         den *= (2 * fs - f.p[i])
     end
 
-    ZeroPoleGain{:z}(z, p, f.k * real(num)/real(den))
+    ZeroPoleGain{:z}(z, p, f.k * real(num) / real(den))
 end
 
 # Pre-warp filter frequencies for digital filtering
@@ -545,12 +545,12 @@ end
 # Get length and alpha for Kaiser window filter with specified
 # transition band width and stopband attenuation in dB
 function kaiserord(transitionwidth::Real, attenuation::Real=60)
-    n = ceil(Int, (attenuation - 7.95)/(π*2.285*transitionwidth))+1
+    n = ceil(Int, (attenuation - 7.95) / (π * 2.285 * transitionwidth)) + 1
 
     if attenuation > 50
-        β = 0.1102*(attenuation - 8.7)
+        β = 0.1102 * (attenuation - 8.7)
     elseif attenuation >= 21
-        β = 0.5842*(attenuation - 21)^0.4 + 0.07886*(attenuation - 21)
+        β = 0.5842 * (attenuation - 21)^0.4 + 0.07886 * (attenuation - 21)
     else
         β = 0.0
     end
@@ -670,50 +670,47 @@ function digitalfilter(ftype::FilterType, proto::FIRWindow; fs::Real=2)
     coefs = firprototype(length(proto.window), ftype, fs)
     @assert length(proto.window) == length(coefs)
     out = (coefs .*= proto.window)
-    proto.scale ? rmul!(out, 1/scalefactor(out, ftype, fs)) : out
+    proto.scale ? rmul!(out, 1 / scalefactor(out, ftype, fs)) : out
 end
 
 
-# Compute FIR coefficients necessary for arbitrary rate resampling
-function resample_filter(rate::AbstractFloat, Nϕ = 32, rel_bw = 1.0, attenuation = 60)
-    f_nyq       = rate >= 1.0 ? 1.0/Nϕ : rate/Nϕ
-    cutoff      = f_nyq * rel_bw
-    trans_width = cutoff * 0.2
-
-    # Determine resampling filter order
-    hLen, α = kaiserord(trans_width, attenuation)
+"""
+    resample_filter(rate::AbstractFloat, Nϕ::Integer=32, rel_bw::Real=1.0, att=60)
 
-    # Round the number of taps up to a multiple of Nϕ.
-    # Otherwise the missing taps will be filled with 0.
-    hLen = Nϕ * ceil(Int, hLen/Nϕ)
+Compute FIR coefficients suitable for arbitrary rate resampling
+with `Nϕ` filters / phases, stop-band attenuation `att`, and relative bandwidth `rel_bw`.
+"""
+function resample_filter(rate::AbstractFloat, Nϕ::Integer=32, rel_bw::Real=1.0, attenuation=60)
+    f_nyq = rate >= 1.0 ? 1.0 / Nϕ : rate / Nϕ
+    return _resample_filter(f_nyq, Nϕ, rel_bw, attenuation)
+end
 
-    # Ensure that the filter is an odd length
-    if (iseven(hLen))
-        hLen += 1
-    end
+"""
+    resample_filter(rate::Union{Integer,Rational}, rel_bw::Real=1.0, attenuation=60)
 
-    # Design filter
-    h = digitalfilter(Lowpass(cutoff), FIRWindow(kaiser(hLen, α)))
-    rmul!(h, Nϕ)
+Compute FIR coefficients suitable for rational rate resampling,
+with stop-band attenuation `att`, and relative bandwidth `rel_bw`.
+"""
+function resample_filter(rate::Union{Integer,Rational}, rel_bw::Real=1.0, attenuation=60)
+    Nϕ         = numerator(rate)
+    decimation = denominator(rate)
+    f_nyq      = min(1 / Nϕ, 1 / decimation)
+    return _resample_filter(f_nyq, Nϕ, rel_bw, attenuation)
 end
 
-# Compute FIR coefficients necessary for rational rate resampling
-function resample_filter(rate::Union{Integer,Rational}, rel_bw = 1.0, attenuation = 60)
-    Nϕ          = numerator(rate)
-    decimation  = denominator(rate)
-    f_nyq       = min(1/Nϕ, 1/decimation)
+function _resample_filter(f_nyq, Nϕ, rel_bw, attenuation)
     cutoff      = f_nyq * rel_bw
     trans_width = cutoff * 0.2
 
     # Determine resampling filter order
     hLen, α = kaiserord(trans_width, attenuation)
 
-    # Round the number of taps up to a multiple of Nϕ (same as interpolation factor).
+    # Round the number of taps up to a multiple of Nϕ.
     # Otherwise the missing taps will be filled with 0.
-    hLen = Nϕ * ceil(Int, hLen/Nϕ)
+    hLen = Nϕ * ceil(Int, hLen / Nϕ)
 
     # Ensure that the filter is an odd length
-    if (iseven(hLen))
+    if iseven(hLen)
         hLen += 1
     end
 

diff --git a/src/Filters/filt.jl b/src/Filters/filt.jl
@@ -340,10 +340,10 @@ end
 # Zero phase digital filtering for second order sections
 function filtfilt(f::SecondOrderSections{:z,T,G}, x::AbstractArray{S}) where {T,G,S}
     zi = filt_stepstate(f)
-    pad_length = 6 * length(f.biquads)
+    pad_length = min(6 * length(f.biquads), size(x, 1) - 1)
     t = promote_type(T, G, S)
     zitmp = similar(zi, t)
-    extrapolated = Vector{t}(undef, size(x, 1)+pad_length*2)
+    extrapolated = Vector{t}(undef, size(x, 1) + pad_length * 2)
     out = similar(x, t)
 
     for col in CartesianIndices(axes(x)[2:end])