Deprecate filt/filt! with si parameter

src/Filters/filt.jl

@@ -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},
@@ -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},
@@ -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)
@@ -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

src/deprecated.jl

@@ -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)
+)

src/dspbase.jl

@@ -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"))
@@ -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
@@ -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
@@ -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
 

test/dsp.jl

@@ -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

test/filt.jl

@@ -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]
@@ -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
 
 #######################################