1 #!/usr/bin/python

     2 # 13.1

     3 

     4 import os, sys

     5 sys.path += [os.getcwdu() + os.sep + ".." + os.sep + "python"]

     6 

     7 import scipy as sp

     8 from polyfuncs import polmul, poladd

     9 

    10 def spec(A, dA, B, dB, rho):

    11     AA = rho * sp.convolve(A, A[::-1])

    12     BB = sp.convolve(B, B[::-1])

    13     diff1 = dA - dB

    14     dBB = 2 * dB

    15     for i in range(diff1):

    16         BB, dBB = polmul(BB, dBB, sp.array([0, 1]), 1)

    17     rbbr, drbbr = poladd(AA, 2*dA, BB, dBB)

    18     rbbr = rbbr.squeeze()

    19     rts = sp.roots(rbbr)

    20     rtsin = rts[dA:2*dA+1]

    21     b = 1

    22     for i in range(dA):

    23         b = sp.convolve(b, sp.array([1, -rtsin[i]]))

    24     bbr = sp.convolve(b, b[::-1])

    25     r = rbbr[0] / bbr[0]

    26     return r,b,rbbr

    27 

    28 if __name__ == "__main__":

    29     A = sp.convolve([-0.5, 1], [-0.9, 1])

    30     dA = 2

    31     B = 0.5*sp.array([-0.9, 1])

    32     dB = 1

    33     rho = 1

    34 

    35     r, beta1, sigma = spec(A, dA, B, dB, rho)

    36     print "sigma", sigma

    37     print "beta1", beta1

    38     print "r", r