5 Infinite Impulse Response (IIR) Filter Design
IIR Filters
The next type of filter we will be designing are called Infinite Impulse Response (IIR) filter. These filters will use a feedback of the output in order to generate the next sample, causing a infinite ripple in the output. These filters are primarily going to be built by translating the analog filters described in Appendix B, into discrete forms.
Section A) Feedback Discrete Systems
Unlike the FIR filters of chapter 4, the filters in this section will employ feedback of previous outputs to generate the next output. This causes the output of the system to ripple for a long time, hence the name of Infinite Impulse Response (IIR) filters. Another way these filters are described as Pole-Zero (PZ) filters, since the feedback creates poles in the z-transform transfer function.
Section A.1) Impulse Response
We have touched on these types of filters in chapter 3, when we were developing the z-transform. The best example would be the simple case of a first order system
[latex]y(n) = x(n) + a y(n-1)[/latex] (5.1)
The signal flow diagram of which is
Figure 5.1 First Order Signal Flow Diagram.
In order to build an image of these systems and their response, consider the output y(n) to an input that is an impulse
[latex]\delta(n) = \{_{0\text{ otherwise } }^{1 \text{ for n = 0}}[/latex]
Since the input is zero for all n < 0, the output is also 0.
[latex]y(n) = 0 \text{ for n } \lt 0[/latex]
At time 0, the input is a one and then
[latex]y(0) = x(0) + a y(-1) = 1 + a 0 = 1[/latex]
[latex]y(1) = a*y(0) + x(1) = a[/latex]
[latex]y(2) = a*y(1) + x(1) = a^2[/latex]
[latex]y(3) = a*y(2) = a^3[/latex]
…
[latex]y(n) = a^n[/latex] (5.2)
As we can see from this result, the output of such a system from an impulse input is and infinite signal. Hence the name of Infinite Impulse Response (IIR) filter. Although interesting, the impulse response of a discrete system such as this is not directly used to design or interpret an IIR filter. Rather we need to look at the frequency response, which is the primary character we will be designing for. A quick review of Section 3, would be helpful in understanding the follow developments.
Section A.2) Frequency Response
As stated in the previous section, chapter 3 addressed the Frequency Response (FR) of discrete systems. And when we were building FIR filters we were able to first simply write out the FR for a generic set of coefficients. Then solve these for a given FR. In the case of IIR filters, we will need to work this from a different angle. Here we will need to look at how the pole locations affect the FR and how we can develop specific patterns to create the FR that we want.
If we look back at chapter 3, we can see the FR of a system with feedback has poles which pull the [latex]H( e^{j \omega} )[/latex] towards infinity. This effect will mean we will need a distinct pattern of poles to get the desired FR. This phenomena is also seen in the design of an analog filter. As shown here, the Butterworth filter has a set of poles that form a circle near the imaginary or jw axis in the analog s-plane. ( A reading of Appendix B might be in order at this time. )
Figure 5.2 Magnitude of Third-Order Butterworth H(a+jw)
Close observation of the previous circle of poles, we can see that the pattern of three poles work together to create a very smooth low-pass filter along the imaginary (jw) axis. In the case of a discrete system, we will see that the pattern is not easily defined in the z-plane. Rather we will work by taking an analog filter and translating to a discrete implementation.
Section B) Bi-Linear Transform
If you have reviewed the analog filters in Appendix B, you should note that the classic design of analog filters involves the transforming (mapping) of low-pass, high-pass and band-pass filters into a normalized low-pass filter. A canonical or standardized filter is then matched to that normalized low-pass, which is then translated back to the desired form. With the IIR filter we will need to find a mapping that takes our digital frequency response and translates it into an analog filter specification. Basically to map between the z and s planes.
We actually had a mapping when we developed the z transform and it was [latex]z = e^{sT}[/latex], however this is not a one-to-one mapping. Also if we take this equation and solve for s, we get [latex]s = \frac{1}{T} ln(z)[/latex], which if we substitute it into H(s), does not produce a valid H(z). Another approach is to use the Bi-Linear (BL) transformation given by
[latex]s = \frac{z-1}{z+1}[/latex] or [latex]z = \frac{1+s}{1-s}[/latex] (5.3)
The BL transformation is a one-to-one transformation and produces a valid H(z). As a way to visualize this transformation, consider the following points.
If [latex]s = 0 \Rightarrow z = 1[/latex]
and [latex]s = +j \Rightarrow z = \frac{1+j}{1-j} = \frac{(1+j)(1+j)}{(1-j)(1+j)} = \frac{1+2j-1}{2} = +j[/latex]
Similarly [latex]s = -j \Rightarrow z = -j[/latex]
If we let [latex]s \to +\infty \Rightarrow z \to \frac{1+\infty}{1-\infty} \to \frac{+\infty}{-\infty} = -1[/latex]
Similarly [latex]s \to -\infty \Rightarrow z \to \frac{1-\infty}{1+\infty} \to \frac{-\infty}{+\infty} = -1[/latex]
A visual representation of this mapping is shown in Figure 5.3, emphasizing the previous five points. An important point is that the [latex]j \omega[/latex] axis maps onto the unit circle. Both of which represents the FR of the systems.
Figure 5.3 Conformal Mapping for BL Transformation
Figure 5.4 shows the mapping in a slightly different form, showing the mapping of points along the frequency axis.
Figure 5.4 BL Mapping Across Frequency Axes.
With this in mind, let’s walk through a process, by which we could create a discreet feedback system with a given FR.
Section B.1 Bi-Linear Filter Design Procedure
The following steps are how most systems design a discrete-feedback or IIR filter to have a given FR.
Design Procedure
1) Set the sampling frequency, [latex]f_s[/latex] and define the discrete frequency response for the range ([latex]\omega = 0 \to pi[/latex], or [latex]f = 0 \to \frac{f_s}{2}[/latex]).
An example would be to assume a sampling rate of 10 KHz and design a high-pass filter as shown here.
Do note that the range is from 0 to [latex]\frac{f_s}{2}[/latex] = 5 kHz which is all that is needed, since all other frequencies are reflections of the response at these frequencies.
2) The critical frequencies should now be translated to the matching analog frequencies. This translation starts by equating the analog and digital frequencies as shown here.
[latex]s = j \omega_a = {\frac{z - 1}{z+1}}|_{z = e^{j \omega T}} = \frac{e^{j \omega T} - 1} {e^{j \omega T} + 1 }[/latex] (5.4)
where [latex]\omega_a[/latex] is the analog radial frequency ([latex]0 \to \infty[/latex]) and [latex]\omega[/latex] is the digital radial frequency ([latex]0 \to \pi[/latex])
[latex]j \omega_a = \frac{e^{j \omega \frac{T}{2}} (e^{j \omega \frac{T}{2}} - e^{-j \omega \frac{T}{2}})}{e^{j \omega \frac{T}{2}} (e^{j \omega \frac{T}{2}} + e^{-j \omega \frac{T}{2}})} = \frac{2 j sin( \omega \frac{T}{2} )}{ cos( \omega \frac{T}{2} ) } = j tan( \omega \frac{T}{2} )[/latex]
Rarely does anyone use radial frequencies, thus we need to translate to plain frequencies using
[latex]f_s = \frac{1}{T}[/latex] and $latex$ \omega = 2 \pi f$
Resulting in the translation
[latex]j \omega_a = j tan( \frac{\pi f }{f_s} ) [/latex] (5.4)
This equation matches with the frequency mapping from the previous section. With f = 0 [latex]\to \omega = 0[/latex] and f = [latex]\frac{f_s}{2} \to \omega = \infty[/latex].
For the example we are working through the critical frequencies are f = 0, 1 kHz, 3 kHz and 5 kHz. Using the (5.4) we have [latex]\omega_a = 0, tan( \frac{\pi}{10} ) \simeq 0.3249, tan( \frac{3 \pi}{10}) \simeq 1.3764[/latex], and $latex\infty$
Figure 5.5 Frequencies Translated by Bi-Linear Transform
3a) Using the process described in Appendix B, design an analog filter that matches this spec. The first step is to translate this high-pass filter into a Normalized Low Pass (NLP) filter via [latex]\hat{s} = \frac{wp}{s}[/latex], where wp = 1.3764. This translates the design into
Figure 5.6 Frequencies Translated to Normalized Low Pass Frequencies.
Note that the upper cutoff (1.3764) has translated to one, while the lower cutoff (0.3249) moves up to 4.236.
3b) Find canonical form that matches this NLP form. The simplest form is the Butterworth filter, where the order of the filter (n) can be computed by the dB per decade or
[latex]-20 n = \frac{20 log(|H(\omega_s)|}{log(\omega_s)} = \frac{-20}{log(4.236)} = - 31.9[/latex]
Rounding up we have n = 2, and the second order Butterworth filter is relatively simple as
[latex]\hat H(s) = \frac{1}{s^2 + \sqrt{2} s + 1}[/latex]
3c) We can now translate this design back to the original frequencies via [latex]\hat s = \frac{1.3764}{s}[/latex], which looks like
[latex]H(s) = \hat H( \hat s ) |_{\hat s = \frac{1.3764}{s}} = \frac{1}{(\frac{1.3764}{s})^2 + \frac{ 1.3764 \sqrt{2}}{s} + 1}[/latex]
[latex]H(s) = \frac{s^2}{s^2 + 1.9465 s + 1.8944}[/latex]
4) We are now ready to create the z-transform transfer function that will define the discrete system. The translation uses the original definition of the BL transform. Which looks like
[latex]H(z) = H(s) |_{s = \frac{z-1}{z+1}} = \frac{\frac{(z-1)^2}{(z+1)^2}}{\frac{(z-1)^2}{(z+1)^2} + 1.9465 \frac{(z-1)}{(z+1)} + 1.8944}[/latex]
multiplying top and bottom by [latex](z+1)^2[/latex] we have
[latex]H(z) = \frac{z^2 - 2 z + 1}{4.8409 z^2 + 1.788 z + 0.9479}[/latex]
To convert this into a difference equation, we begin my multiplying the top and bottom by 0.20657 and z-2 to give
[latex]H(z) = \frac{0.20657 - 0.4131 z^{-1} + 0.20657 z^{-2}}{1 + 0.36957 z^{-1} + 0.1958 z^{-2}}[/latex]
As a check on the FR, we will factor the numerator into +1 ( zeros ) and the denominator into -0.18476 + j 0.402092 (poles). Visually these look like
Figure 5.7 Pole-Zero Plot of Filter to Implement Called for High Pass Filter.
Which indicates a high-pass filter (zeros at 1, f = 0, and poles over closer to -1).
5) Implementation
This transfer function represents the difference equation
[latex]y(n) = -0.36957 y(n-1) - 0.1958 y(n-2) + 0.20657 x(n) -0.4131 x(n-1) + 0.20657 x(n-2)[/latex]
Or the following signal flow graph.
Figure 5.8 Signal Flow Diagram for High Pass Filter.
Note this unit is the basic building block for most IIR digital filters. The next section will show how these filters can be developed with the MATLAB signal processing toolbox.
Section C) MATLAB Designs
The software environment of MATLAB has a collection of functions that allow the user to create filters without having to go through these calculations by hand.
I will start with the basic call to build the filter we worked through in the previous section. As was described in the previous section, the order of the filter we will need to match the specifications, is important. The order can be found using the MATLAB command “buttord”, which will need the passband-corner and the stopband-corner frequencies as well as the passband attenuation tolerance and the stopband attenuation. The call looks something like
n = buttord( 3/5, 1/5, -3, -20 );
The four arguments that are used in command are
- The passband corner frequency (3 k) divided by [latex]\frac{f_s}{2}[/latex] (5 k).
- The stopband corner frequency (1 k) divided by [latex]\frac{f_s}{2}[/latex] (5 k).
- The passband gain tolerance (-3 dB).
- The stopband attenuation (-20 dB).
Included here is the original specification of the filter, in the digital frequency domain.
Figure 5.9 Specification of Desired High-Pass Filter.
Once we know the order needed, saved in the variable “n” found by the previous function, we can generate the filter in question using the butter command, shown here
[B,A] = butter(n,3/5,’high’)
B =
0.2066 -0.4131 0.2066
A =
1.0000 0.3695 0.1958
The arguments are
- The order of the filter ( 2 )
- The passband cutoff frequency ( 3k / 5k )
- Type of filter ( Optional, if non-given a low-pass filter is assumed )
and the outputs are
- The numerator of H(z) for the filter
- The denominator of H(z)
It should be noted that the numbers generated by the command “butter”, match those computed previously.
Section C1) MATLAB Designs
Much like the analog filters, IIR filters can be generated using the Butterworth or Chebyshev configuration and MATLAB supports these types, with similar commands. The following examples demonstrate how these two types of filters can be created in MATLAB.
It should be noted that high-order filters cannot be accurately represented by polynomials. In the next posting, it will be shown to be a case of round off error. Fortunately, MATLAB has an approach that represents the filters in the form of poles and zeros, which represents the filter perfectly.
Section D) Implementation of IIR Filters
When a filter is moved to a lower platform, such as an embedded system, implementing the filter in a language such as C may be required. This begins with the export of the filters designed previously to a header filter that can be called in a C program.
Figure 5.10 MATLAB Script For the Export of a Set of Filters.
// Define Low-Pass Filter
int LP_N = 6;
float LP_b[6] = { // Numerator
0.00642473F, 0.0064248F, 0.0F,
0.00642473F, 0.0128494F, 0.00642466F};
float LP_a[6] = { // Denominator
1.0F, -0.951955F, 0.0F,
1.0F, -1.94516F, 0.952038F};
float LP_Delays[4] = { 0.0F };//Delays for Low Pass
// Define Band-Pass Filter
int BP_N = 18;
float BP_b[18] = { // Numerator
0.258373F, 0.516747F, 0.258373F,
0.258373F, -0.516747F, 0.258372F,
0.258373F, 0.51732F, 0.258948F,
0.258373F, 0.516172F, 0.2578F,
0.258373F, -0.51732F, 0.258948F,
0.258373F, -0.516173F, 0.257801F};
float BP_a[18] = { // Denominator
1.0F, -1.47463F, 0.725089F,
1.0F, -1.19651F, 0.765093F,
1.0F, -1.75324F, 0.834348F,
1.0F, -1.05038F, 0.908539F,
1.0F, -1.89809F, 0.933409F,
1.0F, -1.95829F, 0.982333F};
float BP_Delays[12] = { 0.0F };//Delays for Band Pass
// Define High-Pass Filter
int HP_N = 9;
float HP_b[9] = { // Numerator
0.227282F, -0.454571F, 0.227286F,
0.227282F, -0.455448F, 0.228169F,
0.227282F, -0.453673F, 0.226394F};
float HP_a[9] = {// Denominator
1.0F, 0.970874F, 0.332497F,
1.0F, 0.40708F, 0.578724F,
1.0F, -0.0211381F, 0.857413F};
float HP_Delays[6] = { 0.0F };//Delays for High Pass
Figure 5.11 Header File Created by Previous MATLAB Script.
In order to provide at least a functional test the header file, shown in Figure 5.11, a simple program that uses the three filters to filter a sweeping cosine, or what some call a chirp, in order to visualize the frequency response of the filters.
// ExampleCode.cpp : Defines the entry point for the console application.
//
#include “stdafx.h” // Silly header created by MS Studio
#include <math.h>
#include <stdio.h>
#include “filters.c” // Defines the filters to be run
// Function that implements a second order difference equation
// using a transposed format.
float SecondOrderfilter( float in, float *b, float *a, float *delays )
{
// Compute output
float out = *(b++) * in + delays[0];
// Service delay line.
delays[0] = *(b++) * in – *(a++) * out + delays[1];
delays[1] = *b * in – *a * out;
// return results.
return out;
} // End of SecondOrderfilter
// Test program that generates a 4096 point chirp and filters it with
// the three filters.
int main(int argc, char* argv[])
{
FILE *fout;
float in,outLP, outBP, outHP;
float scale = 4.0F*(float) atan(1.0) / 4096.0F; // equals ( pi / 4096 )
int n, k,m;
// Open output file and test for valid.
fout = fopen( “Output.Txt”, “w” );
if( fout )
{
// Loop through time index of 0 to 4095
for( n = 0; n < 4096; n++ )
{
// Chirp at index n
in = (float) cos( scale*n*n );
// Loop for implementing Low-Pass Filter.
outLP = in;
m = 0;
for( k = 0; k < LP_N; k+=3 )
{
outLP = SecondOrderfilter( outLP, LP_b+k, LP_a+k+1, LP_Delays+m );
m += 2;
}
// Loop for implementing Band-Pass Filter.
outBP = in;
m = 0;
for( k = 0; k < BP_N; k+=3 )
{
outBP = SecondOrderfilter( outBP, BP_b+k, BP_a+k+1, BP_Delays+m );
m += 2;
} // End of filter implementation loop.
// Loop for implementing High-Pass Filter.
outHP = in;
m = 0;
for( k = 0; k < HP_N; k+=3 )
{
outHP = SecondOrderfilter( outHP, HP_b+k, HP_a+k+1, HP_Delays+m );
m += 2;
} // End of filter implementation loop.
// Write out data for this time index.
fprintf( fout, “%f,%f,%f,%f\n”, in, outLP, outBP, outHP ); // Write out
} // End of time loop.
fclose( fout ); // Close file
} // End of valid open test
return 0;
} // end of main
Figure 5.12 Program in C, to Test Filters.
Figure 5.13 Plot of Input and Filtered Outputs.