Ville Jalkanen
Matlab examples
1 Filter description
A filter is appropriately described by the transfer function. It is a ratio between two polynomials
𝑁(𝑠) 𝑏𝑛 𝑠 𝑛 + 𝑏𝑛−1 𝑠 𝑛−1 + ⋯ + 𝑏0 𝑠 0
𝐻(𝑠) = = .
𝐷(𝑠) 𝑎𝑚 𝑠 𝑚 + 𝑎𝑚−1 𝑠 𝑚−1 + ⋯ + 𝑎0 𝑠 0
In Matlab the filter (or system) is described by setting the coefficents in the numerator and
denominator polynomials. The coefficients must be set in the correct order. The following example
demonstrates this for an analog second order high pass filter
𝑠2
𝐻(𝑠) = .
𝑠 2 + √2 ∙ 𝑠 + 1
Clear all The transfer function appears in the command
num = [1 0 0] window:
den = [1 sqrt(2) 1]
tf(num,den) Transfer function:
s^2
-----------------
s^2 + 1.414 s + 1
2 Poles and zeros
If the transfer function is known, then it is possible to find the poles and zeros. In this case there are
two zeros in origo and two complex conjugated poles.
Clear all The pole-zero diagram becomes:
num = [1 0 0]
den = [1 sqrt(2) 1]
[z,p,k] = tf2zp(num,den)
zplane(z,p) 1
The result in the command
window: 0.5
Imaginary Part
z =
0 2
0
0
p =
-0.5
-0.7071 + 0.7071i
-0.7071 - 0.7071i
k = -1
1 -1 -0.5 0 0.5 1
Real Part
1
Ville Jalkanen
If the poles and zeros are known, then it is possible to find the transfer function. Assume that a
system has two zeros, one at -2 and one in the origin, and two poles at -2 ± i.
clear all The pole-zero diagram
k = 1;
z = [-2; 0];
p = [-2+i;-2-i]; 1.5
zplane(z,p)
axis([-3 1 -1.5 1.5])
[num,den] = zp2tf(z,p,k) 1
tf(num,den)
Imaginary Part
0.5
The result in the command
window: 0
num = -0.5
1 2 0
-1
den =
1 4 5
-1.5
Transfer function: -3 -2 -1 0 1
s^2 + 2 s Real Part
-------------
s^2 + 4 s + 5
3 The transfer function in a diagram - The frequency response
The magnitude and phase of the frequency response can be obtained by the following commands.
clear all;
num = [1 2 0]; 0
10
den = [1 4 5];
freqs(num,den)
-1
Magnitude
10
Another way to do this:
-2
H = tf(num,den); 10
bode(H)
-3
10
Note that with the last 10
-2
10
-1
10
0 1
10
command, we obtain the Frequency (rad/s)
magnitude in dB.
100
Phase (degrees)
50
0
-2 -1 0 1
10 10 10 10
Frequency (rad/s)
2
Ville Jalkanen
4 The magnitude response in decibel and different plot-commands
We can also obtain the magnitude in dB by using the freqs-function and controlling the appearance
of the graph.
clear all; Result from semilogx-command:
num = [1 2 0];
den = [1 4 5]; 0
[h,w] = freqs(num,den)
semilogx(w,20*log10(abs(h)))
-10
ylabel('Magnitude (dB)')
xlabel('Frequency (rad/s)')
Magnitude (dB)
-20
Try the following commands in separate
figures (use figure to create a new -30
figure-window):
-40
figure
plot(w,abs(h))
figure -50
-2 -1 0 1
plot(log10(w),20*log10(abs(h))) 10 10 10 10
figure Frequency (rad/s)
loglog(w,abs(h))
figure
semilogy(log10(w),abs(h))
We can also control the frequency axis by
creating frequency vectors to be used in
freqs:
w = linspace(0,100)
or
w = logspace(-2,2)
freqs(num,den,w)
5 The phase in degrees
We can also create a phase diagram using a new frequency vector.
clear all; The resulting figure:
num = [1 2 0];
den = [1 4 5];
w = logspace(-2,2,100)
h = freqs(num,den,w)
semilogx(w,180/pi*angle(h),'.')
ylabel('Phase (degrees)')
xlabel('Frequency (rad/s)')
3
Ville Jalkanen
100
80
Phase (degrees)
60
40
20
0
-2 0 2
10 10 10
Frequency (rad/s)
6 Numerical values
In the command window numerical values can be shown in columns instead of being arranged as
rows.
clear all; The result in the command window:
num = [0 0 1]; q =
den = [1 sqrt(2) 1];
w = logspace(-1,1,10); 0.1000 -0.0004
h = freqs(num,den,w); 0.1668 -0.0034
h_dB = 20*log10(abs(h)); 0.2783 -0.0260
q = [w' h_dB'] 0.4642 -0.1970
0.7743 -1.3334
For a linear frequency scale: 1.2915 -5.7779
2.1544 -13.5304
w = 1:1:10; 3.5938 -22.2482
h = freqs(num,den,w); 5.9948 -31.1145
h_dB = 20*log10(abs(h)); 10.0000 -40.0004
q = [w' h_dB']
For the linear frequency vector:
q =
1.0000 -3.0103
2.0000 -12.3045
3.0000 -19.1381
4.0000 -24.0993
5.0000 -27.9657
6.0000 -31.1294
7.0000 -33.8057
8.0000 -36.1247
9.0000 -38.1704
10.0000 -40.0004
4
Ville Jalkanen
7 Fifth order Butterworth filter
clear all;
[z,p,k] = buttap(5);
figure(1);
zplane(z,p)
title('5:th order Butterworth');
figure(2);
[num,den]=zp2tf(z,p,k);
w=logspace(-1,1,500);
h=freqs(num,den,w);
semilogx(w,20.*log10(abs(h)));
title('5:th order Butterworth');
axis([.1,10,-80,+5]);
8 Fifth order Chebyshev filter, type 1, 3 dB ripple in the passband
clear all
[z,p,k]=cheb1ap(5,3);
figure(1)
zplane(z,p);
title('5:th order Chebyshev-1');
figure(2)
[num,den]=zp2tf(z,p,k);
w=logspace(-1,1,500);
h=freqs(num,den,w);
semilogx(w,20.*log10(abs(h)));
title('5:th order Chebyshev-1');
axis([.1,10,-80,+5]);
9 Fifth order Chebyshev filter, type 2, 40 dB ripple in the stopband
[z,p,k]=cheb2ap(5,40);
figure(3)
zplane(z,p);
title('5:th order Chebyshev-2');
figure(4)
[num,den]=zp2tf(z,p,k);
w=logspace(-1,1,500);
h=freqs(num,den,w);
semilogx(w,20.*log10(abs(h)));
title('5:th order Chebyshev-2');
axis([.1,10,-80,+5]);
10 Fifth order Bessel filter
[z,p,k]=besselap(5);
figure(1)
zplane(z,p);
title('5:th order Bessel filter');
figure(2)
[num,den]=zp2tf(z,p,k);
w=logspace(-1,1,500);
h=freqs(num,den,w);
semilogx(w,20.*log10(abs(h)));
5
Ville Jalkanen
title('5:th order Bessel filter');
axis([.1,10,-80,+5]);
11 Fifth order Cauer/Elliptic filter
[z,p,k]=ellipap(5,3,40);
figure(1)
zplane(z,p);
title('5:th order Cauer filter');
figure(2)
[num,den]=zp2tf(z,p,k);
w=logspace(-1,1,500);
h=freqs(num,den,w);
semilogx(w,20.*log10(abs(h)));
title('5:th order Cauer filter');
axis([.1,10,-80,+5]);
12 Step response
clear all;
[z,p,k]=besselap(5);
[num,den]=zp2tf(z,p,k);
dt=0.01;
t=0:dt:30;
h=step(tf(num,den),t);
plot(t,h,'b');
hold on
[z,p,k]=cheb1ap(5,3);
[num,den]=zp2tf(z,p,k);
dt=0.01;
t=0:dt:30;
h=step(tf(num,den),t);
plot(t,h,'r');
hold off
title('red=chebyshev, blue=bessel');
13 Impulse response
clear all;
[z,p,k]=besselap(5);
[num,den]=zp2tf(z,p,k);
dt=0.01;
t=0:dt:30;
h=impulse(tf(num,den),t);
plot(t,h,'b');
hold on
[z,p,k]=cheb1ap(5,3);
[num,den]=zp2tf(z,p,k);
dt=0.01;
t=0:dt:30;
h=impulse(tf(num,den),t);
plot(t,h,'r');
hold off
title('red=chebyshev, blue=bessel');
6
Ville Jalkanen
14 Output signal from the system
clear all; A short duration sinusoid with frequency 1
dt=0.01;
t=0:dt:30;
rad/s is filtered by Chebyshev and Bessel
uin=[zeros(1,100),sin(t),zeros(1,100)]; filter respectively.
[zc,pc,kc]=cheb1ap(5,3);
[numc,denc]=zp2tf(zc,pc,kc);
hc=impulse(numc,denc,t);
uutc=conv(hc,uin).*dt;
[zb,pb,kb]=besselap(5);
[numb,denb]=zp2tf(zb,pb,kb);
hb=impulse(numb,denb,t);
uutb=conv(hb,uin).*dt;
T=0:dt:(length(uutb)-1.)*dt;
UIN=[zeros(1,100),
sin(t),zeros(1,length(T)-
length(uin)+100)];
plot(T,uutb,'r',T, uutc,'b', T,UIN,'k')
axis([0,60,-1.2,1.2])
15 Low pass-to-low pass (not using lp2lp)
A system has the poles -1 ± 3j. We wish to double the cutoff frequency. The Matlab code to do this is:
clear all;
z0=[]; 5
p0=[-1+3.*j;-1-3.*j];
k0=p0(1).*p0(2);
[num0,den0]=zp2tf(z0,p0,k0); 0
dw=0.1;
w=0:dw:10; -5
H0=freqs(num0,den0,w);
semilogx(w,20.*log10(abs(H0)),'r');
hold on; -10
%transformed system is marked by 1
z1=[];
p1=2.*p0; -15
k1=2.^2.*k0;
[num1,den1]=zp2tf(z1,p1,k1);
H1=freqs(num1,den1,w); -20
-1 0 1
semilogx(w,20.*log10(abs(H1)),'b'); 10 10 10
hold off
16 Low pass-to-high pass
Determine the poles and zeros for a 5:th order high pass Chebyshev type 1 filter with 3 dB ripple and
a cutoff frequency of 1 kHz.
7
Ville Jalkanen
clear all; magnitude of H as a function of frequency
%original system is marked with 0 5
[z0, p0,k0]=cheb1ap(5,3);
[num0,den0]=zp2tf(z0,p0,k0); 0
%We move the cutoff frequency to
% 2*pi*1000 and transform to HP.
% Transformed filter is marked by 1 -5
[num1,den1]=lp2hp(num0,den0,2.*pi.*1000);
dw=100; -10
w=1000:dw:30000;
H1=freqs(num1,den1,w); -15
figure(1)
plot(w./2./pi,20.*log10(abs(H1)),'b');
hold off -20
0 1000 2000 3000 4000 5000
axis([0,5000,-20,+5])
title('magnitude of H as a function of 4
frequency') x 10
[z1,p1,k1]=tf2zp(num1,den1)
figure(2) 1
zplane(z1,p1)
Imaginary Part
In the command window the zeros and poles are: 0
z1 =
-0.0248
-1
0.0230
0.0018
-0.0000 + 0.0004i -3 -2 -1 0
-0.0000 - 0.0004i Real Part 4
x 10
p1 =
1.0e+004 *
-3.5392
-0.2394 + 0.9949i
-0.2394 - 0.9949i
-0.0368 + 0.6484i
-0.0368 - 0.6484i
17 Low pass-to-band pass
Starting from a 5:th order Chebyshev (type 1) filter with 3 dB ripple, we wish to transform it to a
band pass filter with cutoff frequencies 800 Hz and 1200 Hz, which gives a bandwidth of 400 Hz and a
center frequency approximately 1 kHz (980 Hz). We determine the pole and zeros of the system.
clear all; magnitude of H as a function of frequency
%original system is marked with 0 5
[z0, p0,k0]=cheb1ap(5,3);
[num0,den0]=zp2tf(z0,p0,k0); 0
% Transformed filter is marked by 1
[num1,den1]=lp2bp(num0,den0,2.*pi.*980,2.*pi.* -5
400);
dw=10; -10
w=1000:dw:30000;
-15
H1=freqs(num1,den1,w);
figure(1)
-20
plot(w./2./pi,20.*log10(abs(H1)),'b'); 0 500 1000 1500 2000
hold off
axis([0,2000,-20,+5])
8
Ville Jalkanen
title('magnitude of H as a function of
frequency')
[z1,p1,k1]=tf2zp(num1,den1) 5000
Imaginary Part
figure(2)
zplane(z1,p1)
0
In the command window the zeros and poles are:
-5000
z1 =
4.1184
-1 0 1
1.2609 + 3.9130i
1.2609 - 3.9130i Real Part 4
x 10
-3.3201 + 2.4045i
-3.3201 - 2.4045i
p1 =
1.0e+003 *
-0.0823 + 7.4895i
-0.0823 - 7.4895i
-0.2023 + 6.9506i
-0.2023 - 6.9506i
-0.2231 + 6.1535i
-0.2231 - 6.1535i
-0.1586 + 5.4503i
-0.1586 - 5.4503i
-0.0556 + 5.0618i
-0.0556 - 5.0618i
