import pylab as pl

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

    if (omega == None):

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

    num = pl.poly1d(H.num)

    den = pl.poly1d(H.den)

    if sys_type == 'd':

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

    else:

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

    fresp = pl.array(fresp)

    return omega, fresp

def phasemag(fresp):

    mag = abs(fresp)

    mag = 20*pl.log10(mag)

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

    return mag, phase

def bode(sys, omega=None):

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

    mag, phase = phasemag(fresp)

    pl.subplot(211)

    pl.grid(linestyle='--')

    pl.semilogx(omega, mag)

    xmin, xmax = pl.xlim()

    pl.xlim(1e-3, xmax)

    pl.subplot(212)

    pl.grid(linestyle='--')

    pl.semilogx(omega, phase)

    xmin, xmax = pl.xlim()

    pl.xlim(1e-3, xmax)

    return (211, 212)