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)