DiscreteTransferFunction

A discrete transfer function

\[H(z) = \dfrac{B(z)}{A(z)} = \dfrac{b_1 z^{n_b-1} + b_2 z^{n_b-2} + \dots + b_{n_b}}{a_1 z^{n_a-1} + a_2 z^{n_a-2} + \dots + a_{n_a}}\]

with the coefficient vectors b and a given in descending powers of $z$. Their lengths nb and na are structural parameters and fix the shape of the model, while the coefficients themselves are ordinary parameters and may be changed without rebuilding.

The realization is the controller canonical form on an internal signal $w$ satisfying $A(z) w = u$ and $y = B(z) w$:

\[a_1 w(k) = u(k) - \sum_{j=1}^{n_a-1} a_{j+1} w(k-j)\]

\[y(k) = \sum_{i=1}^{n_a} \bar{b}_i w(k-i+1)\]

where $\bar{b}$ is b right-aligned to length $n_a$ by zero-padding in front, matching the convention of tf(num, den, Ts). The leading denominator coefficient a[1] must be nonzero. The state x is the shift register holding $w(k), w(k-1), \dots, w(k-n_a+2)$, so the state dimension is $n_a - 1$. A numerator longer than the denominator is improper and is rejected at model construction. When the numerator is strictly shorter than the denominator, $\bar{b}_1 = 0$ and the component is strictly proper, so the output does not depend on the current input. The sample period is that of the connected clock.

The initial condition is selected through the initialization enum. Because the state is a shift register of a single signal, every variant describes a rest history in which all past values of $w$ share one level:

• InitialOutput(y0) sets $w = y_0 / B(1)$, so the first output equals y0 and, driven by the matching DC input, remains there. This requires nonzero DC transmission $B(1) = \sum_i b_i$; for a zero-DC-gain system such as a differentiator the initial state is not finite, so use InitialState instead.
• InitialState(x0) sets every shift-register entry to x0.
• SteadyState sets $w = u(k) / A(1)$, the equilibrium for the initial input. This requires a finite steady state, $A(1) = \sum_j a_j \neq 0$; for an integrator-like system with a pole at $z = 1$ the initial state is not finite, so use InitialOutput or InitialState instead.

The static-gain case $n_a = 1$ has no state and ignores initialization. The array variants of InitialCondition are not applicable to this scalar component and raise an error.

Usage

DiscreteComponents.DiscreteTransferFunction(a=vcat([1.0], fill(0.0, na - 1)), b=vcat([1.0], fill(0.0, nb - 1)), bb=_dtf_bb(b, na, nb), x_init=_ic_state_level(initialization))

Parameters:

Connectors

• u - This connector represents a real signal as an input to a component (RealInput)
• y - This connector represents a real signal as an output from a component (RealOutput)

Variables

Behavior

using DiscreteComponents #hide
using ModelingToolkit #hide
@variables a #hide
@variables b #hide
@named sys = DiscreteComponents.DiscreteTransferFunction(a=a, b=b) #hide
full_equations(sys) #hide

Source

