//! # Transfer functions for discrete time systems.
//!
//! Specialized struct and methods for discrete time transfer functions
//! * time delay
//! * initial value
//! * static gain
//! * ARMA (autoregressive moving average) time evaluation method
//!
//! This module contains the discretization struct of a continuous time
//! transfer function
//! * forward Euler mehtod
//! * backward Euler method
//! * Tustin (trapezoidal) method
use nalgebra::RealField;
use num_complex::Complex;
use num_traits::{Float, Zero};
use std::{
cmp::Ordering,
collections::VecDeque,
fmt::Debug,
iter::Sum,
ops::{Add, Div, Mul},
};
use crate::{enums::Discrete, plots::Plotter, transfer_function::TfGen};
/// Discrete transfer function
pub type Tfz<T> = TfGen<T, Discrete>;
impl<T: Float> Tfz<T> {
/// Time delay for discrete time transfer function.
/// `y(k) = u(k - h)`
/// `G(z) = z^(-h)
///
/// # Arguments
///
/// * `h` - Time delay
///
/// # Example
/// ```
/// use au::{num_complex::Complex, units::Seconds, Tfz};
/// let d = Tfz::delay(2);
/// assert_eq!(0.010000001, d(Complex::new(0., 10.0_f32)).norm());
/// ```
pub fn delay(k: i32) -> impl Fn(Complex<T>) -> Complex<T> {
move |z| z.powi(-k)
}
/// System inital value response to step input.
/// `y(0) = G(z->infinity)`
///
/// # Example
/// ```
/// use au::{poly, Tfz};
/// let tf = Tfz::new(poly!(4.), poly!(1., 5.));
/// assert_eq!(0., tf.init_value());
/// ```
#[must_use]
pub fn init_value(&self) -> T {
let n = self.num().degree();
let d = self.den().degree();
match n.cmp(&d) {
Ordering::Less => T::zero(),
Ordering::Equal => self.num().leading_coeff() / self.den().leading_coeff(),
Ordering::Greater => T::infinity(),
}
}
}
impl<'a, T: 'a + Add<&'a T, Output = T> + Div<Output = T> + Zero> Tfz<T> {
/// Static gain `G(1)`.
/// Ratio between constant output and constant input.
/// Static gain is defined only for transfer functions of 0 type.
///
/// Example
///
/// ```
/// use au::{poly, Tfz};
/// let tf = Tfz::new(poly!(5., -3.),poly!(2., 5., -6.));
/// assert_eq!(2., tf.static_gain());
/// ```
#[must_use]
pub fn static_gain(&'a self) -> T {
let n = self
.num()
.as_slice()
.iter()
.fold(T::zero(), |acc, c| acc + c);
let d = self
.den()
.as_slice()
.iter()
.fold(T::zero(), |acc, c| acc + c);
n / d
}
}
impl<T: Float + RealField> Tfz<T> {
/// System stability. Checks if all poles are inside the unit circle.
///
/// # Example
///
/// ```
/// use au::{Poly, Tfz};
/// let tfz = Tfz::new(Poly::new_from_coeffs(&[1.]), Poly::new_from_roots(&[0.5, 1.5]));
/// assert!(!tfz.is_stable());
/// ```
#[must_use]
pub fn is_stable(&self) -> bool {
self.complex_poles().iter().all(|p| p.norm() < T::one())
}
}
/// Macro defining the common behaviour when creating the arma iterator.
///
/// # Arguments
///
/// * `self` - `self` parameter keyword
/// * `y_coeffs` - vector containing the coefficients of the output
/// * `u_coeffs` - vector containing the coefficients of the input
/// * `y` - queue containing the calculated outputs
/// * `u` - queue containing the supplied inputs
macro_rules! arma {
($self:ident, $y_coeffs:ident, $u_coeffs:ident, $y:ident, $u:ident) => {{
let g = $self.normalize();
let n_n = g.num().degree().unwrap_or(0);
let n_d = g.den().degree().unwrap_or(0);
let n = n_n.max(n_d);
// The front is the lowest order coefficient.
// The back is the higher order coefficient.
// The higher degree terms are the more recent.
// The last coefficient is always 1, because g is normalized.
// [a0, a1, a2, ..., a(n-1), 1]
let mut output_coefficients = g.den().coeffs();
// Remove the last coefficient by truncating the vector by one.
// This is done because the last coefficient of the denominator corresponds
// to the currently calculated output.
output_coefficients.truncate(n_d);
// [a0, a1, a2, ..., a(n-1)]
$y_coeffs = output_coefficients;
// [b0, b1, b2, ..., bm]
$u_coeffs = g.num().coeffs();
// The coefficients do not need to be extended with zeros,
// when the coffiecients are 'zipped' with the VecDeque, the zip stops at the
// shortest iterator.
let length = n + 1;
// The front is the oldest calculated output.
// [y(k-n), y(k-n+1), ..., y(k-1), y(k)]
$y = VecDeque::from(vec![T::zero(); length]);
// The front is the oldest input.
// [u(k-n), u(k-n+1), ..., u(k-1), u(k)]
$u = VecDeque::from(vec![T::zero(); length]);
}};
}
impl<T: Float + Mul<Output = T> + Sum> Tfz<T> {
/// Autoregressive moving average representation of a discrete transfer function
/// It transforms the transfer function into time domain input-output
/// difference equation.
/// ```text
/// b_n*z^n + b_(n-1)*z^(n-1) + ... + b_1*z + b_0
/// Y(z) = G(z)U(z) = --------------------------------------------- U(z)
/// z^n + a_(n-1)*z^(n-1) + ... + a_1*z + a_0
///
/// y(k) = - a_(n-1)*y(k-1) - ... - a_1*y(k-n+1) - a_0*y(k-n) +
/// + b_n*u(k) + b_(n-1)*u(k-1) + ... + b_1*u(k-n+1) + b_0*u(k-n)
/// ```
///
/// # Arguments
///
/// * `input` - Input function
///
/// # Example
/// ```
/// use au::{poly, signals::discrete, Tfz};
/// let tfz = Tfz::new(poly!(1., 2., 3.), poly!(0., 0., 0., 1.));
/// let mut iter = tfz.arma_fn(discrete::step(1., 0));
/// assert_eq!(Some(0.), iter.next());
/// assert_eq!(Some(3.), iter.next());
/// assert_eq!(Some(5.), iter.next());
/// assert_eq!(Some(6.), iter.next());
/// ```
pub fn arma_fn<F>(&self, input: F) -> ArmaFn<F, T>
where
F: Fn(usize) -> T,
{
let y_coeffs: Vec<_>;
let u_coeffs: Vec<_>;
let y: VecDeque<_>;
let u: VecDeque<_>;
arma!(self, y_coeffs, u_coeffs, y, u);
ArmaFn {
y_coeffs,
u_coeffs,
y,
u,
input,
k: 0,
}
}
/// Autoregressive moving average representation of a discrete transfer function
/// It transforms the transfer function into time domain input-output
/// difference equation.
/// ```text
/// b_n*z^n + b_(n-1)*z^(n-1) + ... + b_1*z + b_0
/// Y(z) = G(z)U(z) = --------------------------------------------- U(z)
/// z^n + a_(n-1)*z^(n-1) + ... + a_1*z + a_0
///
/// y(k) = - a_(n-1)*y(k-1) - ... - a_1*y(k-n+1) - a_0*y(k-n) +
/// + b_n*u(k) + b_(n-1)*u(k-1) + ... + b_1*u(k-n+1) + b_0*u(k-n)
/// ```
///
/// # Arguments
///
/// * `iter` - Iterator supplying the input data to the model
///
/// # Example
/// ```
/// use au::{poly, signals::discrete, Tfz};
/// let tfz = Tfz::new(poly!(1., 2., 3.), poly!(0., 0., 0., 1.));
/// let mut iter = tfz.arma_iter(std::iter::repeat(1.));
/// assert_eq!(Some(0.), iter.next());
/// assert_eq!(Some(3.), iter.next());
/// assert_eq!(Some(5.), iter.next());
/// assert_eq!(Some(6.), iter.next());
/// ```
pub fn arma_iter<I, II>(&self, iter: II) -> ArmaIter<I, T>
where
II: IntoIterator<Item = T, IntoIter = I>,
I: Iterator<Item = T>,
{
let y_coeffs: Vec<_>;
let u_coeffs: Vec<_>;
let y: VecDeque<_>;
let u: VecDeque<_>;
arma!(self, y_coeffs, u_coeffs, y, u);
ArmaIter {
y_coeffs,
u_coeffs,
y,
u,
iter: iter.into_iter(),
}
}
}
/// Iterator for the autoregressive moving average model of a discrete
/// transfer function.
/// The input is supplied through a function.
#[derive(Debug)]
pub struct ArmaFn<F, T>
where
F: Fn(usize) -> T,
{
/// y coefficients
y_coeffs: Vec<T>,
/// u coefficients
u_coeffs: Vec<T>,
/// y queue buffer
y: VecDeque<T>,
/// u queue buffer
u: VecDeque<T>,
/// input function
input: F,
/// step
k: usize,
}
/// Macro containing the common iteration steps of the ARMA model
///
/// # Arguments
///
/// * `self` - `self` keyword parameter
macro_rules! arma_iter {
($self:ident, $current_input:ident) => {{
// Push the current input into the most recent position of the input buffer.
$self.u.push_back($current_input);
// Discard oldest input.
$self.u.pop_front();
let input: T = $self
.u_coeffs
.iter()
.zip(&$self.u)
.map(|(&c, &u)| c * u)
.sum();
// Push zero in the last position shifting output values one step back
// in time, zero suppress last coefficient which shall be the current
// calculated output value.
$self.y.push_back(T::zero());
// Discard oldest output.
$self.y.pop_front();
let old_output: T = $self
.y_coeffs
.iter()
.zip(&$self.y)
.map(|(&c, &y)| c * y)
.sum();
// Calculate the output.
let new_y = input - old_output;
// Put the new calculated value in the last position of the buffer.
// `back_mut` returns None if the Deque is empty, this should never happen.
debug_assert!(!$self.y.is_empty());
*$self.y.back_mut()? = new_y;
Some(new_y)
}};
}
impl<F, T> Iterator for ArmaFn<F, T>
where
F: Fn(usize) -> T,
T: Float + Mul<Output = T> + Sum,
{
type Item = T;
fn next(&mut self) -> Option<Self::Item> {
let current_input = (self.input)(self.k);
self.k += 1;
arma_iter!(self, current_input)
}
}
/// Iterator for the autoregressive moving average model of a discrete
/// transfer function.
/// The input is supplied through an iterator.
#[derive(Debug)]
pub struct ArmaIter<I, T>
where
I: Iterator,
{
/// y coefficients
y_coeffs: Vec<T>,
/// u coefficients
u_coeffs: Vec<T>,
/// y queue buffer
y: VecDeque<T>,
/// u queue buffer
u: VecDeque<T>,
/// input iterator
iter: I,
}
impl<I, T> Iterator for ArmaIter<I, T>
where
I: Iterator<Item = T>,
T: Float + Mul<Output = T> + Sum,
{
type Item = T;
fn next(&mut self) -> Option<Self::Item> {
let current_input = self.iter.next()?;
arma_iter!(self, current_input)
}
}
impl<T: Float> Plotter<T> for Tfz<T> {
/// Evaluate the transfer function at the given value.
///
/// # Arguments
///
/// * `theta` - angle at which the function is evaluated.
/// Evaluation occurs at G(e^(i*theta)).
fn eval_point(&self, theta: T) -> Complex<T> {
self.eval(&Complex::from_polar(T::one(), theta))
}
}
#[cfg(test)]
mod tests {
use super::*;
use crate::{poly, polynomial::Poly, signals::discrete, units::ToDecibel};
use num_complex::Complex64;
#[test]
fn tfz() {
let _ = Tfz::new(poly!(1.), poly!(1., 2., 3.));
}
#[test]
fn delay() {
let d = Tfz::delay(2);
assert_relative_eq!(0.010_000_001, d(Complex::new(0., 10.0_f32)).norm());
}
#[test]
fn initial_value() {
let tf = Tfz::new(poly!(4.), poly!(1., 5.));
assert_relative_eq!(0., tf.init_value());
let tf = Tfz::new(poly!(4., 10.), poly!(1., 5.));
assert_relative_eq!(2., tf.init_value());
let tf = Tfz::new(poly!(4., 1.), poly!(5.));
assert_relative_eq!(std::f32::INFINITY, tf.init_value());
}
#[test]
fn static_gain() {
let tf = Tfz::new(poly!(5., -3.), poly!(2., 5., -6.));
assert_relative_eq!(2., tf.static_gain());
}
#[test]
fn stability() {
let stable_den = Poly::new_from_roots(&[-0.3, 0.5]);
let stable_tf = Tfz::new(poly!(1., 2.), stable_den);
assert!(stable_tf.is_stable());
let unstable_den = Poly::new_from_roots(&[0., -2.]);
let unstable_tf = Tfz::new(poly!(1., 2.), unstable_den);
assert!(!unstable_tf.is_stable());
}
#[test]
fn eval() {
let tf = Tfz::new(
Poly::new_from_coeffs(&[2., 20.]),
Poly::new_from_coeffs(&[1., 0.1]),
);
let z = 0.5 * Complex64::i();
let g = tf.eval(&z);
assert_relative_eq!(20.159, g.norm().to_db(), max_relative = 1e-4);
assert_relative_eq!(75.828, g.arg().to_degrees(), max_relative = 1e-4);
}
#[test]
fn arma() {
let tfz = Tfz::new(poly!(0.5_f32), poly!(-0.5, 1.));
let mut iter = tfz.arma_fn(discrete::impulse(1., 0)).take(6);
assert_eq!(Some(0.), iter.next());
assert_eq!(Some(0.5), iter.next());
assert_eq!(Some(0.25), iter.next());
assert_eq!(Some(0.125), iter.next());
assert_eq!(Some(0.0625), iter.next());
assert_eq!(Some(0.03125), iter.next());
assert_eq!(None, iter.next());
}
#[test]
fn arma_iter() {
use std::iter;
let tfz = Tfz::new(poly!(0.5_f32), poly!(-0.5, 1.));
let mut iter = tfz.arma_iter(iter::once(1.).chain(iter::repeat(0.)).take(6));
assert_eq!(Some(0.), iter.next());
assert_eq!(Some(0.5), iter.next());
assert_eq!(Some(0.25), iter.next());
assert_eq!(Some(0.125), iter.next());
assert_eq!(Some(0.0625), iter.next());
assert_eq!(Some(0.03125), iter.next());
assert_eq!(None, iter.next());
}
}