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