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