--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/python/margin.py	Fri May 27 14:24:59 2011 +0530
@@ -0,0 +1,19 @@
+import pylab as pl
+def p_margin(H, type='d', Ts=1):
+    // |h(e^(i*w*dt))|=1 <-- h(e^(i*w*dt))*h(e^(-i*w*dt))
+    z=poly(0,varn(h.den));
+    den = H.den
+    sm=simp_mode();simp_mode(%f);hh=h*horner(h,1/z)-1;simp_mode(sm)
+    //find the numerator roots
+    z=roots(hh.num);
+    z(abs(abs(z)-1)>eps)=[];// retain only roots with modulus equal to 1
+    w=log(z)/(%i*dt);
+    ws=real(w(abs(imag(w))<eps&real(w)>0)); //frequency points with unitary modulus
+    if ws==[] then phm=%inf,fr=[],return,end
+    f=horner(h,exp(%i*ws*dt));
+  end
+    return phm, fr
+
+def g_margin(H, type='d', Ts=1):
+    
+    return gm, fr