# Converting continuous filters to discrete implementation

Here’s a trick to quickly convert a continuous, strictly proper transfer function into pseudocode for a discrete time implementation.  Let’s start with a first order low pass filter:

$$\frac{Y(s)}{U(s)} = \frac{1}{\tau s + 1}$$

The first step is to cross multiply the terms (for simplicity I’ll omit the Y(s) and U(s) for the rest of this post).

$$\tau sY + Y = U \rightarrow \tau \dot{Y} + Y = U$$

Then solve the transfer function for the highest ordered derivative of $Y$ (meaning the $Y$ term with the most $\cdot$ above it), then integrate both sides of the equation until $Y$ (with no dots) becomes the left hand side.  The integration is with respect to time, but I’ve omitted that for clarity in the expression.

$$\dot{Y} = \frac{U-Y}{\tau} \rightarrow \int{\dot{Y}} = \int{\frac{U-Y}{\tau}} \rightarrow Y = \int{\frac{U-Y}{\tau}}$$

Now replace the $=\int\left(\right)$ with a $+= dt\left(\right)$ to achieve the following pseudocode

$$Y += dt\left( \frac{U-Y}{\tau} \right)$$

Here, $dt$ is the time between filter updates, and $Y$ is a static variable whose value persists between filter updates.  What we did is solve for the derivative of the filter output ($\dot{Y}$) and then integrated the derivative using Euler’s method to get the filter’s software implementation.

We can do the same thing with a second order filter:

$$\frac{Y}{U} = \frac{\omega_n^2}{s^2 + 2\zeta\omega_ns + \omega_n^2}$$

$$\ddot{Y} + 2\zeta\omega_n \dot{Y} + \omega_n^2 Y = \omega_n^2 U \rightarrow \ddot{Y} = -2\zeta\omega_n \dot{Y} + \omega_n^2(U-Y)$$

Here we’ve again solved for the highest ordered derivative of $Y$.  Now what we’ll do is integrate both sides of the equation until we are again left with only $Y$ on the left hand side.

$$\int{\ddot{Y}} = \int{-2\zeta\omega_n \dot{Y} + \omega_n^2(U-Y)} \rightarrow \dot{Y} = -2\zeta\omega_n Y + \int{\omega_n^2(U-Y)}$$

$$\int{\dot{Y}} = \int{\left(-2\zeta\omega_n Y + \int{\omega_n^2(U-Y)}\right)} \rightarrow Y = \int{\left(-2\zeta\omega_n Y + \int{\omega_n^2(U-Y)}\right)}$$

Let’s break the previous expression into two parts so only one $\int()$ appears on the right hand of any equals sign.

$$k = \int{\omega_n^2(U-Y)}$$

$$Y = \int{\left(-2\zeta\omega_n Y + k\right)}$$

Now we pull the same trick we did before by replacing each $=\int()$ with $+= dt()$ to arrive at the software pseudocode.

$$k += dt\left(\omega_n^2(U-Y) \right)$$

$$Y += dt\left(-2\zeta\omega_n Y + k \right)$$

Now we have two static variables, $k$ and $Y$, that each represent a separate state variable of the second order filter.

As an exercise for the reader (oh man, did I really just say that?), show that the second order band pass filter

$$\frac{Y}{U} = \frac{2\zeta\omega_n s}{s^2 + 2\zeta\omega_ns + \omega_n^2}$$

has the following pseudocode implementation

$$k += dt(-\omega_n^2Y)$$

$$Y += dt\left(k + 2\zeta\omega_n(U – Y)\right)$$

# Testing equivalence of rationals using natural encodings

I’ve recently been spending a decent amount of my free time in self-study of Ulrich Kohlenbach’s Applied Proof Theory: Proof Interpretations and their Use in Mathematics.  As someone who has no experience with Constructivism, it is definitely challenging, but there is no better way to learn than hands-on.  Here I explore one of the topics the book mentions – encoding rational numbers into the naturals – using Haskell.

The idea is simple – take operations on $\mathbb{Q}$ and transform them into operations on $\mathbb{N}$.  We’ll use a special encoding to assign a unique $n \in \mathbb{N}$ for each $q \in \mathbb{Q}$ (yes, there are just as many natural numbers as there are rationals!).  We can then perform operations (in this exercise, equivalence tests) in the encoded set ($\mathbb{N}$) and have the results be valid in the set we care about ($\mathbb{Q}$).  Some of the steps below will seem silly, but this is just an exercise in theory.

The gist is available here.

