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) $$