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
- T is the sampling time of your digital filter in seconds (the time between subsequent samples)
- \(\tau\) is the filter time constant, \(f_c = {1 \over \tau} \left[{rad \over s}\right]\)
- \(input[n]\) is the current signal input sample
- \(input[n-1]\) is the previous signal input sample
- \(output[n]\) is the current filter output sample, and
- \(output[n-1]\) is the previous filter output sample
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:
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\).
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.
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.
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:
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.
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.