##// END OF EJS Templates
jro_soft » compressed_sensing
RSS
Implementation of voltages simulation for estimating correlations (some things need to be fixed)
yamaro -
r25:26
parent child
Show More
Show file before | Show file after | Hide comments
@@ -1,527 +1,572
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 No newline at end of file
4 # Original MATLAB code developed by Brian Harding
5 # Rewritten in Python by Yolian Amaro No newline at end of file
5 # Rewritten in Python by Yolian Amaro
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 No newline at end of file
10
11 import time No newline at end of file
11 import time
12 import matplotlib.pyplot as plt No newline at end of file
12 import matplotlib.pyplot as plt
13 from scipy.optimize import root No newline at end of file
13 from scipy.optimize import root
14 from scipy.optimize import newton_krylov No newline at end of file
14 from scipy.stats import nanmean No newline at end of file
15 from scipy.stats import nanmean
15 No newline at end of file
16
16 from y_hysell96 import * No newline at end of file
17 from y_hysell96 import *
17 from deb4_basis import * No newline at end of file
18 from deb4_basis import *
18 from modelf import * No newline at end of file
19 from modelf import *
19 from irls_dn2 import * No newline at end of file
20 from irls_dn2 import *
20 No newline at end of file
21
21 No newline at end of file
22
22 #------------------------------------------------------------------------------------------------- No newline at end of file
23 #-------------------------------------------------------------------------------------------------
23 # Set parameters No newline at end of file
24 # Set parameters
24 #------------------------------------------------------------------------------------------------- No newline at end of file
25 #-------------------------------------------------------------------------------------------------
25 No newline at end of file
26
26 ## Calculate Forward Model No newline at end of file
27 ## Calculate Forward Model
27 lambda1 = 6.0 No newline at end of file
28 lambda1 = 6.0
28 k = 2*np.pi/lambda1 No newline at end of file
29 k = 2*np.pi/lambda1
29 No newline at end of file
30
30 ## Magnetic Declination No newline at end of file
31 ## Magnetic Declination
31 dec = -1.24 No newline at end of file
32 dec = -1.24
32 No newline at end of file
33
33 ## Loads Jicamarca antenna positions No newline at end of file
34 ## Loads Jicamarca antenna positions
34 antpos = np.loadtxt("antpos.txt", comments="#", delimiter=";", unpack=False) No newline at end of file
35 antpos = np.loadtxt("antpos.txt", comments="#", delimiter=";", unpack=False)
35 No newline at end of file
36
36 ## rx and ry -- for plotting purposes No newline at end of file
37 ## rx and ry -- for plotting purposes
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
38 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
39 ry = np.array( [[127.5000], [91.5000], [91.5000], [55.5000], [-19.5000], [-127.5000], [-127.5000], [-322.2940]] )
39 No newline at end of file
40
40 ## Plot of antenna positions No newline at end of file
41 ## Plot of antenna positions
41 plt.figure(1) No newline at end of file
42 plt.figure(1)
42 plt.plot(rx, ry, 'ro') No newline at end of file
43 plt.plot(rx, ry, 'ro')
43 plt.title("Jicamarca Antenna Positions") No newline at end of file
44 plt.title("Jicamarca Antenna Positions")
44 plt.xlabel("rx") No newline at end of file
45 plt.xlabel("rx")
45 plt.ylabel("ry") No newline at end of file
46 plt.ylabel("ry")
46 plt.grid(True) No newline at end of file
47 plt.grid(True)
47 plt.plot(rx, rx, 'b-') No newline at end of file
48 plt.plot(rx, rx, 'b-')
48 plt.draw() No newline at end of file
49 plt.draw()
49 No newline at end of file
50
50 ## Jicamarca is nominally at a 45 degree angle No newline at end of file
51 ## Jicamarca is nominally at a 45 degree angle
51 theta = 45 - dec; No newline at end of file
52 theta = 45 - dec;
52 No newline at end of file
53
53 ## Rotation matrix from antenna coord to magnetic coord (East North) -> 45 degrees No newline at end of file
54 ## Rotation matrix from antenna coord to magnetic coord (East North) -> 45 degrees
54 theta_rad = np.radians(theta) # trig functions take radians as argument No newline at end of file
55 theta_rad = np.radians(theta) # trig functions take radians as argument
55 val1 = float( np.cos(theta_rad) )
56 val1 = float( np.cos(theta_rad) )
No newline at end of file
57 val2 = float( -1*np.sin(theta_rad))
56 val2 = float( np.sin(theta_rad) )
No newline at end of file
No newline at end of file
58 val3 = float( np.sin(theta_rad) ) No newline at end of file
57 val3 = float( -1*np.sin(theta_rad)) No newline at end of file
58 val4 = float( np.cos(theta_rad) ) No newline at end of file
59 val4 = float( np.cos(theta_rad) )
59 No newline at end of file
60
60 # Rotation matrix from antenna coord to magnetic coord (East North)
61 # Rotation matrix from antenna coord to magnetic coord (East North)
No newline at end of file
62 R = np.array( [[val1, val2], [val3, val4]] ) No newline at end of file
61 R = np.array( [[val1, val3], [val2, val4]] ) No newline at end of file
62 No newline at end of file
63
63 # Rotate antenna positions to magnetic coord. No newline at end of file
64 # Rotate antenna positions to magnetic coord.
64 AR = np.dot(R.T, antpos) No newline at end of file
65 AR = np.dot(R.T, antpos)
65
66
No newline at end of file
67 plt.plot(AR[0], 0*AR[1], 'yo') # antenna positions only considering magnetic East-West component (x, 0) No newline at end of file
66 plt.plot(AR[0], 0*AR[1], 'yo') No newline at end of file
67 plt.title("Jicamarca Antenna Positions") No newline at end of file
68 plt.title("Jicamarca Antenna Positions")
68 plt.xlabel("rx") No newline at end of file
69 plt.xlabel("rx")
69 plt.ylabel("ry") No newline at end of file
70 plt.ylabel("ry")
70 plt.grid(True) No newline at end of file
71 plt.grid(True)
71 plt.plot(rx, 0*rx, 'g-') No newline at end of file
72 plt.plot(rx, 0*rx, 'g-')
72 plt.draw() No newline at end of file
73 plt.draw()
73
74
No newline at end of file
74 print antpos
No newline at end of file
75 print "\n"
No newline at end of file
76 print AR No newline at end of file
77 No newline at end of file
75
78 # Only take the East component No newline at end of file
76 # Only take the East component
79 r = AR[0,:] No newline at end of file
77 r = AR[0,:]
80 r.sort()
78 r.sort()
No newline at end of file
81
No newline at end of file
82 # ---NEW---
No newline at end of file
83 r1 = AR[1,:] # longitudes / distances
No newline at end of file
84 d = 40 # separation of antennas in meters
No newline at end of file
85 theta_w = 15 # angle between antenna and point of observation in degrees
No newline at end of file
86 # arg = k*(r + d*np.cos(theta_w)) # this is the exponent argument in the corr integral
No newline at end of file
87 # -------- No newline at end of file
88 No newline at end of file
79
89 No newline at end of file
80
90 # Truth model (high and low resolution) No newline at end of file
81 # Truth model (high and low resolution)
91 Nt = (1024.0)*(16.0) # number of pixels in truth image: high resolution No newline at end of file
82 Nt = (1024.0)*(16.0) # number of pixels in truth image: high resolution
92 thbound = 9.0/180*np.pi # the width of the domain in angle space No newline at end of file
83 thbound = 9.0/180*np.pi # the width of the domain in angle space
93 thetat = np.linspace(-thbound, thbound,Nt) # image domain No newline at end of file
84 thetat = np.linspace(-thbound, thbound,Nt) # image domain
94 thetat = thetat.T # transpose No newline at end of file
85 thetat = thetat.T # transpose
95 Nr = (256.0) # number of pixels in reconstructed image: low res No newline at end of file
86 Nr = (256.0) # number of pixels in reconstructed image: low res
96 thetar = np.linspace(-thbound, thbound,Nr) # reconstruction domain No newline at end of file
87 thetar = np.linspace(-thbound, thbound,Nr) # reconstruction domain
97 thetar = thetar.T # transpose No newline at end of file
88 thetar = thetar.T # transpose
98 No newline at end of file
89
99 No newline at end of file
90
100 #------------------------------------------------------------------------------------------------- No newline at end of file
91 #-------------------------------------------------------------------------------------------------
101 # Model for f: Gaussian(s) with amplitudes a, centers mu, widths sig, and background constant b. No newline at end of file
92 # Model for f: Gaussian(s) with amplitudes a, centers mu, widths sig, and background constant b.
102 #------------------------------------------------------------------------------------------------- No newline at end of file
93 #-------------------------------------------------------------------------------------------------
103 No newline at end of file
94
104 # Triple Gaussian No newline at end of file
95 # Triple Gaussian
105 # a = np.array([3, 5, 2]); No newline at end of file
96 # a = np.array([3, 5, 2]);
106 # mu = np.array([-5.0/180*np.pi, 2.0/180*np.pi, 7.0/180*np.pi]); No newline at end of file
97 # mu = np.array([-5.0/180*np.pi, 2.0/180*np.pi, 7.0/180*np.pi]);
107 # sig = np.array([2.0/180*np.pi, 1.5/180*np.pi, 0.3/180*np.pi]); No newline at end of file
98 # sig = np.array([2.0/180*np.pi, 1.5/180*np.pi, 0.3/180*np.pi]);
108 # b = 0; # background No newline at end of file
99 # b = 0; # background
109 No newline at end of file
100
110 # Double Gaussian No newline at end of file
101 # Double Gaussian
111 # a = np.array([3, 5]); No newline at end of file
102 # a = np.array([3, 5]);
112 # mu = np.array([-5.0/180*np.pi, 2.0/180*np.pi]); No newline at end of file
103 # mu = np.array([-5.0/180*np.pi, 2.0/180*np.pi]);
113 # sig = np.array([2.0/180*np.pi, 1.5/180*np.pi]); No newline at end of file
104 # sig = np.array([2.0/180*np.pi, 1.5/180*np.pi]);
114 # b = 0; # background No newline at end of file
105 # b = 0; # background
115 No newline at end of file
106
116 # Single Gaussian No newline at end of file
107 # Single Gaussian
117 a = np.array( [3] ) # amplitude(s) No newline at end of file
108 a = np.array( [3] ) # amplitude(s)
118 mu = np.array( [-3.0/180*np.pi] ) # center No newline at end of file
109 mu = np.array( [-3.0/180*np.pi] ) # center
119 sig = np.array( [2.0/180*np.pi] ) # width No newline at end of file
110 sig = np.array( [2.0/180*np.pi] ) # width
120 b = 0 # background constant No newline at end of file
111 b = 0 # background constant
121 No newline at end of file
112
122 # Empty matrices for factors No newline at end of file
113 # Empty matrices for factors
123 fact = np.zeros(shape=(Nt,1)) No newline at end of file
114 fact = np.zeros(shape=(Nt,1))
124 factr = np.zeros(shape=(Nr,1)) No newline at end of file
115 factr = np.zeros(shape=(Nr,1))
125
116
No newline at end of file
117 # Constructing Gaussian model - brightness function (?) No newline at end of file
126 # Constructing Gaussian model No newline at end of file
127 for i in range(0, a.size): No newline at end of file
118 for i in range(0, a.size):
128 temp = (-(thetat-mu[i])**2/(sig[i]**2)) No newline at end of file
119 temp = (-(thetat-mu[i])**2/(sig[i]**2))
129 tempr = (-(thetar-mu[i])**2/(sig[i]**2)) No newline at end of file
120 tempr = (-(thetar-mu[i])**2/(sig[i]**2))
130 for j in range(0, temp.size): No newline at end of file
121 for j in range(0, temp.size):
131 fact[j] = fact[j] + a[i]*np.exp(temp[j]); No newline at end of file
122 fact[j] = fact[j] + a[i]*np.exp(temp[j]);
132 for m in range(0, tempr.size): No newline at end of file
123 for m in range(0, tempr.size):
133 factr[m] = factr[m] + a[i]*np.exp(tempr[m]); No newline at end of file
124 factr[m] = factr[m] + a[i]*np.exp(tempr[m]);
134 fact = fact + b; No newline at end of file
125 fact = fact + b;
135 factr = factr + b; No newline at end of file
126 factr = factr + b;
136 No newline at end of file
127
137 No newline at end of file
128
138 # #------------------------------------------------------------------------------------------------- No newline at end of file
129 # #-------------------------------------------------------------------------------------------------
139 # # Model for f: Square pulse No newline at end of file
130 # # Model for f: Square pulse
140 # #------------------------------------------------------------------------------------------------- No newline at end of file
131 # #-------------------------------------------------------------------------------------------------
141 # for j in range(0, fact.size): No newline at end of file
132 # for j in range(0, fact.size):
142 # if (theta > -5.0/180*np.pi and theta < 2.0/180*np.pi): No newline at end of file
133 # if (theta > -5.0/180*np.pi and theta < 2.0/180*np.pi):
143 # fact[j] = 0 No newline at end of file
134 # fact[j] = 0
144 # else: No newline at end of file
135 # else:
145 # fact[j] = 1 No newline at end of file
136 # fact[j] = 1
146 # for k in range(0, factr.size): No newline at end of file
137 # for k in range(0, factr.size):
147 # if (thetar[k] > -5.0/180*np.pi and thetar[k] < 2/180*np.pi): No newline at end of file
138 # if (thetar[k] > -5.0/180*np.pi and thetar[k] < 2/180*np.pi):
148 # factr[k] = 0 No newline at end of file
139 # factr[k] = 0
149 # else: No newline at end of file
140 # else:
150 # factr[k] = 1 No newline at end of file
141 # factr[k] = 1
151 No newline at end of file
142
152 # #------------------------------------------------------------------------------------------------- No newline at end of file
143 # #-------------------------------------------------------------------------------------------------
153 # # Model for f: Triangle pulse No newline at end of file
144 # # Model for f: Triangle pulse
154 # #------------------------------------------------------------------------------------------------- No newline at end of file
145 # #-------------------------------------------------------------------------------------------------
155 # mu = -1.0/180*np.pi; No newline at end of file
146 # mu = -1.0/180*np.pi;
156 # sig = 5.0/180*np.pi; No newline at end of file
147 # sig = 5.0/180*np.pi;
157 # wind1 = theta > mu-sig and theta < mu; No newline at end of file
148 # wind1 = theta > mu-sig and theta < mu;
158 # wind2 = theta < mu+sig and theta > mu; No newline at end of file
149 # wind2 = theta < mu+sig and theta > mu;
159 # fact = wind1 * (theta - (mu - sig)); No newline at end of file
150 # fact = wind1 * (theta - (mu - sig));
160 # factr = wind1 * (thetar - (mu - sig)); No newline at end of file
151 # factr = wind1 * (thetar - (mu - sig));
161 # fact = fact + wind2 * (-(theta-(mu+sig))); No newline at end of file
152 # fact = fact + wind2 * (-(theta-(mu+sig)));
162 # factr = factr + wind2 * (-(thetar-(mu+sig))); No newline at end of file
153 # factr = factr + wind2 * (-(thetar-(mu+sig)));
163 # No newline at end of file
154 #
164 # fact = np.divide(fact, (np.sum(fact)*2*thbound/Nt)); # normalize to integral(f)==1 No newline at end of file
155 # fact = np.divide(fact, (np.sum(fact)*2*thbound/Nt)); # normalize to integral(f)==1
165 No newline at end of file
156
166 No newline at end of file
157
167 I = sum(fact)[0]; No newline at end of file
158 I = sum(fact)[0];
168 fact = fact/I; # normalize to sum(f)==1 No newline at end of file
159 fact = fact/I; # normalize to sum(f)==1
169 factr = factr/I; # normalize to sum(f)==1 No newline at end of file
160 factr = factr/I; # normalize to sum(f)==1
170 No newline at end of file
161
171 # Plot Gaussian pulse(s) No newline at end of file
162 # Plot Gaussian pulse(s)
172 plt.figure(2)
163 plt.figure(2)
No newline at end of file
164 plt.plot(thetat, fact, 'bo') No newline at end of file
173 plt.plot(thetat, fact, 'r--') No newline at end of file
174 plt.plot(thetar, factr, 'ro')
165 plt.plot(thetar, factr, 'ro')
No newline at end of file
166 plt.title("Image and Reconstruction Domain") No newline at end of file
175 plt.title("Truth models") No newline at end of file
176 plt.ylabel("Amplitude")
167 plt.ylabel("Amplitude")
No newline at end of file
168 plt.legend(('High res (img domain)','Low res (reconstruction domain)'), loc='upper right') No newline at end of file
177 plt.legend(('High res','Low res'), loc='upper right') No newline at end of file
178 plt.draw() No newline at end of file
169 plt.draw()
179 No newline at end of file
170
180 No newline at end of file
171
181 #------------------------------------------------------------------------------------------------- No newline at end of file
172 #-------------------------------------------------------------------------------------------------
182 # Control the type and number of inversions with: No newline at end of file
173 # Control the type and number of inversions with:
183 # SNRdBvec: the SNRs that will be used. No newline at end of file
174 # SNRdBvec: the SNRs that will be used.
184 # NN: the number of trials for each SNR No newline at end of file
175 # NN: the number of trials for each SNR
185 #------------------------------------------------------------------------------------------------- No newline at end of file
176 #-------------------------------------------------------------------------------------------------
186 No newline at end of file
177
187 #SNRdBvec = np.linspace(5,20,10); No newline at end of file
178 #SNRdBvec = np.linspace(5,20,10);
188 SNRdBvec = np.array([15]); # 15 dB No newline at end of file
179 SNRdBvec = np.array([15]); # 15 dB
189 NN = 1; # number of trials at each SNR No newline at end of file
180 NN = 1; # number of trials at each SNR
190 No newline at end of file
181
191 # Statistics simulation (correlation, root mean square error) No newline at end of file
182 # Statistics simulation (correlation, root mean square error)
192 corr = np.zeros(shape=(4,SNRdBvec.size,NN)); # (method, SNR, trial) No newline at end of file
183 corr = np.zeros(shape=(4,SNRdBvec.size,NN)); # (method, SNR, trial)
193 corrc = np.zeros(shape=(4,SNRdBvec.size,NN)); # (method, SNR, trial) No newline at end of file
184 corrc = np.zeros(shape=(4,SNRdBvec.size,NN)); # (method, SNR, trial)
194 rmse = np.zeros(shape=(4,SNRdBvec.size,NN)); # (method, SNR, trial) No newline at end of file
185 rmse = np.zeros(shape=(4,SNRdBvec.size,NN)); # (method, SNR, trial)
195 No newline at end of file
186
196 # For each SNR and trial No newline at end of file
187 # For each SNR and trial
197 for snri in range(0, SNRdBvec.size): No newline at end of file
188 for snri in range(0, SNRdBvec.size):
198 SNRdB = SNRdBvec[snri]; No newline at end of file
189 SNRdB = SNRdBvec[snri];
199 SNR = 10**(SNRdB/10.0); No newline at end of file
190 SNR = 10**(SNRdB/10.0);
200 No newline at end of file
191
201 for Ni in range(0, NN): No newline at end of file
192 for Ni in range(0, NN):
193
No newline at end of file
194 # -------------------------------------- Voltages estimation----------------------------------------
No newline at end of file
195 v = np.zeros(shape=(r.size,1), dtype=object); # voltages matrix
No newline at end of file
196 R = np.zeros(shape=(r.size, r.size), dtype='complex128'); # corr matrix
No newline at end of file
197 mu1, sigma1 = 0, 1 # mean and standard deviation
No newline at end of file
198 num = 1
No newline at end of file
199
No newline at end of file
200 for t in range (0, 10):
No newline at end of file
201 num = t+1
No newline at end of file
202 print num
No newline at end of file
203 n1 = np.random.normal(mu1, sigma1, [Nt,1]) # random variable
No newline at end of file
204 n2 = np.random.normal(mu1, sigma1, [Nt,1]) # random variable
No newline at end of file
205
No newline at end of file
206 rho = (n1 + n2*1j)*np.sqrt(fact/2) # rho.size is Nt = 16834
No newline at end of file
207 # fact is size (Nt, 1) -> column vector
No newline at end of file
208 thetat = thetat.reshape(Nt,1) # thetat is now (Nt, 1)
No newline at end of file
209
No newline at end of file
210 for i in range(0, r.size):
No newline at end of file
211 v[i] = sum((rho)*np.exp(-1j*k*((r[i])*np.cos(thetat)))); # v is sum ( (Nt, 1)) = 1 value
No newline at end of file
212
No newline at end of file
213 # Computing correlations
No newline at end of file
214 for i1 in range(0, r.size):
No newline at end of file
215 for i2 in range(0,r.size):
No newline at end of file
216 R[i1,i2] = R[i1,i2] + (v[i1][0]*v[i2][0].conjugate())/ (np.sqrt(np.abs(v[i1][0])**2 * np.abs(v[i2][0])**2));
No newline at end of file
217
No newline at end of file
218 R = R / num
No newline at end of file
219 print "R from voltages"
No newline at end of file
220 print R
No newline at end of file
221
No newline at end of file
222 plt.figure(3);
No newline at end of file
223 plt.pcolor(abs(R));
No newline at end of file
224 plt.colorbar();
No newline at end of file
225 plt.title("Antenna Correlations")
No newline at end of file
226 plt.draw()
No newline at end of file
227 plt.show()
No newline at end of file
228
No newline at end of file
229 #--------------------------------------------------------------------------------------------------
No newline at end of file
230 No newline at end of file
202 # Calculate cross-correlation matrix (Fourier components of image) No newline at end of file
231 # Calculate cross-correlation matrix (Fourier components of image)
203 # This is an inefficient way to do this. No newline at end of file
232 # This is an inefficient way to do this.
204 No newline at end of file
233
234 # MATLAB CORR::: R(i1,i2) = fact.' * exp(j*k*(r(i1)-r(i2))*sin(thetat));
No newline at end of file
235 No newline at end of file
205 R = np.zeros(shape=(r.size, r.size), dtype=object); No newline at end of file
236 R = np.zeros(shape=(r.size, r.size), dtype=object);
206 No newline at end of file
237
238 thetat = thetat.reshape(Nt,) No newline at end of file
207 for i1 in range(0, r.size): No newline at end of file
239 for i1 in range(0, r.size):
208 for i2 in range(0,r.size):
240 for i2 in range(0,r.size):
No newline at end of file
241 R[i1,i2] = sum(np.dot(fact.T, np.exp(1j*k*((r[i1]-r[i2])*np.sin(thetat)))))
209 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
No newline at end of file
242
210 R[i1,i2] = sum(R[i1,i2]) No newline at end of file
No newline at end of file
243 print "R w/o voltages\n", R
No newline at end of file
244 # plt.figure(3);
No newline at end of file
245 # plt.pcolor(abs(R));
No newline at end of file
246 # plt.colorbar();
No newline at end of file
247 # plt.title("Antenna Correlations")
No newline at end of file
248 # plt.xlabel("Antenna")
No newline at end of file
249 # plt.ylabel("Antenna")
No newline at end of file
250 # plt.draw()
No newline at end of file
251 # plt.show()
No newline at end of file
252 No newline at end of file
211 No newline at end of file
253
212 # Add uncertainty No newline at end of file
254 # Add uncertainty
213 # This is an ad-hoc way of adding "noise". It models some combination of No newline at end of file
255 # This is an ad-hoc way of adding "noise". It models some combination of
214 # receiver noise and finite integration times. We could use a more No newline at end of file
256 # receiver noise and finite integration times. We could use a more
215 # advanced model (like in Yu et al 2000) in the future. No newline at end of file
257 # advanced model (like in Yu et al 2000) in the future.
216 No newline at end of file
258
217 # This is a way of adding noise while maintaining the No newline at end of file
259 # This is a way of adding noise while maintaining the
218 # positive-semi-definiteness of the matrix. No newline at end of file
260 # positive-semi-definiteness of the matrix.
219 No newline at end of file
261
220 U = linalg.cholesky(R.astype(complex), lower=False); # U'*U = R No newline at end of file
262 U = linalg.cholesky(R.astype(complex), lower=False); # U'*U = R
221 No newline at end of file
263
222 sigma_noise = (np.linalg.norm(U,'fro')/SNR); No newline at end of file
264 sigma_noise = (np.linalg.norm(U,'fro')/SNR);
223 No newline at end of file
265
224 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)); No newline at end of file
266 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));
225 No newline at end of file
267
226 Unz = U + nz; No newline at end of file
268 Unz = U + nz;
227 Rnz = np.dot(Unz.T.conj(),Unz); # the noisy version of R
269 Rnz = np.dot(Unz.T.conj(),Unz); # the noisy version of R
No newline at end of file
270 # plt.figure(3);
228 plt.figure(3);
No newline at end of file
No newline at end of file
271 # plt.pcolor(abs(R));
229 plt.pcolor(abs(R));
No newline at end of file
No newline at end of file
272 # plt.colorbar(); No newline at end of file
230 plt.colorbar(); No newline at end of file
231 No newline at end of file
273
232 #------------------------------------------------------------------------------------------------- No newline at end of file
274 #-------------------------------------------------------------------------------------------------
233 # Fourier Inversion No newline at end of file
275 # Fourier Inversion
234 #------------------------------------------------------------------------------------------------- No newline at end of file
276 #-------------------------------------------------------------------------------------------------
235 f_fourier = np.zeros(shape=(Nr,1), dtype=complex); No newline at end of file
277 f_fourier = np.zeros(shape=(Nr,1), dtype=complex);
236 No newline at end of file
278
237 for i in range(0, thetar.size): No newline at end of file
279 for i in range(0, thetar.size):
238 th = thetar[i]; No newline at end of file
280 th = thetar[i];
239 w = np.exp(1j*k*np.dot(r,np.sin(th)));
281 w = np.exp(1j*k*np.dot(r,np.sin(th)));
No newline at end of file
282 f_fourier[i] = np.dot(np.dot(w.T.conj(),U) , w); No newline at end of file
240 temp = np.dot(w.T.conj(),U)
No newline at end of file
241 f_fourier[i] = np.dot(temp, w); No newline at end of file
242 No newline at end of file
283
243 f_fourier = f_fourier.real; # get rid of numerical imaginary noise No newline at end of file
284 f_fourier = f_fourier.real; # get rid of numerical imaginary noise
244 No newline at end of file
285
245 No newline at end of file
286
246 #------------------------------------------------------------------------------------------------- No newline at end of file
287 #-------------------------------------------------------------------------------------------------
247 # Capon Inversion No newline at end of file
288 # Capon Inversion
248 #------------------------------------------------------------------------------------------------- No newline at end of file
289 #-------------------------------------------------------------------------------------------------
249 f_capon = np.zeros(shape=(Nr,1)); No newline at end of file
290 f_capon = np.zeros(shape=(Nr,1));
250 No newline at end of file
291
251 tic_capon = time.time(); No newline at end of file
292 tic_capon = time.time();
252 No newline at end of file
293
253 for i in range(0, thetar.size): No newline at end of file
294 for i in range(0, thetar.size):
254 th = thetar[i]; No newline at end of file
295 th = thetar[i];
255 w = np.exp(1j*k*np.dot(r,np.sin(th))); No newline at end of file
296 w = np.exp(1j*k*np.dot(r,np.sin(th)));
256 f_capon[i] = 1/ ( np.dot( w.T.conj(), (linalg.solve(Rnz,w)) ) ).real No newline at end of file
297 f_capon[i] = 1/ ( np.dot( w.T.conj(), (linalg.solve(Rnz,w)) ) ).real
257 No newline at end of file
298
258 toc_capon = time.time() No newline at end of file
299 toc_capon = time.time()
259 elapsed_time_capon = toc_capon - tic_capon; No newline at end of file
300 elapsed_time_capon = toc_capon - tic_capon;
260 No newline at end of file
301
261 f_capon = f_capon.real; # get rid of numerical imaginary noise No newline at end of file
302 f_capon = f_capon.real; # get rid of numerical imaginary noise
262 No newline at end of file
303
263 No newline at end of file
304
264 #------------------------------------------------------------------------------------------------- No newline at end of file
305 #-------------------------------------------------------------------------------------------------
265 # MaxEnt Inversion No newline at end of file
306 # MaxEnt Inversion
266 #------------------------------------------------------------------------------------------------- No newline at end of file
307 #-------------------------------------------------------------------------------------------------
267 No newline at end of file
308
268 # Create the appropriate sensing matrix (split into real and imaginary # parts) No newline at end of file
309 # Create the appropriate sensing matrix (split into real and imaginary # parts)
269 M = (r.size-1)*(r.size); No newline at end of file
310 M = (r.size-1)*(r.size);
270 Ht = np.zeros(shape=(M,Nt)); # "true" sensing matrix No newline at end of file
311 Ht = np.zeros(shape=(M,Nt)); # "true" sensing matrix
271 Hr = np.zeros(shape=(M,Nr)); # approximate sensing matrix for reconstruction No newline at end of file
312 Hr = np.zeros(shape=(M,Nr)); # approximate sensing matrix for reconstruction
272 No newline at end of file
313
273 # Need to re-index our measurements from matrix R into vector g No newline at end of file
314 # Need to re-index our measurements from matrix R into vector g
274 g = np.zeros(shape=(M,1)); No newline at end of file
315 g = np.zeros(shape=(M,1));
275 gnz = np.zeros(shape=(M,1)); # noisy version of g No newline at end of file
316 gnz = np.zeros(shape=(M,1)); # noisy version of g
276 No newline at end of file
317
277 # Triangular indexing to perform this re-indexing No newline at end of file
318 # Triangular indexing to perform this re-indexing
278 T = np.ones(shape=(r.size,r.size)); No newline at end of file
319 T = np.ones(shape=(r.size,r.size));
279 [i1v,i2v] = np.where(np.triu(T,1) > 0); # converts linear to triangular indexing No newline at end of file
320 [i1v,i2v] = np.where(np.triu(T,1) > 0); # converts linear to triangular indexing
280 No newline at end of file
321
281 # Build H No newline at end of file
322 # Build H
282 for i1 in range(0, r.size): No newline at end of file
323 for i1 in range(0, r.size):
283 for i2 in range(i1+1, r.size):
324 for i2 in range(i1+1, r.size):
No newline at end of file
325 idx = np.where(np.logical_and((i1==i1v), (i2==i2v)))[0]; # kind of awkward (1->27) No newline at end of file
284 idx = np.where(np.logical_and((i1==i1v), (i2==i2v)))[0]; # kind of awkward 1->27 No newline at end of file
285 idx1 = 2*idx; No newline at end of file
326 idx1 = 2*idx;
286 idx2 = 2*idx+1; No newline at end of file
327 idx2 = 2*idx+1;
287 Hr[idx1,:] = np.cos(k*(r[i1]-r[i2])*np.sin(thetar)).T.conj(); No newline at end of file
328 Hr[idx1,:] = np.cos(k*(r[i1]-r[i2])*np.sin(thetar)).T.conj();
288 Hr[idx2,:] = np.sin(k*(r[i1]-r[i2])*np.sin(thetar)).T.conj(); No newline at end of file
329 Hr[idx2,:] = np.sin(k*(r[i1]-r[i2])*np.sin(thetar)).T.conj();
289 Ht[idx1,:] = np.cos(k*(r[i1]-r[i2])*np.sin(thetat)).T.conj()*Nr/Nt; No newline at end of file
330 Ht[idx1,:] = np.cos(k*(r[i1]-r[i2])*np.sin(thetat)).T.conj()*Nr/Nt;
290 Ht[idx2,:] = np.sin(k*(r[i1]-r[i2])*np.sin(thetat)).T.conj()*Nr/Nt; No newline at end of file
331 Ht[idx2,:] = np.sin(k*(r[i1]-r[i2])*np.sin(thetat)).T.conj()*Nr/Nt;
291 g[idx1] = (R[i1,i2]).real*Nr/Nt; No newline at end of file
332 g[idx1] = (R[i1,i2]).real*Nr/Nt;
292 g[idx2] = (R[i1,i2]).imag*Nr/Nt; No newline at end of file
333 g[idx2] = (R[i1,i2]).imag*Nr/Nt;
293 gnz[idx1] = (Rnz[i1,i2]).real*Nr/Nt; No newline at end of file
334 gnz[idx1] = (Rnz[i1,i2]).real*Nr/Nt;
294 gnz[idx2] = (Rnz[i1,i2]).imag*Nr/Nt; No newline at end of file
335 gnz[idx2] = (Rnz[i1,i2]).imag*Nr/Nt;
295 No newline at end of file
336
296 # Inversion No newline at end of file
337 # Inversion
297 F = Nr/Nt; # normalization No newline at end of file
338 F = Nr/Nt; # normalization
298 sigma = 1; # set to 1 because the difference is accounted for in G No newline at end of file
339 sigma = 1; # set to 1 because the difference is accounted for in G
299 No newline at end of file
340
300 G = np.linalg.norm(g-gnz)**2 ; # pretend we know in advance the actual value of chi^2 No newline at end of file
341 G = np.linalg.norm(g-gnz)**2 ; # pretend we know in advance the actual value of chi^2
301 lambda0 = 1e-5*np.ones(shape=(M,1)); # initial condition (can be set to anything) No newline at end of file
342 lambda0 = 1e-5*np.ones(shape=(M,1)); # initial condition (can be set to anything)
302 No newline at end of file
343
303 No newline at end of file
344
304 # Whitened solution No newline at end of file
345 # Whitened solution
305 def myfun(lambda1): No newline at end of file
346 def myfun(lambda1):
306 return y_hysell96(lambda1,gnz,sigma,F,G,Hr); No newline at end of file
347 return y_hysell96(lambda1,gnz,sigma,F,G,Hr);
307 No newline at end of file
348
308 No newline at end of file
349
309 tic_maxEnt = time.time(); # start time maxEnt No newline at end of file
350 tic_maxEnt = time.time(); # start time maxEnt
310
351
No newline at end of file
352 lambda1 = root(myfun,lambda0, method='krylov', tol=1e-14); # krylov
311 lambda1 = root(myfun,lambda0, method='krylov', tol=1e-14); No newline at end of file
No newline at end of file
353 # lambda1 = newton_krylov(myfun,lambda0, method='lgmres', x_tol=1e-14) No newline at end of file
312 No newline at end of file
354
313 toc_maxEnt = time.time() No newline at end of file
355 toc_maxEnt = time.time()
314 elapsed_time_maxent = toc_maxEnt - tic_maxEnt; No newline at end of file
356 elapsed_time_maxent = toc_maxEnt - tic_maxEnt;
315 No newline at end of file
357
316 # Solution No newline at end of file
358 # Solution
317 lambda1 = lambda1.x; No newline at end of file
359 lambda1 = lambda1.x;
318 No newline at end of file
360
319 f_maxent = modelf(lambda1, Hr, F); No newline at end of file
361 f_maxent = modelf(lambda1, Hr, F);
320 No newline at end of file
362
321 #------ what is this needed for?-------
363 #------ what is this needed for?-------
No newline at end of file
364 # ystar = myfun(lambda1);
322 ystar = myfun(lambda1);
No newline at end of file
No newline at end of file
365 # Lambda = np.sqrt(sum(lambda1**2*sigma**2)/(4*G));
323 Lambda = np.sqrt(sum(lambda1**2*sigma**2)/(4*G));
No newline at end of file
No newline at end of file
366 # ep = np.multiply(-lambda1,sigma**2)/ (2*Lambda);
324 ep = np.multiply(-lambda1,sigma**2)/ (2*Lambda);
No newline at end of file
No newline at end of file
367 # es = np.dot(Hr, f_maxent) - gnz;
325 es = np.dot(Hr, f_maxent) - gnz;
No newline at end of file
No newline at end of file
368 # chi2 = np.sum((es/sigma)**2); No newline at end of file
326 chi2 = np.sum((es/sigma)**2); No newline at end of file
327 #-------------------------------------- No newline at end of file
369 #--------------------------------------
328 No newline at end of file
370
329 No newline at end of file
371
330 #------------------------------------------------------------------------------------------------- No newline at end of file
372 #-------------------------------------------------------------------------------------------------
331 # CS inversion using Iteratively Reweighted Least Squares (IRLS) No newline at end of file
373 # CS inversion using Iteratively Reweighted Least Squares (IRLS)
332 #------------------------------------------------------------------------------------------------- No newline at end of file
374 #-------------------------------------------------------------------------------------------------
333 # (Use Nr, thetar, gnz, and Hr from MaxEnt above) No newline at end of file
375 # (Use Nr, thetar, gnz, and Hr from MaxEnt above)
334 No newline at end of file
376
335 Psi = deb4_basis(Nr); No newline at end of file
377 Psi = deb4_basis(Nr);
336 No newline at end of file
378
337 # add "sum to 1" constraint No newline at end of file
379 # add "sum to 1" constraint
338 H2 = np.concatenate( (Hr, np.ones(shape=(1,Nr))), axis=0 ); No newline at end of file
380 H2 = np.concatenate( (Hr, np.ones(shape=(1,Nr))), axis=0 );
339 g2 = np.concatenate( (gnz, np.array([[Nr/Nt]])), axis=0 ); No newline at end of file
381 g2 = np.concatenate( (gnz, np.array([[Nr/Nt]])), axis=0 );
340 No newline at end of file
382
341 tic_cs = time.time(); No newline at end of file
383 tic_cs = time.time();
342 No newline at end of file
384
343 # Inversion No newline at end of file
385 # Inversion
344 s = irls_dn2(np.dot(H2,Psi),g2,0.5,G); No newline at end of file
386 s = irls_dn2(np.dot(H2,Psi),g2,0.5,G);
345 No newline at end of file
387
346 toc_cs = time.time() No newline at end of file
388 toc_cs = time.time()
347 elapsed_time_cs = toc_cs - tic_cs; No newline at end of file
389 elapsed_time_cs = toc_cs - tic_cs;
348 No newline at end of file
390
349 # Brightness function No newline at end of file
391 # Brightness function
350 f_cs = np.dot(Psi,s); No newline at end of file
392 f_cs = np.dot(Psi,s);
351 No newline at end of file
393
352 No newline at end of file
394
353 print "shapes:" No newline at end of file
395 print "shapes:"
354 print "Psi ", Psi.shape No newline at end of file
396 print "Psi ", Psi.shape
355 print "s ", s.shape No newline at end of file
397 print "s ", s.shape
356 print "f_CS ", f_cs.shape No newline at end of file
398 print "f_CS ", f_cs.shape
357
399
No newline at end of file
400
358 # Plot
No newline at end of file
No newline at end of file
401 # # Plot
359 plt.figure(4)
No newline at end of file
No newline at end of file
402 # plt.figure(4)
360 plt.plot(thetar,f_cs,'r.-');
No newline at end of file
No newline at end of file
403 # plt.plot(thetar,f_cs,'r.-');
361 plt.plot(thetat,fact,'k-');
No newline at end of file
No newline at end of file
404 # plt.plot(thetat,fact,'k-');
362 plt.title('Reconstruction with Compressed Sensing') No newline at end of file
No newline at end of file
405 # plt.title('Reconstruction with Compressed Sensing') No newline at end of file
363 No newline at end of file
406
364 No newline at end of file
407
365 #------------------------------------------------------------------------------------------------- No newline at end of file
408 #-------------------------------------------------------------------------------------------------
366 # Scaling and shifting No newline at end of file
409 # Scaling and shifting
367 # (Only necessary for Capon solution) No newline at end of file
410 # (Only necessary for Capon solution)
368 #------------------------------------------------------------------------------------------------- No newline at end of file
411 #-------------------------------------------------------------------------------------------------
369 f_capon = f_capon/np.max(f_capon)*np.max(fact); No newline at end of file
412 f_capon = f_capon/np.max(f_capon)*np.max(fact);
370 No newline at end of file
413
371 No newline at end of file
414
372 #------------------------------------------------------------------------------------------------- No newline at end of file
415 #-------------------------------------------------------------------------------------------------
373 # Analyze stuff No newline at end of file
416 # Analyze stuff
374 #------------------------------------------------------------------------------------------------- No newline at end of file
417 #-------------------------------------------------------------------------------------------------
375 No newline at end of file
418
376 # Calculate MSE No newline at end of file
419 # Calculate MSE
377 rmse_fourier = np.sqrt(np.mean((f_fourier - factr)**2)); No newline at end of file
420 rmse_fourier = np.sqrt(np.mean((f_fourier - factr)**2));
378 rmse_capon = np.sqrt(np.mean((f_capon - factr)**2)); No newline at end of file
421 rmse_capon = np.sqrt(np.mean((f_capon - factr)**2));
379 rmse_maxent = np.sqrt(np.mean((f_maxent - factr)**2)); No newline at end of file
422 rmse_maxent = np.sqrt(np.mean((f_maxent - factr)**2));
380 rmse_cs = np.sqrt(np.mean((f_cs - factr)**2)); No newline at end of file
423 rmse_cs = np.sqrt(np.mean((f_cs - factr)**2));
381 No newline at end of file
424
382 No newline at end of file
425
383 relrmse_fourier = rmse_fourier / np.linalg.norm(fact); No newline at end of file
426 relrmse_fourier = rmse_fourier / np.linalg.norm(fact);
384 relrmse_capon = rmse_capon / np.linalg.norm(fact); No newline at end of file
427 relrmse_capon = rmse_capon / np.linalg.norm(fact);
385 relrmse_maxent = rmse_maxent / np.linalg.norm(fact); No newline at end of file
428 relrmse_maxent = rmse_maxent / np.linalg.norm(fact);
386 relrmse_cs = rmse_cs / np.linalg.norm(fact); No newline at end of file
429 relrmse_cs = rmse_cs / np.linalg.norm(fact);
387 No newline at end of file
430
388 No newline at end of file
431
389 # Calculate correlation No newline at end of file
432 # Calculate correlation
390 corr_fourier = np.dot(f_fourier.T.conj(),factr) / (np.linalg.norm(f_fourier)*np.linalg.norm(factr)); No newline at end of file
433 corr_fourier = np.dot(f_fourier.T.conj(),factr) / (np.linalg.norm(f_fourier)*np.linalg.norm(factr));
391 corr_capon = np.dot(f_capon.T.conj(),factr) / (np.linalg.norm(f_capon)*np.linalg.norm(factr)); No newline at end of file
434 corr_capon = np.dot(f_capon.T.conj(),factr) / (np.linalg.norm(f_capon)*np.linalg.norm(factr));
392 corr_maxent = np.dot(f_maxent.T.conj(),factr) / (np.linalg.norm(f_maxent)*np.linalg.norm(factr)); No newline at end of file
435 corr_maxent = np.dot(f_maxent.T.conj(),factr) / (np.linalg.norm(f_maxent)*np.linalg.norm(factr));
393 corr_cs = np.dot(f_cs.T.conj(),factr) / (np.linalg.norm(f_cs)*np.linalg.norm(factr)); No newline at end of file
436 corr_cs = np.dot(f_cs.T.conj(),factr) / (np.linalg.norm(f_cs)*np.linalg.norm(factr));
394 No newline at end of file
437
395 No newline at end of file
438
396 # Calculate centered correlation No newline at end of file
439 # Calculate centered correlation
397 f0 = factr - np.mean(factr); No newline at end of file
440 f0 = factr - np.mean(factr);
398 f1 = f_fourier - np.mean(f_fourier); No newline at end of file
441 f1 = f_fourier - np.mean(f_fourier);
399 No newline at end of file
442 corrc_fourier = np.dot(f0.T.conj(),f1) / (np.linalg.norm(f0)*np.linalg.norm(f1));
400 corrc_fourier = np.dot(f0.T.conj(),f1) / (np.linalg.norm(f0)*np.linalg.norm(f1)); No newline at end of file
443 f1 = f_capon - np.mean(f_capon);
401 f1 = f_capon - np.mean(f_capon); No newline at end of file
444 corrc_capon = np.dot(f0.T.conj(),f1) / (np.linalg.norm(f0)*np.linalg.norm(f1));
402 corrc_capon = np.dot(f0.T.conj(),f1) / (np.linalg.norm(f0)*np.linalg.norm(f1)); No newline at end of file
445 f1 = f_maxent - np.mean(f_maxent);
403 f1 = f_maxent - np.mean(f_maxent); No newline at end of file
446 corrc_maxent = np.dot(f0.T.conj(),f1) / (np.linalg.norm(f0)*np.linalg.norm(f1));
404 corrc_maxent = np.dot(f0.T.conj(),f1) / (np.linalg.norm(f0)*np.linalg.norm(f1)); No newline at end of file
447 f1 = f_cs - np.mean(f_cs);
405 f1 = f_cs - np.mean(f_cs); No newline at end of file
448 corrc_cs = np.dot(f0.T.conj(),f1) / (np.linalg.norm(f0)*np.linalg.norm(f1));
406 corrc_cs = np.dot(f0.T.conj(),f1) / (np.linalg.norm(f0)*np.linalg.norm(f1)); No newline at end of file
449
407 No newline at end of file
450
408 No newline at end of file
451 #-------------------------------------------------------------------------------------------------
409 #------------------------------------------------------------------------------------------------- No newline at end of file
452 # Plot stuff
410 # Plot stuff No newline at end of file
453 #-------------------------------------------------------------------------------------------------
411 #------------------------------------------------------------------------------------------------- No newline at end of file
454
412 No newline at end of file
455 #---- Capon----#
413 #---- Capon----# No newline at end of file
456 plt.figure(5)
414 plt.figure(5) No newline at end of file
457 plt.subplot(3, 1, 1)
415 plt.subplot(3, 1, 1) No newline at end of file
458 plt.plot(180/np.pi*thetar, f_capon, 'r', label='Capon')
416 plt.plot(180/np.pi*thetar, f_capon, 'r', label='Capon') No newline at end of file
459 plt.plot(180/np.pi*thetat,fact, 'k--', label='Truth')
417 plt.plot(180/np.pi*thetat,fact, 'k--', label='Truth') No newline at end of file
460 plt.ylabel('Power (arbitrary units)')
418 plt.ylabel('Power (arbitrary units)') No newline at end of file
461 plt.legend(loc='upper right')
419 plt.legend(loc='upper right') No newline at end of file
462
420 No newline at end of file
463 # formatting y-axis
421 # formatting y-axis No newline at end of file
464 locs,labels = plt.yticks()
422 locs,labels = plt.yticks() No newline at end of file
465 plt.yticks(locs, map(lambda x: "%.1f" % x, locs*1e4))
423 plt.yticks(locs, map(lambda x: "%.1f" % x, locs*1e4)) No newline at end of file
466 plt.text(0.0, 1.01, '1e-4', fontsize=10, transform = plt.gca().transAxes)
424 plt.text(0.0, 1.01, '1e-4', fontsize=10, transform = plt.gca().transAxes) No newline at end of file
467
425 No newline at end of file
468
426 No newline at end of file
469 #---- MaxEnt----#
427 #---- MaxEnt----# No newline at end of file
470 plt.subplot(3, 1, 2)
428 plt.subplot(3, 1, 2) No newline at end of file
471 plt.plot(180/np.pi*thetar, f_maxent, 'r', label='MaxEnt')
429 plt.plot(180/np.pi*thetar, f_maxent, 'r', label='MaxEnt') No newline at end of file
472 plt.plot(180/np.pi*thetat,fact, 'k--', label='Truth')
430 plt.plot(180/np.pi*thetat,fact, 'k--', label='Truth') No newline at end of file
473 plt.ylabel('Power (arbitrary units)')
431 plt.ylabel('Power (arbitrary units)') No newline at end of file
474 plt.legend(loc='upper right')
432 plt.legend(loc='upper right') No newline at end of file
475
433 No newline at end of file
476 # formatting y-axis
434 # formatting y-axis No newline at end of file
477 locs,labels = plt.yticks()
435 locs,labels = plt.yticks() No newline at end of file
478 plt.yticks(locs, map(lambda x: "%.1f" % x, locs*1e4))
436 plt.yticks(locs, map(lambda x: "%.1f" % x, locs*1e4)) No newline at end of file
479 plt.text(0.0, 1.01, '1e-4', fontsize=10, transform = plt.gca().transAxes)
437 plt.text(0.0, 1.01, '1e-4', fontsize=10, transform = plt.gca().transAxes) No newline at end of file
480
438 No newline at end of file
481
439 No newline at end of file
482 #---- Compressed Sensing----#
440 #---- Compressed Sensing----# No newline at end of file
483 plt.subplot(3, 1, 3)
441 plt.subplot(3, 1, 3) No newline at end of file
484 plt.plot(180/np.pi*thetar, f_cs, 'r', label='CS')
442 plt.plot(180/np.pi*thetar, f_cs, 'r', label='CS') No newline at end of file
485 plt.plot(180/np.pi*thetat,fact, 'k--', label='Truth')
443 plt.plot(180/np.pi*thetat,fact, 'k--', label='Truth') No newline at end of file
486 plt.ylabel('Power (arbitrary units)')
444 plt.ylabel('Power (arbitrary units)') No newline at end of file
487 plt.legend(loc='upper right')
445 plt.legend(loc='upper right') No newline at end of file
488
446 No newline at end of file
489 # formatting y-axis
447 # formatting y-axis No newline at end of file
490 locs,labels = plt.yticks()
448 locs,labels = plt.yticks() No newline at end of file
491 plt.yticks(locs, map(lambda x: "%.1f" % x, locs*1e4))
449 plt.yticks(locs, map(lambda x: "%.1f" % x, locs*1e4)) No newline at end of file
492 plt.text(0.0, 1.01, '1e-4', fontsize=10, transform = plt.gca().transAxes)
450 plt.text(0.0, 1.01, '1e-4', fontsize=10, transform = plt.gca().transAxes) No newline at end of file
493
451 No newline at end of file
494 # # Store Results
452 # # Store Results No newline at end of file
495 corr[0, snri, Ni] = corr_fourier;
453 corr[0, snri, Ni] = corr_fourier; No newline at end of file
496 corr[1, snri, Ni] = corr_capon;
454 corr[1, snri, Ni] = corr_capon; No newline at end of file
497 corr[2, snri, Ni] = corr_maxent;
455 corr[2, snri, Ni] = corr_maxent; No newline at end of file
498 corr[3, snri, Ni] = corr_cs;
456 corr[3, snri, Ni] = corr_cs; No newline at end of file
499
457 No newline at end of file
500 rmse[0,snri,Ni] = relrmse_fourier;
458 rmse[0,snri,Ni] = relrmse_fourier; No newline at end of file
501 rmse[1,snri,Ni] = relrmse_capon;
459 rmse[1,snri,Ni] = relrmse_capon; No newline at end of file
502 rmse[2,snri,Ni] = relrmse_maxent;
460 rmse[2,snri,Ni] = relrmse_maxent; No newline at end of file
503 rmse[3,snri,Ni] = relrmse_cs;
461 rmse[3,snri,Ni] = relrmse_cs; No newline at end of file
504
462 No newline at end of file
505 corrc[0,snri,Ni] = corrc_fourier;
463 corrc[0,snri,Ni] = corrc_fourier; No newline at end of file
506 corrc[1,snri,Ni] = corrc_capon;
464 corrc[1,snri,Ni] = corrc_capon; No newline at end of file
507 corrc[2,snri,Ni] = corrc_maxent;
465 corrc[2,snri,Ni] = corrc_maxent; No newline at end of file
508 corrc[3,snri,Ni] = corrc_cs;
466 corrc[3,snri,Ni] = corrc_cs; No newline at end of file
509
467 No newline at end of file
510 print "--------Time performace--------"
468 print "--------Time performace--------" No newline at end of file
511 print 'Capon:\t', elapsed_time_capon, 'sec';
469 print 'Capon:\t', elapsed_time_capon, 'sec'; No newline at end of file
512 print 'Maxent:\t',elapsed_time_maxent, 'sec';
470 print 'Maxent:\t',elapsed_time_maxent, 'sec'; No newline at end of file
513 print 'CS:\t',elapsed_time_cs, 'sec\n';
471 print 'CS:\t',elapsed_time_cs, 'sec\n'; No newline at end of file
514
472 No newline at end of file
515 print (NN*(snri)+Ni+1), '/', (SNRdBvec.size*NN), '\n';
473 print (NN*(snri)+Ni+1), '/', (SNRdBvec.size*NN), '\n'; No newline at end of file
516
474 No newline at end of file
517
475 No newline at end of file
518 #-------------------------------------------------------------------------------------------------
476 #------------------------------------------------------------------------------------------------- No newline at end of file
519 # Analyze and plot statistics
477 # Analyze and plot statistics No newline at end of file
520 #-------------------------------------------------------------------------------------------------
478 #------------------------------------------------------------------------------------------------- No newline at end of file
521
479 No newline at end of file
522 metric = corr; # set this to rmse, corr, or corrc
480 metric = corr; # set this to rmse, corr, or corrc No newline at end of file
523
481 No newline at end of file
524 # Remove outliers (this part was experimental and wasn't used in the paper)
482 # Remove outliers (this part was experimental and wasn't used in the paper) No newline at end of file
525 # nsig = 3;
483 # nsig = 3; No newline at end of file
526 # for i = 1:4
484 # for i = 1:4 No newline at end of file
527 # for snri = 1:length(SNRdBvec)
485 # for snri = 1:length(SNRdBvec) No newline at end of file
528 # av = mean(met(i,snri,:));
486 # av = mean(met(i,snri,:)); No newline at end of file
529 # s = std(met(i,snri,:));
487 # s = std(met(i,snri,:)); No newline at end of file
530 # idx = abs(met(i,snri,:) - av) > nsig*s;
488 # idx = abs(met(i,snri,:) - av) > nsig*s; No newline at end of file
531 # met(i,snri,idx) = nan;
489 # met(i,snri,idx) = nan; No newline at end of file
532 # if sum(idx)>0
490 # if sum(idx)>0 No newline at end of file
533 # fprintf('i=%i, snr=%i, %i/%i pts removed\n',...
491 # fprintf('i=%i, snr=%i, %i/%i pts removed\n',... No newline at end of file
534 # i,round(SNRdBvec(snri)),sum(idx),length(idx));
492 # i,round(SNRdBvec(snri)),sum(idx),length(idx)); No newline at end of file
535 # end
493 # end No newline at end of file
536 # end
494 # end No newline at end of file
537 # end
495 # end No newline at end of file
538
496 No newline at end of file
539
497 No newline at end of file
540 # Avg ignoring NaNs
498 # Avg ignoring NaNs No newline at end of file
541 ave = np.zeros(shape=(4))
499 ave = np.zeros(shape=(4)) No newline at end of file
542
500 No newline at end of file
543 print metric
501 print metric
544 print metric.shape
No newline at end of file
545
502
No newline at end of file
No newline at end of file
546 ave[0] = nanmean(metric[0,0,:]); # Fourier
503 ave[0] = nanmean(metric[0,:,:]); # Fourier
No newline at end of file
No newline at end of file
547 ave[1] = nanmean(metric[1,0,:]); # Capon
504 ave[1] = nanmean(metric[1,:,:]); # Capon
No newline at end of file
No newline at end of file
548 ave[2] = nanmean(metric[2,0,:]); # MaxEnt
505 ave[2] = nanmean(metric[2,:,:]); # MaxEnt
No newline at end of file
No newline at end of file
549 ave[3] = nanmean(metric[3,0,:]); # Compressed Sensing
506 ave[3] = nanmean(metric[3,:,:]); # Compressed Sensing No newline at end of file
No newline at end of file
550
No newline at end of file
551 print ave
No newline at end of file
552 No newline at end of file
507 No newline at end of file
553 # Plot based on chosen metric
508 # Plot based on chosen metric No newline at end of file
554 plt.figure(6);
509 plt.figure(6); No newline at end of file
555 f = plt.scatter(SNRdBvec, ave[0], marker='+', color='b', s=60); # Fourier
510 f = plt.scatter(SNRdBvec, ave[0], marker='+', color='b', s=60); # Fourier No newline at end of file
556 c = plt.scatter(SNRdBvec, ave[1], marker='o', color= 'g', s=60); # Capon
511 c = plt.scatter(SNRdBvec, ave[1], marker='o', color= 'g', s=60); # Capon No newline at end of file
557 me= plt.scatter(SNRdBvec, ave[2], marker='s', color= 'c', s=60); # MaxEnt
512 me= plt.scatter(SNRdBvec, ave[2], marker='s', color= 'c', s=60); # MaxEnt No newline at end of file
558 cs= plt.scatter(SNRdBvec, ave[3], marker='*', color='r', s=60); # Compressed Sensing
513 cs= plt.scatter(SNRdBvec, ave[3], marker='*', color='r', s=60); # Compressed Sensing No newline at end of file
559
514 No newline at end of file
560 plt.legend((f,c,me,cs),('Fourier','Capon', 'MaxEnt', 'Comp. Sens.'),scatterpoints=1, loc='upper right')
515 plt.legend((f,c,me,cs),('Fourier','Capon', 'MaxEnt', 'Comp. Sens.'),scatterpoints=1, loc='upper right') No newline at end of file
561 plt.xlabel('SNR')
516 plt.xlabel('SNR') No newline at end of file
562 plt.ylabel('Correlation with Truth')
517 plt.ylabel('Correlation with Truth') No newline at end of file
563
518 No newline at end of file
564 print "--------Correlations--------"
519 print "--------Correlations--------" No newline at end of file
565 print "Fourier:", ave[0]
520 print "Fourier:", ave[0] No newline at end of file
566 print "Capon:\t", ave[1]
521 print "Capon:\t", ave[1] No newline at end of file
567 print "MaxEnt:\t", ave[2]
522 print "MaxEnt:\t", ave[2] No newline at end of file
568 print "CS:\t", ave[3]
523 print "CS:\t", ave[3] No newline at end of file
569
524 No newline at end of file
570 plt.show()
525 plt.show() No newline at end of file
571
526 No newline at end of file
572
527 No newline at end of file
573
Show file before | Show file after | Hide comments
@@ -1,65 +1,65
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 No newline at end of file
8
9 def dualfilt1(): No newline at end of file
9 def dualfilt1():
10 No newline at end of file
10
11 # Kingsbury Q-filters for the dual-tree complex DWT No newline at end of file
11 # Kingsbury Q-filters for the dual-tree complex DWT
12 # No newline at end of file
12 #
13 # USAGE: No newline at end of file
13 # USAGE:
14 # [af, sf] = dualfilt1 No newline at end of file
14 # [af, sf] = dualfilt1
15 # OUTPUT: No newline at end of file
15 # OUTPUT:
16 # af{i}, i = 1,2 - analysis filters for tree i No newline at end of file
16 # af{i}, i = 1,2 - analysis filters for tree i
17 # sf{i}, i = 1,2 - synthesis filters for tree i No newline at end of file
17 # sf{i}, i = 1,2 - synthesis filters for tree i
18 # note: af{2} is the reverse of af{1} No newline at end of file
18 # note: af{2} is the reverse of af{1}
19 # REFERENCE: No newline at end of file
19 # REFERENCE:
20 # N. G. Kingsbury, "A dual-tree complex wavelet No newline at end of file
20 # N. G. Kingsbury, "A dual-tree complex wavelet
21 # transform with improved orthogonality and symmetry No newline at end of file
21 # transform with improved orthogonality and symmetry
22 # properties", Proceedings of the IEEE Int. Conf. on No newline at end of file
22 # properties", Proceedings of the IEEE Int. Conf. on
23 # Image Proc. (ICIP), 2000 No newline at end of file
23 # Image Proc. (ICIP), 2000
24 # See dualtree No newline at end of file
24 # See dualtree
25 # No newline at end of file
25 #
26 # WAVELET SOFTWARE AT POLYTECHNIC UNIVERSITY, BROOKLYN, NY No newline at end of file
26 # WAVELET SOFTWARE AT POLYTECHNIC UNIVERSITY, BROOKLYN, NY
27 # http://taco.poly.edu/WaveletSoftware/ No newline at end of file
27 # http://taco.poly.edu/WaveletSoftware/
28
28
No newline at end of file
29 # These coefficients are rounded to 8 decimal places. No newline at end of file
29 # These cofficients are rounded to 8 decimal places. No newline at end of file
30 No newline at end of file
30
31 a1 = np.array([ No newline at end of file
31 a1 = np.array([
32 [ 0.03516384000000, 0], No newline at end of file
32 [ 0.03516384000000, 0],
33 [ 0, 0], No newline at end of file
33 [ 0, 0],
34 [-0.08832942000000, -0.11430184000000], No newline at end of file
34 [-0.08832942000000, -0.11430184000000],
35 [ 0.23389032000000, 0], No newline at end of file
35 [ 0.23389032000000, 0],
36 [ 0.76027237000000, 0.58751830000000], No newline at end of file
36 [ 0.76027237000000, 0.58751830000000],
37 [ 0.58751830000000, -0.76027237000000], No newline at end of file
37 [ 0.58751830000000, -0.76027237000000],
38 [ 0, 0.23389032000000], No newline at end of file
38 [ 0, 0.23389032000000],
39 [-0.11430184000000, 0.08832942000000], No newline at end of file
39 [-0.11430184000000, 0.08832942000000],
40 [ 0, 0], No newline at end of file
40 [ 0, 0],
41 [ 0, -0.03516384000000] No newline at end of file
41 [ 0, -0.03516384000000]
42 ]); No newline at end of file
42 ]);
43 No newline at end of file
43
44 a2 = np.array([ No newline at end of file
44 a2 = np.array([
45 [ 0, -0.03516384000000], No newline at end of file
45 [ 0, -0.03516384000000],
46 [ 0, 0], No newline at end of file
46 [ 0, 0],
47 [-0.11430184000000, 0.08832942000000], No newline at end of file
47 [-0.11430184000000, 0.08832942000000],
48 [ 0, 0.23389032000000], No newline at end of file
48 [ 0, 0.23389032000000],
49 [ 0.58751830000000, -0.76027237000000], No newline at end of file
49 [ 0.58751830000000, -0.76027237000000],
50 [ 0.76027237000000, 0.58751830000000], No newline at end of file
50 [ 0.76027237000000, 0.58751830000000],
51 [ 0.23389032000000, 0], No newline at end of file
51 [ 0.23389032000000, 0],
52 [ -0.08832942000000, -0.11430184000000], No newline at end of file
52 [ -0.08832942000000, -0.11430184000000],
53 [ 0, 0], No newline at end of file
53 [ 0, 0],
54 [ 0.03516384000000, 0] No newline at end of file
54 [ 0.03516384000000, 0]
55 ]); No newline at end of file
55 ]);
56 No newline at end of file
56
57 af = np.array([ [a1,a2] ], dtype=object) No newline at end of file
57 af = np.array([ [a1,a2] ], dtype=object)
58 No newline at end of file
58
59 s1 = a1[::-1] No newline at end of file
59 s1 = a1[::-1]
60 s2 = a2[::-1] No newline at end of file
60 s2 = a2[::-1]
61 No newline at end of file
61
62 sf = np.array([ [s1,s2] ], dtype=object) No newline at end of file
62 sf = np.array([ [s1,s2] ], dtype=object)
63 No newline at end of file
63
64 No newline at end of file
64
65 return af, sf No newline at end of file
65 return af, sf
Show file before | Show file after | Hide comments
@@ -1,52 +1,52
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 brentq No newline at end of file
8 from scipy.optimize import brentq
9 No newline at end of file
9
10 def irls_dn2(A,b,p,G): No newline at end of file
10 def irls_dn2(A,b,p,G):
11 No newline at end of file
11
12 # Minimize ||u||_p subject to ||A*u-b||_2^2 <= G (0 < p <= 1) No newline at end of file
12 # Minimize ||u||_p subject to ||A*u-b||_2^2 <= G (0 < p <= 1)
13 No newline at end of file
13
14 # What this function actually does is finds the lambda1 so that the solution No newline at end of file
14 # What this function actually does is finds the lambda1 so that the solution
15 # to the following problem satisfies ||A*u-b||_2^2 <= G: No newline at end of file
15 # to the following problem satisfies ||A*u-b||_2^2 <= G:
16 # Minimize lambda1*||u||_p + ||A*u-b||_2 No newline at end of file
16 # Minimize lambda1*||u||_p + ||A*u-b||_2
17 No newline at end of file
17
18 # Start with a large lambda1, and do a line search until fidelity <= G. No newline at end of file
18 # Start with a large lambda1, and do a line search until fidelity <= G.
19 # (Inversions with large lambda1 are really fast anyway). No newline at end of file
19 # (Inversions with large lambda1 are really fast anyway).
20 No newline at end of file
20
21 # Then spin up fsolve to localize the root even better No newline at end of file
21 # Then spin up fsolve to localize the root even better
22 No newline at end of file
22
23 # Line Search No newline at end of file
23 # Line Search
24 No newline at end of file
24
25 alpha = 2.0; # Line search parameter No newline at end of file
25 alpha = 2.0; # Line search parameter
26 lambda1 = 1e5; # What's a reasonable but safe initial guess? No newline at end of file
26 lambda1 = 1e5; # What's a reasonable but safe initial guess?
27 u = irls_dn(A,b,p,lambda1); No newline at end of file
27 u = irls_dn(A,b,p,lambda1);
28 No newline at end of file
28
29 No newline at end of file
29
30 fid = norm(np.dot(A,u)-b)**2; No newline at end of file
30 fid = norm(np.dot(A,u)-b)**2;
31 No newline at end of file
31
32 print '----------------------------------\n'; No newline at end of file
32 print '----------------------------------\n';
33 No newline at end of file
33
34 while (fid >= G): No newline at end of file
34 while (fid >= G):
35 lambda1 = lambda1 / alpha; # Balance between speed and accuracy No newline at end of file
35 lambda1 = lambda1 / alpha; # Balance between speed and accuracy
36 u = irls_dn(A,b,p,lambda1); No newline at end of file
36 u = irls_dn(A,b,p,lambda1);
37 fid = norm(np.dot(A,u)-b)**2; No newline at end of file
37 fid = norm(np.dot(A,u)-b)**2;
38 print 'lambda = %2e \t' % lambda1, '||A*u-b||^2 = %.1f' % fid; No newline at end of file
38 print 'lambda = %2e \t' % lambda1, '||A*u-b||^2 = %.1f' % fid;
39
39
No newline at end of file
40 # Refinement using brentq (equiv to fzero in matlab) No newline at end of file
40 # Refinement using fzero/ brentq No newline at end of file
41 lambda0 = np.array([lambda1,lambda1*alpha]); # interval with zero-crossing No newline at end of file
41 lambda0 = np.array([lambda1,lambda1*alpha]); # interval with zero-crossing
42 No newline at end of file
42
43 def myfun(lambda1): No newline at end of file
43 def myfun(lambda1):
44 return norm(np.dot(A, irls_dn(A,b,p,lambda1)) - b)**2 - G; No newline at end of file
44 return norm(np.dot(A, irls_dn(A,b,p,lambda1)) - b)**2 - G;
45 No newline at end of file
45
46 # Find zero-crossing at given interval (lambda1, lambda1*alpha) No newline at end of file
46 # Find zero-crossing at given interval (lambda1, lambda1*alpha)
47 lambda1 = brentq(myfun, lambda0[0], lambda0[1], xtol=0.01*lambda1) No newline at end of file
47 lambda1 = brentq(myfun, lambda0[0], lambda0[1], xtol=0.01*lambda1)
48 No newline at end of file
48
49 u = irls_dn(A,b,p,lambda1); No newline at end of file
49 u = irls_dn(A,b,p,lambda1);
50 No newline at end of file
50
51 return u; No newline at end of file
51 return u;
52 No newline at end of file
52
General Comments 0
You need to be logged in to leave comments. Login now