IIR Filter and C Implementation Using Octave GNU Tool

IIR Filter and C Implementation Using Octave GNU Tool

1. Abbreviation

FIR – Finite Impulse Response
IIR – Infinite Impulse Response
SOS – second order section
DFT – Discrete Fourier Transformation
IDFT – Inverse Discrete Fourier Transformation
FFT – Fast Fourier Transformation

2. Introduction

The IIR filter is a linear digital filter with an infinite impulse response and the z-transform:

This corresponds to the implementation of the filter:

or

There is still an alternative representation of the IIR filter:

or

Do not confuse these two variants, because different Tools use different views and there are often implementation errors in this regard.

In the future, I will use the representation (1), (2) and (3).

3. IIR Filter in direct Form

This chapter will briefly show the main direct implementations of the IIR filter. I recommend reading more about this in [1], [2], [3].

Using formula (1), you can implement an IIR filter, as shown in Figure 3-1. This implementation is called direct Form 1 and requires storing the x(n-1),…,x(n-M) prehistory of the input signal, as well as y(n-1), y(n-2), y(n-3), …, y (n-N) – output signal. In total, you need to reserve (N+M) memory cells.

                                   Figure 3-1 IIR Filter. Direct Form 1

It turns out that there is a simpler implementation of the IIR filter. Let’s write (1) in another form:

where

(8) in discrete time:

and

(10) in discrete time:

Using (9) and (11) , we can do the canonic implementation of the filter shown in Figure 3-2

                                Figure 3-2 IIR Filter. Direct Form 2 canonic

Using the signal flow graph theory and direct form 2 we can get the transposed direct form 2. See the Figure 3-3

                  Figure 3-3 IIR Filter. Transposed Direct Form 2 canonic

Notes:

• The last two canonical forms require max (M,N) memory locations to store the filter state or prehistory. Usually, N≥ M, so max (M, N)=N

• Of the two canonical forms, I prefer transposed direct form 2, since the SW implementation is simpler in this case, and more accurate in the case of fixed point arithmetic

• When implementing the IIR filter, it is important to monitor the internal computational overflows of the filter, especially in the case of fixed point arithmetic

4. IIR Filter in cascade Form

The direct form of the filter is not used, since there is a high sensitivity of the nulls and poles of the filter from quantization or rounding of the filter coefficients. To prevent this problem, split the z transform into second-order section (SOS) blocks.

Thus, we get the following filter chain, as shown in Figure 4-1.

                         Figure 4-1: IIR Filter. Cascade Realization.

Notes:

• All my examples will use SOS blocks in cascading form. See the iir_filter.c/h

• In this case, each second-order block is implemented in one of the ways shown in the previous chapter, see the Figure 3-1, Figure 3-2, Figure 3-3.

• Octave GNU / Matlab Tool have some functions for the cascade realization support:
  zp2sos(), tf2sos(), sos2zp(), sos2tf(), sosfilt()

5. Design of the IIR Filter

There are two directions in the development of digital filters:

1. Direct calculation of the z-transform coefficients in (1) that meet the requirements of the frequency response and pulse response of the digital filter. I have never used this method and refer you to [1], where you will find a detailed description

2. In the case is used experience and methods of developing analog filters, i.e. choosing a suitable s-transform of the prototype filter, and then convert the s-plane (Laplace form-continuous time) to the z-plane (z form-discrete time)

Let’s take a look at the 2nd direction. There are three main methods for converting analog filters to IIR digital form:

• Backward Euler.
  Approximation of the derivative using a simple formula:

  
  where Ts – sampling time

• Impulse invariance.
  In this case, the pulse response h(t) of the prototype analog filter is sampled h(nTs)=h(n) with the Ts time and then the z-transform is calculated:

  
  It can be proved that the poles from the s-plane in this case to the z-plane as follows:
  
  where Ts – sampling time

• Bilinear transform. The bilinear method is the most accurate approximation. This transformation maps the left s-half-plane to the unit circle of the z-plane.
  
  where Ts – sampling time

Notes:

1. All three methods are nonlinear. This means that we get a digital filter different from the analog prototype

2. The first two methods (Backward Euler and Impulse invariance) are used only in the simplest applications, for example:

   • The implementation of the integrator:

     

   • First order IIR filter. See the article “Simple Low Pass IIR Filter Implementation using the C Language“

3. I will use bilinear transform as the most accurate method

4. Octave GNU / Matlab Tool is supported the bilinear transform with the function: bilinear()

6. Using Example of the IIR Filter

In my demo SW (iir_filter. m/c/h) as an example, I use Bessel IIR filter from the article “Bessel Filter in Digital Form?“: Low Pass Filter, Fband=1kHz, Tsampling=83.33us (Fsampling≈12kHz) with the z-Transform:

           Figure 6-1: Low Pass Bessel IIR Filter 4th Order. Magnitude Response.

7. Octave GNU file IIR_Filter.m

In this demo SW, where I analyze the resulting Bessel IIR filter after quantizing and rounding filter coefficients, generate test signals and pass them through the filter. To do this, the following functions are defined in this file:

% Plotting the Magnitude Response for the 4th order digital filter: (b11+b12*z^-1+b13*z^-2)(b21+b22*z^-1+b23*z^-2)/(a11+a12*z^-1+a13*z^-2)(a21+a22*z^-1+a23*z^-2)
% Input parameters:
% ZB1 – numerator 1 of the second order section
% ZA1 – denominator 1 of the second order section
% ZB2 – numerator 2 of the second order section
% ZA2 – denominator 2 of the second order section
% Tsample – sampling time in seconds
function plot_DigitalMagnitudeResponse(ZB1, ZA1, ZB2, ZA2, Tsample)

