I'm writing this guide because I haven't found anywhere all of this information is collected in one place in a step by step manner. If you just need to get a digital filter implemented without worrying too much about the underlying theory, this is your guide!

Step 1: Find the Transfer Function of Your Filter, \(H(s)\)

The starting point is the continuous time transfer function for your filter. For example, a first order low-pass filter with a DC gain of 1 has the transfer function:

$$\large{H(s) = {1 \over \tau s + 1}}$$

Where \(\tau\) is the time constant of the filter. The cutoff frequency, \(f_c\), is equal to \({1 \over \tau} \left[{rad \over s}\right] \).

Step 2: Perform Tustin's Method (Bilinear Transform)

The next step is to convert from a continuous time representation of the system in the s-domain to a discrete time representation in the z-domain using Tustin's Method. This can be achieved with the following substitution:

$$\large{s \leftarrow {2 \over T}{(z-1) \over (z+1)}}$$

That is to say

$$\large{H(z) = H(s) \bigg| _{s = {2 \over T}{(z-1) \over (z+1)} }}$$

where \(T\) is the sampling period of your system. This is an approximation and the faster (smaller) your sampling time, \(T\), the closer the final representation will mimic the analog filter's behaviour. In general, your sampling time should be at least 10 times faster than your time constant:

$$\large{T \leq {\tau \over 10}}$$

Our first order filter then becomes:

$$\large{H(z) = {1 \over {\tau {2 \over T}{{z-1}\over{z+1}}} + 1}}$$

Step 3: Simplify the Numerator and Denominator to Polynomials

In order to represent this filter as a difference equation and implement it in software, we first need to simplify the transfer function until the numerator and denominator are both polynomials in z.

Here are the steps for our first order low-pass filter example.

$$\large{ \begin{align} H(z) & = {1 \over {\tau {2 \over T}{{z-1}\over{z+1}}} + 1} \\[10pt] & = {1 \over {{2 \tau (z-1) \over T(z+1)}} + {T(z+1) \over T(z+1)}} \\[10pt] & = {T(z+1) \over {2 \tau (z-1) + T(z+1)}} \\[10pt] & = {{Tz + T} \over {2 \tau z - 2 \tau + Tz + T}} \\[10pt] & = {{Tz + T} \over {(2 \tau + T)z + (T - 2 \tau)}} \\[10pt] \end{align} } $$

Step 4: Make the Transfer Function Causal

In order to make our transfer function causal, we must multiply the numerator and denominator by \(z^{-n}\) where \(n\) is the highest order of \(z\) found in \(H(z)\). In our first order low-pass filter example, we will multiply the numerator and denominator by \(z^{-1}\).

$$\large{H(z) = {{Tz + T} \over {(2 \tau + T)z + (T- 2 \tau)}} \times {z^{-1} \over z^{-1}} = {{T + Tz^{-1}} \over {(2 \tau + T) + (T - 2 \tau)z^{-1}}}}$$

Step 5: Write the Discrete Transfer Function as a Difference Equation Using Direct Form I

To write \(H(z)\) as a difference equation using Direct Form I, we will first rewrite \(H(z)\) as \(Output(z) \over Input(z)\) and cross multiply.

$$ \large{ \begin{align} H(z) = {Output(z) \over Input(z)} & = {{T + Tz^{-1}} \over {(2 \tau + T) + (T - 2 \tau)z^{-1}}} \\[10pt] Output(z) \left[(2 \tau + T) + (T - 2 \tau)z^{-1}\right] & = \left(T + Tz^{-1}\right)Input(z) \\[10pt] (2 \tau + T) \cdot Output(z) + (T - 2 \tau) \cdot Output(z) \cdot z^{-1} & = T \cdot Input(z) + T \cdot Input(z) \cdot z^{-1} \end{align} } $$

Finally, each multiplication by \(z^{-1}\) is equal to a digital delay of one sample. So the difference equation becomes:

$$\large{(2 \tau + T) \cdot output[n] + (T - 2 \tau) \cdot output[n-1] = T \cdot input[n] + T \cdot input[n-1]}$$

Isolating \(output[n]\) we get:

$$\large{output[n] = {{T \cdot input[n] + T \cdot input[n-1] - (T - 2 \tau) \cdot output[n-1]} \over {2 \tau + T}}}$$

where

This equation can then be implemented in code. For example, if your application is running in real time in an embedded system, it might look something like this:

