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