Filter

block Filter "Continuous low pass, high pass, band pass or band stop IIR-filter of type CriticalDamping, Bessel, Butterworth or ChebyshevI"
    import Modelica.Blocks.Continuous.Internal;

    extends Modelica.Blocks.Interfaces.SISO;

    parameter Modelica.Blocks.Types.AnalogFilter analogFilter = Modelica.Blocks.Types.AnalogFilter.CriticalDamping "Analog filter characteristics (CriticalDamping/Bessel/Butterworth/ChebyshevI)";
    parameter Modelica.Blocks.Types.FilterType filterType = Modelica.Blocks.Types.FilterType.LowPass "Type of filter (LowPass/HighPass/BandPass/BandStop)";
    parameter Integer order(min = 1) = 2 "Order of filter";
    parameter Modelica.SIunits.Frequency f_cut "Cut-off frequency";
    parameter Real gain = 1 "Gain (= amplitude of frequency response at zero frequency)";
    parameter Real A_ripple(unit = "dB") = 0.5 "Pass band ripple for Chebyshev filter (otherwise not used); > 0 required"
        annotation (Dialog(enable = analogFilter == Modelica.Blocks.Types.AnalogFilter.ChebyshevI));
    parameter Modelica.SIunits.Frequency f_min = 0 "Band of band pass/stop filter is f_min (A=-3db*gain) .. f_cut (A=-3db*gain)"
        annotation (Dialog(enable = filterType == Modelica.Blocks.Types.FilterType.BandPass or filterType == Modelica.Blocks.Types.FilterType.BandStop));
    parameter Boolean normalized = true "= true, if amplitude at f_cut = -3db, otherwise unmodified filter";
    parameter Modelica.Blocks.Types.Init init = Modelica.Blocks.Types.Init.SteadyState "Type of initialization (no init/steady state/initial state/initial output)"
        annotation (
            Evaluate = true,
            Dialog(tab = "Advanced"));
    final parameter Integer nx = if filterType == Modelica.Blocks.Types.FilterType.LowPass or filterType == Modelica.Blocks.Types.FilterType.HighPass then order else 2 * order;
    parameter Real x_start[nx] = zeros(nx) "Initial or guess values of states"
        annotation (Dialog(tab = "Advanced"));
    parameter Real y_start = 0 "Initial value of output"
        annotation (Dialog(tab = "Advanced"));
    parameter Real u_nominal = 1 "Nominal value of input (used for scaling the states)"
        annotation (Dialog(tab = "Advanced"));
    Modelica.Blocks.Interfaces.RealOutput x[nx] "Filter states";
protected
    parameter Integer ncr = if analogFilter == Modelica.Blocks.Types.AnalogFilter.CriticalDamping then order else mod(order, 2);
    parameter Integer nc0 = if analogFilter == Modelica.Blocks.Types.AnalogFilter.CriticalDamping then 0 else integer(0.5 * order);
    parameter Integer na = if filterType == Modelica.Blocks.Types.FilterType.BandPass or filterType == Modelica.Blocks.Types.FilterType.BandStop then order else if analogFilter == Modelica.Blocks.Types.AnalogFilter.CriticalDamping then 0 else integer(0.5 * order);
    parameter Integer nr = if filterType == Modelica.Blocks.Types.FilterType.BandPass or filterType == Modelica.Blocks.Types.FilterType.BandStop then 0 else if analogFilter == Modelica.Blocks.Types.AnalogFilter.CriticalDamping then order else mod(order, 2);
    parameter Real cr[ncr](each fixed = false);
    parameter Real c0[nc0](each fixed = false);
    parameter Real c1[nc0](each fixed = false);
    parameter Real r[nr](each fixed = false);
    parameter Real a[na](each fixed = false);
    parameter Real b[na](each fixed = false);
    parameter Real ku[na](each fixed = false);
    parameter Real k1[if filterType == Modelica.Blocks.Types.FilterType.LowPass then 0 else na](each fixed = false);
    parameter Real k2[if filterType == Modelica.Blocks.Types.FilterType.LowPass then 0 else na](each fixed = false);
    Real uu[na + nr + 1];
