1 import pylab as pl

     2 

     3 def freqresp(H, omega, sys_type='c'):

     4     if (omega == None):

     5         omega = pl.logspace(-2, 2);

     6 

     7     num = pl.poly1d(H.num)

     8     den = pl.poly1d(H.den)

     9     if sys_type == 'd':

    10         fresp = map(lambda w: num(w) / den(w), pl.exp(2*pl.pi*omega*1j))

    11     else:

    12         fresp = map(lambda w: num(w) / den(w), 2*pl.pi*omega*1j)

    13     fresp = pl.array(fresp)

    14     return omega, fresp

    15 

    16 def phasemag(fresp):

    17     mag = abs(fresp)

    18     mag = 20*pl.log10(mag)

    19 

    20     phase = pl.angle(fresp)*180/pl.pi

    21 

    22     return mag, phase

    23     

    24 def bode(sys, omega=None):

    25     omega, fresp  = freqresp(sys, omega, 'd')

    26     mag, phase = phasemag(fresp)

    27 

    28     pl.subplot(211)

    29     pl.grid(linestyle='--')

    30     pl.semilogx(omega, mag)

    31     xmin, xmax = pl.xlim()

    32     pl.xlim(1e-3, xmax)

    33 

    34     pl.subplot(212)

    35     pl.grid(linestyle='--')

    36     pl.semilogx(omega, phase)

    37     xmin, xmax = pl.xlim()

    38     pl.xlim(1e-3, xmax)

    39 

    40     return (211, 212)

    41