1 import pylab as pl

     2 def p_margin(H, type='d', Ts=1):

     3     // |h(e^(i*w*dt))|=1 <-- h(e^(i*w*dt))*h(e^(-i*w*dt))

     4     z=poly(0,varn(h.den));

     5     den = H.den

     6     sm=simp_mode();simp_mode(%f);hh=h*horner(h,1/z)-1;simp_mode(sm)

     7     //find the numerator roots

     8     z=roots(hh.num);

     9     z(abs(abs(z)-1)>eps)=[];// retain only roots with modulus equal to 1

    10     w=log(z)/(%i*dt);

    11     ws=real(w(abs(imag(w))<eps&real(w)>0)); //frequency points with unitary modulus

    12     if ws==[] then phm=%inf,fr=[],return,end

    13     f=horner(h,exp(%i*ws*dt));

    14   end

    15     return phm, fr

    16 

    17 def g_margin(H, type='d', Ts=1):

    18     

    19     return gm, fr