initial equation
    if init == Modelica.Blocks.Types.Init.InitialState then 
        x = x_start;
    elseif init == Modelica.Blocks.Types.Init.SteadyState then 
        der(x) = zeros(nx);
    elseif init == Modelica.Blocks.Types.Init.InitialOutput then 
        y = y_start;
        if 1 < nx then 
            der(x[1:nx - 1]) = zeros(nx - 1);
        end if;
    end if;
    if analogFilter == Modelica.Blocks.Types.AnalogFilter.CriticalDamping then 
        cr = Internal.Filter.base.CriticalDamping(order, normalized);
    elseif analogFilter == Modelica.Blocks.Types.AnalogFilter.Bessel then 
        (cr, c0, c1) = Internal.Filter.base.Bessel(order, normalized);
    elseif analogFilter == Modelica.Blocks.Types.AnalogFilter.Butterworth then 
        (cr, c0, c1) = Internal.Filter.base.Butterworth(order, normalized);
    elseif analogFilter == Modelica.Blocks.Types.AnalogFilter.ChebyshevI then 
        (cr, c0, c1) = Internal.Filter.base.ChebyshevI(order, A_ripple, normalized);
    end if;
    if filterType == Modelica.Blocks.Types.FilterType.LowPass then 
        (r, a, b, ku) = Internal.Filter.roots.lowPass(cr, c0, c1, f_cut);
    elseif filterType == Modelica.Blocks.Types.FilterType.HighPass then 
        (r, a, b, ku, k1, k2) = Internal.Filter.roots.highPass(cr, c0, c1, f_cut);
    elseif filterType == Modelica.Blocks.Types.FilterType.BandPass then 
        (a, b, ku, k1, k2) = Internal.Filter.roots.bandPass(cr, c0, c1, f_min, f_cut);
    elseif filterType == Modelica.Blocks.Types.FilterType.BandStop then 
        (a, b, ku, k1, k2) = Internal.Filter.roots.bandStop(cr, c0, c1, f_min, f_cut);
    end if;