Before we jump in, we’ll establish some prerequisite machinery for our exploration.  First, we’ll take the Nat instance from the exploration on infinite search and expand it a bit using the nomenclature in Kohlenbach.  Next, we’ll borrow the least function from Seemingly Impossible Functional Programs, slightly modified:

least p = if p 0 then 0 else 1 + least(\n -> p (1 + n))

Lastly, we’ll build the one-point compactification of $\mathbb{N}$ for use in the infinite search monad.  The premise of this monad is that by defining an infinite set, an algorithm for exhaustively searching that set is automatically generated.

natural :: Set Nat
natural = union (singleton Z) (S <$> natural) Now we’re ready for math. Rationals can be thought of as pairs of integers, one being the numerator$num$and the other the denominator$den$. There are two challenges in trying to encode this information into$\mathbb{N}$. 1. The encoding has two inputs and one output. 2. Rationals have sign, but the naturals have no sign. We’ll tackle 2 first. If$q$is positive, then its numerator$num$is mapped to the even natural$2|num|$and the odd natural$2|num|-1$otherwise. In this way, we can encode the sign information in$q$into the parity of a natural$n_{num}$. The denominator$den$, on the other hand, maps to separate natural$n_{den}$through$|den| – 1$. Now that we have$num$and$den$mapping to two naturals, we’ll map the pair of naturals$(n_{num}, n_{den})$to a single natural$q_c \in \mathbb{N}$, with$q_c$meaning an encoded rational. Using the Cantor pairing function $$j(x^0, y^0) := \left\{ \begin{array}{lr} \mbox{min }u \le (x+y)^2 + 3x + y [2u = (x+y)^2 +3x + y] \mbox{ if it exists} \\ 0\mbox{ otherwise} \end{array} \right.$$ The function numden2QCode below performs this encoding, returning a natural number of type QCode representing a rational from the rational’s numerator and denominator. numden2QCode :: Integer -> Integer -> QCode numden2QCode num den | (num >= 0 && den > 0) || (num <= 0 && den < 0) = j (fromInteger$ 2*abs(num)) (fromInteger $abs(den)-1) | otherwise = j (fromInteger$ 2*abs(num)-1) (fromInteger $abs(den)-1) j :: Nat -> Nat -> QCode j x y = if k <= t then QCode k else QCode 0 where k = least (\u -> 2*u == t) t = (x+y)*(x+y) + 3*x + y Testing for equivalency of rationals using QCodes might seem like a straightforward process of comparing the QCodes themselves, but one characteristic of the numden2QCode encoding we’ve developed is that$(num, den)$(1, 2) has a different QCode than (2,4), even though these rationals are equivalent. To test for equivalence, we’ll project the QCode back into the$n_{num}$and$n_{den}$pairs . $$j_1(z) := \mbox{min }x \le z[\exists y \le z (j(x, y) = z)]$$ $$j_2(z) := \mbox{min }y \le z[\exists x \le z (j(x, y) = z)]$$ Here,$z$is our rational encoding$q_c$,$j_1$returns$n_{num}$and$j_2$returns$n_{den}$. For example, >j1$ j 2 3
S (S Z)

>j2 $j 3 4 S (S (S (S Z))) Lastly, to test equivalence, we’ll extend QCode to be an instance of the Eq typeclass. The trick here is that we need to handle odd and even$n_{num}$cases separately, since that indicates the rational is positive or negative. Once we back out the$n_{num}$and$n_{den}$, we convert back to rational num/den form and equate using cross multiplication. instance Eq QCode where (==) n1 n2 | evenNat j1n1 && evenNat j1n2 = (j1n1/2)*(j2n2 + 1) == (j1n2/2)*(j2n1 + 1) | oddNat j1n1 && oddNat j1n2 = ((j1n1+1)/2)*(j2n2 + 1) == ((j1n2+1)/2)*(j2n1 + 1) | otherwise = False where j1n1 = j1 n1 j1n2 = j1 n2 j2n1 = j2 n1 j2n2 = j2 n2 Finally, we can equate some rationals! print$ numden2QCode (-1) 1 == numden2QCode (-2) 2 -- True
print $numden2QCode (-1) 1 == numden2QCode (-3) 3 -- True print$ numden2QCode (-1) 2 == numden2QCode 2 (-4) -- True
print $numden2QCode 1 2 == numden2QCode 2 4 -- True print$ numden2QCode 1 2    == numden2QCode 1 3    -- False

Note that you’ll probably want to compile these into an executable like I did because they are slow otherwise.  In executable form, they are done in about 1.5 seconds.  In a future post, I’ll explore why a few changes to the Num instance for Nat reduced execution time from 88 seconds down to 1.5.