equal
deleted
inserted
replaced
|
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 |