Skip to content

Commit

Permalink
Deprecate filt/filt! with si parameter (#603)
Browse files Browse the repository at this point in the history 
Co-authored-by: wheeheee <104880306+wheeheee@users.noreply.github.com>
  • Loading branch information
@martinholters @wheeheee
martinholters and wheeheee authored Dec 6, 2024
1 parent ade83dd commit 16f1ec0
Show file tree
Hide file tree
Showing 5 changed files with 146 additions and 69 deletions.
52 changes: 17 additions & 35 deletions src/Filters/filt.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -6,37 +6,30 @@
using ..DSP: _filt_fir!, _filt_iir!

## PolynomialRatio
_zerosi(f::PolynomialRatio{:z,T}, ::AbstractArray{S}) where {T,S} =
zeros(promote_type(T, S), max(-firstindex(f.a), -firstindex(f.b)))

"""
filt!(out, f, x[, si])
filt!(out, f, x)
Same as [`filt()`](@ref) but writes the result into the `out`
argument. Output array `out` may not be an alias of `x`, i.e. filtering may
not be done in place.
"""
filt!(out, f::PolynomialRatio{:z}, x::AbstractArray, si=_zerosi(f, x)) =
filt!(out, coefb(f), coefa(f), x, si)
filt!(out, f::PolynomialRatio{:z}, x::AbstractArray) = filt!(out, coefb(f), coefa(f), x)

"""
filt(f::FilterCoefficients{:z}, x::AbstractArray[, si])
filt(f::FilterCoefficients{:z}, x::AbstractArray)
Apply filter or filter coefficients `f` along the first dimension
of array `x`. If `f` is a filter coefficient object, `si`
is an optional array representing the initial filter state (defaults
to zeros). If `f` is a `PolynomialRatio`, `Biquad`, or
of array `x`. If `f` is a `PolynomialRatio`, `Biquad`, or
`SecondOrderSections`, filtering is implemented directly. If
`f` is a `ZeroPoleGain` object, it is first converted to a
`SecondOrderSections` object. If `f` is a Vector, it is
interpreted as an FIR filter, and a naïve or FFT-based algorithm is
selected based on the data and filter length.
"""
filt(f::PolynomialRatio{:z}, x, si=_zerosi(f, x)) = filt(coefb(f), coefa(f), x, si)
filt(f::PolynomialRatio{:z}, x) = filt(coefb(f), coefa(f), x)

## SecondOrderSections
_zerosi(f::SecondOrderSections{:z,T,G}, ::AbstractArray{S}) where {T,G,S} =
zeros(promote_type(T, G, S), 2, length(f.biquads))

# filt! algorithm (no checking, returns si)
function _filt!(out::AbstractArray, si::AbstractArray{S,N}, f::SecondOrderSections{:z},
Expand All @@ -58,29 +51,21 @@ function _filt!(out::AbstractArray, si::AbstractArray{S,N}, f::SecondOrderSectio
si
end

function filt!(out::AbstractArray, f::SecondOrderSections{:z}, x::AbstractArray,
si::AbstractArray{S,N}=_zerosi(f, x)) where {S,N}
biquads = f.biquads

function filt!(out::AbstractArray, f::SecondOrderSections{:z,T,G}, x::AbstractArray{S}) where {T,G,S}
size(x) != size(out) && throw(DimensionMismatch("out size must match x"))
(size(si, 1) != 2 || size(si, 2) != length(biquads) || (N > 2 && size(si)[3:end] != size(x)[2:end])) &&
throw(ArgumentError("si must be 2 x nbiquads or 2 x nbiquads x nsignals"))

initial_si = si
si = similar(si, axes(si)[1:2])
si = Matrix{promote_type(T, G, S)}(undef, 2, length(f.biquads))
for col in CartesianIndices(axes(x)[2:end])
copyto!(si, view(initial_si, :, :, N > 2 ? col : CartesianIndex()))
fill!(si, zero(eltype(si)))
_filt!(out, si, f, x, col)
end
out
end

filt(f::SecondOrderSections{:z,T,G}, x::AbstractArray{S}, si=_zerosi(f, x)) where {T,G,S<:Number} =
filt!(similar(x, promote_type(T, G, S)), f, x, si)
filt(f::SecondOrderSections{:z,T,G}, x::AbstractArray{S}) where {T,G,S<:Number} =
filt!(similar(x, promote_type(T, G, S)), f, x)

## Biquad
_zerosi(::Biquad{:z,T}, ::AbstractArray{S}) where {T,S} =
zeros(promote_type(T, S), 2)

# filt! algorithm (no checking, returns si)
function _filt!(out::AbstractArray, si1::Number, si2::Number, f::Biquad{:z},
Expand All @@ -95,21 +80,17 @@ function _filt!(out::AbstractArray, si1::Number, si2::Number, f::Biquad{:z},
(si1, si2)
end

# filt! variant that preserves si
function filt!(out::AbstractArray, f::Biquad{:z}, x::AbstractArray,
si::AbstractArray{S,N}=_zerosi(f, x)) where {S,N}
function filt!(out::AbstractArray, f::Biquad{:z,T}, x::AbstractArray{S}) where {T,S}
size(x) != size(out) && throw(DimensionMismatch("out size must match x"))
(size(si, 1) != 2 || (N > 1 && size(si)[2:end] != size(x)[2:end])) &&
throw(ArgumentError("si must have two rows and 1 or nsignals columns"))

for col in CartesianIndices(axes(x)[2:end])
_filt!(out, si[1, N > 1 ? col : CartesianIndex()], si[2, N > 1 ? col : CartesianIndex()], f, x, col)
_filt!(out, zero(promote_type(T, S)), zero(promote_type(T, S)), f, x, col)
end
out
end

filt(f::Biquad{:z,T}, x::AbstractArray{S}, si=_zerosi(f, x)) where {T,S<:Number} =
filt!(similar(x, promote_type(T, S)), f, x, si)
filt(f::Biquad{:z,T}, x::AbstractArray{S}) where {T,S<:Number} =
filt!(similar(x, promote_type(T, S)), f, x)

## For arbitrary filters, convert to SecondOrderSections
filt(f::FilterCoefficients{:z}, x) = filt(convert(SecondOrderSections, f), x)
Expand Down Expand Up @@ -375,8 +356,9 @@ function filtfilt(f::SecondOrderSections{:z,T,G}, x::AbstractArray{S}) where {T,
istart = 1
for i = 1:Base.trailingsize(x, 2)
extrapolate_signal!(extrapolated, 1, x, istart, size(x, 1), pad_length)
reverse!(filt!(extrapolated, f, extrapolated, mul!(zitmp, zi, extrapolated[1])))
filt!(extrapolated, f, extrapolated, mul!(zitmp, zi, extrapolated[1]))
_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]
end
Expand Down
97 changes: 97 additions & 0 deletions src/deprecated.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -3,3 +3,100 @@ import .Util.nextfastfft
@deprecate nextfastfft(ns...) nextfastfft.(ns) false

@deprecate (conv(u::AbstractVector{T}, v::AbstractVector{T}, A::AbstractMatrix{T}) where T) conv(u, transpose(v), A)

@deprecate(
filt!(out, f::PolynomialRatio{:z}, x::AbstractVector, si::AbstractVector),
filt!(out, DF2TFilter(f, copy(si)), x)
)
@deprecate(
filt!(out, f::PolynomialRatio{:z}, x::AbstractArray{<:Any, N}, si::AbstractArray{<:Any, N}) where {N},
filt!(out, DF2TFilter(f, copy(si)), x)
)
@deprecate(
filt!(out, f::PolynomialRatio{:z}, x::AbstractArray, si::AbstractVector),
filt!(out, DF2TFilter(f, repeat(si; outer=(1, size(x)[2:end]...))), x)
)
@deprecate(
filt(f::PolynomialRatio{:z}, x::AbstractVector, si::AbstractVector),
filt(DF2TFilter(f, copy(si)), x)
)
@deprecate(
filt(f::PolynomialRatio{:z}, x::AbstractArray{<:Any, N}, si::AbstractArray{<:Any, N}) where {N},
filt(DF2TFilter(f, copy(si)), x)
)
@deprecate(
filt(f::PolynomialRatio{:z}, x::AbstractArray, si::AbstractVector),
filt(DF2TFilter(f, repeat(si; outer=(1, size(x)[2:end]...))), x)
)
@deprecate(
filt!(out, f::SecondOrderSections{:z}, x::AbstractVector, si::AbstractVecOrMat),
filt!(out, DF2TFilter(f, copy(si)), x)
)
@deprecate(
filt!(out, f::SecondOrderSections{:z}, x::AbstractArray, si::AbstractArray),
filt!(out, DF2TFilter(f, copy(si)), x)
)
@deprecate(
filt!(out, f::SecondOrderSections{:z}, x::AbstractArray, si::AbstractVecOrMat),
filt!(out, DF2TFilter(f, repeat(si; outer=(1, 1, size(x)[2:end]...))), x)
)
@deprecate(
filt(f::SecondOrderSections{:z}, x::AbstractVector, si::AbstractVecOrMat),
filt(DF2TFilter(f, copy(si)), x)
)
@deprecate(
filt(f::SecondOrderSections{:z}, x::AbstractArray, si::AbstractArray),
filt(DF2TFilter(f, copy(si)), x)
)
@deprecate(
filt(f::SecondOrderSections{:z}, x::AbstractArray, si::AbstractVecOrMat),
filt(DF2TFilter(f, repeat(si; outer=(1, 1, size(x)[2:end]...))), x)
)
@deprecate(
filt!(out, f::Biquad{:z}, x::AbstractVector, si::AbstractVector),
filt!(out, DF2TFilter(f, copy(si)), x)
)
@deprecate(
filt!(out, f::Biquad{:z}, x::AbstractArray{<:Any, N}, si::AbstractArray{<:Any, N}) where {N},
filt!(out, DF2TFilter(f, copy(si)), x)
)
@deprecate(
filt!(out, f::Biquad{:z}, x::AbstractArray, si::AbstractVector),
filt!(out, DF2TFilter(f, repeat(si; outer=(1, size(x)[2:end]...))), x)
)
@deprecate(
filt(f::Biquad{:z}, x::AbstractVector, si::AbstractVector),
filt(DF2TFilter(f, copy(si)), x)
)
@deprecate(
filt(f::Biquad{:z}, x::AbstractArray{<:Any, N}, si::AbstractArray{<:Any, N}) where {N},
filt(DF2TFilter(f, copy(si)), x)
)
@deprecate(
filt(f::Biquad{:z}, x::AbstractArray, si::AbstractVector),
filt(DF2TFilter(f, repeat(si; outer=(1, size(x)[2:end]...))), x)
)
@deprecate(
filt!(out, b::Union{AbstractVector, Number}, a::Union{AbstractVector, Number}, x::AbstractVector, si::AbstractVector),
filt!(out, DF2TFilter(PolynomialRatio(b, a), copy(si)), x)
)
@deprecate(
filt!(out, b::Union{AbstractVector, Number}, a::Union{AbstractVector, Number}, x::AbstractArray{<:Any, N}, si::AbstractArray{<:Any, N}) where {N},
filt!(out, DF2TFilter(PolynomialRatio(b, a), copy(si)), x)
)
@deprecate(
filt!(out, b::Union{AbstractVector, Number}, a::Union{AbstractVector, Number}, x::AbstractArray, si::AbstractVector),
filt!(out, DF2TFilter(PolynomialRatio(b, a), repeat(si; outer=(1, size(x)[2:end]...))), x)
)
@deprecate(
filt(b::Union{AbstractVector, Number}, a::Union{AbstractVector, Number}, x::AbstractVector, si::AbstractVector),
filt(DF2TFilter(PolynomialRatio(b, a), copy(si)), x)
)
@deprecate(
filt(b::Union{AbstractVector, Number}, a::Union{AbstractVector, Number}, x::AbstractArray{<:Any, N}, si::AbstractArray{<:Any, N}) where {N},
filt(DF2TFilter(PolynomialRatio(b, a), copy(si)), x)
)
@deprecate(
filt(b::Union{AbstractVector, Number}, a::Union{AbstractVector, Number}, x::AbstractArray, si::AbstractVector),
filt(DF2TFilter(PolynomialRatio(b, a), repeat(si; outer=(1, size(x)[2:end]...))), x)
)
37 changes: 12 additions & 25 deletions src/dspbase.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2,35 +2,29 @@

const SMALL_FILT_CUTOFF = 66

_zerosi(b,a,T) = zeros(promote_type(eltype(b), eltype(a), T), max(length(a), length(b))-1)

"""
filt(b::Union{AbstractVector,Number},
a::Union{AbstractVector,Number},
x::AbstractArray,
[si::AbstractArray])
x::AbstractArray)
Apply filter described by vectors `a` and `b` to vector `x`, with an optional initial filter
state vector `si` (defaults to zeros).
Apply filter described by vectors `a` and `b` to vector `x`.
Inputs that are `Number`s are treated as one-element `Vector`s.
"""
function filt(b::Union{AbstractVector, Number}, a::Union{AbstractVector, Number},
x::AbstractArray{T}, si::AbstractArray{S} = _zerosi(b,a,T)) where {T,S}
filt!(similar(x, promote_type(eltype(b), eltype(a), T, S)), b, a, x, si)
end
filt(b::Union{AbstractVector, Number}, a::Union{AbstractVector, Number}, x::AbstractArray{T}) where {T} =
filt!(similar(x, promote_type(eltype(b), eltype(a), T)), b, a, x)

# in-place filtering: returns results in the out argument, which may shadow x
# (and does so by default)

"""
filt!(out, b, a, x, [si])
filt!(out, b, a, x)
Same as [`filt`](@ref) but writes the result into the `out` argument, which may
alias the input `x` to modify it in-place.
"""
function filt!(out::AbstractArray, b::Union{AbstractVector, Number}, a::Union{AbstractVector, Number},
x::AbstractArray{T}, si::AbstractArray{S,N} = _zerosi(b,a,T)) where {T,S,N}
x::AbstractArray{T}) where {T}
isempty(b) && throw(ArgumentError("filter vector b must be non-empty"))
isempty(a) && throw(ArgumentError("filter vector a must be non-empty"))
a[1] == 0 && throw(ArgumentError("filter vector a[1] must be nonzero"))
Expand All @@ -41,13 +35,6 @@ function filt!(out::AbstractArray, b::Union{AbstractVector, Number}, a::Union{Ab
as = length(a)
bs = length(b)
sz = max(as, bs)
silen = sz - 1

if size(si, 1) != silen
throw(ArgumentError("initial state vector si must have max(length(a),length(b))-1 rows"))
elseif N > 1 && size(si)[2:end] != size(x)[2:end]
throw(ArgumentError("initial state si must be a vector or have the same number of columns as x"))
end

iszero(size(x, 1)) && return out
isone(sz) && return (k = b[1] / a[1]; @noinline mul!(out, x, k)) # Simple scaling without memory
Expand All @@ -62,14 +49,15 @@ function filt!(out::AbstractArray, b::Union{AbstractVector, Number}, a::Union{Ab
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)

if as == 1 && bs <= SMALL_FILT_CUTOFF
fill!(si, zero(eltype(si)))
_small_filt_fir!(out, b, x, si, Val(bs))
else
initial_si = si
si = similar(si, axes(si, 1))
for col in CartesianIndices(axes(x)[2:end])
# Reset the filter state
copyto!(si, view(initial_si, :, N > 1 ? col : CartesianIndex()))
fill!(si, zero(eltype(si)))
if as > 1
_filt_iir!(out, b, a, x, si, col)
else
Expand Down Expand Up @@ -146,14 +134,13 @@ end
# Convert array filter tap input to tuple for small-filtering
function _small_filt_fir!(
out::AbstractArray, h::AbstractVector, x::AbstractArray,
si::AbstractArray{S,N}, ::Val{bs}) where {S,N,bs}
si::AbstractVector, ::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))
for col in CartesianIndices(axes(x)[2:end])
v_si = N > 1 ? view(si, :, col) : si
_filt_fir!(out, b, x, v_si, col)
_filt_fir!(out, b, x, si, col)
end
end

Expand Down
6 changes: 3 additions & 3 deletions test/dsp.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -20,16 +20,16 @@ using DSP: filt, filt!, deconv, conv, xcorr,
@test filt(b, 1., [x 1.0:8.0]) == [filt(b, 1., x) filt(b, 1., 1.0:8.0)]
@test filt(b, [1., -0.5], [x 1.0:8.0]) == [filt(b, [1., -0.5], x) filt(b, [1., -0.5], 1.0:8.0)]
si = zeros(3)
@test filt(b, 1., [x 1.0:8.0], si) == [filt(b, 1., x, si) filt(b, 1., 1.0:8.0, si)]
@test @test_deprecated(filt(b, 1., [x 1.0:8.0], si) == [filt(b, 1., x, si) filt(b, 1., 1.0:8.0, si)])
@test si == zeros(3) # Will likely fail if/when arrayviews are implemented
si = [zeros(3) ones(3)]
@test filt(b, 1., [x 1.0:8.0], si) == [filt(b, 1., x, zeros(3)) filt(b, 1., 1.0:8.0, ones(3))]
@test @test_deprecated(filt(b, 1., [x 1.0:8.0], si) == [filt(b, 1., x, zeros(3)) filt(b, 1., 1.0:8.0, ones(3))])
# With initial conditions: a lowpass 5-pole butterworth filter with W_n = 0.25,
# and a stable initial filter condition matched to the initial value
b = [0.003279216306360201,0.016396081531801006,0.03279216306360201,0.03279216306360201,0.016396081531801006,0.003279216306360201]
a = [1.0,-2.4744161749781606,2.8110063119115782,-1.703772240915465,0.5444326948885326,-0.07231566910295834]
si = [0.9967207836936347,-1.4940914728163142,1.2841226760316475,-0.4524417279474106,0.07559488540931815]
@test filt(b, a, ones(10), si) ≈ ones(10) # Shouldn't affect DC offset
@test @test_deprecated(filt(b, a, ones(10), si)) ≈ ones(10) # Shouldn't affect DC offset

@test_throws ArgumentError filt!([1, 2], [1], [1], [1])
end
Expand Down
23 changes: 17 additions & 6 deletions test/filt.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -80,10 +80,10 @@ end
@test all(col -> col ≈ y_ref, eachslice(filt(Biquad(PolynomialRatio(b, a)), x); dims=slicedims))
@test all(col -> col ≈ y_ref, eachslice(filt(SecondOrderSections(PolynomialRatio(b, a)), x); dims=slicedims))
# with si given
@test all(col -> col ≈ y_ref, eachslice(filt(b, a, x, zeros(1, sz[2:end]...)); dims=slicedims))
@test all(col -> col ≈ y_ref, eachslice(filt(PolynomialRatio(b, a), x, zeros(1, sz[2:end]...)); dims=slicedims))
@test all(col -> col ≈ y_ref, eachslice(filt(Biquad(PolynomialRatio(b, a)), x, zeros(2, sz[2:end]...)); dims=slicedims))
@test all(col -> col ≈ y_ref, eachslice(filt(SecondOrderSections(PolynomialRatio(b, a)), x, zeros(2, 1, sz[2:end]...)); dims=slicedims))
@test all(col -> col ≈ y_ref, eachslice(@test_deprecated(filt(b, a, x, zeros(1, sz[2:end]...))); dims=slicedims))
@test all(col -> col ≈ y_ref, eachslice(@test_deprecated(filt(PolynomialRatio(b, a), x, zeros(1, sz[2:end]...))); dims=slicedims))
@test all(col -> col ≈ y_ref, eachslice(@test_deprecated(filt(Biquad(PolynomialRatio(b, a)), x, zeros(2, sz[2:end]...))); dims=slicedims))
@test all(col -> col ≈ y_ref, eachslice(@test_deprecated(filt(SecondOrderSections(PolynomialRatio(b, a)), x, zeros(2, 1, sz[2:end]...))); dims=slicedims))
# use _small_filt_fir!
b = [0.1, 0.1]
a = [1.0]
Expand Down Expand Up @@ -186,10 +186,21 @@ end
a = [0.9, 0.6]
b = [0.4, 1]
z = [0.4750]
x = readdlm(joinpath(dirname(@__FILE__), "data", "spectrogram_x.txt"),'\t')
filt!(vec(x), b, a, vec(x), z)
x = vec(readdlm(joinpath(dirname(@__FILE__), "data", "spectrogram_x.txt"),'\t'))
@test_deprecated(filt!(x, b, a, x, z))

@test matlab_filt ≈ x

x = vec(readdlm(joinpath(dirname(@__FILE__), "data", "spectrogram_x.txt"),'\t'))
filt!(x, DF2TFilter(PolynomialRatio(b, a), z), x)

@test matlab_filt ≈ x

# With initial conditions: a lowpass 5-pole butterworth filter with W_n = 0.25,
# and a stable initial filter condition matched to the initial value
zpg = digitalfilter(Lowpass(0.25), Butterworth(5))
si = [0.9967207836936347, -1.4940914728163142, 1.2841226760316475, -0.4524417279474106, 0.07559488540931815]
@test filt(DF2TFilter(PolynomialRatio(zpg), si), ones(10)) ≈ ones(10) # Shouldn't affect DC offset
end

#######################################
Expand Down

2 comments on commit 16f1ec0

@martinholters
Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@JuliaRegistrator register

@JuliaRegistrator
Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Registration pull request created: JuliaRegistries/General/120797

Tip: Release Notes

Did you know you can add release notes too? Just add markdown formatted text underneath the comment after the text
"Release notes:" and it will be added to the registry PR, and if TagBot is installed it will also be added to the
release that TagBot creates. i.e.

@JuliaRegistrator register

Release notes:

## Breaking changes

- blah

To add them here just re-invoke and the PR will be updated.

Tagging

After the above pull request is merged, it is recommended that a tag is created on this repository for the registered package version.

This will be done automatically if the Julia TagBot GitHub Action is installed, or can be done manually through the github interface, or via:

git tag -a v0.8.0 -m "<description of version>" 16f1ec071b1044700e8d110cefb9387c6237a146
git push origin v0.8.0

Please sign in to comment.