% Plotting the Phase Response for the 4th order digital filter: (b11+b12*z^-1+b13*z^-2)(b21+b22*z^-1+b23*z^-2)/(a11+a12*z^-1+a13*z^-2)(a21+a22*z^-1+a23*z^-2)
% Input parameters:
% ZB1 – numerator 1 of the second order section
% ZA1 – denominator 1 of the second order section
% ZB2 – numerator 2 of the second order section
% ZA2 – denominator 2 of the second order section
% Tsample – sampling time in seconds
function plot_DigitalPhaseResponse(ZB1, ZA1, ZB2, ZA2, Tsample)

% Plotting the Signal Response for the 4th order digital filter: (b11+b12*z^-1+b13*z^-2)(b21+b22*z^-1+b23*z^-2)/(a11+a12*z^-1+a13*z^-2)(a21+a22*z^-1+a23*z^-2)
% Input parameters:
% ZB1 – numerator 1 of the second order section
% ZA1 – denominator 1 of the second order section
% ZB2 – numerator 2 of the second order section
% ZA2 – denominator 2 of the second order section
% SampleNumber – Sample number of the step function
function plot_SignalResponseDigitalFilter(ZB1, ZA1, ZB2, ZA2, x_input, SampleNumber)

% IIR Filter Coefficients Quantization
% Input parameters:
% ZA, ZB – IIR coefficients
% ScaleFactor – integer scale
% Output parameters:
% ZA_quantization – integer format
% ZB_quantization – integer format
function [ZA_quantization ZB_quantization] = IirFilterQuantization(ZA, ZB, ScaleFactor)

% Integer Test Signal to File
% Support C-Code
% Input parameters:
% SignalAmpl – amplitude of the test signal
% SignalFreqInHz – frequency of the test signal in Hz
% SignalLength – table length
% Tsample – sampling time in s
% FileNameString – output file name
% Output parameters:
% TestSignal – test signal vector
function TestSignal = IntegerTestSignals2file(SignalAmpl, SignalFreqInHz, SignalLength, Tsample, FileNameString)

% Integer Table to File
% Support C-Code
% Input parameters:
% ParameterVector – data array
% FileNameString – output file name
function IntegerParamVector2file(ParameterVector, FileNameString)

% Float Table to File
% Support C-Code
% Input parameters:
% ParameterVector – data array
% FileNameString – output file name
function FloatParamVector2file(ParameterVector, FileNameString)

Test signals are generated in the IIR_Filter.m:

• Step signal

• Sine curve 500Hz in pass band of the filter

• Sine curve 5kHz in stop band of the filter

• Sum of the two sine curves 500Hz and 5kHz

             Figure 7-1: Low Pass Bessel IIR Filter 4th Order. Step Response.
                                  Blue – input signal, red – output signal

            Figure 7-2: Low Pass Bessel IIR Filter 4th Order. 500Hz Response.
                                 Blue – input signal, red – output signal

       

           Figure 7-3: Low Pass Bessel IIR Filter 4th Order. 500Hz and 5kHz
                                Response. Blue – input signal, red – output signal

8. C Language iir_filter.c/h

The demo SW is realized the cascade SOS-es with the form I / form II / transposed form II. There are two demo SW:

• iir_filter_integer.c/h – realization of the IIR filter with the fixed point arithmetic

• iir_filter_float.c/h – realization of the IIR filter with the float arithmetic

Functions:
– Form I. See the Figure 3-1.
void iir_form_I_initialization(void); – clear the filter history (working arrays)
sint16 iir_form_I(sint16 x_input); – calculate the IIR filter: x_input – input, return value is filter output

– Form II. See the Figure 3-2.
void iir_form_II_initialization(void); – clear the filter history (working arrays)
sint16 iir_form_II(sint16 x_input); – calculate the IIR filter: x_input – input, return value is filter output

– Transposed Form II. See the Figure 3-3.
void iir_transposed_form_II_initialization(void); – clear the filter history (working arrays)
sint16 iir_transposed_form_II(sint16 x_input); – calculate the IIR filter: x_input – input, return value is filter output

Data:
#define IIR_SOS_NMB 2 – number of the SOS-es
iir_coef_integer/_float[ ][ ] – IIR coefficient arrays of the SOS-es

The iir_filter_integer.c/h has a scale factor still:
/* IIR coefficient factor: 2^15=32768 */
#define COEFF_SHIFT_FACTOR 15    – 2^COEFF_SHIFT_FACTOR for the IIR coefficients

9. Download the IIR_Filter.m and iir_filter.c/h

You can download the files:
IIR_Filter.m
iir_filter_integer.c/h
iir_filter_float.c/h

with the button:

10. Literature / References

[1] L.Rabiner, B.Gold “Theory and application of digital signal processing“, Prentice-Hall, Ing Englewood Cliffs, New Jercy 1975

[2] A.Oppenheim, R.Schafer “Discrete-Time Signal Processing“, published by Pearson Education, Inc, publishing as Prentice Hall, 1999

[3] John G. Proakis, Dimitris K. Manolakis, “Digital Signal Processing”, Fourth Edition, Pearson Education Limited 2014, ISBN 10: 1-292-02573-5, ISBN 13: 978-1-292-02573-5