#!/usr/bin/env python
# coding: utf-8

# # Characterization of Discrete Systems in the Time Domain
# 
# *This Jupyter notebook is part of a [collection of notebooks](../index.ipynb) in the bachelors module Signals and Systems, Comunications Engineering, Universität Rostock. Please direct questions and suggestions to [Sascha.Spors@uni-rostock.de](mailto:Sascha.Spors@uni-rostock.de).*

# ## Difference Equations
# 
# The relation between input $x[k]$ and output $y[k] = \mathcal{H} \{ x[k] \}$ of a discrete system can be described by relating actual and past samples of the input and output signal to each other. When restricting the relation to a linear superposition of these, this results in a [linear difference equation](https://en.wikipedia.org/wiki/Linear_difference_equation). The interconnection of difference equations to [linear time-invariant (LTI) systems](https://en.wikipedia.org/wiki/LTI_system_theory) is discussed in the following.

# ### Linear Difference Equations with Constant Coefficients
# 
# The relation between the input $x[k] \in \mathbb{C}$ with $k \in \mathbb{Z}$ and output $y[k] \in \mathbb{C}$ of a system in terms of a linear difference equation with constant coefficients $a_n, b_m \in \mathbb{C}$ reads
# 
# \begin{equation}
# \sum_{n = 0}^{N} a_n \; y[k - n] = \sum_{m = 0}^{M} b_m \; x[k -m]
# \end{equation}
# 
# with $a_N \neq 0$. The order of the difference equation is given by $N$. For a unique solution of the difference equation, $N$ initial conditions are required. Under the assumption that the input signal of the system $x[k] = 0$ for $k < 0$ these are given by $y[0], y[-1], \dots, y[-N + 1]$. Other cases can be handled accordingly.

# ### Linear Time-Invariant Systems
# 
# LTI systems are an important class of systems. In order to show that linear difference equations with constant coefficients represent LTI systems, the linearity and time-invariance property has to be checked. A system is linear when the superposition principle holds
# 
# \begin{equation}
# \mathcal{H} \{ A x_1[k] + B x_2[k] \} = A \cdot \mathcal{H} \{ x_1[k] \} + B \cdot \mathcal{H} \{ x_2[k] \}
# \end{equation}
# 
# for arbitraty signals $x_1[k], x_2[k] \in \mathbb{C}$ and coefficients $A, B \in \mathbb{C}$. Above condition is fulfilled by a linear difference equation, since the [summation](https://en.wikipedia.org/wiki/Summation) on left- and right-hand side is a linear operation. A system is time-invariant when
# 
# \begin{equation}
# y[k - \kappa] = \mathcal{H} \{ x[k - \kappa] \}
# \end{equation}
# 
# with $y[k] = \mathcal{H} \{ x[k] \}$ holds for $\kappa \in \mathbb{Z}$. Substituting $k$ by $k -\kappa$ in above linear difference equation shows that the time-invariance is fulfilled.

# ### Computation of the Output Signal
# 
# The solution of above difference equation is derived by extracting $y[k]$ from the sum on the left-hand-side and rearranging terms. This results in
# 
# \begin{equation}
# y[k] = \frac{1}{a_0} \left( \sum_{m=0}^{M} b_m \; x[k-m] - \sum_{n=1}^{N} a_n \; y[k-n] \right)
# \end{equation}
# 
# for $a_0 \neq 0$. The output signal is given by a weighted superposition of the actual and past samples of the input signal $x[k]$ and past samples of the output signal $y[k]$. The second sum is known as *recursive* part since it can be interpreted as a feedback of past output samples. The dependency on the input samples (the first sum) is known as *non-recursive* part.
# 
# It can be concluded that linear difference equations with constant coefficients represent LTI systems. The coefficients together with the initial conditions characterize an LTI system fully. However, not all LTI systems can be described by linear difference equations of finite order $N$. Idealized systems may result in a linear difference equation of infinite order, for instance.

# ### Recursive and Non-Recursive Systems
# 
# Above analysis of the solution of a linear difference equation revealed two contributions to the output signal $y[k]$: a non-recursive part which depends only on the current and past input samples $x[k]$ and a recursive part which depends on past output samples $y[k]$. A system is said to be
# 
# * non-recursive if its output depends only on the input signal $x[k]$. This implies that $a_n = 0$ for $n > 0$.
# * recursive if its output depends on the input signal $x[k]$ and past output samples $y[k]$. This implies that at least one $a_n \neq 0$ for $n > 0$.
# 
# The classification into non-recursive and recursive systems has far reaching implications on the theory and practical realization of recursive systems. The feedback of past output samples in recursive systems for instance, may lead to instabilities and numerical issues. The former are discussed in detail later.

# ### Examples
# 
# Some examples for discrete systems defined by their input/output relation formulated in terms of a linear difference equation are given in the following.

