filtfilt with second order sections

[gnadt]

Attempting to use second order sections to filter measurements. Getting a mul! error.

sos = SecondOrderSections{Float64,Array{Float64,1}}
biquads –> Vector{Biquad{Float64}} with 21 elements
g –> Vector{Float64} with 21 elements
meas = Vector{Float64} with 15500 elements

filtfilt(sos,meas)

ERROR: DimensionMismatch("result C has dimensions (2, 21), needs (2,1)")
Stacktrace:
[1] _generic_matmatmul!(::Array{Any,2}, ::Char, ::Char, ::Array{Float64,2}, ::Array{Float64,1}, ::LinearAlgebra.MulAddMul{true,true,Bool,Bool}) at /Users/vagrant/worker/juliapro-release-osx1011-0_6/build/tmp_julia/Julia-1.4.app/Contents/Resources/julia/share/julia/stdlib/v1.4/LinearAlgebra/src/matmul.jl:739
[2] generic_matmatmul!(::Array{Any,2}, ::Char, ::Char, ::Array{Float64,2}, ::Array{Float64,1}, ::LinearAlgebra.MulAddMul{true,true,Bool,Bool}) at /Users/vagrant/worker/juliapro-release-osx1011-0_6/build/tmp_julia/Julia-1.4.app/Contents/Resources/julia/share/julia/stdlib/v1.4/LinearAlgebra/src/matmul.jl:727
[3] mul!(::Array{Any,2}, ::Array{Float64,2}, ::Array{Float64,1}, ::Bool, ::Bool) at /Users/vagrant/worker/juliapro-release-osx1011-0_6/build/tmp_julia/Julia-1.4.app/Contents/Resources/julia/share/julia/stdlib/v1.4/LinearAlgebra/src/matmul.jl:235
[4] mul! at /Users/vagrant/worker/juliapro-release-osx1011-0_6/build/tmp_julia/Julia-1.4.app/Contents/Resources/julia/share/julia/stdlib/v1.4/LinearAlgebra/src/matmul.jl:208 [inlined]
[5] filtfilt(::SecondOrderSections{Float64,Array{Float64,1}}, ::Array{Float64,1}) at /Users/gnadt/.juliapro/JuliaPro_v1.4.0-1/packages/DSP/KGj8O/src/Filters/filt.jl:325
[6] top-level scope at none:0

[martinholters]

Who ever gets around to tackling this could also look into why the output array is apparently allocated with eltype Any here.

[wheeheee]

Who ever gets around to tackling this could also look into why the output array is apparently allocated with eltype Any here.

Would this be because of incorrect type promotion? -> promote_type(T,G,S) should be promote_type(T,eltype(G),S)
However, it also looks like the package is missing an implementation of _filt! for vector g in a SOS, like documented in Matlab? The same page also says that g should be one longer than biquads, I think.

[martinholters]

Ah, no, the issue is that we expect g to be a scalar.

julia> filtfilt(SecondOrderSections([Biquad(1.0, 0.0, 0.0, 0.0, 0.0) for _ in 1:21], ones(21)), rand(15500))
ERROR: DimensionMismatch("result C has dimensions (2, 21), needs (2,1)")
[...]

julia> filtfilt(SecondOrderSections([Biquad(1.0, 0.0, 0.0, 0.0, 0.0) for _ in 1:21], 1.0), rand(15500))
15500-element Vector{Float64}:
[...]

We should probably restrict the type in the definition of SecondOrderSections accordingly. Given that our Biquads are not normalized (e.g. b0 == 1) in any way, one might even wonder why we need g at all. The only reason I see is the case of zero biquads, and indeed it's nice that

julia> h = convert(SecondOrderSections, PolynomialRatio([5], [1]))
SecondOrderSections{:z, Float64, Float64}(Biquad{:z, Float64}[], 5.0)

julia> filt(h, ones(3))
3-element Vector{Float64}:
 5.0
 5.0
 5.0

just works without any special-casing.

So from my perspective, the starting point of a fix here would be

--- a/src/Filters/coefficients.jl
+++ b/src/Filters/coefficients.jl
@@ -276,7 +276,7 @@ Filter representation in terms of a cascade of second-order
 sections and gain. `biquads` must be specified as a vector of
 `Biquads`.
 """
-struct SecondOrderSections{Domain,T,G} <: FilterCoefficients{Domain}
+struct SecondOrderSections{Domain,T,G<:Number} <: FilterCoefficients{Domain}
     biquads::Vector{Biquad{Domain,T}}
     g::G
 end

[wheeheee]

Resolved by #522

[martinholters]

#552