"""
A discrete transfer function
```math
H(z) = \dfrac{B(z)}{A(z)} = \dfrac{b_1 z^{n_b-1} + b_2 z^{n_b-2} + \dots + b_{n_b}}{a_1 z^{n_a-1} + a_2 z^{n_a-2} + \dots + a_{n_a}}
```
with the coefficient vectors `b` and `a` given in descending powers of ``z``. Their
lengths `nb` and `na` are structural parameters and fix the shape of the model, while the
coefficients themselves are ordinary parameters and may be changed without rebuilding.

The realization is the controller canonical form on an internal signal ``w`` satisfying
``A(z) w = u`` and ``y = B(z) w``:
```math
a_1 w(k) = u(k) - \sum_{j=1}^{n_a-1} a_{j+1} w(k-j)
```
```math
y(k) = \sum_{i=1}^{n_a} \bar{b}_i w(k-i+1)
```
where ``\bar{b}`` is `b` right-aligned to length ``n_a`` by zero-padding in front, matching
the convention of `tf(num, den, Ts)`. The leading denominator coefficient `a[1]` must be
nonzero. The state `x` is the shift register holding
``w(k), w(k-1), \dots, w(k-n_a+2)``, so the state dimension is ``n_a - 1``. A numerator
longer than the denominator is improper and is rejected at model construction. When the
numerator is strictly shorter than the denominator, ``\bar{b}_1 = 0`` and the component is
strictly proper, so the output does not depend on the current input. The sample period is
that of the connected clock.

The initial condition is selected through the `initialization` enum. Because the state is
a shift register of a single signal, every variant describes a rest history in which all
past values of ``w`` share one level:
- `InitialOutput(y0)` sets ``w = y_0 / B(1)``, so the first output equals `y0` and, driven
  by the matching DC input, remains there. This requires nonzero DC transmission
  ``B(1) = \sum_i b_i``; for a zero-DC-gain system such as a differentiator the initial
  state is not finite, so use `InitialState` instead.
- `InitialState(x0)` sets every shift-register entry to `x0`.
- `SteadyState` sets ``w = u(k) / A(1)``, the equilibrium for the initial input. This
  requires a finite steady state, ``A(1) = \sum_j a_j \neq 0``; for an integrator-like
  system with a pole at ``z = 1`` the initial state is not finite, so use `InitialOutput`
  or `InitialState` instead.

The static-gain case ``n_a = 1`` has no state and ignores `initialization`. The array
variants of `InitialCondition` are not applicable to this scalar component and raise an
error.
"""
component DiscreteTransferFunction@[input clk extends Discrete]
  "The input signal"
  u = RealInput@[clk]() {
    "Dyad": {
      "placement": {
        "diagram": {"iconName": "default", "x1": -50, "y1": 450, "x2": 50, "y2": 550, "rot": 0}
      },
      "tags": []
    }
  }
  "The output signal"
  y = RealOutput@[clk]() {
    "Dyad": {
      "placement": {
        "diagram": {"iconName": "default", "x1": 960, "y1": 450, "x2": 1060, "y2": 550, "rot": 0}
      },
      "tags": []
    }
  }
  "Number of denominator coefficients, one more than the state dimension"
  structural parameter na::Integer = 1
  "Number of numerator coefficients; must not exceed na"
  structural parameter nb::Integer = 1
  "Denominator coefficients (descending powers of z)"
  parameter a::Real[na] = vcat([1.0], fill(0.0, na - 1))
  "Numerator coefficients (descending powers of z)"
  parameter b::Real[nb] = vcat([1.0], fill(0.0, nb - 1))
  "Initial-condition specification"
  structural parameter initialization::InitialCondition = DiscreteComponents.InitialCondition.InitialOutput(y0 = 0.0)
  "State dimension"
  final structural parameter nx::Integer = na - 1
  "Numerator right-aligned to length na by zero-padding in front"
  final parameter bb::Real[na] = _dtf_bb(b, na, nb)
  "Common initial level of the shift register (only used by the InitialState variant)"
  final parameter x_init::Real = _ic_state_level(initialization)
  "Shift register: x[i] is the canonical-form signal delayed by i-1 samples"
  variable x::Real[nx]
  "Denominator products `a[j+1]*x[j](k-1)`"
  variable ax::Real[nx]
  "Numerator products `bb[j+1]*x[j](k-1)`"
  variable bx::Real[nx]
