Skip to content

Commit

Permalink
Let _filt_iir! handle a, b of different length
Browse files Browse the repository at this point in the history
  • Loading branch information
martinholters committed Dec 16, 2024
1 parent 3e79857 commit 878281a
Show file tree
Hide file tree
Showing 2 changed files with 15 additions and 13 deletions.
9 changes: 1 addition & 8 deletions src/Filters/filt.jl
Original file line number Diff line number Diff line change
Expand Up @@ -162,14 +162,7 @@ function filt!(out::AbstractArray{<:Any,N}, f::DF2TFilter{<:PolynomialRatio,Arra
mul!(out, x, b[1])
else
a = coefa(f.coef)
as = length(a)
bs = length(b)
if as != 1
if as < n
append!(a, zero(eltype(a)) for _ in 1:(n-as))
elseif bs < n
append!(b, zero(eltype(b)) for _ in 1:(n-bs))
end
if length(a) != 1
for col in CartesianIndices(axes(x)[2:end])
_filt_iir!(out, b, a, x, view(si, :, col), col)
end
Expand Down
19 changes: 14 additions & 5 deletions src/dspbase.jl
Original file line number Diff line number Diff line change
Expand Up @@ -45,9 +45,6 @@ function filt!(out::AbstractArray, b::Union{AbstractVector, Number}, a::Union{Ab
a = @noinline broadcast(/, a, norml)
b = @noinline broadcast(/, b, norml)
end
# Pad the coefficients with zeros if needed
bs<sz && (b = copyto!(zeros(eltype(b), sz), b))
1<as<sz && (a = copyto!(zeros(eltype(a), sz), a))

si = Vector{promote_type(eltype(b), eltype(a), T)}(undef, sz - 1)

Expand Down Expand Up @@ -75,10 +72,22 @@ function _filt_iir!(out, b, a, x, si, col)
xi = x[i, col]
val = muladd(xi, b[1], si[1])
out[i, col] = val
for j=1:(silen-1)
for j in 1:min(length(a), length(b), silen) - 1
si[j] = muladd(val, -a[j+1], muladd(xi, b[j+1], si[j+1]))
end
si[silen] = muladd(xi, b[silen+1], -a[silen+1]*val)
if length(a) == length(b)
si[silen] = muladd(xi, b[silen+1], -a[silen+1]*val)
elseif length(a) > length(b)
for j in length(b):silen-1
si[j] = muladd(val, -a[j+1], si[j+1])
end
si[silen] = -a[silen+1]*val
else
for j in length(a):silen-1
si[j] = muladd(xi, b[j+1], si[j+1])
end
si[silen] = xi*b[silen+1]
end
end
end

Expand Down

0 comments on commit 878281a

Please sign in to comment.