float tau = 10; float T = 0.1; previousInput = 0; previousOutput = 0; ... if(sampleReady) { sampleReady = 0; input = readInput(); //update your input from whatever source your data comes from output = (T*input + T*previousInput - (T - 2*tau)*previousOutput)/(2*tau + T); previousOutput = output; //keep track of previous values for next sample previousInput = input; }

The previous example is a direct implementation of the difference equation. However, it is more computationally efficient to pre-calculate the difference equation coefficients. The input coefficients are called \(b_n\) and the output coefficients are called \(a_n\).

float tau = 10; float T = 0.1; float a0 = 2*tau + T; float a1 = (T - 2*tau); float b0 = T; float b1 = T; previousInput = 0; previousOutput = 0; ... if(sampleReady) { sampleReady = 0; input = readInput(); //update your input from whatever source your data comes from output = (b0*input + b1*previousInput - a1*previousOutput)/a0; previousOutput = output; //keep track of previous values for next sample previousInput = input; }

An even more efficient way to implement the filter is to normalize your coefficients so that \(a_0\) equals \(1\). You can do that simply by dividing all of your coefficients by \(a_0\). This will eliminate the final division by \(a_0\) seen in the previous example.

float tau = 10; float T = 0.1; float a0 = 2*tau + T; float a1 = (T - 2*tau)/a0; float b0 = T/a0; float b1 = T/a0; a0 = 1; // a0 is now 1 since all other coefficients have been normalized. previousInput = 0; previousOutput = 0; ... if(sampleReady) { sampleReady = 0; input = readInput(); //update your input from whatever source your data comes from output = b0*input + b1*previousInput - a1*previousOutput; previousOutput = output; //keep track of previous values for next sample previousInput = input; }

The difference equation is now reduced to three multiplications, one addition, and one subtraction. If you are using a system that has hardware multiply instructions, this will run very fast. In the second order low-pass filter example below, you will see that pre-calculating and normalizing the coefficients makes even more of a difference on the number of instructions required to carry out the difference equation.

For a first order filter, you require one sample of memory for both the input and output. You can use Direct Form II to realise your difference equation and cut the memory requirement in half.

In order to see accurate results, the value you set for \(T\) in your code must be close to the actual time interval between samples!

Here is a sample output using the above difference equation in Microsoft Excel on a simulated unit step input.

First Order Unit Step Response

Second Order Low-Pass Filter Example

Here I will go through the steps again with a second order low-pass filter. Once again, DC gain was set to 1.

1.

$$\large{H(s) = {\omega_n^2 \over {s^2 + 2 \zeta \omega_n s + \omega_n^2}}}$$

2.

$$\large{H(z) = H(s) \bigg| _{s = {2 \over T}{(z-1) \over (z+1)} } = { \omega_n^2 \over {\left[ {2 \over T}{(z-1) \over (z+1)} \right]^2 + 2 \zeta \omega_n \left[{2 \over T}{(z-1) \over (z+1)}\right] + \omega_n^2}}}$$

3.

$$\large{H(z) = { {\omega_n^2 T^2 z^2 + 2 \omega_n^2 T^2 z + \omega_n^2 T^2 } \over {(4 + 4 \zeta \omega_n T + \omega_n^2 T^2)z^2 + (2 \omega_n^2 T^2 - 8) z + (4 - 4 \zeta \omega_n T + \omega_n^2 T^2)} } }$$

4.

$$\large{H(z) = { {\omega_n^2 T^2 + 2 \omega_n^2 T^2 z^{-1} + \omega_n^2 T^2 z^{-2} } \over {(4 + 4 \zeta \omega_n T + \omega_n^2 T^2) + (2 \omega_n^2 T^2 - 8) z^{-1} + (4 - 4 \zeta \omega_n T + \omega_n^2 T^2) z^{-2}} } }$$

5.

$$\large{output[n] = { {\omega_n^2 T^2 \cdot input[n] + 2 \omega_n^2 T^2 \cdot input[n-1] + \omega_n^2 T^2 \cdot input[n-2] - (2 \omega_n^2 T^2 - 8) \cdot output[n-1] - (4 - 4 \zeta \omega_n T + \omega_n^2 T^2) \cdot output[n-2]} \over {4 + 4 \zeta \omega_n T + \omega_n^2 T^2} } }$$

Of course, we would not want to compute this version of the difference equation multiple times as there are over 40 multiplications and dozens of additions/subtractions. We will instead compute the coefficients and normalize them to \(a_0\).

$$ \large{ \begin{align} a_0 & = 1; \\[10pt] a_1 & = {{2 \omega_n^2 T^2 - 8} \over {4 + 4 \zeta \omega_n T + \omega_n^2 T^2}} \\[10pt] a_2 & = {{4 - 4 \zeta \omega_n T + \omega_n^2 T^2} \over {4 + 4 \zeta \omega_n T + \omega_n^2 T^2}} \\[10pt] b_0 & = {{\omega_n^2 T^2} \over {4 + 4 \zeta \omega_n T + \omega_n^2 T^2}} \\[10pt] b_1 & = {{2 \omega_n^2 T^2} \over {4 + 4 \zeta \omega_n T + \omega_n^2 T^2}} \\[10pt] b_2 & = {{\omega_n^2 T^2} \over {4 + 4 \zeta \omega_n T + \omega_n^2 T^2}} \end{align} } $$

Finally, your code will look something like the following:

float omega_n = 0.2; float zeta = 1; float T = 0.1; float a0 = 4 + 4*zeta*omega_n*T + omega_n*omega_n*T*T; float a1 = (2*omega_n*omega_n*T*T - 8)/a0; float a2 = (4 - 4*zeta*omega_n*T + omega_n*omega_n*T*T)/a0; float b0 = (omega_n*omega_n*T*T)/a0; float b1 = (2*omega_n*omega_n*T*T)/a0; // This could also be b1 = 2*b0 float b2 = (omega_n*omega_n*T*T)/a0; // and this could be b2 = b0 a0 = 1; // a0 is now 1 since all other coefficients have been normalized. input_n1 = 0; // input[n-1] input_n2 = 0; // input[n-2] ... etc output_n1 = 0; output_n2 = 0; ... if(sampleReady) { sampleReady = 0; input = readInput(); //update your input from whatever source your data comes from output = b0*input + b1*input_n1 + b2*input_n2 - a1*output_n1 - a2*output_n2; output_n2 = output_n1; //keep track of previous values for next sample output_n1 = output; input_n2 = input_n1; input_n1 = input; }

I've shown a few optimizations in the code comments. There are others you can do as well if you wish. For example, \(\omega_n^2 T^2\) is computed several times and it may be worth the extra memory to compute it once and reuse the result.

Here are two sample outputs using the above difference equation in Microsoft Excel on a unit step input.

First Order Unit Step Response

First Order Unit Step Response

The second order filter is not for the faint of heart! The intermediate steps have been left out. It would be good practice to see if you can follow this guide and get to the same result yourself!

If you found this guide useful, please let me know!

This guide was made possible by MathJax.