general cleanups

docs/src/util.md

@@ -19,6 +19,7 @@ rms
 rmsfft
 meanfreq
 finddelay
+shiftin!
 shiftsignal
 shiftsignal!
 alignsignals

src/Filters/Filters.jl

@@ -7,8 +7,9 @@ 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 ..DSP: optimalfftfiltlength, os_fft_complexity, SMALL_FILT_CUTOFF
+import ..DSP: filt, filt!
+using Compat: Compat
 using FFTW
 
 include("coefficients.jl")

src/Filters/coefficients.jl

@@ -56,11 +56,9 @@ end
 # Transfer function form
 #
 
-function shiftpoly(p::LaurentPolynomial, i)
-    if i > 0
-        return p * LaurentPolynomial([one(eltype(p))], 1, indeterminate(p))^i
-    elseif i < 0
-        return p * LaurentPolynomial([one(eltype(p))], -1, indeterminate(p))^-i
+function shiftpoly(p::LaurentPolynomial{T}, i::Integer) where {T<:Number}
+    if !iszero(i)
+        return p * LaurentPolynomial([one(T)], i, indeterminate(p))
     end
     return p
 end
@@ -104,7 +102,7 @@ Filter with:
 ```math
 H(z) = \\frac{\\verb!b[1]! + \\ldots + \\verb!b[m]! z^{-m+1}}{\\verb!a[1]! + \\ldots + \\verb!a[n]! z^{-n+1}}
 ```
