Add interface to DSP

docs/src/lib/constructors.md

@@ -20,4 +20,5 @@ zpk
 delay
 pade
 ssdata
+ControlSystemsBase.seriesform
 ```

docs/src/man/creating_systems.md

@@ -433,4 +433,19 @@ uF  │       │ │  │       ├──────►         │ yC │  uP
                              └───────────────────────────────────┘
 ```
 
-See code example [complicated_feedback.jl](https://github.com/JuliaControl/RobustAndOptimalControl.jl/blob/master/examples/complicated_feedback.jl).
+See code example [complicated_feedback.jl](https://github.com/JuliaControl/RobustAndOptimalControl.jl/blob/master/examples/complicated_feedback.jl).
+
+## Filter design
+Filters can be designed using [DSP.jl](https://docs.juliadsp.org/stable/filters/). This results in filter objects with types from the DSP package, which can be converted to transfer functions using [`tf`](@ref) from ControlSystemsBase.
+
+```@example FilterDesign
+using DSP, ControlSystemsBase, Plots
+
+fs = 100
+df = digitalfilter(Bandpass(5, 10; fs), Butterworth(2))
+G = tf(df, 1/fs) # Sample time must be provided in the conversion to get the correct frequency scale in the Bode plot
+bodeplot(G, xscale=:identity, yscale=:identity, hz=true)
+```
+
+See also
+- [`ControlSystemsBase.seriesform`](@ref)

lib/ControlSystemsBase/src/ControlSystemsBase.jl

@@ -129,6 +129,7 @@ import LinearAlgebra: BlasFloat
 export lyap # Make sure LinearAlgebra.lyap is available
 import Printf
 import Printf: @printf, @sprintf
+import DSP
 import DSP: conv
 using ForwardDiff
 import MatrixPencils
@@ -198,6 +199,7 @@ include("nonlinear_components.jl")
 include("types/staticsystems.jl")
 
 include("plotting.jl")
+include("dsp.jl")
 
 @deprecate pole poles
 @deprecate tzero tzeros

lib/ControlSystemsBase/src/dsp.jl

@@ -0,0 +1,39 @@
+
+tf(p::DSP.PolynomialRatio{:z}, h::Real = 1) = tf(DSP.coefb(p), DSP.coefa(p), h)
+tf(p::DSP.PolynomialRatio{:s}) = tf(DSP.coefb(p), DSP.coefa(p))
+tf(p::DSP.ZeroPoleGain{:z}, h::Real = 1) = tf(DSP.PolynomialRatio(p), h)
+tf(p::DSP.ZeroPoleGain{:s}) = tf(DSP.PolynomialRatio(p))
+
+function DSP.PolynomialRatio(G::TransferFunction{<:Discrete})
+    DSP.PolynomialRatio(
+        numvec(G)[],
+        denvec(G)[],
+    )
+end
+
+function TransferFunction(b::DSP.Biquad, h::Real = 1)
+    b0, b1, b2, a1, a2 = b.b0, b.b1, b.b2, b.a1, b.a2
+    tf([b0, b1, b2], [1, a1, a2], h)
+end
+
+
+"""
+    Gs, k = seriesform(G::TransferFunction{Discrete})
+
+Convert a transfer function `G` to a vector of second-order transfer functions and a scalar gain `k`, the product of which equals `G`.
+"""
+function seriesform(G::TransferFunction{<:Discrete})
+    Gs = DSP.SecondOrderSections(DSP.PolynomialRatio(G))
+    bqs = TransferFunction.(Gs.biquads, G.Ts)
+    bqs, Gs.g
+end
+
+
+function DSP.SecondOrderSections(G::TransferFunction{<:Discrete})
+    DSP.SecondOrderSections(DSP.PolynomialRatio(G))
+end
+
+function zpk(p::DSP.ZeroPoleGain, h::Real)
+    z,p,k = p.z, p.p, p.k
+    zpk(z, p, k, h)
+end

lib/ControlSystemsBase/test/runtests.jl

@@ -39,7 +39,8 @@ my_tests = [
             "test_hammerstein_wiener",
             "test_demo_systems",
             "test_autovec",
-            "test_plots"
+            "test_plots",
+            "test_dsp",
             ]
 
 @testset "All Tests" begin

lib/ControlSystemsBase/test/test_dsp.jl

@@ -0,0 +1,30 @@
+@testset "DSP interoperability" begin
+    @info "Testing DSP interoperability"
+    import DSP 
+    G = DemoSystems.resonant()*DemoSystems.resonant(ω0=2) |> tf
+    Gd0 = c2d(DemoSystems.resonant(), 0.1) |> tf
+    Gd = c2d(G, 0.1)
+    Gs, k = ControlSystemsBase.seriesform(Gd)
+    @test k*prod(Gs) ≈ Gd atol=1e-3
+    Gds = DSP.SecondOrderSections(Gd)
+
+    u = randn(100)
+    uf = DSP.filt(Gds, u, zeros(2,2))
+    uls = lsim(Gd, u').y'
+    @test uf[1:end-1] ≈ uls[2:end]
+
+    z,p,k = randn(3), randn(3), randn()
+
+    f = DSP.ZeroPoleGain(z,p,k)
+    fcs = zpk(f, 1)
+
+    @test fcs.matrix[1].z ≈ z
+    @test fcs.matrix[1].p ≈ p
+    @test fcs.matrix[1].k ≈ k
+
+    u = randn(10)
+    uf = DSP.filt(f, u)
+    uls = lsim(fcs, u').y'
+    @test uf ≈ uls
+
+end

src/dsp.jl

@@ -0,0 +1,31 @@
+ControlSystems.tf(p::DSP.PolynomialRatio, h::Real = 1) = tf(DSP.coefb(p), DSP.coefa(p), h)
+ControlSystems.tf(p::DSP.ZeroPoleGain, h::Real = 1) = tf(DSP.PolynomialRatio(p), h)
+
+DSP.PolynomialRatio(G::TransferFunction{<:Discrete}) = DSP.PolynomialRatio(numpoly(G)[], denpoly(G)[])
+
+function ControlSystems.TransferFunction(b::DSP.Biquad, h::Real = 1)
+    b0, b1, b2, a1, a2 = b.b0, b.b1, b.b2, b.a1, b.a2
+    tf([b0, b1, b2], [1, a1, a2], h)
+end
+
+
+"""
+    Gs, k = seriesform(G::TransferFunction)
+
+Convert a transfer function `G` to a vector of second-order transfer functions, the product of which equals `G`. `k` is a scalar gain.
+"""
+function seriesform(G::TransferFunction{<:Discrete})
+    Gs = DSP.SecondOrderSections(DSP.PolynomialRatio(G))
+    bqs = TransferFunction.(Gs.biquads, G.Ts)
+    bqs, Gs.g
+end
+
+
+function DSP.SecondOrderSections(G::TransferFunction{<:Discrete})
+    DSP.SecondOrderSections(DSP.PolynomialRatio(G))
+end
+
+function ControlSystems.zpk(p::DSP.ZeroPoleGain, h::Real)
+    z,p,k = p.z, p.p, p.z
+    zpk(z,p,k,h)
+end

test/test_dsp.jl

@@ -0,0 +1,10 @@
+@testset "DSP interoperability" begin
+    @info "Testing DSP interoperability"
+    import DSP 
+    G = DemoSystems.resonant()*DemoSystems.resonant(ω0=2) |> tf
+    Gd = c2d(G, 0.1)
+    Gs, k = seriesform(Gd)
+    @test k*prod(Gs) ≈ Gd
+    Gds = DSP.SecondOrderSections(Gd)
+
+end