Simplify extrapolate_signal!, use CartesianIndices in filtfilt,…

… housekeeping (#611)

* use `findfirst!`, type stability in `_freqrange`

also `setindex!`, pushfirst!

* LazyString for errors

* remove `UnitRange`s

* rearrangements

* exp10, one

typos

* zero once

= not +=

* multiline string

* while to for loops

* remove inbounds

* one line loop

* CartesianIndex() for iir_filtfilt

* Fill once

* simplify `extrapolate_signal!`

* restrict `dsp_reverse`

* reuse extrapolated

* ensure newb is mutable

src/Filters/design.jl

@@ -38,7 +38,7 @@ Butterworth(n::Integer) = Butterworth(Float64, n)
 
 function chebyshev_poles(::Type{T}, n::Integer, ε::Real) where {T<:Real}
     p = Vector{Complex{T}}(undef, n)
-    μ = asinh(convert(T, 1)/ε)/n
+    μ = asinh(one(T) / ε) / n
     b = -sinh(μ)
     c = cosh(μ)
     for i = 1:n÷2
@@ -60,7 +60,7 @@ function Chebyshev1(::Type{T}, n::Integer, ripple::Real) where {T<:Real}
     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)
+    ε = sqrt(exp10(convert(T, ripple) / 10) - 1)
     p = chebyshev_poles(T, n, ε)
     k = one(T)
     for i = 1:n÷2
@@ -86,7 +86,7 @@ function Chebyshev2(::Type{T}, n::Integer, ripple::Real) where {T<:Real}
     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)
+    ε = 1 / sqrt(exp10(convert(T, ripple) / 10) - 1)
     p = chebyshev_poles(T, n, ε)
     map!(inv, p, p)
 
@@ -163,8 +163,8 @@ function Elliptic(::Type{T}, n::Integer, rp::Real, rs::Real) where {T<:Real}
     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)
+    εp = sqrt(exp10(convert(T, rp) / 10) - 1)
+    εs = sqrt(exp10(convert(T, rs) / 10) - 1)
 
     # Eq. (3)
     k1 = εp/εs
@@ -178,22 +178,22 @@ function Elliptic(::Type{T}, n::Integer, rp::Real, rs::Real) where {T<:Real}
     # Eq. (47)
     k′ = one(T)
     for i = 1:n÷2
-        k′ *= sne(convert(T, 2i-1)/n, k1′_landen)
+        k′ *= sne(convert(T, 2i - 1) / n, k1′_landen)
     end
-    k′ = k1′²^(convert(T, n)/2)*k′^4
+    k′ = k1′²^(convert(T, n) / 2) * k′^4
 
     k = sqrt(1 - abs2(k′))
     k_landen = landen(k)
 
     # Eq. (65)
-    v0 = -im/convert(T, n)*asne(im/εp, k1)
+    v0 = -im / convert(T, n) * asne(im / εp, k1)
 
     z = Vector{Complex{T}}(undef, 2*(n÷2))
     p = Vector{Complex{T}}(undef, n)
     gain = one(T)
     for i = 1:n÷2
         # Eq. (43)
-        w = convert(T, 2i-1)/n
+        w = convert(T, 2i - 1) / n
 
         # Eq. (62)
         ze = complex(zero(T), -inv(k*cde(w, k_landen)))
@@ -213,7 +213,7 @@ function Elliptic(::Type{T}, n::Integer, rp::Real, rs::Real) where {T<:Real}
         p[end] = pole
         gain *= abs(pole)
     else
-        gain *= 10^(-convert(T, rp)/20)
+        gain *= exp10(-convert(T, rp) / 20)
     end
 
     ZeroPoleGain{:s}(z, p, gain)
@@ -235,12 +235,12 @@ Elliptic(n::Integer, rp::Real, rs::Real) = Elliptic(Float64, n, rp, rs)
 function normalize_freq(w::Real, fs::Real)
     w <= 0 && throw(DomainError(w, "frequencies must be positive"))
     f = 2 * w / fs
-    f >= 1 && throw(DomainError(f, "frequencies must be less than the Nyquist frequency $(fs/2)"))
+    f >= 1 && throw(DomainError(f, lazy"frequencies must be less than the Nyquist frequency $(fs/2)"))
     f
 end
 function normalize_complex_freq(w::Real, fs::Real)
     f = 2 * w / fs
-    f >= 2 && throw(DomainError(f, "frequencies must be less than the sampling frequency $(fs)"))
+    f >= 2 && throw(DomainError(f, lazy"frequencies must be less than the sampling frequency $(fs)"))
     f
 end
 
@@ -407,7 +407,7 @@ function transform_prototype(ftype::Bandstop, proto::ZeroPoleGain{:s})
         newp[2i] = b + pm
     end
 