relations
  if nx > 0
    # Dyad has no sum over a generator, so each dot product against the register is
    # materialized as an array of products and reduced with `sum`.
    for j in 1:nx
      ax[j] = a[j + 1] * x[j]@(clk-1)
      bx[j] = bb[j + 1] * x[j]@(clk-1)
    end
    "Register head: `a[1]*x[1](k) = u(k) - a[2:na]*x(k-1)`"
    x[1]@clk = (u@clk - sum(ax)) / a[1]
    for j in 2:nx
      x[j]@clk = x[j - 1]@(clk-1)
    end
    if nb == na
      y@clk = b[1] * x[1]@clk + sum(bx)
    else
      # A right-aligned numerator shorter than the denominator has bb[1] = 0, so the
      # output does not depend on u(k). Omitting the term instead of multiplying by a
      # parameter that happens to be zero keeps the component structurally strictly
      # proper, which is what breaks an algebraic loop around the filter.
      y@clk = sum(bx)
    end
    switch initialization
      case InitialOutput
        initial x@(clk-1) = fill(initialization.y0 / sum(b), nx)
      case SteadyState
        initial x@(clk-1) = fill(u@clk / sum(a), nx)
      case InitialState
        initial x@(clk-1) = fill(x_init, nx)
      case InitialOutputArray
        initial x@(clk-1) = error("InitialOutputArray is not supported for DiscreteTransferFunction; use InitialOutput, InitialState, or SteadyState")
      case InitialStateArray
        initial x@(clk-1) = error("InitialStateArray is not supported for DiscreteTransferFunction; use InitialOutput, InitialState, or SteadyState")
    end
  else
    # A(z) is a constant, so H(z) = b[1]/a[1] is a static gain with no state.
    y = b[1] / a[1] * u
  end
metadata {
  "Dyad": {"icons": {"default": "dyad://DiscreteComponents/DiscreteTransferFunction.svg"}}
}
end

"""
A discrete transfer function
```math
H(z) = \dfrac{B(z)}{A(z)} = \dfrac{b_1 z^{n_b-1} + b_2 z^{n_b-2} + \dots + b_{n_b}}{a_1 z^{n_a-1} + a_2 z^{n_a-2} + \dots + a_{n_a}}
```
with the coefficient vectors `b` and `a` given in descending powers of ``z``. Their
lengths `nb` and `na` are structural parameters and fix the shape of the model, while the
coefficients themselves are ordinary parameters and may be changed without rebuilding.

The realization is the controller canonical form on an internal signal ``w`` satisfying
``A(z) w = u`` and ``y = B(z) w``:
```math
a_1 w(k) = u(k) - \sum_{j=1}^{n_a-1} a_{j+1} w(k-j)
```
```math
y(k) = \sum_{i=1}^{n_a} \bar{b}_i w(k-i+1)
```
where ``\bar{b}`` is `b` right-aligned to length ``n_a`` by zero-padding in front, matching
the convention of `tf(num, den, Ts)`. The leading denominator coefficient `a[1]` must be
nonzero. The state `x` is the shift register holding
``w(k), w(k-1), \dots, w(k-n_a+2)``, so the state dimension is ``n_a - 1``. A numerator
longer than the denominator is improper and is rejected at model construction. When the
numerator is strictly shorter than the denominator, ``\bar{b}_1 = 0`` and the component is
strictly proper, so the output does not depend on the current input. The sample period is
that of the connected clock.

The initial condition is selected through the `initialization` enum. Because the state is
a shift register of a single signal, every variant describes a rest history in which all
past values of ``w`` share one level:
- `InitialOutput(y0)` sets ``w = y_0 / B(1)``, so the first output equals `y0` and, driven
  by the matching DC input, remains there. This requires nonzero DC transmission
  ``B(1) = \sum_i b_i``; for a zero-DC-gain system such as a differentiator the initial
  state is not finite, so use `InitialState` instead.
- `InitialState(x0)` sets every shift-register entry to `x0`.
- `SteadyState` sets ``w = u(k) / A(1)``, the equilibrium for the initial input. This
  requires a finite steady state, ``A(1) = \sum_j a_j \neq 0``; for an integrator-like
  system with a pole at ``z = 1`` the initial state is not finite, so use `InitialOutput`
  or `InitialState` instead.

The static-gain case ``n_a = 1`` has no state and ignores `initialization`. The array
variants of `InitialCondition` are not applicable to this scalar component and raise an
error.
"""
component DiscreteTransferFunction
  "The input signal"
  u = RealInput@[clk]() {
    "Dyad": {
      "placement": {
        "diagram": {"iconName": "default", "x1": -50, "y1": 450, "x2": 50, "y2": 550, "rot": 0}
      },
      "tags": []
    }
  }
  "The output signal"
  y = RealOutput@[clk]() {
    "Dyad": {
      "placement": {
        "diagram": {"iconName": "default", "x1": 960, "y1": 450, "x2": 1060, "y2": 550, "rot": 0}
      },
      "tags": []
    }
  }
  "Number of denominator coefficients, one more than the state dimension"
  structural parameter na::Integer = 1
  "Number of numerator coefficients; must not exceed na"
  structural parameter nb::Integer = 1
  "Denominator coefficients (descending powers of z)"
  parameter a::Real[na] = vcat([1.0], fill(0.0, na - 1))
  "Numerator coefficients (descending powers of z)"
  parameter b::Real[nb] = vcat([1.0], fill(0.0, nb - 1))
  "Initial-condition specification"
  structural parameter initialization::InitialCondition = DiscreteComponents.InitialCondition.InitialOutput(y0 = 0.0)
  "State dimension"
  final structural parameter nx::Integer = na - 1
  "Numerator right-aligned to length na by zero-padding in front"
  final parameter bb::Real[na] = _dtf_bb(b, na, nb)
  "Common initial level of the shift register (only used by the InitialState variant)"
  final parameter x_init::Real = _ic_state_level(initialization)
  "Shift register: x[i] is the canonical-form signal delayed by i-1 samples"
  variable x::Real[nx]
  "Denominator products `a[j+1]*x[j](k-1)`"
  variable ax::Real[nx]
  "Numerator products `bb[j+1]*x[j](k-1)`"
  variable bx::Real[nx]
