Matlab filter design example
Example session
Design a 1000 Hz 2nd order Butterworth lowpass filter at a 4000 Hz sample rate using all of the
above methods.
Note: type "help command" to get help on matlab commands, i.e., "help butter."
New:
You may also wish to experiment with the new Matlab "Filter design and analysis tool." To run
this, type "fdatool" at the Matlab prompt. Use Analysis->FullView to rescale Y-axis.
matlab
Commands to get started: intro, demo, help help
Commands for more information: help, whatsnew, info, subscribe
>> [z,p,k]=butter(2,6283,'s')
z = (NOTE: no finite zeroes in butterworth)
[]
p = (NOTE:this is location of the 2 poles in s-
plane)
1.0e+03 *
-4.4428 + 4.4428i
-4.4428 - 4.4428i
k =
39476089 (NOTE: this is a constant multiplying the
result)
(==========================)
(<<< BEWARE MATLAB BUG!!!)
(==========================)
Warning: There is an error
in the [Link] version
of matlab. There should be no zeroes
in
the returned value of z. To fix this,
you must suppress the zeroes by
entering the command "z=[]" after
In the current version you will get
z=[ 0 0 0.2500 ] which is wrong!
For the new (broken) version, type the following to correct the
problem:
>> z=[]
z =
[]
>> [n,d]=zp2tf(z,p,k)
n = (NOTE: numerator is this constant )
0 0 39476089
d =
1.0e+07 * (==========================)
0.0000 0.0009 3.9476 (<<< BEWARE MATLAB BUG!!!)
(==========================)
Warning: There is a bug in some versions
of matlab. Although it looks like
the values are zero, they are not.
The print format is suppressing
the significant digits. To fix this
use:
>> h=tf(n,d)
Transfer function:
3.948e07
-----------------------
s^2 + 8886 s + 3.948e07
>> pzmap(h)
>> format short e
(Now you can see the denominator.)
>> d
d = (NOTE: denominator is s^2 + 8886 s +
3.9476e+07)
1.0000e+00 8.8855e+03 3.9476e+07
For high order filters, you need even greater
numerical precision to avoid errors. So,
in this case use:
>> format long e
>> d
d =
1.000000000000000e+00 8.885503812390156e+03
3.947608900000000e+07
>> bode(n,d) You should get this plot
>> print [Link]
>> [nd,dd]=impinvar(n,d,4000)
nd = (NOTE: numerator = 0 + 2622 z^-1 )
(NOTE: matlab does NOT premultiply by Ts)
1.0e+03 *
0 2.6220
dd = (NOTE: denominator = 1 -0.29 z^-1 + 0.10 z^-2)
1.0000 -0.2925 0.1085
Warning: In the [Link] version
of matlab, they now seem to be
premultiplying by Ts, and so the
result now is:
>> [nd,dd]=impinvar(n,d,4000)
nd =
0 6.5549e-01 0
dd =
1.0000e+00 -2.9248e-01 1.0846e-01
>> freqz(nd,dd,256) You should get this plot
Warning: In the current ([Link]) version
of matlab, they now seem to be
premultiplying by Ts, and so
You should get this plot
>> print [Link]
>>[top,bot]=residuez(nd,dd) (NOTE: use this to find partial fraction form)
top =
0- 4.4428e+03i
0+ 4.4428e+03i
bot =
1.4624e-01+ 2.9508e-01i
1.4624e-01- 2.9508e-01i
(NOTE: Then partial frac expansion is:
-4442 i + 4442 i
--------------------- ---------------------
1- (.146 + .29 i)z^-1 1- (.146 - .29 i)z^-1
Warning: In the [Link] version
of matlab, (will they ever get
their act together)
they now seem to be
premultiplying by Ts, and so the
result now is:
>> [top,bot]=residuez(nd,dd)
top =
0 + 1.1107e+00i
0 - 1.1107e+00i
bot =
1.4624e-01 - 2.9508e-01i
1.4624e-01 + 2.9508e-01i
>> [nd,dd]=bilinear(n,d,4000)
nd =
0.2261 0.4523 0.2261
dd =
1.0000 -0.2810 0.1856
>> freqz(nd,dd,256) You should get this plot
>> print [Link]
>> [nd,dd]=bilinear(n,d,4000,1000)
nd =
0.2929 0.5858 0.2929
dd =
1.0000 -0.0000 0.1716
>> freqz(nd,dd,256) You should get this plot
>> print [Link]
>> h=tf(nd,dd)
Transfer function:
0.2929 s^2 + 0.5858 s + 0.2929
------------------------------
s^2 - 3.455e-05 s + 0.1716
>> pzmap(h)
>> zgrid
>> zgrid([],[])
>> quit
127007 flops.
You may also want to try:<
[ppp,zzz]=pzmap[nd,dd]
pzmap[nd,dd]
