LowpassButterworth

block LowpassButterworth "Output the input signal filtered with a low pass Butterworth filter of any order"
    import Modelica.Blocks.Types.Init;
    import Modelica.Constants.pi;

    extends Modelica.Blocks.Interfaces.SISO;

    parameter Integer n(min = 1) = 2 "Order of filter";
    parameter SI.Frequency f(start = 1) "Cut-off frequency";
    parameter Modelica.Blocks.Types.Init initType = Modelica.Blocks.Types.Init.NoInit "Type of initialization (1: no init, 2: steady state, 3: initial state, 4: initial output)"
        annotation (
            Evaluate = true,
            Dialog(group = "Initialization"));
    parameter Real x1_start[m] = zeros(m) "Initial or guess values of states 1 (der(x1)=x2)"
        annotation (Dialog(group = "Initialization"));
    parameter Real x2_start[m] = zeros(m) "Initial or guess values of states 2"
        annotation (Dialog(group = "Initialization"));
    parameter Real xr_start = 0 "Initial or guess value of real pole for uneven order otherwise dummy"
        annotation (Dialog(group = "Initialization"));
    parameter Real y_start = 0 "Initial value of output (states are initialized in steady state if possible)"
        annotation (Dialog(
            enable = initType == Init.InitialOutput,
            group = "Initialization"));
    output Real x1[m](start = x1_start) "states 1 of second order filters (der(x1) = x2)";
    output Real x2[m](start = x2_start) "states 2 of second order filters";
    output Real xr(start = xr_start) "state of real pole for uneven order otherwise dummy";
protected
    parameter Integer m = integer(0.5 * n);
    parameter Boolean evenOrder = 2 * m == n;
    parameter Real w = 2 * pi * f;
    Real z[m + 1];
    Real polereal[m];
    Real poleimag[m];
    Real realpol;
    Real k2[m];
    Real D[m];
    Real w0[m];
    Real k1;
    Real T;
initial equation
    if initType == Init.SteadyState then 
        der(x1) = zeros(m);
        der(x2) = zeros(m);
        if not evenOrder then 
            der(xr) = 0;
        end if;
    elseif initType == Init.InitialState then 
        x1 = x1_start;
        x2 = x2_start;
        if not evenOrder then 
            xr = xr_start;
        end if;
    elseif initType == Init.InitialOutput then 
        y = y_start;
        der(x1) = zeros(m);
        if evenOrder then 
            if 1 < m then 
                der(x2[1:m - 1]) = zeros(m - 1);
            end if;
        else 
            der(x1) = zeros(m);
        end if;
    end if;
equation
    for i in 1:m loop
        der(x1[i]) = x2[i];
        der(x2[i]) = k2[i] * w0[i] ^ 2 * z[i] - 2 * D[i] * w0[i] * x2[i] - w0[i] ^ 2 * x1[i];
        z[i + 1] = x1[i];
    end for;
    for i in 1:m loop
        polereal[i] = Modelica.Math.cos(0.5 * pi + pi / n * (i - 0.5));
        poleimag[i] = Modelica.Math.sin(0.5 * pi + pi / n * (i - 0.5));
        w0[i] = (polereal[i] ^ 2 + poleimag[i] ^ 2) * w;
        D[i] = -polereal[i] / w0[i] * w;
    end for;
    if evenOrder then 
        xr = 0;
        y = z[m + 1];
    else 
        der(xr) = (k1 * z[m + 1] - xr) / T;
        y = xr;
    end if;
    T = realpol ^ (-1);
    k1 = 1;
    k2 = ones(m);
    realpol = 1 * w;
    z[1] = u;

    annotation (
        Icon(
            coordinateSystem(
                preserveAspectRatio = true,
                extent = {
                    {-100, -100}, 
                    {100, 100}}),
            graphics = {
                Line(
                    points = {
                        {-80, 78}, 
                        {-80, -90}},
                    color = {192, 192, 192}), 
                Polygon(
                    lineColor = {192, 192, 192},
                    fillColor = {192, 192, 192},
                    fillPattern = FillPattern.Solid,
                    points = {
                        {-79.5584, 91.817}, 
                        {-87.5584, 69.817}, 
                        {-71.5584, 69.817}, 
                        {-79.5584, 91.817}}), 
                Line(
                    origin = {-1.939, -1.816},
                    points = {
                        {81.939, 36.056}, 
                        {65.362, 36.056}, 
                        {14.39, -26.199}, 
                        {-29.966, 113.485}, 
                        {-65.374, -61.217}, 
                        {-78.061, -78.184}},
                    color = {0, 0, 127},
                    smooth = Smooth.Bezier), 
                Line(
                    points = {
                        {-90.9779, -80.7697}, 
                        {81.0221, -80.7697}},
                    color = {192, 192, 192}), 
                Polygon(
                    lineColor = {192, 192, 192},
                    fillColor = {192, 192, 192},
                    fillPattern = FillPattern.Solid,
                    points = {
                        {91.3375, -79.8233}, 
                        {69.3375, -71.8233}, 
                        {69.3375, -87.8233}, 
                        {91.3375, -79.8233}}), 
                Text(
                    lineColor = {192, 192, 192},
                    extent = {
                        {-45.1735, -68}, 
                        {92, -11.47}},
                    textString = "LowpassButterworthFilter"), 
                Text(
                    extent = {
                        {8, -146}, 
                        {8, -106}},
                    textString = "f=%f"), 
                Text(
                    lineColor = {192, 192, 192},
                    extent = {
                        {-2, 48}, 
                        {94, 94}},
                    textString = "%n")}),
        Diagram(
            coordinateSystem(
                preserveAspectRatio = true,
                extent = {
                    {-100, -100}, 
                    {100, 100}}),
            graphics = {
                Line(points = {
                    {40, 0}, 
                    {-40, 0}}), 
                Text(
                    extent = {
                        {-55, 55}, 
                        {55, 5}},
                    textString = "1"), 
                Text(
                    extent = {
                        {-55, -5}, 
                        {55, -55}},
                    textString = "a(s)"), 
                Rectangle(
                    extent = {
                        {-60, 60}, 
                        {60, -60}},
                    lineColor = {0, 0, 255}), 
                Line(
                    points = {
                        {-100, 0}, 
                        {-60, 0}},
                    color = {0, 0, 255}), 
                Line(
                    points = {
                        {60, 0}, 
                        {100, 0}},
                    color = {0, 0, 255})}),
        Documentation(info = "<html>\n<p>\nThis block defines the transfer function between the input u\nand the output y as an n-th order low pass filter with <em>Butterworth</em>\ncharacteristics and cut-off frequency f. It is implemented as\na series of second order filters and a first order filter.\nButterworth filters have the feature that the amplitude at the\ncut-off frequency f is 1/sqrt(2) (= 3 dB), i.e., they are\nalways \"normalized\". Step responses of the Butterworth filter of\ndifferent orders are shown in the next figure:\n</p>\n\n<p>\n<img src=\"modelica://Modelica/Resources/Images/Blocks/Butterworth.png\"\n     alt=\"Butterworth.png\">\n</p>\n\n<p>\nIf transients at the simulation start shall be avoided, the filter\nshould be initialized in steady state (e.g., using option\ninitType=Modelica.Blocks.Types.Init.SteadyState).\n</p>\n\n</html>"));
end LowpassButterworth;