relations
  if nx > 0
    # Dyad has no sum over a generator, so each dot product against the register is
    # materialized as an array of products and reduced with `sum`.
    for j in 1:nx
      ax[j] = a[j + 1] * x[j]@(clk-1)
      bx[j] = bb[j + 1] * x[j]@(clk-1)
    end
    "Register head: `a[1]*x[1](k) = u(k) - a[2:na]*x(k-1)`"
    x[1]@clk = (u@clk - sum(ax)) / a[1]
    for j in 2:nx
      x[j]@clk = x[j - 1]@(clk-1)
    end
    if nb == na
      y@clk = b[1] * x[1]@clk + sum(bx)
    else
      # A right-aligned numerator shorter than the denominator has bb[1] = 0, so the
      # output does not depend on u(k). Omitting the term instead of multiplying by a
      # parameter that happens to be zero keeps the component structurally strictly
      # proper, which is what breaks an algebraic loop around the filter.
      y@clk = sum(bx)
    end
    switch initialization
      case InitialOutput
        initial x@(clk-1) = fill(initialization.y0 / sum(b), nx)
      case SteadyState
        initial x@(clk-1) = fill(u@clk / sum(a), nx)
      case InitialState
        initial x@(clk-1) = fill(x_init, nx)
      case InitialOutputArray
        initial x@(clk-1) = error("InitialOutputArray is not supported for DiscreteTransferFunction; use InitialOutput, InitialState, or SteadyState")
      case InitialStateArray
        initial x@(clk-1) = error("InitialStateArray is not supported for DiscreteTransferFunction; use InitialOutput, InitialState, or SteadyState")
    end
  else
    # A(z) is a constant, so H(z) = b[1]/a[1] is a static gain with no state.
    y = b[1] / a[1] * u
  end
metadata {
  "Dyad": {"icons": {"default": "dyad://DiscreteComponents/DiscreteTransferFunction.svg"}}
}
end

Test Cases

No test cases defined.

Related

• Examples
• Experiments
• Analyses