-    # Any emaining poles/zeros are real and not cancelled
+    # Any remaining poles/zeros are real and not cancelled
     npm = sqrt(-complex(ftype.w2 * ftype.w1))
     for (n, newc) in ((np, newp), (nz, newz))
         for i = n+1:npairs
@@ -476,8 +476,8 @@ function bilinear(f::ZeroPoleGain{:s,Z,P,K}, fs::Real) where {Z,P,K}
     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{typeof(zero(P)/fs)}(undef, length(f.p))
+    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)
@@ -534,7 +534,8 @@ function iirnotch(w::Real, bandwidth::Real; fs=2)
     b = 1 / (1 + tanpi(bandwidth / 2))
     # Eq. 8.2.22
     cosw0 = cospi(w)
-    Biquad(b, -2b*cosw0, b, -2b*cosw0, 2b-1)
+    b1 = -2b * cosw0
+    Biquad(b, b1, b, b1, 2b-1)
 end
 
 #

src/Filters/filt.jl

@@ -36,17 +36,16 @@ function _filt!(out::AbstractArray, si::AbstractArray{S,N}, f::SecondOrderSectio
                 x::AbstractArray, col::Union{Int,CartesianIndex}) where {S,N}
     g = f.g
     biquads = f.biquads
-    n = length(biquads)
     @inbounds for i in axes(x, 1)
         yi = x[i, col]
-        for fi = 1:n
+        for fi in eachindex(biquads)
             biquad = biquads[fi]
             xi = yi
             yi = muladd(biquad.b0, xi, si[1, fi])
             si[1, fi] = muladd(biquad.a1, -yi, muladd(biquad.b1, xi, si[2, fi]))
             si[2, fi] = muladd(biquad.b2, xi, -biquad.a2 * yi)
         end
-        out[i, col] = yi*g
+        out[i, col] = g * yi
     end
     si
 end
