Skip to content

Commit 73d58fb

Browse files
Let _filt_iir! handle a, b of different length (#614)
1 parent 3e79857 commit 73d58fb

File tree

2 files changed

+15
-13
lines changed

2 files changed

+15
-13
lines changed

src/Filters/filt.jl

Lines changed: 1 addition & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -162,14 +162,7 @@ function filt!(out::AbstractArray{<:Any,N}, f::DF2TFilter{<:PolynomialRatio,Arra
162162
mul!(out, x, b[1])
163163
else
164164
a = coefa(f.coef)
165-
as = length(a)
166-
bs = length(b)
167-
if as != 1
168-
if as < n
169-
append!(a, zero(eltype(a)) for _ in 1:(n-as))
170-
elseif bs < n
171-
append!(b, zero(eltype(b)) for _ in 1:(n-bs))
172-
end
165+
if length(a) != 1
173166
for col in CartesianIndices(axes(x)[2:end])
174167
_filt_iir!(out, b, a, x, view(si, :, col), col)
175168
end

src/dspbase.jl

Lines changed: 14 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -45,9 +45,6 @@ function filt!(out::AbstractArray, b::Union{AbstractVector, Number}, a::Union{Ab
4545
a = @noinline broadcast(/, a, norml)
4646
b = @noinline broadcast(/, b, norml)
4747
end
48-
# Pad the coefficients with zeros if needed
49-
bs<sz && (b = copyto!(zeros(eltype(b), sz), b))
50-
1<as<sz && (a = copyto!(zeros(eltype(a), sz), a))
5148

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

@@ -75,10 +72,22 @@ function _filt_iir!(out, b, a, x, si, col)
7572
xi = x[i, col]
7673
val = muladd(xi, b[1], si[1])
7774
out[i, col] = val
78-
for j=1:(silen-1)
75+
for j in 1:min(length(a), length(b), silen) - 1
7976
si[j] = muladd(val, -a[j+1], muladd(xi, b[j+1], si[j+1]))
8077
end
81-
si[silen] = muladd(xi, b[silen+1], -a[silen+1]*val)
78+
if length(a) == length(b)
79+
si[silen] = muladd(xi, b[silen+1], -a[silen+1]*val)
80+
elseif length(a) > length(b)
81+
for j in length(b):silen-1
82+
si[j] = muladd(val, -a[j+1], si[j+1])
83+
end
84+
si[silen] = -a[silen+1]*val
85+
else
86+
for j in length(a):silen-1
87+
si[j] = muladd(xi, b[j+1], si[j+1])
88+
end
89+
si[silen] = xi*b[silen+1]
90+
end
8291
end
8392
end
8493

0 commit comments

Comments
 (0)