diff --git a/src/Filters/filt.jl b/src/Filters/filt.jl index 4ec76e003..2b3b9523a 100644 --- a/src/Filters/filt.jl +++ b/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 diff --git a/src/deprecated.jl b/src/deprecated.jl index da3573722..3675be8d2 100644 --- a/src/deprecated.jl +++ b/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) +) diff --git a/src/dspbase.jl b/src/dspbase.jl index c2135e83a..afb7c326b 100644 --- a/src/dspbase.jl +++ b/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 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 diff --git a/test/dsp.jl b/test/dsp.jl index 610b44bc4..897023a15 100644 --- a/test/dsp.jl +++ b/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 diff --git a/test/filt.jl b/test/filt.jl index 79473b500..570a4ec42 100644 --- a/test/filt.jl +++ b/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 #######################################