# How to implement a Bessel filter with ModelingToolkit?

**URL:** <https://discourse.julialang.org/t/how-to-implement-a-bessel-filter-with-modelingtoolkit/101220>\
**Category:** Modelling & Simulations\
**Tags:** question, package, modelingtoolkit\
**Created:** [July 5, 2023, 5:27pm UTC](https://discourse.julialang.org/t/how-to-implement-a-bessel-filter-with-modelingtoolkit/101220 "2023-07-05T17:27:31Z")\
**Posts on this page:** 14\
**Page:** 1

<div class="post-metadata">

**Author:** ![ufechner7](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ufechner7/32/51363_2.png) [@ufechner7](https://discourse.julialang.org/u/ufechner7)\
**Post date:** [July 5, 2023, 5:27pm UTC](https://discourse.julialang.org/t/how-to-implement-a-bessel-filter-with-modelingtoolkit/101220/1 "2023-07-05T17:27:31Z")

</div>

How can I implement the following Bessel filter with ModelingToolkit?

\frac{w\_1²}{0.618 s²~+~1.3617 w\_1 s~+~w\_1²}

I know there is this filter component in the standard library: [Basic Blocks · ModelingToolkitStandardLibrary.jl](https://docs.sciml.ai/ModelingToolkitStandardLibrary/stable/API/blocks/#ModelingToolkitStandardLibrary.Blocks.SecondOrder)

But I don’t think I can use it directly…

There are a lot of blocks in the standard library, but I could not find a generic transfer function block…

UPDATE:  
I found this tutorial how to convert a transfer function into a differential equation: [Single Diff Eq → Transfer Function](https://lpsa.swarthmore.edu/Representations/SysRepTransformations/TF2SDE.html)

I will try this approach…

---

<div class="post-metadata">

**Author:** ![baggepinnen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/baggepinnen/32/693_2.png) [@baggepinnen](https://discourse.julialang.org/u/baggepinnen)\
**Post date:** [July 5, 2023, 7:31pm UTC](https://discourse.julialang.org/t/how-to-implement-a-bessel-filter-with-modelingtoolkit/101220/2 "2023-07-05T19:31:23Z")

</div>

I’m sure you’re aware that it’s straightforward to convert a transfer function to a linear statespace system. A statespace system is a system of linear differential equations. You should be able to take it from there.

---

<div class="post-metadata">

**Author:** ![ufechner7](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ufechner7/32/51363_2.png) [@ufechner7](https://discourse.julialang.org/u/ufechner7)\
**Post date:** [July 5, 2023, 7:32pm UTC](https://discourse.julialang.org/t/how-to-implement-a-bessel-filter-with-modelingtoolkit/101220/3 "2023-07-05T19:32:51Z")

</div>

Well, I will manage, but it is no so very user friendly for people who try to convert Simulink models to ModelingToolkit…

---

<div class="post-metadata">

**Author:** ![baggepinnen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/baggepinnen/32/693_2.png) [@baggepinnen](https://discourse.julialang.org/u/baggepinnen)\
**Post date:** [July 5, 2023, 7:33pm UTC](https://discourse.julialang.org/t/how-to-implement-a-bessel-filter-with-modelingtoolkit/101220/4 "2023-07-05T19:33:48Z")

</div>

> [@ufechner7](#):
>
> But I don’t think I can use it directly

Sure you can, just scale both numerator and denominator so that the denominator is monic (this does not change the input-output relationship) and then identify the parameters with those of the SecondOrder type

---

<div class="post-metadata">

**Author:** ![baggepinnen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/baggepinnen/32/693_2.png) [@baggepinnen](https://discourse.julialang.org/u/baggepinnen)\
**Post date:** [July 5, 2023, 7:40pm UTC](https://discourse.julialang.org/t/how-to-implement-a-bessel-filter-with-modelingtoolkit/101220/5 "2023-07-05T19:40:20Z")

</div>

This answer

> [@Designing a band pass filter with modeling toolkit](https://discourse.julialang.org/t/designing-a-band-pass-filter-with-modeling-toolkit/99118/6):
>
> Before doing analog filter design with DSP, beware of [analogfilter design seems to be broken · Issue #341 · JuliaDSP/DSP.jl · GitHub](https://github.com/JuliaDSP/DSP.jl/issues/341) TLDR; analog filter design is quite broken in DSP but the problem can be worked around. Having said that, here’s how you can convert a filter to an MTK system with the PR [#843](https://github.com/JuliaControl/ControlSystems.jl/pull/843) using DSP, ControlSystemsBase # Digital filter fs = 100 df = digitalfilter(Bandpass(5, 10; fs), Butterworth(2)) G = tf(df, 1/fs) bodeplot(G, xscale=:identity, yscale=:identity, hz=true) …

demonstrates one particular way of obtaining an MTK model for a transfer function.

---

<div class="post-metadata">

**Author:** ![ufechner7](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ufechner7/32/51363_2.png) [@ufechner7](https://discourse.julialang.org/u/ufechner7)\
**Post date:** [July 5, 2023, 7:48pm UTC](https://discourse.julialang.org/t/how-to-implement-a-bessel-filter-with-modelingtoolkit/101220/6 "2023-07-05T19:48:54Z")

</div>

Thanks!

---

<div class="post-metadata">

**Author:** ![ufechner7](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ufechner7/32/51363_2.png) [@ufechner7](https://discourse.julialang.org/u/ufechner7)\
**Post date:** [July 6, 2023, 8:59am UTC](https://discourse.julialang.org/t/how-to-implement-a-bessel-filter-with-modelingtoolkit/101220/7 "2023-07-06T08:59:25Z")

</div>

Using the link in the first post I found the following solution:

Step one: convert the transfer function into a differential equation:  
0.618 ~ \ddot y~+~1.3617 ~w\_1~ \dot y ~+~ w\_1^2~y~=~w\_1^2~u

Step two: implement the differential equation with ModelingToolkit:

```julia
# Bessel filter with step signal as input
function model_bessel(step_time, step_size; w1=3.14, y0=0.0, yd0=0.0)
    println("Defining the model...")

    @variables t u(t)=0 y(t)=y0 yd(t)=yd0
    @parameters w=w1
    
    D = Differential(t)
    eqs = [D(y) ~ yd,
           D(yd) ~ 1/0.618 * (w*w * u - 1.3617w * yd - w*w*y),
           u ~ if_else(t > step_time, step_size, 0.0)]
    
    @named sys = ODESystem(eqs, t)
    println("Model defined...")
    
    sys = structural_simplify(sys)
    println("Model simplified!")
    sys, u, y
end

```

This works fine and has the following step response:

 ![2nd_order_Bessel_filter](https://global.discourse-cdn.com/julialang/original/3X/e/b/eb8ca1065642a52565990fe8d94a7d1f6aa71687.png)

In contrast to a Butterworth filter a Bessel filter has no overshoot. 😀

---

<div class="post-metadata">

**Author:** ![baggepinnen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/baggepinnen/32/693_2.png) [@baggepinnen](https://discourse.julialang.org/u/baggepinnen)\
**Post date:** [July 6, 2023, 9:33am UTC](https://discourse.julialang.org/t/how-to-implement-a-bessel-filter-with-modelingtoolkit/101220/8 "2023-07-06T09:33:47Z")

</div>

Simulating high-order transfer functions this way is prone to numerical difficulties, the method I the post I linked mitigates this by performing numerical balancing of the statespace system.

---

<div class="post-metadata">

**Author:** ![ufechner7](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ufechner7/32/51363_2.png) [@ufechner7](https://discourse.julialang.org/u/ufechner7)\
**Post date:** [July 6, 2023, 10:00am UTC](https://discourse.julialang.org/t/how-to-implement-a-bessel-filter-with-modelingtoolkit/101220/9 "2023-07-06T10:00:54Z")

</div>

Well, second order is not high order, or would you see that differently?

---

<div class="post-metadata">

**Author:** ![baggepinnen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/baggepinnen/32/693_2.png) [@baggepinnen](https://discourse.julialang.org/u/baggepinnen)\
**Post date:** [July 6, 2023, 10:34am UTC](https://discourse.julialang.org/t/how-to-implement-a-bessel-filter-with-modelingtoolkit/101220/10 "2023-07-06T10:34:02Z")

</div>

No, second order is not really all that high. You may still experience numerical difficulties due to poor scaling though, even for this simple system. Consider what happens if you increase `w1`, for `w1=1e6`, the solver exits with

```julia
┌ Warning: dt(8.881784197001252e-16) <= dtmin(8.881784197001252e-16) at t=1.9999999999999998, and step error estimate = 8.930759781475992. Aborting. There is either an error in your model specification or the true solution is unstable.
└ @ SciMLBase ~/.julia/packages/SciMLBase/s9wrq/src/integrator_interface.jl:599
retcode: DtLessThanMin

```

Already for `w1=1000`, you can see part of the problem in the solution:

```julia
model, u, y = model_bessel(2, 100.0, w1=1e3)
prob = ODEProblem(model, [], tspan)
sol = solve(prob, Tsit5())
plot(sol)

```

 ![image](https://global.discourse-cdn.com/julialang/original/3X/3/2/322784f298ff7fddd126dd847fe018bbbe7e48a2.png)

The derivative state becomes enormous, much larger that the position state.

Compare this to a balanced representation:

```julia
using ControlSystemsMTK, ControlSystemsBase

function balanced_model_bessel(step_time, step_size; w1=3.14, y0=0.0, yd0=0.0)
    println("Defining the model...")

    @variables t u(t)=0 y(t)=y0 yd(t)=yd0
    @parameters w=w1

    H = tf(w1^2, [0.618, 1.3617w1, w1^2])
    Hss = ss(H) # This performs balancing
    @named filter = ODESystem(Hss)
    
    eqs = [filter.input.u ~ u
           u ~ ifelse(t > step_time, step_size, 0.0)]
    
    @named sys = ODESystem(eqs, t, systems=[filter])
    println("Model defined...")
    
    sys = structural_simplify(sys)
    println("Model simplified!")
    sys, u
end

model, u = balanced_model_bessel(2, 100.0, w1=1e3)
prob = ODEProblem(model, [], tspan)
sol = solve(prob, Tsit5())
@show length(sol.t)
plot(sol)

```

 ![image](https://global.discourse-cdn.com/julialang/original/3X/9/a/9af11eff13a9457340dc43b5a30e0cad11dcfd74.png)

Here, the two state variables are of comparable magnitude, making the job for the solver much easier.

To plot the correctly scaled filter output in this case, you’d have to plot `plot(sol, idxs=model.filter₊output₊u)`, since the output is scaled w.r.t. the state.

---

<div class="post-metadata">

**Author:** ![ufechner7](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ufechner7/32/51363_2.png) [@ufechner7](https://discourse.julialang.org/u/ufechner7)\
**Post date:** [July 6, 2023, 11:22am UTC](https://discourse.julialang.org/t/how-to-implement-a-bessel-filter-with-modelingtoolkit/101220/11 "2023-07-06T11:22:19Z")

</div>

Thank you very much for explaining the balanced solution, but I am hesitant to add more packages to my project unless I really need them. I am using already way too many packages… And wind turbines don’t operate in the MHz range… 🙂

---

<div class="post-metadata">

**Author:** ![Oscar\_Smith](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/oscar_smith/32/25343_2.png) [@Oscar\_Smith](https://discourse.julialang.org/u/Oscar_Smith)\
**Post date:** [July 6, 2023, 11:52am UTC](https://discourse.julialang.org/t/how-to-implement-a-bessel-filter-with-modelingtoolkit/101220/12 "2023-07-06T11:52:46Z")

</div>

is that a challenge?

---

<div class="post-metadata">

**Author:** ![zdenek\_hurak](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/zdenek_hurak/32/53118_2.png) [@zdenek\_hurak](https://discourse.julialang.org/u/zdenek_hurak)\
**Post date:** [July 7, 2023, 10:10pm UTC](https://discourse.julialang.org/t/how-to-implement-a-bessel-filter-with-modelingtoolkit/101220/13 "2023-07-07T22:10:16Z")

</div>

I am sorry for the digression from the topic of this thread but reading the code

```julia
plot(sol, idxs=model.filter₊output₊u)

```

I have to share my confusion or perhaps even frustration regarding the separation symbols in MTK. So, dots or sub\_pluses? `.` or `₊`? In the above code they are combined and yet the REPL output (as well as the generated plot legend) reads:

```plaintext
julia> model.filter₊output₊u
sys₊filter₊output₊u(t)

```

But clearly the sub\_plus cannot be used everywhere:

```julia
julia> model₊filter₊output₊u
ERROR: UndefVarError: `model₊filter₊output₊u` not defined

```

If I remember well, you were also involved in [the discussion of the choice of the namespace separation char](https://github.com/SciML/ModelingToolkit.jl/issues/907). The conclusion seemed to go back to visualising the separation char as `.`, but now it does not even work when writing the code.

```julia
julia> model.filter.output.u
ERROR: ArgumentError: System sys: variable filter does not exist
Stacktrace:
 [1] getvar(sys::ODESystem, name::Symbol; namespace::Bool)
   @ ModelingToolkit ~/.julia/packages/ModelingToolkit/GFEA2/src/systems/abstractsystem.jl:358
 [2] getvar
   @ ~/.julia/packages/ModelingToolkit/GFEA2/src/systems/abstractsystem.jl:311 [inlined]
 [3] getproperty(sys::ODESystem, name::Symbol; namespace::Bool)
   @ ModelingToolkit ~/.julia/packages/ModelingToolkit/GFEA2/src/systems/abstractsystem.jl:309
 [4] getproperty(sys::ODESystem, name::Symbol)
   @ ModelingToolkit ~/.julia/packages/ModelingToolkit/GFEA2/src/systems/abstractsystem.jl:308
 [5] top-level scope
   @ REPL[10]:1

```

---

<div class="post-metadata">

**Author:** ![baggepinnen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/baggepinnen/32/693_2.png) [@baggepinnen](https://discourse.julialang.org/u/baggepinnen)\
**Post date:** [July 8, 2023, 3:47am UTC](https://discourse.julialang.org/t/how-to-implement-a-bessel-filter-with-modelingtoolkit/101220/14 "2023-07-08T03:47:46Z")

</div>

Yes, I don’t really like this either. In this case, the plus was used since the function written by Uwe returned the simplified system rather than the original system. The intended usage is to access names in the original, unsimplified system, in which case you never need the subplus.

Unfortunately, you still occasionally need to call `complete` on the unsimplified system for this to work properly 😕