IIR filter without Matlab inbuilt function

Digital signal processing is an essential part of modern communication and audio processing systems. One important aspect of DSP is the use of digital filters to remove unwanted frequencies from signals. Here, we will discuss the implementation of an Infinite Impulse Response (IIR) filter using MATLAB without using built-in functions.

An IIR filter is a type of digital filter that uses feedback to achieve the desired frequency response. It is a recursive filter that uses the current and past input samples as well as past output samples to calculate the current output sample. The IIR filter is given by the following equation called difference equation:

\( y(n)=-\sum_{k=0}^{N-1}a_{k} y(n-k) + \sum_{k=0}^{M-1}b_{k} x(n-k) \) ---->(1)

where,  \(b_0,b_1,b_2...b_N\) and, \(a_0,a_1,a_2...a_N\) are the filter coefficients and N is the filter order.

The equation (1) can be rewritten in terms of filter transfer equation which are of three types:Z-transform, DTFT and DTF.

The Z-Transform of IIR filter is,

\( H(z) = \frac{\sum_{k=0}^{M}b_k z^{-k}}{\sum_{k=0}^{N}a_k z^{-k}} = \frac{b_0 + b_1 z^{-1} + b_2 z^{-2} + ... + b_n z^{-n}}{a_0 + a_1 z^{-1} + a_2 z^{-2} + ... + a_n z^{-n}}\) ---->(2)

See derivation of Z-transform equation.

The DTFT of IIR filter is,

\( H(e^{jw}) = \sum_{0}^{\infty} h(n) e^{-jwn} \) ---->(3)

The DFT of IIR filter is given by,

\( H(k) = \sum_{0}^{N-1} x(n) e^{-j \frac{2 \pi}{N}kn} \) ---->(4)

We can use any of the above four equation to solve IIR filter problem depending upon the given requirements. Here it is shown how to use the feedforward and feedback coefficients and using a loop to calculate the output samples of IIR filter without using Matlab inbuilt functions. That is the difference equation (1) above i used.

Code Explanation

The following code implements an IIR filter in MATLAB without using built-in functions:

clear

% Define the filter coefficients
% LPF with fc=1kHz, fs=44.1kHz, Butterworth
b = [2.15e-05,8.60e-05,1.29e-04,8.60e-05,2.15e-05];
a = [1,-3.62,4.95,-3.01,0.68]; % Denominator coefficients

% Define the input signal
N = 100; % Number of samples
fs = 44100; % Sampling frequency (Hz)
f1 = 500; % Frequency of the first sine wave (Hz)
f2 = 2000; % Frequency of the second sine wave (Hz)
f3 = 3000; % Frequency of the third sine wave (Hz)
t = 0:1/fs:(N-1)/fs; % Time vector
x = sin(2*pi*f1*t) + sin(2*pi*f2*t) + sin(2*pi*f3*t); % Composite sine signal

% Initialize the output signal
y = zeros(1, length(x));

% Apply the filter to the input signal
for i = 5:length(x)
    y(i) = b(1)*x(i) + b(2)*x(i-1) + b(3)*x(i-2) + b(4)*x(i-3) + b(5)*x(i-4) ...
           - a(2)*y(i-1) - a(3)*y(i-2) - a(4)*y(i-3) - a(5)*y(i-4);
end

% Plot the input and output signals
figure(1)
clf
n = 1:length(x);
subplot(2,1,1)
plot(n, x)
title('Input signal(x)')
xlabel('Sample')
ylabel('Amplitude')
subplot(2,1,2)
plot(n, y)
title('Output signal(y)')
xlabel('Sample')
ylabel('Amplitude')

% Compute the frequency response
[h, w] = freqz(b, a);

% Convert from normalized frequency to Hz for axis
f = w/pi*fs/2; 

% Plot the magnitude and phase responses
figure(2)
clf
subplot(2,1,1);
plot(f, 20*log10(abs(h)));
xlabel('Frequency (Hz)');
ylabel('Magnitude (dB)');
title('Magnitude Response');

subplot(2,1,2);
plot(f, unwrap(angle(h)));
xlabel('Frequency (Hz)');
ylabel('Phase (rad)');
title('Phase Response');

This code implements a Butterworth low-pass filter with a cutoff frequency of 1 kHz, and applies it to a composite sine signal composed of three different sine waves of frequencies 500 Hz, 2000 Hz, and 3000 Hz.

The filter coefficients for the Butterworth filter are defined in the following lines:

b = [2.15e-05,8.60e-05,1.29e-04,8.60e-05,2.15e-05];
a = [1,-3.62,4.95,-3.01,0.68];

See how to calculate IIR filter coefficients in Matlab. These coefficients are used in the difference equation of the filter, which is implemented in the loop:

for i = 5:length(x)
    y(i) = b(1)*x(i) + b(2)*x(i-1) + b(3)*x(i-2) + b(4)*x(i-3) + b(5)*x(i-4) ...
           - a(2)*y(i-1) - a(3)*y(i-2) - a(4)*y(i-3) - a(5)*y(i-4);
end

This loop iterates over each sample in the input signal x, and applies the filter to the current sample and the previous four samples to produce the output sample y(i). The output signal is initialized as a vector of zeros of the same length as the input signal:

y = zeros(1, length(x));

The input signal x is defined as a composite sine wave signal, sampled at 44.1 kHz, with three different frequencies:

N = 100; % Number of samples
fs = 44100; % Sampling frequency (Hz)
f1 = 500; % Frequency of the first sine wave (Hz)
f2 = 2000; % Frequency of the second sine wave (Hz)
f3 = 3000; % Frequency of the third sine wave (Hz)
t = 0:1/fs:(N-1)/fs; % Time vector
x = sin(2*pi*f1*t) + sin(2*pi*f2*t) + sin(2*pi*f3*t); % Composite sine signal

The input and output signals are then plotted using the plot function, and the frequency response of the filter is computed using the freqz function(see difference between filter() and freqz() matlab functions).

figure(1)
clf
n = 1:length(x);
subplot(2,1,1)
plot(n, x)
title('Input signal(x)')
xlabel('Sample')
ylabel('Amplitude')
subplot(2,1,2)
plot(n, y)
title('Output signal(y)')
xlabel('Sample')
ylabel('Amplitude')

The frequency response is then plotted using the plot function, with the magnitude and phase responses displayed in separate subplots.

The freqz function computes the frequency response of the IIR filter, which is returned as a complex vector h. The frequency vector w is also returned, which represents the normalized frequency ranging from 0 to π. The freqz function is called with the filter coefficients b and a:

[h, w] = freqz(b, a);

The frequency vector w is then converted from normalized frequency to Hz for plotting on the x-axis of the frequency response plots:

f = w/pi*fs/2;

The magnitude and phase responses of the filter are plotted using the plot function, with the unwrap function used to unwrap the phase response:

figure(2)
clf
subplot(2,1,1);
plot(f, 20*log10(abs(h)));
xlabel('Frequency (Hz)');
ylabel('Magnitude (dB)');
title('Magnitude Response');

subplot(2,1,2);
plot(f, unwrap(angle(h)));
xlabel('Frequency (Hz)');
ylabel('Phase (rad)');
title('Phase Response');

The first subplot displays the magnitude response of the filter in dB, and the second subplot displays the phase response of the filter in radians. The unwrap function is used to remove any phase jumps that may occur due to the periodic nature of the phase response.