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