# #### Delay
# 
# The output signal $y[k] = \mathcal{H} \{ x[k] \}$ of a discrete system which delays the input $x[k]$ by $\kappa \in \mathbb{N}$ samples is given as
# 
# \begin{equation}
# y[k] = x[k - \kappa]
# \end{equation}
# 
# The order of the system is $N=0$. The coefficients of the linear difference equation read
# 
# \begin{align}
# a_0 &= 1 \\
# b_m &= \begin{cases}
# 1 & \text{for } m=\kappa \\
# 0 & \text{otherwise}
# \end{cases}
# \end{align}
# 
# The system is non-recursive as can be deduced from its coefficients.
# 
# The output signal $y[k]$ for a cosine signal at the input $x[k] = \cos[\Omega k]$ is computed by solving the difference equation for zero initial conditions and $\kappa = 5$. The Python module [`scipy.signal`](http://docs.scipy.org/doc/scipy/reference/signal.html) offers various functionality for continuous and discrete systems. The function [`scipy.signal.filter`](http://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.lfilter.html#scipy.signal.lfilter) implements above formula for the computation of the output signal $y[k]$ of a discrete LTI system given its coefficients and initial conditions. The coefficients $a_n$ and $b_m$ of the system are stored in the vectors `a` and `b`.

# In[1]:


import numpy as np
from scipy import signal
import matplotlib.pyplot as plt
get_ipython().run_line_magic('matplotlib', 'inline')

a = 1.0
b = [0.0, 0.0, 0.0, 0.0, 1.0]

k = np.arange(30)
x = np.cos(2*np.pi/15 * k)
y = signal.lfilter(b, a, x)


# The input $x[k]$ and output signal $y[k]$ are plotted for illustration.

# In[2]:


plt.figure(figsize=(10, 4))
plt.subplot(121)
plt.stem(k, x)
plt.xlabel('$k$')
plt.ylabel(r'$x[k]$')
plt.ylim([-1.1, 1.1])

plt.subplot(122)
plt.stem(k, y)
plt.xlabel('$k$')
plt.ylabel('$y[k]$')
plt.ylim([-1.1, 1.1])
plt.tight_layout()


# **Exercise**
# 
# * Change the delay $\kappa$ of the system and rerun the example.

# #### Moving Average
# 
# The [moving average](https://en.wikipedia.org/wiki/Moving_average) of a time-series is frequently used to smooth out short-term fluctuations in order to emphasize longer-term trends. A system which calculates the moving average of its input signal $x[k]$ is represented by
# 
# \begin{equation}
# y[k] = \frac{1}{M} \sum_{\kappa = 0}^{M-1} x[k - \kappa]
# \end{equation}
# 
# The order of the system is $N=0$. The coefficients of the linear difference equation read
# 
# \begin{align}
# a_0 &= 1 \\
# b_m &= \begin{cases}
# \frac{1}{M} & \text{for } 0 \leq m < M \\
# 0 & \text{otherwise}
# \end{cases}
# \end{align}
# 
# The system is non-recursive as can be deduced from its coefficients. The output of the system for a cosine signal at the input superimposed with noise is computed as illustration of the functionality of the moving average.

# In[3]:


np.random.seed(seed=0)

M = 10
a = [1.0]
b = 1/M * np.ones(M)

k = np.arange(0, 60)
x = np.cos(2*np.pi/30 * k) + .2 * np.random.normal(size=(len(k)))
y = signal.lfilter(b, a, x)

plt.figure(figsize=(10, 4))
plt.subplot(121)
plt.stem(k, x)
plt.xlabel('$k$')
plt.ylabel(r'$x[k]$')
plt.ylim([-1.5, 1.5])

plt.subplot(122)
plt.stem(k, y)
plt.xlabel('$k$')
plt.ylabel('$y[k]$')
plt.ylim([-1.5, 1.5])
plt.tight_layout()


# **Exercise**
# 
# * Change the length `M` of the moving average. How does the output of the system change?

# #### Second-Order System
# 
# The recursive LTI system with the following linear difference equation is given
# 
# \begin{equation}
# y[k] - y[k-1] + \frac{1}{2} y[k-2] = x[k]
# \end{equation}
# 
# The order of the system is $N=2$. The coefficients of the system can be derived from its difference equation as $a_0 = 1$, $a_1 = -1$, $a_2 = \frac{1}{2}$ and $b_0 = 1$. The output signal for a rectangular signal at the input $x[k] = \text{rect}_{20}[k]$ is computed and plotted for zero initial conditions.

# In[4]:


def rect(k, N):
    return np.where((k >= 0) & (k < N), 1.0, 0.0)


a = [1.0, -1.0, 1/2]
b = [1.0]

k = np.arange(0, 40)
x = rect(k, 20)
y = signal.lfilter(b, a, x)

plt.figure(figsize=(10, 4))
plt.subplot(121)
plt.stem(k, x)
plt.xlabel('$k$')
plt.ylabel(r'$x[k]$')
plt.ylim([-0.7, 2.6])

plt.subplot(122)
plt.stem(k, y)
plt.xlabel('$k$')
plt.ylabel('$y[k]$')
plt.ylim([-0.7, 2.6])
plt.tight_layout()


# **Copyright**
# 
# The notebooks are provided as [Open Educational Resource](https://de.wikipedia.org/wiki/Open_Educational_Resources). Feel free to use the notebooks for your own educational purposes. The text is licensed under [Creative Commons Attribution 4.0](https://creativecommons.org/licenses/by/4.0/), the code of the IPython examples under the [MIT license](https://opensource.org/licenses/MIT). Please attribute the work as follows: *Lecture Notes on Signals and Systems* by Sascha Spors.