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

                            < M A T L A B (R) >
                (c) Copyright 1984-94 The MathWorks, Inc.
                            All Rights Reserved
                               Version 4.2c
                                Dec 31 1994

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 5.1.0.421 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 analogbut2.eps
>> [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 5.1.0.421 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  (5.1.0.421) version
				of matlab, they now seem to be 
				premultiplying by Ts, and so  				
				 You should get this plot 
>> print impinvar.eps 

>>[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  5.1.0.421 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 bilin.eps
>> [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 bilinprewarp.eps

>> 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]