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 ztransform:
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 31. This implementation is called direct Form 1 and requires storing the x(n1),…,x(nM) prehistory of the input signal, as well as y(n1), y(n2), y(n3), …, y (nN) – output signal. In total, you need to reserve (N+M) memory cells.
Figure 31 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 32
Figure 32 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 33
Figure 33 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 secondorder section (SOS) blocks.
Thus, we get the following filter chain, as shown in Figure 41.
Figure 41: 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 secondorder block is implemented in one of the ways shown in the previous chapter, see the Figure 31, Figure 32, Figure 33.

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:

Direct calculation of the ztransform 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

In the case is used experience and methods of developing analog filters, i.e. choosing a suitable stransform of the prototype filter, and then convert the splane (Laplace formcontinuous time) to the zplane (z formdiscrete 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 ztransform is calculated:
It can be proved that the poles from the splane in this case to the zplane as follows:
where Ts – sampling time 
Bilinear transform. The bilinear method is the most accurate approximation. This transformation maps the left shalfplane to the unit circle of the zplane.
where Ts – sampling time
Notes:

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

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“


I will use bilinear transform as the most accurate method

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 zTransform:
Figure 61: 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 CCode
% 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 CCode
% Input parameters:
% ParameterVector – data array
% FileNameString – output file name
function IntegerParamVector2file(ParameterVector, FileNameString)
% Float Table to File
% Support CCode
% 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 71: Low Pass Bessel IIR Filter 4th Order. Step Response.
Blue – input signal, red – output signal
Figure 72: Low Pass Bessel IIR Filter 4th Order. 500Hz Response.
Blue – input signal, red – output signal
Figure 73: 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 SOSes 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