{
"cells": [
{
"cell_type": "markdown",
"id": "f2c59caa",
"metadata": {},
"source": [
"# How to Implement a First-Order Low-Pass Filter in Discrete Time\n",
"\n",
"Copyright © 2021 [Joshua Marshall](https://www.ece.queensu.ca/people/j-marshall/)\n",
"\n",
"We often teach or learn about filters in continuous time, but then need to implement them in discrete time (e.g., in code) on data acquired at discrete sample times. This notebook shows one way to design and implement a simple first-order low-pass filter in discrete time."
]
},
{
"attachments": {},
"cell_type": "markdown",
"id": "9fd2f6bd",
"metadata": {},
"source": [
"## Filter Model Derivation\n",
"\n",
"Let's start with an integrator, which can be approximated in discrete time by any number of methods. For this example, we use a [trapezoidal approximation](https://en.wikipedia.org/wiki/Trapezoidal_rule)\n",
"\n",
"$$\\int_{(k-1)T}^{kT}f(\\tau)d\\tau \\approx \\frac{T}{2}\\left(f_k+f_{k-1}\\right),$$\n",
"\n",
"where $T>0$ be the step tiem (i.e., sample period), $k=1,2,\\ldots$ is the time index, and $f_k:=f(kT)$. Thus, the integral from $t=0$ to the current time $kT$ is \n",
"\n",
"$$i_k=\\int_{0}^{kT}f(\\tau)d\\tau \\approx i_{k-1} + \\frac{T}{2}\\left(f_k+f_{k-1}\\right).$$\n",
"\n",
"Now, let the [$z$-transform](https://en.wikipedia.org/wiki/Z-transform) of $i(t)$ and $f(t)$ be $I(z)$ and $F(z)$, respectively, where $z^{-1}$ denotes a unit time delay of period $T$. Therefore, taking the $z$-transform of the integral approximation above yields\n",
"\n",
"$$G_I(z) = \\frac{T}{2}\\left(\\frac{1+z^{-1}}{1-z^{-1}}\\right) =\\frac{I(z)}{F(z)}.$$\n",
"\n",
"This formula is often called [Tustin's method](https://en.wikipedia.org/wiki/Bilinear_transform). For those more familar with the Laplace domain, Tustin's method uses a first-order approximation for the integral $1/s$ such that\n",
"\n",
"$$s\\approx \\frac{2}{T}\\left(\\frac{1-z^{-1}}{1+z^{-1}}\\right).$$\n",
"\n",
"We now have the machinery to build a first-order low-pass filter. For simplicity, let $Y(z)$ be the out put of our filter with input $X(z)$. In continuious time, the [frequency response of a low-pass filter](https://en.wikipedia.org/wiki/Low-pass_filter#Frequency_response) is given by \n",
"\n",
"$$G_{\\rm LPF}(s) = \\frac{\\omega_c}{s+\\omega_c},$$\n",
"\n",
"where $\\omega_c~\\textup{[rad/s]}$ is the filter's [cut-off frequency](https://en.wikipedia.org/wiki/Cutoff_frequency). Substituting our above approximtion for $s$ into $G_{\\rm LPF}(s)$ we obtain\n",
"\n",
"$$G_{\\rm LPF}(z) = \\frac{T\\omega_c(1+z^{-1})}{(T\\omega_c+2)+(T\\omega_c-2)z^{-1}} = \\frac{Y(z)}{X(z)},$$\n",
"\n",
"where $Y(z)$ is the filtered version of $X(z)$.\n",
"\n",
"Therefore, using our expression for the filter transfer function $G_{\\rm LPF}(z)$ we can write\n",
"\n",
"$$ Y(z)\\left(T\\omega_c+2 + (T\\omega_c-2)z^{-1}\\right) = T\\omega_c(1+z^{-1})X(z).$$ \n",
"\n",
"Thus, in the time domain, we obtain\n",
"\n",
"$$(T\\omega_c + 2)y_k + (T\\omega_c-2)y_{k-1} = T\\omega_c\\left(x_k+x_{k-1}\\right),$$\n",
"\n",
"which can be rearranged to obtain\n",
"\n",
"$$y_k = \\left(\\frac{2-T\\omega_c}{2+T\\omega_c}\\right)y_{k-1} + \\left(\\frac{T\\omega_c}{2+T\\omega_c}\\right)\\left(x_k+x_{k-1}\\right),$$ \n",
"\n",
"where $x_k$ is the input at time $t=kT$ and $y_k$ is the filtered output at time $t=kT$."
]
},
{
"cell_type": "markdown",
"id": "9302491a",
"metadata": {},
"source": [
"## Example Implementation\n",
"\n",
"Let's try out our design on some contrived data.\n",
"\n",
"### Create Test Input Data\n",
"\n",
"Let's create some input data with two frequencies that are sufficiently far apart."
]
},
{
"cell_type": "code",
"execution_count": 1,
"id": "306f534c",
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n"
]
},
{
"cell_type": "markdown",
"id": "2ec264ba",
"metadata": {},
"source": [
"Set the simulation time as `SIM_TIME` \\[s\\] and the sample period $T$ as `T` \\[s\\]."
]
},
{
"cell_type": "code",
"execution_count": 2,
"id": "d3621dfa",
"metadata": {},
"outputs": [],
"source": [
"SIM_TIME = 15\n",
"T = 0.01\n"
]
},
{
"cell_type": "code",
"execution_count": 3,
"id": "5ea46ebc",
"metadata": {},
"outputs": [],
"source": [
"t = np.arange(0.0, SIM_TIME, T)\n",
"N = np.size(t)\n"
]
},
{
"cell_type": "markdown",
"id": "aa5ed675",
"metadata": {},
"source": [
"Create two contrived periodic signals $x_1$ and $x_2$ with frequencies $f_1$ and $f_2$ \\[Hz\\], respectively."
]
},
{
"cell_type": "code",
"execution_count": 4,
"id": "a9dbd990",
"metadata": {},
"outputs": [],
"source": [
"f_1 = 0.25\n",
"f_2 = 5.0\n",
"x_1 = np.sin(2.0 * np.pi * f_1 * t)\n",
"x_2 = np.sin(2.0 * np.pi * f_2 * t)\n"
]
},
{
"cell_type": "markdown",
"id": "9cb708eb",
"metadata": {},
"source": [
"### Plot the Input Signal\n",
"\n",
"Now plow the signals."
]
},
{
"cell_type": "code",
"execution_count": 5,
"id": "e1fcbdb0",
"metadata": {},
"outputs": [],
"source": [
"import matplotlib.pyplot as plt"
]
},
{
"cell_type": "code",
"execution_count": 6,
"id": "75462c8f",
"metadata": {},
"outputs": [],
"source": [
"%config InlineBackend.figure_format = 'svg'\n",
"plt.rc('text', usetex=True)\n",
"plt.rc('text.latex', preamble=r'\\usepackage{amsmath}')"
]
},
{
"cell_type": "code",
"execution_count": 7,
"id": "3279cb7b",
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": "\n\n\n",
"text/plain": [
"