@@ -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 |
@@ -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 |
@@ -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 |
@@ -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 |
@@ -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 |
@@ -9,20 +9,17 | |||
|
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 |
General Comments 0
You need to be logged in to leave comments.
Login now