equation
    for i in 1:na loop
        der(x[nr + 2 * i - 1]) = a[i] * x[nr + 2 * i - 1] - b[i] * x[nr + 2 * i] + ku[i] * uu[nr + i];
        der(x[nr + 2 * i]) = b[i] * x[nr + 2 * i - 1] + a[i] * x[nr + 2 * i];
    end for;
    for i in 1:nr loop
        der(x[i]) = r[i] * (x[i] - uu[i]);
    end for;
    if filterType == Modelica.Blocks.Types.FilterType.LowPass then 
        for i in 1:nr loop
            uu[i + 1] = x[i];
        end for;
        for i in 1:na loop
            uu[nr + i + 1] = x[nr + 2 * i];
        end for;
    elseif filterType == Modelica.Blocks.Types.FilterType.HighPass then 
        for i in 1:nr loop
            uu[i + 1] = -x[i] + uu[i];
        end for;
        for i in 1:na loop
            uu[nr + i + 1] = k1[i] * x[nr + 2 * i - 1] + k2[i] * x[nr + 2 * i] + uu[nr + i];
        end for;
    elseif filterType == Modelica.Blocks.Types.FilterType.BandPass then 
        for i in 1:na loop
            uu[nr + i + 1] = k1[i] * x[nr + 2 * i - 1] + k2[i] * x[nr + 2 * i];
        end for;
    elseif filterType == Modelica.Blocks.Types.FilterType.BandStop then 
        for i in 1:na loop
            uu[nr + i + 1] = k1[i] * x[nr + 2 * i - 1] + k2[i] * x[nr + 2 * i] + uu[nr + i];
        end for;
    else 
        assert(false, "filterType (= " + String(filterType) + ") is unknown");
        uu = zeros(na + nr + 1);
    end if;
    assert(filterType == Modelica.Blocks.Types.FilterType.LowPass or filterType == Modelica.Blocks.Types.FilterType.HighPass or 0 < f_min, "f_min > 0 required for band pass and band stop filter");
    assert(0 < A_ripple, "A_ripple > 0 required");
    assert(0 < f_cut, "f_cut > 0 required");
    assert(0 < u_nominal, "u_nominal > 0 required");
    y = gain * u_nominal * uu[nr + na + 1];
    uu[1] = u / u_nominal;

    annotation (
        Icon(
            coordinateSystem(
                preserveAspectRatio = true,
                extent = {
                    {-100, -100}, 
                    {100, 100}}),
            graphics = {
                Line(
                    points = {
                        {-80, 80}, 
                        {-80, -88}},
                    color = {192, 192, 192}), 
                Polygon(
                    lineColor = {192, 192, 192},
                    fillColor = {192, 192, 192},
                    fillPattern = FillPattern.Solid,
                    points = {
                        {-80, 92}, 
                        {-88, 70}, 
                        {-72, 70}, 
                        {-80, 92}}), 
                Line(
                    points = {
                        {-90, -78}, 
                        {82, -78}},
                    color = {192, 192, 192}), 
                Polygon(
                    lineColor = {192, 192, 192},
                    fillColor = {192, 192, 192},
                    fillPattern = FillPattern.Solid,
                    points = {
                        {90, -78}, 
                        {68, -70}, 
                        {68, -86}, 
                        {90, -78}}), 
                Text(
                    lineColor = {192, 192, 192},
                    extent = {
                        {-66, 52}, 
                        {88, 90}},
                    textString = "%order"), 
                Text(
                    extent = {
                        {-138, -140}, 
                        {162, -110}},
                    textString = "f_cut=%f_cut"), 
                Rectangle(
                    lineColor = {160, 160, 164},
                    fillColor = {255, 255, 255},
                    fillPattern = FillPattern.Backward,
                    extent = {
                        {-80, -78}, 
                        {22, 10}}), 
                Line(
                    origin = {3.333, -6.667},
                    points = {
                        {-83.333, 34.667}, 
                        {24.667, 34.667}, 
                        {42.667, -71.333}},
                    color = {0, 0, 127},
                    smooth = Smooth.Bezier)}),
        Documentation(
            info = "<html>\n\n<p>\nThis blocks models various types of filters:\n</p>\n\n<blockquote>\n<strong>low pass, high pass, band pass, and band stop filters</strong>\n</blockquote>\n\n<p>\nusing various filter characteristics:\n</p>\n\n<blockquote>\n<strong>CriticalDamping, Bessel, Butterworth, Chebyshev Type I filters</strong>\n</blockquote>\n\n<p>\nBy default, a filter block is initialized in <strong>steady-state</strong>, in order to\navoid unwanted oscillations at the beginning. In special cases, it might be\nuseful to select one of the other initialization options under tab\n\"Advanced\".\n</p>\n\n<p>\nTypical frequency responses for the 4 supported low pass filter types\nare shown in the next figure:\n</p>\n\n<blockquote>\n<img src=\"modelica://Modelica/Resources/Images/Blocks/LowPassOrder4Filters.png\"\n     alt=\"LowPassOrder4Filters.png\">\n</blockquote>\n\n<p>\nThe step responses of the same low pass filters are shown in the next figure,\nstarting from a steady state initial filter with initial input = 0.2:\n</p>\n\n<blockquote>\n<img src=\"modelica://Modelica/Resources/Images/Blocks/LowPassOrder4FiltersStepResponse.png\"\n     alt=\"LowPassOrder4FiltersStepResponse.png\">\n</blockquote>\n\n<p>\nObviously, the frequency responses give a somewhat wrong impression\nof the filter characteristics: Although Butterworth and Chebyshev\nfilters have a significantly steeper magnitude as the\nCriticalDamping and Bessel filters, the step responses of\nthe latter ones are much better. This means for example, that\na CriticalDamping or a Bessel filter should be selected,\nif a filter is mainly used to make a non-linear inverse model\nrealizable.\n</p>\n\n<p>\nTypical frequency responses for the 4 supported high pass filter types\nare shown in the next figure:\n</p>\n\n<blockquote>\n<img src=\"modelica://Modelica/Resources/Images/Blocks/HighPassOrder4Filters.png\"\n     alt=\"HighPassOrder4Filters.png\">\n</blockquote>\n\n<p>\nThe corresponding step responses of these high pass filters are\nshown in the next figure:\n</p>\n<blockquote>\n<img src=\"modelica://Modelica/Resources/Images/Blocks/HighPassOrder4FiltersStepResponse.png\"\n     alt=\"HighPassOrder4FiltersStepResponse.png\">\n</blockquote>\n\n<p>\nAll filters are available in <strong>normalized</strong> (default) and non-normalized form.\nIn the normalized form, the amplitude of the filter transfer function\nat the cut-off frequency f_cut is -3 dB (= 10^(-3/20) = 0.70794..).\nNote, when comparing the filters of this function with other software systems,\nthe setting of \"normalized\" has to be selected appropriately. For example, the signal processing\ntoolbox of MATLAB provides the filters in non-normalized form and\ntherefore a comparison makes only sense, if normalized = <strong>false</strong>\nis set. A normalized filter is usually better suited for applications,\nsince filters of different orders are \"comparable\",\nwhereas non-normalized filters usually require to adapt the\ncut-off frequency, when the order of the filter is changed.\nSee a comparison of \"normalized\" and \"non-normalized\" filters at hand of\nCriticalDamping filters of order 1,2,3:\n</p>\n\n<blockquote>\n<img src=\"modelica://Modelica/Resources/Images/Blocks/CriticalDampingNormalized.png\"\n     alt=\"CriticalDampingNormalized.png\">\n</blockquote>\n\n<blockquote>\n<img src=\"modelica://Modelica/Resources/Images/Blocks/CriticalDampingNonNormalized.png\"\n     alt=\"CriticalDampingNonNormalized.png\">\n</blockquote>\n\n<h4>Implementation</h4>\n\n<p>\nThe filters are implemented in the following, reliable way:\n</p>\n\n<ol>\n<li> A prototype low pass filter with a cut-off angular frequency of 1 rad/s is constructed\n     from the desired analogFilter and the desired normalization.</li>\n\n<li> This prototype low pass filter is transformed to the desired filterType and the\n     desired cut-off frequency f_cut using a transformation on the Laplace variable \"s\".</li>\n\n<li> The resulting first and second order transfer functions are implemented in\n     state space form, using the \"eigen value\" representation of a transfer function:\n     <pre>\n\n  // second order block with eigen values: a +/- jb\n  <strong>der</strong>(x1) = a*x1 - b*x2 + (a^2 + b^2)/b*u;\n  <strong>der</strong>(x2) = b*x1 + a*x2;\n       y  = x2;\n     </pre>\n     The dc-gain from the input to the output of this block is one and the selected\n     states are in the order of the input (if \"u\" is in the order of \"one\", then the\n     states are also in the order of \"one\"). In the \"Advanced\" tab, a \"nominal\" value for\n     the input \"u\" can be given. If appropriately selected, the states are in the order of \"one\" and\n     then step-size control is always appropriate.</li>\n</ol>\n\n<h4>References</h4>\n\n<dl>\n<dt>Tietze U., and Schenk C. (2002):</dt>\n<dd> <strong>Halbleiter-Schaltungstechnik</strong>.\n     Springer Verlag, 12. Auflage, pp. 815-852.</dd>\n</dl>\n\n</html>",
            revisions = "<html>\n<dl>\n  <dt><strong>Main Author:</strong></dt>\n  <dd><a href=\"http://www.robotic.dlr.de/Martin.Otter/\">Martin Otter</a>,\n      DLR Oberpfaffenhofen.</dd>\n</dl>\n\n<h4>Acknowledgement</h4>\n\n<p>\nThe development of this block was partially funded by BMBF within the\n     <a href=\"http://www.eurosyslib.com/\">ITEA2 EUROSYSLIB</a>\n      project.\n</p>\n\n</html>"));
end Filter;