##// 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 '''
2 Created on May 26, 2014
3
4 @author: Yolian Amaro
5 '''
6
7 import pywt
8 import numpy as np
9
10 def FSfarras():
11 #function [af, sf] = FSfarras
12
13 # Farras filters organized for the dual-tree
14 # complex DWT.
15 #
16 # USAGE:
17 # [af, sf] = FSfarras
18 # OUTPUT:
19 # af{i}, i = 1,2 - analysis filters for tree i
20 # sf{i}, i = 1,2 - synthesis filters for tree i
21 # See farras, dualtree, dualfilt1.
22 #
23 # WAVELET SOFTWARE AT POLYTECHNIC UNIVERSITY, BROOKLYN, NY
24 # http://taco.poly.edu/WaveletSoftware/
25 #
26 # Translated to Python by Yolian Amaro
27
28
29
30 a1 = np.array( [
31 [ 0, 0],
32 [-0.08838834764832, -0.01122679215254],
33 [ 0.08838834764832, 0.01122679215254],
34 [ 0.69587998903400, 0.08838834764832],
35 [ 0.69587998903400, 0.08838834764832],
36 [ 0.08838834764832, -0.69587998903400],
37 [-0.08838834764832, 0.69587998903400],
38 [ 0.01122679215254, -0.08838834764832],
39 [ 0.01122679215254, -0.08838834764832],
40 [0, 0]
41 ] );
42
43 a2 = np.array([
44 [ 0.01122679215254, 0],
45 [ 0.01122679215254, 0],
46 [-0.08838834764832, -0.08838834764832],
47 [ 0.08838834764832, -0.08838834764832],
48 [ 0.69587998903400, 0.69587998903400],
49 [ 0.69587998903400, -0.69587998903400],
50 [ 0.08838834764832, 0.08838834764832],
51 [-0.08838834764832, 0.08838834764832],
52 [ 0, 0.01122679215254],
53 [ 0, -0.01122679215254]
54 ]);
55
56 #print a2.shape
57
58 af = np.array([ [a1,a2] ], dtype=object)
59
60 s1 = a1[::-1]
61 s2 = a2[::-1]
62
63 sf = np.array([ [s1,s2] ], dtype=object)
64
65 return af, sf
Show file before | Show file after | Hide comments
@@ -0,0 +1,38
1 '''
2 Created on May 26, 2014
3
4 @author: Yolian Amaro
5 '''
6
7 import numpy as np
8 import FSfarras
9 import dualfilt1
10
11 def deb4_basis(N):
12
13 Psi = np.zeros(shape=(N,2*N+1));
14 idx = 1;
15
16 J = 4;
17 [Faf, Fsf] = FSfarras;
18 [af, sf] = dualfilt1;
19 # #
20 # # # compute transform of zero vector
21 # # x = zeros(1,N);
22 # # w = dualtree(x, J, Faf, af);
23 # #
24 # # # Uses both real and imaginary wavelets
25 # # for i in range (1, J+1):
26 # # for j in range (1, 2):
27 # # for k in range (1, (w[i][j]).size):
28 # # w[i][j](k) = 1;
29 # # y = idualtree(w, J, Fsf, sf);
30 # # w[i][j](k) = 0;
31 # # # store it
32 # # Psi(:,idx) = y.T.conj();
33 # # idx = idx + 1;
34 # #
35 # # # Add uniform vector (seems to be useful if there's a background
36 # # Psi(:,2*N+1) = 1/np.sqrt(N);
37 #
38 # return Psi No newline at end of file
Show file before | Show file after | Hide comments
@@ -0,0 +1,65
1 '''
2 Created on May 29, 2014
3
4 @author: Yolian Amaro
5 '''
6
7 import numpy as np
8
9 def dualfilt1():
10
11 # Kingsbury Q-filters for the dual-tree complex DWT
12 #
13 # USAGE:
14 # [af, sf] = dualfilt1
15 # OUTPUT:
16 # af{i}, i = 1,2 - analysis filters for tree i
17 # sf{i}, i = 1,2 - synthesis filters for tree i
18 # note: af{2} is the reverse of af{1}
19 # REFERENCE:
20 # N. G. Kingsbury, "A dual-tree complex wavelet
21 # transform with improved orthogonality and symmetry
22 # properties", Proceedings of the IEEE Int. Conf. on
23 # Image Proc. (ICIP), 2000
24 # See dualtree
25 #
26 # WAVELET SOFTWARE AT POLYTECHNIC UNIVERSITY, BROOKLYN, NY
27 # http://taco.poly.edu/WaveletSoftware/
28
29 # These cofficients are rounded to 8 decimal places.
30
31 a1 = np.array([
32 [ 0.03516384000000, 0],
33 [ 0, 0],
34 [-0.08832942000000, -0.11430184000000],
35 [ 0.23389032000000, 0],
36 [ 0.76027237000000, 0.58751830000000],
37 [ 0.58751830000000, -0.76027237000000],
38 [ 0, 0.23389032000000],
39 [-0.11430184000000, 0.08832942000000],
40 [ 0, 0],
41 [ 0, -0.03516384000000]
42 ]);
43
44 a2 = np.array([
45 [ 0, -0.03516384000000],
46 [ 0, 0],
47 [-0.11430184000000, 0.08832942000000],
48 [ 0, 0.23389032000000],
49 [ 0.58751830000000, -0.76027237000000],
50 [ 0.76027237000000, 0.58751830000000],
51 [ 0.23389032000000, 0],
52 [ -0.08832942000000, -0.11430184000000],
53 [ 0, 0],
54 [ 0.03516384000000, 0]
55 ]);
56
57 af = np.array([ [a1,a2] ], dtype=object)
58
59 s1 = a1[::-1]
60 s2 = a2[::-1]
61
62 sf = np.array([ [s1,s2] ], dtype=object)
63
64
65 return af, sf
Show file before | Show file after | Hide comments
@@ -0,0 +1,80
1 '''
2 Created on May 27, 2014
3
4 @author: Yolian Amaro
5 '''
6
7 #from scipy.sparse import eye
8 from scipy import linalg
9 import scipy.sparse as sps
10 import numpy as np
11
12 def irls_dn(A,b,p,lambda1):
13
14
15 # Minimize lambda*||u||_p + ||A*u-b||_2, 0 < p <= 1
16 # using Iterative Reweighted Least Squares
17 # (see http://math.lanl.gov/Research/Publications/Docs/chartrand-2008-iteratively.pdf
18 # and http://web.eecs.umich.edu/~aey/sparse/sparse11.pdf)
19
20 # Note to self: I found that "warm-starting" didn't really help too much.
21
22 [M,N] = A.shape;
23 # Initialize and precompute:
24 eps = 1e-2; # damping parameter
25 [Q,R] = linalg.qr(A.T.conj(),0);
26 print A.shape
27 print R.shape
28 print b.shape
29 c = linalg.solve(R.T.conj(),b); # will be used later also
30 u = Q*c; # minimum 2-norm solution
31 I = sps.eye(M);
32
33 # Spacing of floating point numbers
34 eps = np.spacing(1)
35
36 # Loop until damping parameter is small enough
37 while (eps > 1e-7):
38 epschange = 0;
39 # Loop until it's time to change eps
40 while (~epschange):
41 # main loop
42 # u_n = W*A'*(A*W*A'+ lambda*I)^-1 * b
43 # where W = diag(1/w)
44 # where w = (u.^2 + eps).^(p/2-1)
45
46 # Update
47 w = (u**2 + eps)**(1-p/2);
48
49 # Empty temporary N x N matrix
50 temp = np.zeros(shape=(N,N))
51
52 # Sparse matrix
53 for i in range (1, N):
54 for j in range (1,N):
55 if(i==j):
56 temp[i,j] = w
57
58 # Compressed Sparse Matrix
59 W = sps.csr_matrix(temp); #Compressed Sparse Row matrix
60
61
62 WAT = W*A.T.conj();
63 u_new = WAT * ( linalg.solve (A*WAT + lambda1*I), b);
64
65 # See if this subproblem is converging
66 delu = np.linalg.norm(u_new-u)/np.linalg.norm(u);
67 epschange = delu < (np.sqrt(eps)/100);
68
69 # Make update
70 u = u_new;
71
72
73 eps = eps/10; # decrease eps
74 # Print info
75 print 'eps =',eps;
76
77 return u
78
79
80
Show file before | Show file after | Hide comments
@@ -0,0 +1,50
1 '''
2 Created on May 27, 2014
3
4 @author: Yolian Amaro
5 '''
6
7 from irls_dn import *
8
9 def irls_dn2(A,b,p,G):
10
11 # Minimize ||u||_p subject to ||A*u-b||_2^2 <= G (0 < p <= 1)
12
13 # What this function actually does is finds the lambda1 so that the solution
14 # to the following problem satisfies ||A*u-b||_2^2 <= G:
15 # Minimize lambda1*||u||_p + ||A*u-b||_2
16
17 # Start with a large lambda1, and do a line search until fidelity <= G.
18 # (Inversions with large lambda1 are really fast anyway).
19
20 # Then spin up fzero to localize the root even better
21
22 # Line Search
23
24 alpha = 2; # Line search parameter
25 lambda1 = 1e5; # What's a reasonable but safe initial guess?
26 u = irls_dn(A,b,p,lambda1);
27 # fid = np.norm(A*u-b)^2;
28 #
29 # print '----------------------------------\n';
30 #
31 # while (fid >= G)
32 # lambda1 = lambda1 / alpha; # Balance between speed and accuracy
33 # u = irls_dn(A,b,p,lambda1);
34 # fid = np.norm(A*u-b)^2;
35 # print 'lambda1 = #2e \t ||A*u-b||^2 = #.1f\n',lambda1,fid);
36 #
37 # # Refinement using fzero
38 # lambda10 = [lambda1 lambda1*alpha]; # interval with zero-crossing
39 # f = @(lambda1) np.norm(A*irls_dn(A,b,p,lambda1) - b)^2 - G;
40 # opts = optimset('fzero');
41 # # opts.Display = 'iter';
42 # opts.Display = 'none';
43 # opts.TolX = 0.01*lambda1;
44 # lambda1 = fzero(f,lambda10,opts);
45 #
46 # u = irls_dn(A,b,p,lambda1);
47 #
48 #
49 # return u;
50
Show file before | Show file after | Hide comments
@@ -9,20 +9,17
9 #----------------------------------------------------------
9 #----------------------------------------------------------
10
10
11 import math
11 import math
12 #import cmath
13 #import scipy
14 #import matplotlib
15 import numpy as np
12 import numpy as np
16 #from numpy import linalg
17 import matplotlib.pyplot as plt
13 import matplotlib.pyplot as plt
18 from scipy import linalg
14 from scipy import linalg
19 import time
15 import time
20 from y_hysell96 import*
16 from y_hysell96 import*
21 from deb4_basis import *
17 from deb4_basis import *
22 from modelf import *
18 from modelf import *
23 from scipy.optimize import fsolve
19 #from scipy.optimize import fsolve
24 from scipy.optimize import root
20 from scipy.optimize import root
25 import pywt
21 import pywt
22 from irls_dn2 import *
26
23
27
24
28 ## Calculate Forward Model
25 ## Calculate Forward Model
@@ -138,8 +135,9
138 I = sum(fact)[0];
135 I = sum(fact)[0];
139 fact = fact/I; # normalize to sum(f)==1
136 fact = fact/I; # normalize to sum(f)==1
140 factr = factr/I; # normalize to sum(f)==1
137 factr = factr/I; # normalize to sum(f)==1
141 #plot(thetat,fact,'r'); hold on;
138 #plt.figure()
142 #plot(thetar,factr,'k.'); hold off;
139 #plt.plot(thetat,fact,'r');
140 #plt.plot(thetar,factr,'k.');
143 #xlim([min(thetat) max(thetat)]);
141 #xlim([min(thetat) max(thetat)]);
144
142
145 #x = np.linspace(thetat.min(), thetat.max) ????
143 #x = np.linspace(thetat.min(), thetat.max) ????
@@ -197,13 +195,14
197 temp4 = (sigma_noise * (temp1 + temp2))/np.sqrt(2.0)
195 temp4 = (sigma_noise * (temp1 + temp2))/np.sqrt(2.0)
198
196
199 nz = np.multiply(temp4, temp3)
197 nz = np.multiply(temp4, temp3)
200
198
199 #---------------------- Eliminar esto:------------------------------------------
200 #nz = ((abs(np.multiply(temp4, temp3)) > 0).astype(int))
201 #nz = ((abs(np.dot(temp4, temp3)) > 0).astype(int))
201 #nz = ((abs(np.dot(temp4, temp3)) > 0).astype(int))
202
203
204 #nz = np.dot(np.dot(sigma_noise, (temp1 + temp2)/math.sqrt(2), temp3 ));
202 #nz = np.dot(np.dot(sigma_noise, (temp1 + temp2)/math.sqrt(2), temp3 ));
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))));
203 #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))));
206
204 #--------------------------------------------------------------------------------
205
207 Unz = U + nz;
206 Unz = U + nz;
208 Rnz = np.dot(Unz.T.conj(),Unz); # the noisy version of R
207 Rnz = np.dot(Unz.T.conj(),Unz); # the noisy version of R
209 plt.figure(3);
208 plt.figure(3);
@@ -230,16 +229,17
230
229
231 f_capon = np.zeros(shape=(Nr,1));
230 f_capon = np.zeros(shape=(Nr,1));
232
231
233 #tic_capon = time.time();
232 tic_capon = time.time();
234
233
235 for i in range(0, thetar.size):
234 for i in range(0, thetar.size):
236 th = thetar[i];
235 th = thetar[i];
237 w = np.exp(1j*k*np.dot(r,np.sin(th)));
236 w = np.exp(1j*k*np.dot(r,np.sin(th)));
238 f_capon[i] = np.divide(1, ( np.dot( w.T.conj(), (linalg.solve(Rnz,w)) ) ).real)
237 f_capon[i] = np.divide(1, ( np.dot( w.T.conj(), (linalg.solve(Rnz,w)) ) ).real)
239
238
240 #toc_capon = time.time()
239
241
240 toc_capon = time.time()
242 # elapsed_time_capon = toc_capon - tic_capon;
241
242 elapsed_time_capon = toc_capon - tic_capon;
243
243
244 f_capon = f_capon.real; # get rid of numerical imaginary noise
244 f_capon = f_capon.real; # get rid of numerical imaginary noise
245
245
@@ -278,10 +278,15
278 sigma = 1; # set to 1 because the difference is accounted for in G
278 sigma = 1; # set to 1 because the difference is accounted for in G
279
279
280 ##### ADD *10 for consistency with old model, NEED TO VERIFY THIS!!!!? line below
280 ##### ADD *10 for consistency with old model, NEED TO VERIFY THIS!!!!? line below
281 G = np.linalg.norm(g-gnz)**2; # pretend we know in advance the actual value of chi^2
281 G = np.linalg.norm(g-gnz)**2 ; # pretend we know in advance the actual value of chi^2
282
283 tic_maxent = time.time();
282
284
283 lambda0 = 1e-5*np.ones(shape=(M,1)); # initial condition (can be set to anything)
285 lambda0 = 1e-5*np.ones(shape=(M,1)); # initial condition (can be set to anything)
284
286
287 toc_maxent = time.time()
288 elapsed_time_maxent = toc_maxent - tic_maxent;
289
285 # Whitened solution
290 # Whitened solution
286 def myfun(lambda1):
291 def myfun(lambda1):
287 return y_hysell96(lambda1,gnz,sigma,F,G,Hr);
292 return y_hysell96(lambda1,gnz,sigma,F,G,Hr);
@@ -304,24 +309,28
304 es = np.dot(Hr, f_maxent) - gnz; # should be same as ep
309 es = np.dot(Hr, f_maxent) - gnz; # should be same as ep
305 chi2 = np.sum((es/sigma)**2);
310 chi2 = np.sum((es/sigma)**2);
306
311
307
312
308
309 # CS inversion using irls ########################
313 # CS inversion using irls ########################
310
314
311 # (Use Nr, thetar, gnz, and Hr from MaxEnt above)
315 # (Use Nr, thetar, gnz, and Hr from MaxEnt above)
312
316
313 #Psi = deb4_basis(Nr); ###### REPLACED BY LINE BELOW (?)
317 Psi = deb4_basis(Nr); ###### REPLACED BY LINE BELOW (?)
314
318
315 wavelet1 = pywt.Wavelet('db4')
319 # REMOVE THIS?--------------------------------
316 Phi, Psi, x = wavelet1.wavefun(level=3)
320 #wavelet1 = pywt.Wavelet('db4')
321 #Phi, Psi, x = wavelet1.wavefun(level=3)
322 # --------------------------------------------
317
323
318 # add "sum to 1" constraint
324 # add "sum to 1" constraint
319 H2 = np.concatenate( (Hr, np.ones(shape=(1,Nr))), axis=0 );
325 H2 = np.concatenate( (Hr, np.ones(shape=(1,Nr))), axis=0 );
320 N_temp = np.array([[Nr/Nt]]);
326 N_temp = np.array([[Nr/Nt]]);
321 g2 = np.concatenate( (gnz, N_temp), axis=0 );
327 g2 = np.concatenate( (gnz, N_temp), axis=0 );
322
328 H2 = H2.T.conj();
323
329
324 #s = irls_dn2(H2*Psi,g2,0.5,G);
330 print 'H2 shape', H2.shape
331 print 'Psi shape', Psi.shape
332
333 s = irls_dn2(H2*Psi,g2,0.5,G);
325 # f_cs = Psi*s;
334 # f_cs = Psi*s;
326 #
335 #
327 # # plot
336 # # plot
@@ -331,100 +340,129
331 # hold off;
340 # hold off;
332
341
333
342
334 # # # Scaling and shifting
343 # # # Scaling and shifting
335 # # # Only necessary for capon solution
344 # # # Only necessary for capon solution
336
337
338 f_capon = f_capon/np.max(f_capon)*np.max(fact);
339
340
341 ### analyze stuff ######################
342 # calculate MSE
343 rmse_fourier = np.sqrt(np.mean((f_fourier - factr)**2));
344 rmse_capon = np.sqrt(np.mean((f_capon - factr)**2));
345 rmse_maxent = np.sqrt(np.mean((f_maxent - factr)**2));
346 #rmse_cs = np.sqrt(np.mean((f_cs - factr).^2));
347
345
348 relrmse_fourier = rmse_fourier / np.linalg.norm(fact);
346
349 relrmse_capon = rmse_capon / np.linalg.norm(fact);
347 f_capon = f_capon/np.max(f_capon)*np.max(fact);
350 relrmse_maxent = rmse_maxent / np.linalg.norm(fact);
348
351 #relrmse_cs = rmse_cs / np.norm(fact);
349
352
350 ### analyze stuff ######################
353 factr = factr.T.conj()
351 # calculate MSE
354
352 rmse_fourier = np.sqrt(np.mean((f_fourier - factr)**2));
355 # calculate correlation
353 rmse_capon = np.sqrt(np.mean((f_capon - factr)**2));
356 corr_fourier = np.dot(f_fourier,factr) / (np.linalg.norm(f_fourier)*np.linalg.norm(factr));
354 rmse_maxent = np.sqrt(np.mean((f_maxent - factr)**2));
357 corr_capon = np.dot(f_capon,factr) / (np.linalg.norm(f_capon)*np.linalg.norm(factr));
355 #rmse_cs = np.sqrt(np.mean((f_cs - factr).^2));
358 corr_maxent = np.dot(f_maxent,factr) / (np.linalg.norm(f_maxent)*np.linalg.norm(factr));
356
359 #corr_cs = np.dot(f_cs,factr) / (norm(f_cs)*norm(factr));
357
360
358 relrmse_fourier = rmse_fourier / np.linalg.norm(fact);
361 # calculate centered correlation
359 relrmse_capon = rmse_capon / np.linalg.norm(fact);
362 f0 = factr - np.mean(factr);
360 relrmse_maxent = rmse_maxent / np.linalg.norm(fact);
363 f1 = f_fourier - np.mean(f_fourier);
361 #relrmse_cs = rmse_cs / np.norm(fact);
364 corrc_fourier = np.dot(f0,f1) / (np.linalg.norm(f0)*np.linalg.norm(f1));
362
365 f1 = f_capon - np.mean(f_capon);
363 # To be able to perform dot product (align matrices) done below within the dot calculations
366 corrc_capon = np.dot(f0,f1) / (np.linalg.norm(f0)*np.linalg.norm(f1));
364
367 f1 = f_maxent - np.mean(f_maxent);
365
368 corrc_maxent = np.dot(f0,f1) / (np.linalg.norm(f0)*np.linalg.norm(f1));
366 #f_fourier = f_fourier.T.conj()
369 #f1 = f_cs - mean(f_cs);
367 #f_capon = f_capon.T.conj()
370 #corrc_cs = dot(f0,f1) / (norm(f0)*norm(f1));
368 #f_maxent = f_maxent.T.conj()
371
369
372
370 #factr = factr.T.conj()
373
371
374 # # # plot stuff #########################
372 # calculate correlation
375
373
376 #---- Capon----
374 corr_fourier = np.dot(f_fourier.T.conj(),factr) / (np.linalg.norm(f_fourier)*np.linalg.norm(factr));
377 plt.figure(4)
375 corr_capon = np.dot(f_capon.T.conj(),factr) / (np.linalg.norm(f_capon)*np.linalg.norm(factr));
378 plt.subplot(2, 1, 1)
376 corr_maxent = np.dot(f_maxent.T.conj(),factr) / (np.linalg.norm(f_maxent)*np.linalg.norm(factr));
379 plt.plot(180/math.pi*thetar, f_capon, 'r', label='Capon')
377 #corr_cs = np.dot(f_cs,factr) / (norm(f_cs)*norm(factr));
380 plt.plot(180/math.pi*thetat,fact, 'k--', label='Truth')
378
381 plt.ylabel('Power (arbitrary units)')
379
382 plt.legend(loc='upper right')
380 # calculate centered correlation
383
381 f0 = factr - np.mean(factr);
384 # formatting y-axis
382 f1 = f_fourier - np.mean(f_fourier);
385 locs,labels = plt.yticks()
383
386 plt.yticks(locs, map(lambda x: "%.1f" % x, locs*1e4))
384 corrc_fourier = np.dot(f0.T.conj(),f1) / (np.linalg.norm(f0)*np.linalg.norm(f1));
387 plt.text(0.0, 1.01, '1e-4', fontsize=10, transform = plt.gca().transAxes)
385 f1 = f_capon - np.mean(f_capon);
388
386 corrc_capon = np.dot(f0.T.conj(),f1) / (np.linalg.norm(f0)*np.linalg.norm(f1));
389
387 f1 = f_maxent - np.mean(f_maxent);
390 #---- MaxEnt----
388 corrc_maxent = np.dot(f0.T.conj(),f1) / (np.linalg.norm(f0)*np.linalg.norm(f1));
391 plt.subplot(2, 1, 2)
389 #f1 = f_cs - mean(f_cs);
392 plt.plot(180/math.pi*thetar, f_maxent, 'r', label='MaxEnt')
390 #corrc_cs = dot(f0,f1) / (norm(f0)*norm(f1));
393 plt.plot(180/math.pi*thetat,fact, 'k--', label='Truth')
391
394 plt.ylabel('Power (arbitrary units)')
392
395 plt.legend(loc='upper right')
393
396
394 # # # plot stuff #########################
397 # formatting y-axis
395
398 locs,labels = plt.yticks()
396 #---- Capon----
399 plt.yticks(locs, map(lambda x: "%.1f" % x, locs*1e4))
397 plt.figure(4)
400 plt.text(0.0, 1.01, '1e-4', fontsize=10, transform = plt.gca().transAxes)
398 plt.subplot(2, 1, 1)
401
399 plt.plot(180/math.pi*thetar, f_capon, 'r', label='Capon')
402 plt.show()
400 plt.plot(180/math.pi*thetat,fact, 'k--', label='Truth')
403
401 plt.ylabel('Power (arbitrary units)')
404
402 plt.legend(loc='upper right')
405 # # subplot(3,1,2);
403
406 # # plot(180/pi*thetar,f_maxent,'r-');
404 # formatting y-axis
407 # # hold on;
405 locs,labels = plt.yticks()
408 # # plot(180/pi*thetat,fact,'k--');
406 plt.yticks(locs, map(lambda x: "%.1f" % x, locs*1e4))
409 # # hold off;
407 plt.text(0.0, 1.01, '1e-4', fontsize=10, transform = plt.gca().transAxes)
410 # # ylim([min(f_cs) 1.1*max(fact)]);
408
411 # # # title(sprintf('rel. RMSE: #.2e\tCorr: #.3f Corrc: #.3f', relrmse_maxent, corr_maxent, corrc_maxent));
409
412 # # ylabel({'Power';'(arbitrary units)'})
410 #---- MaxEnt----
413 # # # title 'Maximum Entropy Method'
411 plt.subplot(2, 1, 2)
414 # # legend('MaxEnt','Truth');
412 plt.plot(180/math.pi*thetar, f_maxent, 'r', label='MaxEnt')
415 # #
413 plt.plot(180/math.pi*thetat,fact, 'k--', label='Truth')
416 # # subplot(3,1,3);
414 plt.ylabel('Power (arbitrary units)')
417 # # plot(180/pi*thetar,f_cs,'r-');
415 plt.legend(loc='upper right')
418 # # hold on;
416
419 # # plot(180/pi*thetat,fact,'k--');
417 # formatting y-axis
420 # # hold off;
418 locs,labels = plt.yticks()
421 # # ylim([min(f_cs) 1.1*max(fact)]);
419 plt.yticks(locs, map(lambda x: "%.1f" % x, locs*1e4))
422 # # # title(sprintf('rel. RMSE: #.2e\tCorr: #.3f Corrc: #.3f', relrmse_cs, corr_cs, corrc_cs));
420 plt.text(0.0, 1.01, '1e-4', fontsize=10, transform = plt.gca().transAxes)
423 # # # title 'Compressed Sensing - Debauchies Wavelets'
421
424 # # xlabel 'Degrees'
422 plt.show()
425 # # ylabel({'Power';'(arbitrary units)'})
423
426 # # legend('Comp. Sens.','Truth');
424
427 # #
425 # # PLOT PARA COMPRESSED SENSING
428 # # # set(gcf,'Position',[749 143 528 881]); # CSL
426 # #
429 # # # set(gcf,'Position',[885 -21 528 673]); # macbook
427 # # subplot(3,1,3);
430 # # pause(0.01); No newline at end of file
428 # # plot(180/pi*thetar,f_cs,'r-');
429 # # hold on;
430 # # plot(180/pi*thetat,fact,'k--');
431 # # hold off;
432 # # ylim([min(f_cs) 1.1*max(fact)]);
433 # # # title(sprintf('rel. RMSE: #.2e\tCorr: #.3f Corrc: #.3f', relrmse_cs, corr_cs, corrc_cs));
434 # # # title 'Compressed Sensing - Debauchies Wavelets'
435 # # xlabel 'Degrees'
436 # # ylabel({'Power';'(arbitrary units)'})
437 # # legend('Comp. Sens.','Truth');
438 # #
439 # # # set(gcf,'Position',[749 143 528 881]); # CSL
440 # # # set(gcf,'Position',[885 -21 528 673]); # macbook
441 # # pause(0.01);
442
443
444 # # Store Results
445 corr[0, snri, Ni] = corr_fourier;
446 corr[1, snri, Ni] = corr_capon;
447 corr[2, snri, Ni] = corr_maxent;
448 #corr[3, snri, Ni] = corr_cs;
449
450 rmse[0,snri,Ni] = relrmse_fourier;
451 rmse[1,snri,Ni] = relrmse_capon;
452 rmse[2,snri,Ni] = relrmse_maxent;
453 #rmse[3,snri,Ni] = relrmse_cs;
454
455 corrc[0,snri,Ni] = corrc_fourier;
456 corrc[1,snri,Ni] = corrc_capon;
457 corrc[2,snri,Ni] = corrc_maxent;
458 #corrc[3,snri,Ni] = corrc_cs;
459
460
461 print 'Capon:\t', elapsed_time_capon, 'sec';
462 print 'Maxent:\t',elapsed_time_maxent, 'sec';
463 #print 'CS:\t%3.3f sec\n',elapsed_time_cs;
464
465 print (NN*(snri+1) + Ni), '/', (SNRdBvec.size*NN);
466
467 print corr
468 No newline at end of file
General Comments 0
You need to be logged in to leave comments. Login now