diff --git a/src/Filters/filt.jl b/src/Filters/filt.jl index aa214650..e3e0238d 100644 --- a/src/Filters/filt.jl +++ b/src/Filters/filt.jl @@ -32,8 +32,8 @@ filt(f::PolynomialRatio{:z}, x) = filt(coefb(f), coefa(f), x) ## SecondOrderSections # filt! algorithm (no checking, returns si) -function _filt!(out::AbstractArray, si::AbstractArray{S,N}, f::SecondOrderSections{:z}, - x::AbstractArray, col::Union{Int,CartesianIndex}) where {S,N} +function _filt!(out::AbstractArray, si::AbstractMatrix{S}, f::SecondOrderSections{:z}, + x::AbstractArray, col::CartesianIndex) where {S} g = f.g biquads = f.biquads @inbounds for i in axes(x, 1) @@ -68,7 +68,7 @@ filt(f::SecondOrderSections{:z,T,G}, x::AbstractArray{S}) where {T,G,S<:Number} # filt! algorithm (no checking, returns si) function _filt!(out::AbstractArray, si1::Number, si2::Number, f::Biquad{:z}, - x::AbstractArray, col::Union{Int,CartesianIndex}) + x::AbstractArray, col::CartesianIndex) @inbounds for i in axes(x, 1) xi = x[i, col] yi = muladd(f.b0, xi, si1) @@ -162,23 +162,40 @@ function filt!(out::AbstractArray{<:Any,N}, f::DF2TFilter{<:PolynomialRatio,Arra mul!(out, x, b[1]) else a = coefa(f.coef) + outer_dims = CartesianIndices(axes(x)[2:end]) + if length(a) != 1 - for col in CartesianIndices(axes(x)[2:end]) + for col in outer_dims _filt_iir!(out, b, a, x, view(si, :, col), col) end elseif n <= SMALL_FILT_CUTOFF - vtup = ntuple(j -> b[j], Val(length(b))) - for col in CartesianIndices(axes(x)[2:end]) - _filt_fir!(out, vtup, x, view(si, :, col), col, Val(true)) - end + _small_filt_fir_storesi!(out, b, x, si, Val(length(b))) else - for col in CartesianIndices(axes(x)[2:end]) + for col in outer_dims _filt_fir!(out, b, x, view(si, :, col), col) end end end return out end +function _small_filt_fir_storesi!( + out::AbstractArray, h::AbstractVector, x::AbstractArray, + si::Array, ::Val{bs}) where {bs} + + bs < 2 && throw(ArgumentError("invalid tuple size")) + length(h) != bs && throw(ArgumentError("length(h) does not match bs")) + b = ntuple(j -> h[j], Val(bs)) + col_len = size(si, 1) + LI = LinearIndices(si) + for col in CartesianIndices(axes(x)[2:end]) + @static if VERSION > v"1.12-" + si_vec = Base.wrap(Array, memoryref(si.ref, LI[1, col]), col_len) + else + si_vec = view(si, :, col) + end + _filt_fir!(out, b, x, si_vec, col, Val(true)) + end +end ## SecondOrderSections DF2TFilter(coef::SecondOrderSections{:z,T}, coldims::Tuple{Vararg{Integer}}) where {T} = DF2TFilter(coef, T, coldims) diff --git a/src/dspbase.jl b/src/dspbase.jl index 7f97a96d..7ac711e2 100644 --- a/src/dspbase.jl +++ b/src/dspbase.jl @@ -115,21 +115,23 @@ end const SMALL_FILT_VECT_CUTOFF = 19 # Transposed direct form II -@generated function _filt_fir!(out, b::NTuple{N,T}, x, siarr, col, ::Val{StoreSI}=Val(false)) where {N,T,StoreSI} +@generated function _filt_fir!(out, b::NTuple{N,T}, x, siarr, col, ::Val{StoreSI}) where {N,T,StoreSI} silen = N - 1 si_end = Symbol(:si_, silen) quote N <= SMALL_FILT_VECT_CUTOFF && checkbounds(siarr, $silen) Base.@nextract $silen si siarr + VECTORIZE_LARGER = N > SMALL_FILT_VECT_CUTOFF for i in axes(x, 1) xi = x[i, col] val = muladd(xi, b[1], si_1) + if VECTORIZE_LARGER + @inbounds out[i, col] = val + end Base.@nexprs $(silen - 1) j -> (si_j = muladd(xi, b[j+1], si_{j + 1})) $si_end = xi * b[N] - if N > SMALL_FILT_VECT_CUTOFF - @inbounds out[i, col] = val - else + if !VECTORIZE_LARGER out[i, col] = val end end @@ -149,7 +151,7 @@ function _small_filt_fir!( length(h) != bs && throw(ArgumentError("length(h) does not match bs")) b = ntuple(j -> h[j], Val(bs)) for col in CartesianIndices(axes(x)[2:end]) - _filt_fir!(out, b, x, si, col) + _filt_fir!(out, b, x, si, col, Val(false)) end end diff --git a/src/periodograms.jl b/src/periodograms.jl index 6bbce435..9c1d1b5c 100644 --- a/src/periodograms.jl +++ b/src/periodograms.jl @@ -482,7 +482,7 @@ function periodogram(s::AbstractMatrix{T}; ptype = 0 elseif radialsum ptype = 1 - elseif radialavg + else # if radialavg ptype = 2 end norm2 = length(s)