-returns `PolynomialRatio` object with `a[1] = 1` and other specified coefficients divided by `a[1]`. 
+returns `PolynomialRatio` object with `a[1] = 1` and other specified coefficients divided by `a[1]`.
 ```jldoctest
 julia> PolynomialRatio([1,1],[1,2])
 PolynomialRatio{:z, Float64}(LaurentPolynomial(1.0*z⁻¹ + 1.0), LaurentPolynomial(2.0*z⁻¹ + 1.0))
@@ -179,14 +177,17 @@ function Base.:^(f::PolynomialRatio{D,T}, e::Integer) where {D,T}
     end
 end
 
+coef_s(p::LaurentPolynomial{T}) where T = (n = firstindex(p);  append!(reverse!(coeffs(p)), zero(T) for _ in 1:n))
+coef_z(p::LaurentPolynomial{T}) where T = (n = -lastindex(p); prepend!(reverse!(coeffs(p)), zero(T) for _ in 1:n))
+
 """
     coefb(f)
 
 Coefficients of the numerator of a PolynomialRatio object, highest power
 first, i.e., the `b` passed to `filt()`
 """
-coefb(f::PolynomialRatio{:s}) = reverse([zeros(firstindex(f.b)); coeffs(f.b)])
-coefb(f::PolynomialRatio{:z}) = reverse([coeffs(f.b); zeros(-lastindex(f.b))])
+coefb(f::PolynomialRatio{:s}) = coef_s(f.b)
+coefb(f::PolynomialRatio{:z}) = coef_z(f.b)
 coefb(f::FilterCoefficients) = coefb(PolynomialRatio(f))
 
 """
@@ -195,8 +196,8 @@ coefb(f::FilterCoefficients) = coefb(PolynomialRatio(f))
 Coefficients of the denominator of a PolynomialRatio object, highest power
 first, i.e., the `a` passed to `filt()`
 """
-coefa(f::PolynomialRatio{:s}) = reverse([zeros(firstindex(f.a)); coeffs(f.a)])
-coefa(f::PolynomialRatio{:z}) = reverse([coeffs(f.a); zeros(-lastindex(f.a))])
+coefa(f::PolynomialRatio{:s}) = coef_s(f.a)
+coefa(f::PolynomialRatio{:z}) = coef_z(f.a)
 coefa(f::FilterCoefficients) = coefa(PolynomialRatio(f))
 
 #
@@ -324,25 +325,16 @@ function groupzp(z, p)
     groupedz = similar(z, n)
     i = 1
     while i <= n
-        closest_zero_idx = 0
-        closest_zero_val = Inf
-        for j = 1:length(z)
-            val = abs(z[j] - p[i])
-            if val < closest_zero_val
-                closest_zero_idx = j
-                closest_zero_val = val
-            end
-        end
+        p_i = p[i]
+        _, closest_zero_idx = findmin(x -> abs(x - p_i), z)
         groupedz[i] = splice!(z, closest_zero_idx)
         if !isreal(groupedz[i])
             i += 1
             groupedz[i] = splice!(z, closest_zero_idx)
         end
         i += 1
     end
-    ret = (groupedz, p[1:n])
-    splice!(p, 1:n)
-    ret
+    return (groupedz, splice!(p, 1:n))
 end
 
 # Sort zeros or poles lexicographically (so that poles are adjacent to
@@ -359,10 +351,10 @@ function split_real_complex(x::Vector{T}; sortby=nothing) where T
     end
 
     c = T[]
-    r = typeof(real(zero(T)))[]
+    r = real(T)[]
     ks = collect(keys(d))
     if sortby !== nothing
-        sort!(ks, by=sortby)
+        sort!(ks; by=sortby)
     end
     for k in ks
         if imag(k) != 0

src/Filters/design.jl

@@ -9,12 +9,13 @@ abstract type FilterType end
 #
 
 function Butterworth(::Type{T}, n::Integer) where {T<:Real}
-    n > 0 || error("n must be positive")
+    n > 0 || throw(DomainError(n, "n must be positive"))
 
-    poles = zeros(Complex{T}, n)
-    for i = 1:div(n, 2)
-        w = convert(T, 2i-1)/2n
-        pole = complex(-sinpi(w), cospi(w))
+    poles = Vector{Complex{T}}(undef, n)
+    for i = 1:n÷2
+        w = convert(T, 2i - 1) / 2n
+        sinpi_w, cospi_w = Compat.@inline sincospi(w)
+        pole = complex(-sinpi_w, cospi_w)
         poles[2i-1] = pole
         poles[2i] = conj(pole)
     end
@@ -36,32 +37,33 @@ Butterworth(n::Integer) = Butterworth(Float64, n)
 #
 
 function chebyshev_poles(::Type{T}, n::Integer, ε::Real) where {T<:Real}
-    p = zeros(Complex{T}, n)
+    p = Vector{Complex{T}}(undef, n)
     μ = asinh(convert(T, 1)/ε)/n
     b = -sinh(μ)
     c = cosh(μ)
-    for i = 1:div(n, 2)
-        w = convert(T, 2i-1)/2n
-        pole = complex(b*sinpi(w), c*cospi(w))
+    for i = 1:n÷2
+        w = convert(T, 2i - 1) / 2n
+        sinpi_w, cospi_w = Compat.@inline sincospi(w)
+        pole = complex(b * sinpi_w, c * cospi_w)
         p[2i-1] = pole
         p[2i] = conj(pole)
     end
     if isodd(n)
-        w = convert(T, 2*div(n, 2)+1)/2n
-        pole = b*sinpi(w)
+        w = convert(T, 2 * (n ÷ 2) + 1) / 2n
+        pole = b * sinpi(w)
         p[end] = pole
     end
     p
 end
 
 function Chebyshev1(::Type{T}, n::Integer, ripple::Real) where {T<:Real}
-    n > 0 || error("n must be positive")
-    ripple >= 0 || error("ripple must be non-negative")
+    n > 0 || throw(DomainError(n, "n must be positive"))
+    ripple >= 0 || throw(DomainError(ripple, "ripple must be non-negative"))
 
     ε = sqrt(10^(convert(T, ripple)/10)-1)
     p = chebyshev_poles(T, n, ε)
     k = one(T)
-    for i = 1:div(n, 2)
+    for i = 1:n÷2
         k *= abs2(p[2i])
     end
     if iseven(n)
@@ -81,21 +83,21 @@ the passband.
 Chebyshev1(n::Integer, ripple::Real) = Chebyshev1(Float64, n, ripple)
 
 function Chebyshev2(::Type{T}, n::Integer, ripple::Real) where {T<:Real}
-    n > 0 || error("n must be positive")
-    ripple >= 0 || error("ripple must be non-negative")
+    n > 0 || throw(DomainError(n, "n must be positive"))
+    ripple >= 0 || throw(DomainError(ripple, "ripple must be non-negative"))
 
     ε = 1/sqrt(10^(convert(T, ripple)/10)-1)
     p = chebyshev_poles(T, n, ε)
     map!(inv, p, p)
 
-    z = zeros(Complex{T}, n-isodd(n))
+    z = Vector{Complex{T}}(undef, n - isodd(n))
     k = one(T)
-    for i = 1:div(n, 2)
-        w = convert(T, 2i-1)/2n
-        ze = complex(zero(T), -inv(cospi(w)))
+    for i = 1:n÷2
+        w = convert(T, 2i - 1) / 2n
+        ze = Compat.@inline complex(zero(T), -inv(cospi(w)))
         z[2i-1] = ze
         z[2i] = conj(ze)
-        k *= abs2(p[2i])/abs2(ze)
+        k *= abs2(p[2i]) / abs2(ze)
     end
     isodd(n) && (k *= -real(p[end]))
 
@@ -130,19 +132,16 @@ end
 # cde computes cd(u*K(k), k)
 # sne computes sn(u*K(k), k)
 # Both accept the Landen sequence as generated by landen above
-for (fn, init) in ((:cde, :(cospi(u/2))), (:sne, :(sinpi(u/2))))
-    @eval begin
-        function $fn(u::Number, landen::Vector{T}) where T<:Real
-            winv = inv($init)
-            # Eq. (55)
-            for i = length(landen):-1:1
-                oldwinv = winv
-                winv = 1/(1+landen[i])*(winv+landen[i]/winv)
-            end
-            w = inv(winv)
-        end
+function _ellip(init::Number, landen::Vector{<:Real})
+    winv = inv(init)
+    # Eq. (55)
+    for x in Iterators.reverse(landen)
+        winv = 1 / (1 + x) * (winv + x / winv)
     end
+    w = inv(winv)
 end
+@inline cde(u::Number, landen::Vector{<:Real}) = Compat.@inline _ellip(cospi(u / 2), landen)
+@inline sne(u::Number, landen::Vector{<:Real}) = Compat.@inline _ellip(sinpi(u / 2), landen)
 
 # sne inverse
 function asne(w::Number, k::Real)
@@ -159,17 +158,17 @@ function asne(w::Number, k::Real)
 end
 
 function Elliptic(::Type{T}, n::Integer, rp::Real, rs::Real) where {T<:Real}
-    n > 0 || error("n must be positive")
-    rp > 0 || error("rp must be positive")
-    rp < rs || error("rp must be less than rs")
+    n > 0 || throw(DomainError(n, "n must be positive"))
+    rp > 0 || throw(DomainError(rp, "rp must be positive"))
+    rp < rs || throw(DomainError(rp, "rp must be less than rs"))
 
     # Eq. (2)
     εp = sqrt(10^(convert(T, rp)/10)-1)
     εs = sqrt(10^(convert(T, rs)/10)-1)
 
     # Eq. (3)
     k1 = εp/εs
-    k1 >= 1 && error("filter order is too high for parameters")
+    k1 >= 1 && throw(ArgumentError("filter order is too high for parameters"))
 
     # Eq. (20)
     k1′² = 1 - abs2(k1)
@@ -178,7 +177,7 @@ function Elliptic(::Type{T}, n::Integer, rp::Real, rs::Real) where {T<:Real}
 
     # Eq. (47)
     k′ = one(T)
-    for i = 1:div(n, 2)
+    for i = 1:n÷2
         k′ *= sne(convert(T, 2i-1)/n, k1′_landen)
     end
     k′ = k1′²^(convert(T, n)/2)*k′^4
@@ -189,10 +188,10 @@ function Elliptic(::Type{T}, n::Integer, rp::Real, rs::Real) where {T<:Real}
     # Eq. (65)
     v0 = -im/convert(T, n)*asne(im/εp, k1)
 
-    z = Vector{Complex{T}}(undef, 2*div(n, 2))
+    z = Vector{Complex{T}}(undef, 2*(n÷2))
     p = Vector{Complex{T}}(undef, n)
     gain = one(T)
-    for i = 1:div(n, 2)
+    for i = 1:n÷2
         # Eq. (43)
         w = convert(T, 2i-1)/n
 
@@ -203,8 +202,8 @@ function Elliptic(::Type{T}, n::Integer, rp::Real, rs::Real) where {T<:Real}
 
         # Eq. (64)
         pole = im*cde(w - im*v0, k_landen)
-        p[2i] = pole
         p[2i-1] = conj(pole)
+        p[2i] = pole
 
         gain *= abs2(pole)/abs2(ze)
     end
@@ -234,9 +233,9 @@ Elliptic(n::Integer, rp::Real, rs::Real) = Elliptic(Float64, n, rp, rs)
 
 # returns frequency in half-cycles per sample ∈ (0, 1)
 function normalize_freq(w::Real, fs::Real)
-    w <= 0 && error("frequencies must be positive")
+    w <= 0 && throw(DomainError(w, "frequencies must be positive"))
     f = 2*w/fs
-    f >= 1 && error("frequencies must be less than the Nyquist frequency $(fs/2)")
+    f >= 1 && throw(DomainError(f, "frequencies must be less than the Nyquist frequency $(fs/2)"))
     f
 end
 
@@ -273,7 +272,7 @@ end
 Band pass filter with pass band frequencies (`Wn1`, `Wn2`).
 """
 function Bandpass(w1::Real, w2::Real)
-    w1 < w2 || error("w1 must be less than w2")
+    w1 < w2 || throw(ArgumentError("w1 must be less than w2"))
     Bandpass{Base.promote_typeof(w1/1, w2/1)}(w1, w2)
 end
 
@@ -288,7 +287,7 @@ end
 Band stop filter with stop band frequencies (`Wn1`, `Wn2`).
 """
 function Bandstop(w1::Real, w2::Real)
-    w1 < w2 || error("w1 must be less than w2")
+    w1 < w2 || throw(ArgumentError("w1 must be less than w2"))
     Bandstop{Base.promote_typeof(w1/1, w2/1)}(w1, w2)
 end
 
@@ -573,22 +572,22 @@ FIRWindow(; transitionwidth::Real=throw(ArgumentError("must specify transitionwi
 function firprototype(n::Integer, ftype::Lowpass, fs::Real)
     w = normalize_freq(ftype.w, fs)
 
-    [w*sinc(w*(k-(n-1)/2)) for k = 0:(n-1)]
+    [w*sinc(w*(k-(n+1)/2)) for k = 1:n]
 end
 
 function firprototype(n::Integer, ftype::Bandpass, fs::Real)
     w1 = normalize_freq(ftype.w1, fs)
     w2 = normalize_freq(ftype.w2, fs)
 
-    [w2*sinc(w2*(k-(n-1)/2)) - w1*sinc(w1*(k-(n-1)/2)) for k = 0:(n-1)]
+    [w2*sinc(w2*(k-(n+1)/2)) - w1*sinc(w1*(k-(n+1)/2)) for k = 1:n]
 end
 
 function firprototype(n::Integer, ftype::Highpass, fs::Real)
     w = normalize_freq(ftype.w, fs)
     isodd(n) || throw(ArgumentError("FIRWindow highpass filters must have an odd number of coefficients"))
 
-    out = [-w*sinc(w*(k-(n-1)/2)) for k = 0:(n-1)]
-    out[div(n, 2)+1] += 1
+    out = [-w*sinc(w*(k-(n+1)/2)) for k = 1:n]
+    out[n÷2+1] += 1
     out
 end
 
@@ -597,25 +596,25 @@ function firprototype(n::Integer, ftype::Bandstop, fs::Real)
     w2 = normalize_freq(ftype.w2, fs)
     isodd(n) || throw(ArgumentError("FIRWindow bandstop filters must have an odd number of coefficients"))
 
-    out = [w1*sinc(w1*(k-(n-1)/2)) - w2*sinc(w2*(k-(n-1)/2)) for k = 0:(n-1)]
-    out[div(n, 2)+1] += 1
+    out = [w1*sinc(w1*(k-(n+1)/2)) - w2*sinc(w2*(k-(n+1)/2)) for k = 1:n]
+    out[n÷2+1] += 1
     out
 end
 
 scalefactor(coefs::Vector, ::Union{Lowpass, Bandstop}, fs::Real) = sum(coefs)
-function scalefactor(coefs::Vector, ::Highpass, fs::Real)
-    c = zero(coefs[1])
+function scalefactor(coefs::Vector{T}, ::Highpass, fs::Real) where T
+    c = zero(T)
     for k = 1:length(coefs)
         c += ifelse(isodd(k), coefs[k], -coefs[k])
     end
     c
 end
-function scalefactor(coefs::Vector, ftype::Bandpass, fs::Real)
+function scalefactor(coefs::Vector{T}, ftype::Bandpass, fs::Real) where T
     n = length(coefs)
     freq = normalize_freq(middle(ftype.w1, ftype.w2), fs)
-    c = zero(coefs[1])
-    for k = 0:n-1
-        c += coefs[k+1]*cospi(freq*(k-(n-1)/2))
+    c = zero(T)
+    for k = 1:n
+        c = muladd(coefs[k], cospi(freq * (k - (n + 1) / 2)), c)
     end
     c
 end