Discrete-time (digital) IIR filter, Matlab example

Summary

Matlab example program that designs a discrete-time IIR filter and plots its frequency response.

Example program

The program generates an example filter, calculates the frequency response on a logarithmic axis, and plots it (figure 1, 2).
It also plots

for frequencies far beyond the sampling frequency, showing that the frequency response is periodic with the sampling frequency (figure 3).

Figure 1: Amplitude response

Figure 2: Phase response

Figure 3: Periodic frequency response

Download

% ***********************************************************************
% Discrete-time filter example
% © Markus Nentwig 2007
% This program is provided without any warranty.
% ***********************************************************************
close all; clear all;

% Filter design parameters
order=5;
fCutoff=2000;

fSample=44100;
fCutoff=fCutoff/(fSample/2);

% ***********************************************************************
% Design a filter
% ***********************************************************************
if 0
  [b,a]=butter(order, fCutoff);
elseif 1
  ripple_dB=0.5;
  [b,a]=cheby1(order, ripple_dB, fCutoff);
elseif 0
  ripple_dB=0.5;
  stopband_dB=60;
  [b,a]=ellip(order, ripple_dB, stopband_dB, fCutoff);
end

% ***********************************************************************
% Build a logarithmic frequency axis
% ***********************************************************************
f1=10;
f2=fSample/2;
nPoints=3000;
fbase=exp(linspace(log(f1), log(f2), nPoints));
fbase=fbase/(fSample);

% ***********************************************************************
% Numerator and denominator are polynomials in z=exp(i*omega)
% ***********************************************************************
H_num=polyval(b, exp(i*2*pi*fbase));
H_denom=polyval(a, exp(i*2*pi*fbase));
% H is the complex gain of the filter at any given frequency
H=H_num ./ H_denom;

% ***********************************************************************
% - Convert to dB. H could be a voltage or current => use 20 log |H|
% - Obtain phase in degrees
% - plot
% ***********************************************************************
gain_dB=20*log10(abs(H)+1e-15); % don't show below -300 dB
phase_deg=unwrap(arg(H))/pi*180;
figure(); 
semilogx(fbase*fSample, gain_dB, '+'); xlabel('f/Hz'); ylabel('|H|/dB'); grid on;
title('amplitude response');

figure();
semilogx(fbase*fSample, phase_deg, '+'); xlabel('f/Hz'); ylabel('phase/deg'); grid on;
title('phase response');

% ***********************************************************************
% Build a linear frequency axis
% ***********************************************************************
fbase=linspace(0, 4*fSample, nPoints);
fbase=fbase/fSample;

% ***********************************************************************
% Evaluate frequency response as before
% ***********************************************************************
H_num=polyval(b, exp(i*2*pi*fbase));
H_denom=polyval(a, exp(i*2*pi*fbase));
H=H_num ./ H_denom;

% ***********************************************************************
% Plot as before, but on linear frequency axis
% ***********************************************************************
gain_dB=20*log10(abs(H)+1e-15); % don't show below -300 dB
figure(); plot(fbase*fSample, gain_dB, '+'); xlabel('f/Hz'); ylabel('|H|/dB'); grid on;
title('amplitude response');

© Markus Nentwig 2007-2008
The content of this page is provided without any warranty and may not be reproduced without permission.

Comments? Questions?

Please send me a mail! mnentwig@elisanet.fi