@@ -250,16 +249,17 @@ DF2TFilter(coef::FilterCoefficients{:z}, ::Type{V}, coldims::Tuple{Vararg{Intege
 #
 # in place in output. The istart and n parameters determine the portion
 # of the input signal x to extrapolate.
-function extrapolate_signal!(out, ostart, sig, istart, n, pad_length)
+function extrapolate_signal!(out, sig, pad_length)
+    n = length(sig)
     length(out) >= n + 2 * pad_length || throw(ArgumentError("output is incorrectly sized"))
-    x = 2 * sig[istart]
-    for i = 0:pad_length-1
-        out[ostart+i] = x - sig[istart+pad_length-i]
+    x = 2 * sig[1]
+    for i = 1:pad_length
+        out[i] = x - sig[2+pad_length-i]
     end
-    copyto!(out, ostart+pad_length, sig, istart, n)
-    x = 2 * sig[istart+n-1]
+    copyto!(out, 1+pad_length, sig, 1, n)
+    x = 2 * sig[n]
     for i = 1:pad_length
-        out[ostart+n+pad_length+i-1] = x - sig[istart+n-1-i]
+        out[n+pad_length+i] = x - sig[n-i]
     end
     out
 end
@@ -273,14 +273,13 @@ function iir_filtfilt(b::AbstractVector, a::AbstractVector, x::AbstractArray)
     extrapolated = Vector{t}(undef, size(x, 1) + 2 * pad_length)
     out = similar(x, t)
 
-    for i = 1:Base.trailingsize(x, 2)
-        istart = 1 + (i - 1) * size(x, 1)
-        extrapolate_signal!(extrapolated, 1, x, istart, size(x, 1), pad_length)
-        _filt_iir!(extrapolated, bn, an, extrapolated, mul!(zitmp, zi, extrapolated[1]), 1)
+    for col in CartesianIndices(axes(x)[2:end])
+        extrapolate_signal!(extrapolated, view(x, :, col), pad_length)
+        _filt_iir!(extrapolated, bn, an, extrapolated, mul!(zitmp, zi, extrapolated[1]), CartesianIndex())
         reverse!(extrapolated)
-        _filt_iir!(extrapolated, bn, an, extrapolated, mul!(zitmp, zi, extrapolated[1]), 1)
-        for j = 1:size(x, 1)
-            out[j, i] = extrapolated[end-pad_length+1-j]
+        _filt_iir!(extrapolated, bn, an, extrapolated, mul!(zitmp, zi, extrapolated[1]), CartesianIndex())
+        for j in axes(x, 1)
+            out[j, col] = extrapolated[end-pad_length+1-j]
         end
     end
 
@@ -310,25 +309,26 @@ function filtfilt(b::AbstractVector, x::AbstractArray)
     nb = length(b)
     # Only need as much padding as the order of the filter
     T = Base.promote_eltype(b, x)
-    extrapolated = similar(x, T, size(x, 1)+2nb-2, Base.trailingsize(x, 2))
+    extrapolated = similar(x, T, (size(x, 1)+2*(nb-1), size(x)[2:end]...))
 
     # Convolve b with its reverse
-    newb = filt(b, reverse(b))
+    newb = reverse(b, 1)
+    filt!(newb, b, newb)
     resize!(newb, 2nb-1)
     for i = 1:nb-1
         newb[nb+i] = newb[nb-i]
     end
 
     # Extrapolate each column
-    for i = 1:size(extrapolated, 2)
-        extrapolate_signal!(extrapolated, (i-1)*size(extrapolated, 1)+1, x, (i-1)*size(x, 1)+1, size(x, 1), nb-1)
+    for col in CartesianIndices(axes(x)[2:end])
+        @views extrapolate_signal!(extrapolated[:, col], x[:, col], nb-1)
     end
 
     # Filter
-    out = filt(newb, extrapolated)
+    filt!(extrapolated, newb, extrapolated)
 
     # Drop garbage at start
-    reshape(out[2nb-1:end, :], size(x))
+    return extrapolated[2nb-1:end, axes(extrapolated)[2:end]...]
 end
 
 # Choose whether to use FIR or iir_filtfilt depending on
@@ -353,16 +353,14 @@ function filtfilt(f::SecondOrderSections{:z,T,G}, x::AbstractArray{S}) where {T,
     extrapolated = Vector{t}(undef, size(x, 1)+pad_length*2)
     out = similar(x, t)
 
-    istart = 1
-    for i = 1:Base.trailingsize(x, 2)
-        extrapolate_signal!(extrapolated, 1, x, istart, size(x, 1), pad_length)
+    for col in CartesianIndices(axes(x)[2:end])
+        extrapolate_signal!(extrapolated, view(x, :, col), pad_length)
         _filt!(extrapolated, mul!(zitmp, zi, extrapolated[1]), f, extrapolated, CartesianIndex())
         reverse!(extrapolated)
         _filt!(extrapolated, mul!(zitmp, zi, extrapolated[1]), f, extrapolated, CartesianIndex())
-        for j = 1:size(x, 1)
-            @inbounds out[j, i] = extrapolated[end-pad_length+1-j]
+        for j in axes(x, 1)
+            out[j, col] = extrapolated[end-pad_length+1-j]
         end
-        istart += size(x, 1)
     end
 
     out
@@ -413,7 +411,7 @@ function filt_stepstate(f::SecondOrderSections{:z,T}) where T
     biquads = f.biquads
     si = Matrix{T}(undef, 2, length(biquads))
     y = one(T)
-    for i = 1:length(biquads)
+    for i in eachindex(biquads)
         biquad = biquads[i]
         a1, a2, b0, b1, b2 = biquad.a1, biquad.a2, biquad.b0, biquad.b1, biquad.b2
 
@@ -515,10 +513,9 @@ function _fftfilt!(
         xstart = off - nb + npadbefore + 1
         n = min(nfft - npadbefore, nx - xstart + 1)
 
-        tmp1[1:npadbefore] .= zero(W)
-        tmp1[npadbefore+n+1:end] .= zero(W)
-
+        fill!(tmp1, zero(W))
         copyto!(tmp1, npadbefore+1, x, colstart+xstart, n)
+
         mul!(tmp2, p1, tmp1)
         broadcast!(*, tmp2, tmp2, filterft)
         mul!(tmp1, p2, tmp2)

src/Filters/filt_order.jl

@@ -122,16 +122,14 @@ function brent(f, x1::T, x2::T) where {T<:AbstractFloat}
     m1, m2 = m, m
     fm = f(m)
     fm1, fm2 = fm, fm
-    iter = 0
 
-    while (iter < 1_000)
+    for _ in 1:1000
         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)

src/Filters/remez_fir.jl

@@ -100,11 +100,9 @@ CALCULATE THE LAGRANGE INTERPOLATION COEFFICIENTS
 function lagrange_interp(k::Integer, n::Integer, m::Integer, x::AbstractVector)
     retval = 1.0
     q = x[k]
-    for l = 1 : m
-        for j = l : m : n
-            if j != k
-                retval *= 2.0 * (q - x[j])
-            end
+    for l = 1:m, j = l:m:n
+        if j != k
+            retval *= 2.0 * (q - x[j])
         end
     end
     1.0 / retval

src/Filters/response.jl

@@ -125,7 +125,7 @@ end
 Impulse response of a digital `filter` with `n` points.
 """
 function impresp(filter::FilterCoefficients{:z}, n=100)
-  i = [1; zeros(n-1)]
+  i = setindex!(zeros(n), 1, 1)
   filt(filter, i)
 end
 
@@ -160,16 +160,16 @@ function _freqrange(filter::FilterCoefficients{:s})
     filter = convert(ZeroPoleGain, filter)
     w_interesting = sort!(Float64.(abs.([filter.p; filter.z])))
     include_zero = !isempty(w_interesting) && iszero(w_interesting[1])
-    w_interesting = collect(Iterators.dropwhile(iszero, w_interesting))
-    if isempty(w_interesting) # no non-zero poles or zeros
-        if !include_zero || !isfinite(1/filter.k)
-            return [0; 10 .^ (0:6)] # fallback
+    first_nonzero = findfirst(!iszero, w_interesting)
+    if isnothing(first_nonzero) # no non-zero poles or zeros
+        if !include_zero || !isfinite(1 / filter.k)
+            return setindex!(10.0 .^ (-1:6), 0.0, 1) # fallback
         end
         # include the point where |H|=1 (if any) and go further by factor 10
-        return range(0.0, stop=10 * Float64(max(filter.k, 1/filter.k)), length=200)
+        return collect(range(0.0, stop=10 * Float64(max(filter.k, 1 / filter.k)), length=200))
     end
-    # normal case: go from smalles to largest pole/zero, extended by factor 10
-    w_min, w_max = w_interesting[[1,end]]
-    w = 10 .^ range(log10(w_min)-1, stop=log10(w_max)+1, length=200)
-    return include_zero ? [0.0; w] : w
+    # normal case: go from smallest to largest pole/zero, extended by factor 10
+    w_min, w_max = w_interesting[first_nonzero], w_interesting[end]
+    w = 10 .^ range(log10(w_min) - 1, stop=log10(w_max) + 1, length=200)
+    return include_zero ? pushfirst!(w, 0.0) : w
 end

src/dspbase.jl

@@ -615,7 +615,7 @@ function _conv_kern_fft!(out::AbstractArray{T, N},
     copyto!(out,
             output_indices,
             raw_out,
-            CartesianIndices(UnitRange.(1, outsize)))
+            CartesianIndices(outsize))
 end
 function _conv_kern_fft!(out::AbstractArray{T}, output_indices, u, v) where {T}
     outsize = size(output_indices)
@@ -631,7 +631,7 @@ function _conv_kern_fft!(out::AbstractArray{T}, output_indices, u, v) where {T}
     copyto!(out,
             output_indices,
             upad,
-            CartesianIndices(UnitRange.(1, outsize)))
+            CartesianIndices(outsize))
 end
 
 function _conv_td!(out, output_indices, u::AbstractArray{<:Number, N}, v::AbstractArray{<:Number, N}) where {N}
@@ -652,7 +652,7 @@ end
 
 # whether the given axis are to be considered to carry an offset for `conv!` and `conv`
 conv_axis_with_offset(::Base.OneTo) = false
-conv_axis_with_offset(a::Any) = throw(ArgumentError("unsupported axis type $(typeof(a))"))
+conv_axis_with_offset(a::Any) = throw(ArgumentError(lazy"unsupported axis type $(typeof(a))"))
 
 function conv_axes_with_offset(as::Tuple...)
     with_offset = ((map(a -> map(conv_axis_with_offset, a), as)...)...,)
@@ -808,8 +808,8 @@ function conv(u::AbstractVector{T}, v::Transpose{T,<:AbstractVector}, A::Abstrac
 end
 
 
-dsp_reverse(v, ::NTuple{<:Any, Base.OneTo{Int}}) = reverse(v, dims = 1)
-function dsp_reverse(v, vaxes)
+dsp_reverse(v::AbstractVector, ::Tuple{Base.OneTo}) = reverse(v)
+function dsp_reverse(v::AbstractVector, vaxes::Tuple{AbstractUnitRange})
     vsize = length(v)
     reflected_start = - first(only(vaxes)) - vsize + 1
     reflected_axes = (reflected_start : reflected_start + vsize - 1,)

src/estimation.jl

@@ -166,11 +166,11 @@ function quinn(x::Vector{<:Real}, f0::Real, Fs::Real; tol=1e-6, maxiters::Intege
 
     # Initialize algorithm
     α = 2cos(ω̂)
-
-    # iteration
-    ξ = zeros(eltype(x), N)
     β = zero(eltype(x))
+    ξ = zeros(eltype(x), N)
     ξ[1] = x[1]
+
+    # iteration
     iter = 0
     for outer iter = 1:maxiters
         ξ[2] = muladd(α, ξ[1], x[2])
@@ -194,11 +194,12 @@ function quinn(x::Vector{<:Complex}, f0::Real, Fs::Real; tol=1e-6, maxiters::Int
 
     # Remove any DC term in s
     x = x .- mean(x)
-    N = length(x)
 
-    # iteration
+    N = length(x)
     ξ = zeros(eltype(x), N)
     ξ[1] = x[1]
+
+    # iteration
     iter = 0
     for outer iter = 1:maxiters
         S = zero(eltype(x))