Skip to content
Open
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
35 changes: 26 additions & 9 deletions src/Filters/filt.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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)
Expand Down
12 changes: 7 additions & 5 deletions src/dspbase.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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

Expand Down
2 changes: 1 addition & 1 deletion src/periodograms.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down
Loading