--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/lqg/spec.py	Fri May 27 14:24:59 2011 +0530
@@ -0,0 +1,38 @@
+#!/usr/bin/python
+# 13.1
+
+import os, sys
+sys.path += [os.getcwdu() + os.sep + ".." + os.sep + "python"]
+
+import scipy as sp
+from polyfuncs import polmul, poladd
+
+def spec(A, dA, B, dB, rho):
+    AA = rho * sp.convolve(A, A[::-1])
+    BB = sp.convolve(B, B[::-1])
+    diff1 = dA - dB
+    dBB = 2 * dB
+    for i in range(diff1):
+        BB, dBB = polmul(BB, dBB, sp.array([0, 1]), 1)
+    rbbr, drbbr = poladd(AA, 2*dA, BB, dBB)
+    rbbr = rbbr.squeeze()
+    rts = sp.roots(rbbr)
+    rtsin = rts[dA:2*dA+1]
+    b = 1
+    for i in range(dA):
+        b = sp.convolve(b, sp.array([1, -rtsin[i]]))
+    bbr = sp.convolve(b, b[::-1])
+    r = rbbr[0] / bbr[0]
+    return r,b,rbbr
+
+if __name__ == "__main__":
+    A = sp.convolve([-0.5, 1], [-0.9, 1])
+    dA = 2
+    B = 0.5*sp.array([-0.9, 1])
+    dB = 1
+    rho = 1
+
+    r, beta1, sigma = spec(A, dA, B, dB, rho)
+    print "sigma", sigma
+    print "beta1", beta1
+    print "r", r