##// END OF EJS Templates
jro_soft » compressed_sensing
RSS
Partial implementation of wavelet transform. Implemented functions: dualfilt1, FSfarras, irls_dn (partial), irls_dn2 (partial), deb4_basis (partial).
yamaro -
r17:18
parent child
Show More
Show file before | Show file after | Hide comments
@@ -0,0 +1,65
1 '''
No newline at end of file
2 Created on May 26, 2014
No newline at end of file
3
No newline at end of file
4 @author: Yolian Amaro
No newline at end of file
5 '''
No newline at end of file
6
No newline at end of file
7 import pywt
No newline at end of file
8 import numpy as np
No newline at end of file
9
No newline at end of file
10 def FSfarras():
No newline at end of file
11 #function [af, sf] = FSfarras
No newline at end of file
12
No newline at end of file
13 # Farras filters organized for the dual-tree
No newline at end of file
14 # complex DWT.
No newline at end of file
15 #
No newline at end of file
16 # USAGE:
No newline at end of file
17 # [af, sf] = FSfarras
No newline at end of file
18 # OUTPUT:
No newline at end of file
19 # af{i}, i = 1,2 - analysis filters for tree i
No newline at end of file
20 # sf{i}, i = 1,2 - synthesis filters for tree i
No newline at end of file
21 # See farras, dualtree, dualfilt1.
No newline at end of file
22 #
No newline at end of file
23 # WAVELET SOFTWARE AT POLYTECHNIC UNIVERSITY, BROOKLYN, NY
No newline at end of file
24 # http://taco.poly.edu/WaveletSoftware/
No newline at end of file
25 #
No newline at end of file
26 # Translated to Python by Yolian Amaro
No newline at end of file
27
No newline at end of file
28
No newline at end of file
29
No newline at end of file
30 a1 = np.array( [
No newline at end of file
31 [ 0, 0],
No newline at end of file
32 [-0.08838834764832, -0.01122679215254],
No newline at end of file
33 [ 0.08838834764832, 0.01122679215254],
No newline at end of file
34 [ 0.69587998903400, 0.08838834764832],
No newline at end of file
35 [ 0.69587998903400, 0.08838834764832],
No newline at end of file
36 [ 0.08838834764832, -0.69587998903400],
No newline at end of file
37 [-0.08838834764832, 0.69587998903400],
No newline at end of file
38 [ 0.01122679215254, -0.08838834764832],
No newline at end of file
39 [ 0.01122679215254, -0.08838834764832],
No newline at end of file
40 [0, 0]
No newline at end of file
41 ] );
No newline at end of file
42
No newline at end of file
43 a2 = np.array([
No newline at end of file
44 [ 0.01122679215254, 0],
No newline at end of file
45 [ 0.01122679215254, 0],
No newline at end of file
46 [-0.08838834764832, -0.08838834764832],
No newline at end of file
47 [ 0.08838834764832, -0.08838834764832],
No newline at end of file
48 [ 0.69587998903400, 0.69587998903400],
No newline at end of file
49 [ 0.69587998903400, -0.69587998903400],
No newline at end of file
50 [ 0.08838834764832, 0.08838834764832],
No newline at end of file
51 [-0.08838834764832, 0.08838834764832],
No newline at end of file
52 [ 0, 0.01122679215254],
No newline at end of file
53 [ 0, -0.01122679215254]
No newline at end of file
54 ]);
No newline at end of file
55
No newline at end of file
56 #print a2.shape
No newline at end of file
57
No newline at end of file
58 af = np.array([ [a1,a2] ], dtype=object)
No newline at end of file
59
No newline at end of file
60 s1 = a1[::-1]
No newline at end of file
61 s2 = a2[::-1]
No newline at end of file
62
No newline at end of file
63 sf = np.array([ [s1,s2] ], dtype=object)
No newline at end of file
64
No newline at end of file
65 return af, sf No newline at end of file
Show file before | Show file after | Hide comments
@@ -0,0 +1,38
1 '''
No newline at end of file
2 Created on May 26, 2014
No newline at end of file
3
No newline at end of file
4 @author: Yolian Amaro
No newline at end of file
5 '''
No newline at end of file
6
No newline at end of file
7 import numpy as np
No newline at end of file
8 import FSfarras
No newline at end of file
9 import dualfilt1
No newline at end of file
10
No newline at end of file
11 def deb4_basis(N):
No newline at end of file
12
No newline at end of file
13 Psi = np.zeros(shape=(N,2*N+1));
No newline at end of file
14 idx = 1;
No newline at end of file
15
No newline at end of file
16 J = 4;
No newline at end of file
17 [Faf, Fsf] = FSfarras;
No newline at end of file
18 [af, sf] = dualfilt1;
No newline at end of file
19 # #
No newline at end of file
20 # # # compute transform of zero vector
No newline at end of file
21 # # x = zeros(1,N);
No newline at end of file
22 # # w = dualtree(x, J, Faf, af);
No newline at end of file
23 # #
No newline at end of file
24 # # # Uses both real and imaginary wavelets
No newline at end of file
25 # # for i in range (1, J+1):
No newline at end of file
26 # # for j in range (1, 2):
No newline at end of file
27 # # for k in range (1, (w[i][j]).size):
No newline at end of file
28 # # w[i][j](k) = 1;
No newline at end of file
29 # # y = idualtree(w, J, Fsf, sf);
No newline at end of file
30 # # w[i][j](k) = 0;
No newline at end of file
31 # # # store it
No newline at end of file
32 # # Psi(:,idx) = y.T.conj();
No newline at end of file
33 # # idx = idx + 1;
No newline at end of file
34 # #
No newline at end of file
35 # # # Add uniform vector (seems to be useful if there's a background
No newline at end of file
36 # # Psi(:,2*N+1) = 1/np.sqrt(N);
No newline at end of file
37 #
No newline at end of file
38 # return Psi No newline at end of file
Show file before | Show file after | Hide comments
@@ -0,0 +1,65
1 '''
No newline at end of file
2 Created on May 29, 2014
No newline at end of file
3
No newline at end of file
4 @author: Yolian Amaro
No newline at end of file
5 '''
No newline at end of file
6
No newline at end of file
7 import numpy as np
No newline at end of file
8
No newline at end of file
9 def dualfilt1():
No newline at end of file
10
No newline at end of file
11 # Kingsbury Q-filters for the dual-tree complex DWT
No newline at end of file
12 #
No newline at end of file
13 # USAGE:
No newline at end of file
14 # [af, sf] = dualfilt1
No newline at end of file
15 # OUTPUT:
No newline at end of file
16 # af{i}, i = 1,2 - analysis filters for tree i
No newline at end of file
17 # sf{i}, i = 1,2 - synthesis filters for tree i
No newline at end of file
18 # note: af{2} is the reverse of af{1}
No newline at end of file
19 # REFERENCE:
No newline at end of file
20 # N. G. Kingsbury, "A dual-tree complex wavelet
No newline at end of file
21 # transform with improved orthogonality and symmetry
No newline at end of file
22 # properties", Proceedings of the IEEE Int. Conf. on
No newline at end of file
23 # Image Proc. (ICIP), 2000
No newline at end of file
24 # See dualtree
No newline at end of file
25 #
No newline at end of file
26 # WAVELET SOFTWARE AT POLYTECHNIC UNIVERSITY, BROOKLYN, NY
No newline at end of file
27 # http://taco.poly.edu/WaveletSoftware/
No newline at end of file
28
No newline at end of file
29 # These cofficients are rounded to 8 decimal places.
No newline at end of file
30
No newline at end of file
31 a1 = np.array([
No newline at end of file
32 [ 0.03516384000000, 0],
No newline at end of file
33 [ 0, 0],
No newline at end of file
34 [-0.08832942000000, -0.11430184000000],
No newline at end of file
35 [ 0.23389032000000, 0],
No newline at end of file
36 [ 0.76027237000000, 0.58751830000000],
No newline at end of file
37 [ 0.58751830000000, -0.76027237000000],
No newline at end of file
38 [ 0, 0.23389032000000],
No newline at end of file
39 [-0.11430184000000, 0.08832942000000],
No newline at end of file
40 [ 0, 0],
No newline at end of file
41 [ 0, -0.03516384000000]
No newline at end of file
42 ]);
No newline at end of file
43
No newline at end of file
44 a2 = np.array([
No newline at end of file
45 [ 0, -0.03516384000000],
No newline at end of file
46 [ 0, 0],
No newline at end of file
47 [-0.11430184000000, 0.08832942000000],
No newline at end of file
48 [ 0, 0.23389032000000],
No newline at end of file
49 [ 0.58751830000000, -0.76027237000000],
No newline at end of file
50 [ 0.76027237000000, 0.58751830000000],
No newline at end of file
51 [ 0.23389032000000, 0],
No newline at end of file
52 [ -0.08832942000000, -0.11430184000000],
No newline at end of file
53 [ 0, 0],
No newline at end of file
54 [ 0.03516384000000, 0]
No newline at end of file
55 ]);
No newline at end of file
56
No newline at end of file
57 af = np.array([ [a1,a2] ], dtype=object)
No newline at end of file
58
No newline at end of file
59 s1 = a1[::-1]
No newline at end of file
60 s2 = a2[::-1]
No newline at end of file
61
No newline at end of file
62 sf = np.array([ [s1,s2] ], dtype=object)
No newline at end of file
63
No newline at end of file
64
No newline at end of file
65 return af, sf No newline at end of file
Show file before | Show file after | Hide comments
@@ -0,0 +1,80
1 '''
No newline at end of file
2 Created on May 27, 2014
No newline at end of file
3
No newline at end of file
4 @author: Yolian Amaro
No newline at end of file
5 '''
No newline at end of file
6
No newline at end of file
7 #from scipy.sparse import eye
No newline at end of file
8 from scipy import linalg
No newline at end of file
9 import scipy.sparse as sps
No newline at end of file
10 import numpy as np
No newline at end of file
11
No newline at end of file
12 def irls_dn(A,b,p,lambda1):
No newline at end of file
13
No newline at end of file
14
No newline at end of file
15 # Minimize lambda*||u||_p + ||A*u-b||_2, 0 < p <= 1
No newline at end of file
16 # using Iterative Reweighted Least Squares
No newline at end of file
17 # (see http://math.lanl.gov/Research/Publications/Docs/chartrand-2008-iteratively.pdf
No newline at end of file
18 # and http://web.eecs.umich.edu/~aey/sparse/sparse11.pdf)
No newline at end of file
19
No newline at end of file
20 # Note to self: I found that "warm-starting" didn't really help too much.
No newline at end of file
21
No newline at end of file
22 [M,N] = A.shape;
No newline at end of file
23 # Initialize and precompute:
No newline at end of file
24 eps = 1e-2; # damping parameter
No newline at end of file
25 [Q,R] = linalg.qr(A.T.conj(),0);
No newline at end of file
26 print A.shape
No newline at end of file
27 print R.shape
No newline at end of file
28 print b.shape
No newline at end of file
29 c = linalg.solve(R.T.conj(),b); # will be used later also
No newline at end of file
30 u = Q*c; # minimum 2-norm solution
No newline at end of file
31 I = sps.eye(M);
No newline at end of file
32
No newline at end of file
33 # Spacing of floating point numbers
No newline at end of file
34 eps = np.spacing(1)
No newline at end of file
35
No newline at end of file
36 # Loop until damping parameter is small enough
No newline at end of file
37 while (eps > 1e-7):
No newline at end of file
38 epschange = 0;
No newline at end of file
39 # Loop until it's time to change eps
No newline at end of file
40 while (~epschange):
No newline at end of file
41 # main loop
No newline at end of file
42 # u_n = W*A'*(A*W*A'+ lambda*I)^-1 * b
No newline at end of file
43 # where W = diag(1/w)
No newline at end of file
44 # where w = (u.^2 + eps).^(p/2-1)
No newline at end of file
45
No newline at end of file
46 # Update
No newline at end of file
47 w = (u**2 + eps)**(1-p/2);
No newline at end of file
48
No newline at end of file
49 # Empty temporary N x N matrix
No newline at end of file
50 temp = np.zeros(shape=(N,N))
No newline at end of file
51
No newline at end of file
52 # Sparse matrix
No newline at end of file
53 for i in range (1, N):
No newline at end of file
54 for j in range (1,N):
No newline at end of file
55 if(i==j):
No newline at end of file
56 temp[i,j] = w
No newline at end of file
57
No newline at end of file
58 # Compressed Sparse Matrix
No newline at end of file
59 W = sps.csr_matrix(temp); #Compressed Sparse Row matrix
No newline at end of file
60
No newline at end of file
61
No newline at end of file
62 WAT = W*A.T.conj();
No newline at end of file
63 u_new = WAT * ( linalg.solve (A*WAT + lambda1*I), b);
No newline at end of file
64
No newline at end of file
65 # See if this subproblem is converging
No newline at end of file
66 delu = np.linalg.norm(u_new-u)/np.linalg.norm(u);
No newline at end of file
67 epschange = delu < (np.sqrt(eps)/100);
No newline at end of file
68
No newline at end of file
69 # Make update
No newline at end of file
70 u = u_new;
No newline at end of file
71
No newline at end of file
72
No newline at end of file
73 eps = eps/10; # decrease eps
No newline at end of file
74 # Print info
No newline at end of file
75 print 'eps =',eps;
No newline at end of file
76
No newline at end of file
77 return u
No newline at end of file
78
No newline at end of file
79
No newline at end of file
80 No newline at end of file
Show file before | Show file after | Hide comments
@@ -0,0 +1,50
1 '''
No newline at end of file
2 Created on May 27, 2014
No newline at end of file
3
No newline at end of file
4 @author: Yolian Amaro
No newline at end of file
5 '''
No newline at end of file
6
No newline at end of file
7 from irls_dn import *
No newline at end of file
8
No newline at end of file
9 def irls_dn2(A,b,p,G):
No newline at end of file
10
No newline at end of file
11 # Minimize ||u||_p subject to ||A*u-b||_2^2 <= G (0 < p <= 1)
No newline at end of file
12
No newline at end of file
13 # What this function actually does is finds the lambda1 so that the solution
No newline at end of file
14 # to the following problem satisfies ||A*u-b||_2^2 <= G:
No newline at end of file
15 # Minimize lambda1*||u||_p + ||A*u-b||_2
No newline at end of file
16
No newline at end of file
17 # Start with a large lambda1, and do a line search until fidelity <= G.
No newline at end of file
18 # (Inversions with large lambda1 are really fast anyway).
No newline at end of file
19
No newline at end of file
20 # Then spin up fzero to localize the root even better
No newline at end of file
21
No newline at end of file
22 # Line Search
No newline at end of file
23
No newline at end of file
24 alpha = 2; # Line search parameter
No newline at end of file
25 lambda1 = 1e5; # What's a reasonable but safe initial guess?
No newline at end of file
26 u = irls_dn(A,b,p,lambda1);
No newline at end of file
27 # fid = np.norm(A*u-b)^2;
No newline at end of file
28 #
No newline at end of file
29 # print '----------------------------------\n';
No newline at end of file
30 #
No newline at end of file
31 # while (fid >= G)
No newline at end of file
32 # lambda1 = lambda1 / alpha; # Balance between speed and accuracy
No newline at end of file
33 # u = irls_dn(A,b,p,lambda1);
No newline at end of file
34 # fid = np.norm(A*u-b)^2;
No newline at end of file
35 # print 'lambda1 = #2e \t ||A*u-b||^2 = #.1f\n',lambda1,fid);
No newline at end of file
36 #
No newline at end of file
37 # # Refinement using fzero
No newline at end of file
38 # lambda10 = [lambda1 lambda1*alpha]; # interval with zero-crossing
No newline at end of file
39 # f = @(lambda1) np.norm(A*irls_dn(A,b,p,lambda1) - b)^2 - G;
No newline at end of file
40 # opts = optimset('fzero');
No newline at end of file
41 # # opts.Display = 'iter';
No newline at end of file
42 # opts.Display = 'none';
No newline at end of file
43 # opts.TolX = 0.01*lambda1;
No newline at end of file
44 # lambda1 = fzero(f,lambda10,opts);
No newline at end of file
45 #
No newline at end of file
46 # u = irls_dn(A,b,p,lambda1);
No newline at end of file
47 #
No newline at end of file
48 #
No newline at end of file
49 # return u;
No newline at end of file
50 No newline at end of file
Show file before | Show file after | Hide comments
@@ -1,430 +1,468
1 1 #!/usr/bin/env python No newline at end of file
2 2 No newline at end of file
3 3 #---------------------------------------------------------- No newline at end of file
4 4 # Original MATLAB code developed by Brian Harding No newline at end of file
5 5 # Rewritten in python by Yolian Amaro No newline at end of file
6 6 # Python version 2.7 No newline at end of file
7 7 # May 15, 2014 No newline at end of file
8 8 # Jicamarca Radio Observatory No newline at end of file
9 9 #---------------------------------------------------------- No newline at end of file
10 10 No newline at end of file
11 11 import math
No newline at end of file
12 #import cmath
No newline at end of file
13 #import scipy
No newline at end of file
14 #import matplotlib No newline at end of file
15 12 import numpy as np
No newline at end of file
16 #from numpy import linalg No newline at end of file
17 13 import matplotlib.pyplot as plt No newline at end of file
18 14 from scipy import linalg No newline at end of file
19 15 import time No newline at end of file
20 16 from y_hysell96 import* No newline at end of file
21 17 from deb4_basis import * No newline at end of file
22 18 from modelf import *
19 No newline at end of file
23 from scipy.optimize import fsolve No newline at end of file
24 20 from scipy.optimize import root No newline at end of file
25 21 import pywt No newline at end of file
22 from irls_dn2 import * No newline at end of file
26 23 No newline at end of file
27 24 No newline at end of file
28 25 ## Calculate Forward Model No newline at end of file
29 26 lambda1 = 6.0 No newline at end of file
30 27 k = 2*math.pi/lambda1 No newline at end of file
31 28 No newline at end of file
32 29 ## Calculate Magnetic Declination No newline at end of file
33 30 No newline at end of file
34 31 # [~,~,dec] = igrf11magm(350e3, -11-56/60, -76-52/60, 2012); check this No newline at end of file
35 32 No newline at end of file
36 33 # or calculate it with the above function No newline at end of file
37 34 dec = -1.24 No newline at end of file
38 35 No newline at end of file
39 36 # loads rx, ry (Jicamarca antenna positions) #this can be done with numpy.loadtxt() No newline at end of file
40 37 rx = np.array( [[127.5000], [91.5000], [127.5000], [19.5000], [91.5000], [-127.5000], [-55.5000], [-220.8240]] ) No newline at end of file
41 38 ry = np.array( [[127.5000], [91.5000], [91.5000], [55.5000], [-19.5000], [-127.5000], [-127.5000], [-322.2940]] ) No newline at end of file
42 39 No newline at end of file
43 40 antpos = np.array( [[127.5000, 91.5000, 127.5000, 19.5000, 91.5000, -127.5000, -55.5000, -220.8240], No newline at end of file
44 41 [127.5000, 91.5000, 91.5000, 55.5000, -19.5000, -127.5000, -127.5000, -322.2940]] ) No newline at end of file
45 42 No newline at end of file
46 43 plt.figure(1) No newline at end of file
47 44 plt.plot(rx, ry, 'ro') No newline at end of file
48 45 plt.draw() No newline at end of file
49 46 No newline at end of file
50 47 # Jicamarca is nominally at a 45 degree angle No newline at end of file
51 48 theta = 45 - dec; No newline at end of file
52 49 No newline at end of file
53 50 # Rotation matrix from antenna coord to magnetic coord (East North) No newline at end of file
54 51 theta_rad = math.radians(theta) # trig functions take radians as argument No newline at end of file
55 52 val1 = float( math.cos(theta_rad) ) No newline at end of file
56 53 val2 = float( math.sin(theta_rad) ) No newline at end of file
57 54 val3 = float( -1*math.sin(theta_rad)) No newline at end of file
58 55 val4 = float( math.cos(theta_rad) ) No newline at end of file
59 56 No newline at end of file
60 57 # Rotation matrix from antenna coord to magnetic coord (East North) No newline at end of file
61 58 R = np.array( [[val1, val3], [val2, val4]] ); No newline at end of file
62 59 No newline at end of file
63 60 # Rotate antenna positions to magnetic coord. No newline at end of file
64 61 AR = np.dot(R.T, antpos); No newline at end of file
65 62 No newline at end of file
66 63 # Only take the East component No newline at end of file
67 64 r = AR[0,:] No newline at end of file
68 65 r.sort() # ROW VECTOR? No newline at end of file
69 66 No newline at end of file
70 67 # Truth model (high and low resolution) No newline at end of file
71 68 Nt = (1024.0)*(16.0); # number of pixels in truth image: high resolution No newline at end of file
72 69 thbound = 9.0/180*math.pi; # the width of the domain in angle space No newline at end of file
73 70 thetat = np.linspace(-thbound, thbound,Nt) # image domain No newline at end of file
74 71 thetat = np.transpose(thetat) # transpose # FUNCIONA?????????????????????????????? No newline at end of file
75 72 Nr = (256.0); # number of pixels in reconstructed image: low res No newline at end of file
76 73 thetar = np.linspace(-thbound, thbound,Nr) # reconstruction domain No newline at end of file
77 74 thetar = np.transpose(thetar) #transpose # FUNCIONA????????????????????????????? No newline at end of file
78 75 No newline at end of file
79 76 # Model for f: Gaussian(s) with amplitudes a, centers mu, widths sig, and No newline at end of file
80 77 # background constant b. No newline at end of file
81 78 No newline at end of file
82 79 # Triple Gaussian No newline at end of file
83 80 # a = np.array([3, 5, 2]); No newline at end of file
84 81 # mu = np.array([-5.0/180*math.pi, 2.0/180*math.pi, 7.0/180*math.pi]); No newline at end of file
85 82 # sig = np.array([2.0/180*math.pi, 1.5/180*math.pi, 0.3/180*math.pi]); No newline at end of file
86 83 # b = 0; # background No newline at end of file
87 84 No newline at end of file
88 85 # Double Gaussian No newline at end of file
89 86 # a = np.array([3, 5]); No newline at end of file
90 87 # mu = np.array([-5.0/180*math.pi, 2.0/180*math.pi]); No newline at end of file
91 88 # sig = np.array([2.0/180*math.pi, 1.5/180*math.pi]); No newline at end of file
92 89 # b = 0; # background No newline at end of file
93 90 No newline at end of file
94 91 # Single Gaussian No newline at end of file
95 92 a = np.array( [3] ); No newline at end of file
96 93 mu = np.array( [-3.0/180*math.pi] ) No newline at end of file
97 94 sig = np.array( [2.0/180*math.pi] ) No newline at end of file
98 95 b = 0; No newline at end of file
99 96 No newline at end of file
100 97 fact = np.zeros(shape=(Nt,1)); No newline at end of file
101 98 factr = np.zeros(shape=(Nr,1)); No newline at end of file
102 99 No newline at end of file
103 100 for i in range(0, a.size): No newline at end of file
104 101 temp = (-(thetat-mu[i])**2/(sig[i]**2)) No newline at end of file
105 102 tempr = (-(thetar-mu[i])**2/(sig[i]**2)) No newline at end of file
106 103 for j in range(0, temp.size): No newline at end of file
107 104 fact[j] = fact[j] + a[i]*math.exp(temp[j]); No newline at end of file
108 105 for m in range(0, tempr.size): No newline at end of file
109 106 factr[m] = factr[m] + a[i]*math.exp(tempr[m]); No newline at end of file
110 107 fact = fact + b; No newline at end of file
111 108 factr = factr + b; No newline at end of file
112 109 No newline at end of file
113 110 # # model for f: Square pulse No newline at end of file
114 111 # for j in range(0, fact.size): No newline at end of file
115 112 # if (theta > -5.0/180*math.pi and theta < 2.0/180*math.pi): No newline at end of file
116 113 # fact[j] = 0 No newline at end of file
117 114 # else: No newline at end of file
118 115 # fact[j] = 1 No newline at end of file
119 116 # for k in range(0, factr.size): No newline at end of file
120 117 # if (thetar[k] > -5.0/180*math.pi and thetar[k] < 2/180*math.pi): No newline at end of file
121 118 # fact[k] = 0 No newline at end of file
122 119 # else: No newline at end of file
123 120 # fact[k] = 1 No newline at end of file
124 121 # No newline at end of file
125 122 # No newline at end of file
126 123 # # model for f: triangle pulse No newline at end of file
127 124 # mu = -1.0/180*math.pi; No newline at end of file
128 125 # sig = 5.0/180*math.pi; No newline at end of file
129 126 # wind1 = theta > mu-sig and theta < mu; No newline at end of file
130 127 # wind2 = theta < mu+sig and theta > mu; No newline at end of file
131 128 # fact = wind1 * (theta - (mu - sig)); No newline at end of file
132 129 # factr = wind1 * (thetar - (mu - sig)); No newline at end of file
133 130 # fact = fact + wind2 * (-(theta-(mu+sig))); No newline at end of file
134 131 # factr = factr + wind2 * (-(thetar-(mu+sig))); No newline at end of file
135 132 No newline at end of file
136 133 No newline at end of file
137 134 # fact = fact/(sum(fact)[0]*2*thbound/Nt); # normalize to integral(f)==1 No newline at end of file
138 135 I = sum(fact)[0]; No newline at end of file
139 136 fact = fact/I; # normalize to sum(f)==1 No newline at end of file
140 137 factr = factr/I; # normalize to sum(f)==1
138 No newline at end of file
141 #plot(thetat,fact,'r'); hold on;
No newline at end of file
139 No newline at end of file
142 #plot(thetar,factr,'k.'); hold off; No newline at end of file
No newline at end of file
140 #plt.plot(thetar,factr,'k.'); No newline at end of file
143 141 #xlim([min(thetat) max(thetat)]); No newline at end of file
144 142 No newline at end of file
145 143 #x = np.linspace(thetat.min(), thetat.max) ???? No newline at end of file
146 144 #for i in range(0, thetat.size): No newline at end of file
147 145 plt.figure(2) No newline at end of file
148 146 plt.plot(thetat, fact, 'r--') No newline at end of file
149 147 plt.plot(thetar, factr, 'ro') No newline at end of file
150 148 plt.draw() No newline at end of file
151 149 # xlim([min(thetat) max(thetat)]); FALTA ARREGLAR ESTO No newline at end of file
152 150 No newline at end of file
153 151 No newline at end of file
154 152 ## No newline at end of file
155 153 # Control the type and number of inversions with: No newline at end of file
156 154 # SNRdBvec: the SNRs that will be used. No newline at end of file
157 155 # NN: the number of trials for each SNR No newline at end of file
158 156 No newline at end of file
159 157 #SNRdBvec = np.linspace(5,20,10); No newline at end of file
160 158 SNRdBvec = np.array([15]); No newline at end of file
161 159 NN = 1; # number of trial at each SNR No newline at end of file
162 160 No newline at end of file
163 161 # if using vector arguments should be: (4,SNRdBvec.size,NN) No newline at end of file
164 162 corr = np.zeros(shape=(4,SNRdBvec.size,NN)); # (method, SNR, trial) No newline at end of file
165 163 corrc = np.zeros(shape=(4,SNRdBvec.size,NN)); # (method, SNR, trial) No newline at end of file
166 164 rmse = np.zeros(shape=(4,SNRdBvec.size,NN)); # (method, SNR, trial) No newline at end of file
167 165 No newline at end of file
168 166 for snri in range(0, SNRdBvec.size): # change 1 for SNRdBvec.size when using SNRdBvec as vector No newline at end of file
169 167 for Ni in range(0, NN): No newline at end of file
170 168 SNRdB = SNRdBvec[snri]; No newline at end of file
171 169 SNR = 10**(SNRdB/10.0); No newline at end of file
172 170 No newline at end of file
173 171 # Calculate cross-correlation matrix (Fourier components of image) No newline at end of file
174 172 # This is an inefficient way to do this. No newline at end of file
175 173 R = np.zeros(shape=(r.size, r.size), dtype=object); No newline at end of file
176 174 No newline at end of file
177 175 for i1 in range(0, r.size): No newline at end of file
178 176 for i2 in range(0,r.size): No newline at end of file
179 177 R[i1,i2] = np.dot(fact.T, np.exp(1j*k*np.dot((r[i1]-r[i2]),np.sin(thetat)))) No newline at end of file
180 178 R[i1,i2] = sum(R[i1,i2]) No newline at end of file
181 179 No newline at end of file
182 180 # Add uncertainty No newline at end of file
183 181 # This is an ad-hoc way of adding "noise". It models some combination of No newline at end of file
184 182 # receiver noise and finite integration times. We could use a more No newline at end of file
185 183 # advanced model (like in Yu et al 2000) in the future. No newline at end of file
186 184 No newline at end of file
187 185 # This is a way of adding noise while maintaining the No newline at end of file
188 186 # positive-semi-definiteness of the matrix. No newline at end of file
189 187 No newline at end of file
190 188 U = linalg.cholesky(R.astype(complex), lower=False); # U'*U = R No newline at end of file
191 189 No newline at end of file
192 190 sigma_noise = (np.linalg.norm(U,'fro')/SNR); No newline at end of file
193 191 No newline at end of file
194 192 temp1 = (-1*np.random.rand(U.shape[0], U.shape[1]) + 0.5) No newline at end of file
195 193 temp2 = 1j*(-1*np.random.rand(U.shape[0], U.shape[1]) + 0.5) No newline at end of file
196 194 temp3 = ((abs(U) > 0).astype(float)) # upper triangle of 1's No newline at end of file
197 195 temp4 = (sigma_noise * (temp1 + temp2))/np.sqrt(2.0) No newline at end of file
198 196 No newline at end of file
199 197 nz = np.multiply(temp4, temp3) No newline at end of file
200 198 No newline at end of file
199 #---------------------- Eliminar esto:------------------------------------------
No newline at end of file
200 #nz = ((abs(np.multiply(temp4, temp3)) > 0).astype(int)) No newline at end of file
201 201 #nz = ((abs(np.dot(temp4, temp3)) > 0).astype(int)) No newline at end of file
202 202 No newline at end of file
203 203 No newline at end of file
204 204 #nz = np.dot(np.dot(sigma_noise, (temp1 + temp2)/math.sqrt(2), temp3 )); No newline at end of file
205 205 #nz = np.dot(sigma_noise, (np.dot((np.random.rand(8,8) + j*np.random.rand(8,8))/math.sqrt(2.0) , (abs(U) > 0).astype(int)))); No newline at end of file
206 Unz = U + nz; No newline at end of file
206 207 No newline at end of file
207 208 Unz = U + nz; No newline at end of file
208 209 Rnz = np.dot(Unz.T.conj(),Unz); # the noisy version of R No newline at end of file
209 210 plt.figure(3); No newline at end of file
210 211 plt.pcolor(abs(Rnz)); No newline at end of file
211 212 plt.colorbar(); No newline at end of file
212 213 No newline at end of file
213 214 # Fourier Inversion ################### No newline at end of file
214 215 f_fourier = np.zeros(shape=(Nr,1), dtype=complex); No newline at end of file
215 216 No newline at end of file
216 217 for i in range(0, thetar.size): No newline at end of file
217 218 th = thetar[i]; No newline at end of file
218 219 w = np.exp(1j*k*np.dot(r,np.sin(th))); No newline at end of file
219 220 No newline at end of file
220 221 temp = np.dot(w.T.conj(),U) No newline at end of file
221 222 No newline at end of file
222 223 f_fourier[i] = np.dot(temp, w); No newline at end of file
223 224 No newline at end of file
224 225 f_fourier = f_fourier.real; # get rid of numerical imaginary noise No newline at end of file
225 226 No newline at end of file
226 227 #print f_fourier No newline at end of file
227 228 No newline at end of file
228 229 No newline at end of file
229 230 # Capon Inversion ###################### No newline at end of file
230 231 No newline at end of file
231 232 f_capon = np.zeros(shape=(Nr,1)); No newline at end of file
232 233
234 No newline at end of file
233 #tic_capon = time.time(); No newline at end of file
234 235 No newline at end of file
235 236 for i in range(0, thetar.size): No newline at end of file
236 237 th = thetar[i]; No newline at end of file
237 238 w = np.exp(1j*k*np.dot(r,np.sin(th))); No newline at end of file
238 239 f_capon[i] = np.divide(1, ( np.dot( w.T.conj(), (linalg.solve(Rnz,w)) ) ).real) No newline at end of file
239 240
241 No newline at end of file
240 #toc_capon = time.time()
No newline at end of file
242 No newline at end of file
241
No newline at end of file
243 No newline at end of file
242 # elapsed_time_capon = toc_capon - tic_capon; No newline at end of file
No newline at end of file
244 f_capon = f_capon.real; # get rid of numerical imaginary noise No newline at end of file
243 245 No newline at end of file
244 246 f_capon = f_capon.real; # get rid of numerical imaginary noise No newline at end of file
245 247 No newline at end of file
246 248 # MaxEnt Inversion ##################### No newline at end of file
247 249 No newline at end of file
248 250 # create the appropriate sensing matrix (split into real and imaginary # parts) No newline at end of file
249 251 M = (r.size-1)*(r.size); No newline at end of file
250 252 Ht = np.zeros(shape=(M,Nt)); # "true" sensing matrix No newline at end of file
251 253 Hr = np.zeros(shape=(M,Nr)); # approximate sensing matrix for reconstruction No newline at end of file
252 254 No newline at end of file
253 255 # need to re-index our measurements from matrix R into vector g No newline at end of file
254 256 g = np.zeros(shape=(M,1)); No newline at end of file
255 257 gnz = np.zeros(shape=(M,1)); # noisy version of g No newline at end of file
256 258 No newline at end of file
257 259 # triangular indexing to perform this re-indexing No newline at end of file
258 260 T = np.ones(shape=(r.size,r.size)); No newline at end of file
259 261 [i1v,i2v] = np.where(np.triu(T,1) > 0); # converts linear to triangular indexing No newline at end of file
260 262 No newline at end of file
261 263 # build H No newline at end of file
262 264 for i1 in range(0, r.size): No newline at end of file
263 265 for i2 in range(i1+1, r.size): No newline at end of file
264 266 idx = np.where(np.logical_and((i1==i1v), (i2==i2v)))[0]; # kind of awkward No newline at end of file
265 267 idx1 = 2*idx; # because index starts at 0 No newline at end of file
266 268 idx2 = 2*idx+1; No newline at end of file
267 269 Hr[idx1,:] = np.cos(k*(r[i1]-r[i2])*np.sin(thetar)).T; No newline at end of file
268 270 Hr[idx2,:] = np.sin(k*(r[i1]-r[i2])*np.sin(thetar)).T; No newline at end of file
269 271 Ht[idx1,:] = np.cos(k*(r[i1]-r[i2])*np.sin(thetat)).T*Nr/Nt; No newline at end of file
270 272 Ht[idx2,:] = np.sin(k*(r[i1]-r[i2])*np.sin(thetat)).T*Nr/Nt; No newline at end of file
271 273 g[idx1] = (R[i1,i2]).real*Nr/Nt; # check this again later No newline at end of file
272 274 g[idx2] = (R[i1,i2]).imag*Nr/Nt; # check again No newline at end of file
273 275 gnz[idx1] = (Rnz[i1,i2]).real*Nr/Nt; No newline at end of file
274 276 gnz[idx2] = (Rnz[i1,i2]).imag*Nr/Nt; No newline at end of file
275 277 No newline at end of file
276 278 # inversion No newline at end of file
277 279 F = Nr/Nt; # normalization No newline at end of file
278 280 sigma = 1; # set to 1 because the difference is accounted for in G No newline at end of file
279 281 No newline at end of file
280 282 ##### ADD *10 for consistency with old model, NEED TO VERIFY THIS!!!!? line below
283 No newline at end of file
281 G = np.linalg.norm(g-gnz)**2; # pretend we know in advance the actual value of chi^2 No newline at end of file
No newline at end of file
284
No newline at end of file
285 lambda0 = 1e-5*np.ones(shape=(M,1)); # initial condition (can be set to anything) No newline at end of file
282 286 No newline at end of file
283 287 lambda0 = 1e-5*np.ones(shape=(M,1)); # initial condition (can be set to anything) No newline at end of file
288 elapsed_time_maxent = toc_maxent - tic_maxent;
No newline at end of file
289
No newline at end of file
290 # Whitened solution No newline at end of file
284 291 No newline at end of file
285 292 # Whitened solution No newline at end of file
286 293 def myfun(lambda1): No newline at end of file
287 294 return y_hysell96(lambda1,gnz,sigma,F,G,Hr); No newline at end of file
288 295 No newline at end of file
289 296 tic_maxEnt = time.time(); No newline at end of file
290 297 No newline at end of file
291 298 #sol1 = fsolve(myfun,lambda0.ravel(), args=(), xtol=1e-14, maxfev=100000); No newline at end of file
292 299 lambda1 = root(myfun,lambda0, method='krylov', tol=1e-14); No newline at end of file
293 300 No newline at end of file
294 301 #print lambda1 No newline at end of file
295 302 #print lambda1.x No newline at end of file
296 303 No newline at end of file
297 304 lambda1 = lambda1.x; No newline at end of file
298 305 No newline at end of file
299 306 toc_maxEnt = time.time(); No newline at end of file
300 307 f_maxent = modelf(lambda1, Hr, F); No newline at end of file
301 308 ystar = myfun(lambda1); No newline at end of file
302 309 Lambda = np.sqrt(sum(lambda1**2.*sigma**2)/(4*G)); No newline at end of file
303 310 ep = np.multiply(-lambda1,sigma**2)/ (2*Lambda); No newline at end of file
304 311 es = np.dot(Hr, f_maxent) - gnz; # should be same as ep No newline at end of file
305 312 chi2 = np.sum((es/sigma)**2); No newline at end of file
306 313 No newline at end of file
307 314 No newline at end of file
308 315 No newline at end of file
309 316 # CS inversion using irls ######################## No newline at end of file
310 317 No newline at end of file
311 318 # (Use Nr, thetar, gnz, and Hr from MaxEnt above) No newline at end of file
312 319
320 No newline at end of file
313 #Psi = deb4_basis(Nr); ###### REPLACED BY LINE BELOW (?)
No newline at end of file
321 No newline at end of file
314
No newline at end of file
322 No newline at end of file
315 wavelet1 = pywt.Wavelet('db4')
No newline at end of file
323 No newline at end of file
316 Phi, Psi, x = wavelet1.wavefun(level=3) No newline at end of file
No newline at end of file
324 # add "sum to 1" constraint
No newline at end of file
325 H2 = np.concatenate( (Hr, np.ones(shape=(1,Nr))), axis=0 ); No newline at end of file
317 326 No newline at end of file
318 327 # add "sum to 1" constraint No newline at end of file
319 328 H2 = np.concatenate( (Hr, np.ones(shape=(1,Nr))), axis=0 ); No newline at end of file
320 329 N_temp = np.array([[Nr/Nt]]); No newline at end of file
321 330 g2 = np.concatenate( (gnz, N_temp), axis=0 );
331 No newline at end of file
322
No newline at end of file
332 No newline at end of file
323
No newline at end of file
333 No newline at end of file
324 #s = irls_dn2(H2*Psi,g2,0.5,G); No newline at end of file
No newline at end of file
334 # f_cs = Psi*s;
No newline at end of file
335 #
No newline at end of file
336 # # plot No newline at end of file
325 337 # f_cs = Psi*s; No newline at end of file
326 338 # No newline at end of file
327 339 # # plot No newline at end of file
328 340 # plot(thetar,f_cs,'r.-'); No newline at end of file
329 341 # hold on; No newline at end of file
330 342 # plot(thetat,fact,'k-'); No newline at end of file
331 343 # hold off; No newline at end of file
332 344 No newline at end of file
333 345 No newline at end of file
334 346 # # # Scaling and shifting No newline at end of file
335 347 # # # Only necessary for capon solution No newline at end of file
336 348 No newline at end of file
337 349 No newline at end of file
338 350 f_capon = f_capon/np.max(f_capon)*np.max(fact); No newline at end of file
339 351 No newline at end of file
340 352 No newline at end of file
341 353 ### analyze stuff ###################### No newline at end of file
342 354 # calculate MSE No newline at end of file
343 355 rmse_fourier = np.sqrt(np.mean((f_fourier - factr)**2)); No newline at end of file
344 356 rmse_capon = np.sqrt(np.mean((f_capon - factr)**2)); No newline at end of file
345 357 rmse_maxent = np.sqrt(np.mean((f_maxent - factr)**2)); No newline at end of file
346 358 #rmse_cs = np.sqrt(np.mean((f_cs - factr).^2)); No newline at end of file
347 359 No newline at end of file
348 360 relrmse_fourier = rmse_fourier / np.linalg.norm(fact); No newline at end of file
349 361 relrmse_capon = rmse_capon / np.linalg.norm(fact); No newline at end of file
350 362 relrmse_maxent = rmse_maxent / np.linalg.norm(fact); No newline at end of file
351 363 #relrmse_cs = rmse_cs / np.norm(fact); No newline at end of file
352 364
365 No newline at end of file
353 factr = factr.T.conj() No newline at end of file
No newline at end of file
366 #f_fourier = f_fourier.T.conj()
No newline at end of file
367 #f_capon = f_capon.T.conj()
No newline at end of file
368 #f_maxent = f_maxent.T.conj()
No newline at end of file
369
No newline at end of file
370 #factr = factr.T.conj()
No newline at end of file
371
No newline at end of file
372 # calculate correlation No newline at end of file
354 373 No newline at end of file
355 374 # calculate correlation
375 No newline at end of file
356 corr_fourier = np.dot(f_fourier,factr) / (np.linalg.norm(f_fourier)*np.linalg.norm(factr));
No newline at end of file
376 No newline at end of file
357 corr_capon = np.dot(f_capon,factr) / (np.linalg.norm(f_capon)*np.linalg.norm(factr));
No newline at end of file
377 No newline at end of file
358 corr_maxent = np.dot(f_maxent,factr) / (np.linalg.norm(f_maxent)*np.linalg.norm(factr)); No newline at end of file
No newline at end of file
378 No newline at end of file
359 379 #corr_cs = np.dot(f_cs,factr) / (norm(f_cs)*norm(factr)); No newline at end of file
360 380 No newline at end of file
361 381 # calculate centered correlation No newline at end of file
362 382 f0 = factr - np.mean(factr); No newline at end of file
363 383 f1 = f_fourier - np.mean(f_fourier);
384 No newline at end of file
364 corrc_fourier = np.dot(f0,f1) / (np.linalg.norm(f0)*np.linalg.norm(f1)); No newline at end of file
No newline at end of file
385 f1 = f_capon - np.mean(f_capon); No newline at end of file
365 386 f1 = f_capon - np.mean(f_capon);
387 No newline at end of file
366 corrc_capon = np.dot(f0,f1) / (np.linalg.norm(f0)*np.linalg.norm(f1)); No newline at end of file
367 388 f1 = f_maxent - np.mean(f_maxent);
389 No newline at end of file
368 corrc_maxent = np.dot(f0,f1) / (np.linalg.norm(f0)*np.linalg.norm(f1)); No newline at end of file
369 390 #f1 = f_cs - mean(f_cs); No newline at end of file
370 391 #corrc_cs = dot(f0,f1) / (norm(f0)*norm(f1)); No newline at end of file
371 392 No newline at end of file
372 393 No newline at end of file
373 394 No newline at end of file
374 395 # # # plot stuff ######################### No newline at end of file
375 396 No newline at end of file
376 397 #---- Capon---- No newline at end of file
377 398 plt.figure(4) No newline at end of file
378 399 plt.subplot(2, 1, 1) No newline at end of file
379 400 plt.plot(180/math.pi*thetar, f_capon, 'r', label='Capon') No newline at end of file
380 401 plt.plot(180/math.pi*thetat,fact, 'k--', label='Truth') No newline at end of file
381 402 plt.ylabel('Power (arbitrary units)') No newline at end of file
382 403 plt.legend(loc='upper right') No newline at end of file
383 404 No newline at end of file
384 405 # formatting y-axis No newline at end of file
385 406 locs,labels = plt.yticks() No newline at end of file
386 407 plt.yticks(locs, map(lambda x: "%.1f" % x, locs*1e4)) No newline at end of file
387 408 plt.text(0.0, 1.01, '1e-4', fontsize=10, transform = plt.gca().transAxes) No newline at end of file
388 409 No newline at end of file
389 410 No newline at end of file
390 411 #---- MaxEnt---- No newline at end of file
391 412 plt.subplot(2, 1, 2) No newline at end of file
392 413 plt.plot(180/math.pi*thetar, f_maxent, 'r', label='MaxEnt') No newline at end of file
393 414 plt.plot(180/math.pi*thetat,fact, 'k--', label='Truth') No newline at end of file
394 415 plt.ylabel('Power (arbitrary units)') No newline at end of file
395 416 plt.legend(loc='upper right') No newline at end of file
396 417 No newline at end of file
397 418 # formatting y-axis No newline at end of file
398 419 locs,labels = plt.yticks() No newline at end of file
399 420 plt.yticks(locs, map(lambda x: "%.1f" % x, locs*1e4)) No newline at end of file
400 421 plt.text(0.0, 1.01, '1e-4', fontsize=10, transform = plt.gca().transAxes) No newline at end of file
401 422 No newline at end of file
402 423 plt.show() No newline at end of file
403 424 No newline at end of file
404 425
426 No newline at end of file
405 # # subplot(3,1,2);
No newline at end of file
406 # # plot(180/pi*thetar,f_maxent,'r-');
No newline at end of file
407 # # hold on;
No newline at end of file
408 # # plot(180/pi*thetat,fact,'k--');
No newline at end of file
409 # # hold off;
No newline at end of file
410 # # ylim([min(f_cs) 1.1*max(fact)]);
No newline at end of file
411 # # # title(sprintf('rel. RMSE: #.2e\tCorr: #.3f Corrc: #.3f', relrmse_maxent, corr_maxent, corrc_maxent));
No newline at end of file
412 # # ylabel({'Power';'(arbitrary units)'})
No newline at end of file
413 # # # title 'Maximum Entropy Method'
No newline at end of file
414 # # legend('MaxEnt','Truth'); No newline at end of file
415 427 # # No newline at end of file
416 428 # # subplot(3,1,3); No newline at end of file
417 429 # # plot(180/pi*thetar,f_cs,'r-'); No newline at end of file
418 430 # # hold on; No newline at end of file
419 431 # # plot(180/pi*thetat,fact,'k--'); No newline at end of file
420 432 # # hold off; No newline at end of file
421 433 # # ylim([min(f_cs) 1.1*max(fact)]); No newline at end of file
422 434 # # # title(sprintf('rel. RMSE: #.2e\tCorr: #.3f Corrc: #.3f', relrmse_cs, corr_cs, corrc_cs)); No newline at end of file
423 435 # # # title 'Compressed Sensing - Debauchies Wavelets' No newline at end of file
424 436 # # xlabel 'Degrees' No newline at end of file
425 437 # # ylabel({'Power';'(arbitrary units)'}) No newline at end of file
426 438 # # legend('Comp. Sens.','Truth'); No newline at end of file
427 439 # # No newline at end of file
428 440 # # # set(gcf,'Position',[749 143 528 881]); # CSL No newline at end of file
429 441 # # # set(gcf,'Position',[885 -21 528 673]); # macbook No newline at end of file
430 442 # # pause(0.01); No newline at end of file
443
No newline at end of file
444 # # Store Results
No newline at end of file
445 corr[0, snri, Ni] = corr_fourier;
No newline at end of file
446 corr[1, snri, Ni] = corr_capon;
No newline at end of file
447 corr[2, snri, Ni] = corr_maxent;
No newline at end of file
448 #corr[3, snri, Ni] = corr_cs;
No newline at end of file
449
No newline at end of file
450 rmse[0,snri,Ni] = relrmse_fourier;
No newline at end of file
451 rmse[1,snri,Ni] = relrmse_capon;
No newline at end of file
452 rmse[2,snri,Ni] = relrmse_maxent;
No newline at end of file
453 #rmse[3,snri,Ni] = relrmse_cs;
No newline at end of file
454
No newline at end of file
455 corrc[0,snri,Ni] = corrc_fourier;
No newline at end of file
456 corrc[1,snri,Ni] = corrc_capon;
No newline at end of file
457 corrc[2,snri,Ni] = corrc_maxent;
No newline at end of file
458 #corrc[3,snri,Ni] = corrc_cs;
No newline at end of file
459
No newline at end of file
460
No newline at end of file
461 print 'Capon:\t', elapsed_time_capon, 'sec';
No newline at end of file
462 print 'Maxent:\t',elapsed_time_maxent, 'sec';
No newline at end of file
463 #print 'CS:\t%3.3f sec\n',elapsed_time_cs;
No newline at end of file
464
No newline at end of file
465 print (NN*(snri+1) + Ni), '/', (SNRdBvec.size*NN);
No newline at end of file
466
No newline at end of file
467 print corr
No newline at end of file
468 No newline at end of file
General Comments 0
You need to be logged in to leave comments. Login now