##// END OF EJS Templates
jro_soft » compressed_sensing
RSS
Fixed errors, implemented sfb and idualtree classes
yamaro -
r19:20
parent child
Show More
Show file before | Show file after | Hide comments
@@ -0,0 +1,46
1 '''
2 Created on Jun 5, 2014
3
4 @author: Yolian Amaro
5 '''
6
7 from sfb import *
8
9 def idualtree(w, J, Fsf, sf):
10
11 # Inverse Dual-tree Complex DWT
12 #
13 # USAGE:
14 # y = idualtree(w, J, Fsf, sf)
15 # INPUT:
16 # w - DWT coefficients
17 # J - number of stages
18 # Fsf - synthesis filters for the last stage
19 # sf - synthesis filters for preceeding stages
20 # OUTUT:
21 # y - output signal
22 # See dualtree
23 #
24 # WAVELET SOFTWARE AT POLYTECHNIC UNIVERSITY, BROOKLYN, NY
25 # http://taco.poly.edu/WaveletSoftware/
26
27 # Tree 1
28 y1 = w[J][0];
29
30 for j in range (J-1, 0, -1):
31 y1 = sfb(y1, w[j][0], sf[0,0]);
32
33 y1 = sfb(y1, w[0][0], Fsf[0,0]);
34
35 # Tree 2
36 y2 = w[J][1];
37
38 for j in range (J-1, 0, -1):
39 y2 = sfb(y2, w[j][2], sf[0,1]);
40
41 y2 = sfb(y2, w[0][1], Fsf[0,1]);
42
43 # normalization
44 y = (y1 + y2)/np.sqrt(2);
45
46 return y
Show file before | Show file after | Hide comments
@@ -0,0 +1,68
1 '''
2 Created on Jun 5, 2014
3
4 @author: Yolian Amaro
5 '''
6
7 from multirate import *
8 import numpy as np
9 from cshift import *
10
11 def sfb(lo, hi, sf):
12
13 # Synthesis filter bank
14 #
15 # USAGE:
16 # y = sfb(lo, hi, sf)
17 # INPUT:
18 # lo - low frqeuency input
19 # hi - high frequency input
20 # sf - synthesis filters
21 # sf(:, 1) - lowpass filter (even length)
22 # sf(:, 2) - highpass filter (even length)
23 # OUTPUT:
24 # y - output signal
25 # See also afb
26 #
27 # WAVELET SOFTWARE AT POLYTECHNIC UNIVERSITY, BROOKLYN, NY
28 # http://taco.poly.edu/WaveletSoftware/
29
30 N = 2*lo.size;
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
42 # Need to change format for upfirdn funct:
43 lo = lo.T.conj()
44 lo = lo.reshape(lo.size)
45
46 print 'sfb hi', hi
47
48 # Need to change format for upfirdn funct:
49 hi = hi.T.conj()
50 hi = hi.reshape(hi.size)
51
52 #hi = hi.reshape(1, hi.size)
53
54 lo = upfirdn(lo, sf[:,0], 2, 1);
55 hi = upfirdn(hi, sf[:,1], 2, 1);
56 y = lo + hi;
57 y[0:L-1] = y[0:L-1] + y[N+ np.arange(0,L-1)]; #CHECK IF ARANGE IS CORRECT
58 y = y[0:N];
59
60 print 'y en sbf\n', y.shape
61
62 y = y.reshape(1, y.size)
63 print 'y en sbf\n', y.shape
64
65 y = cshift(y, 1-L/2);
66
67 return y;
68
Show file before | Show file after | Hide comments
@@ -1,65 +1,63
1 1 '''
2 2 Created on May 26, 2014
3 3
4 4 @author: Yolian Amaro
5 5 '''
6 6
7 7 import pywt
8 8 import numpy as np
9 9
10 10 def FSfarras():
11 11 #function [af, sf] = FSfarras
12 12
13 13 # Farras filters organized for the dual-tree
14 14 # complex DWT.
15 15 #
16 16 # USAGE:
17 17 # [af, sf] = FSfarras
18 18 # OUTPUT:
19 19 # af{i}, i = 1,2 - analysis filters for tree i
20 20 # sf{i}, i = 1,2 - synthesis filters for tree i
21 21 # See farras, dualtree, dualfilt1.
22 22 #
23 23 # WAVELET SOFTWARE AT POLYTECHNIC UNIVERSITY, BROOKLYN, NY
24 24 # http://taco.poly.edu/WaveletSoftware/
25 25 #
26 26 # Translated to Python by Yolian Amaro
27 27
28 28
29
30 29 a1 = np.array( [
31 30 [ 0, 0],
32 31 [-0.08838834764832, -0.01122679215254],
33 32 [ 0.08838834764832, 0.01122679215254],
34 33 [ 0.69587998903400, 0.08838834764832],
35 34 [ 0.69587998903400, 0.08838834764832],
36 35 [ 0.08838834764832, -0.69587998903400],
37 36 [-0.08838834764832, 0.69587998903400],
38 37 [ 0.01122679215254, -0.08838834764832],
39 38 [ 0.01122679215254, -0.08838834764832],
40 39 [0, 0]
41 40 ] );
42 41
43 42 a2 = np.array([
44 43 [ 0.01122679215254, 0],
45 44 [ 0.01122679215254, 0],
46 45 [-0.08838834764832, -0.08838834764832],
47 46 [ 0.08838834764832, -0.08838834764832],
48 47 [ 0.69587998903400, 0.69587998903400],
49 48 [ 0.69587998903400, -0.69587998903400],
50 49 [ 0.08838834764832, 0.08838834764832],
51 50 [-0.08838834764832, 0.08838834764832],
52 51 [ 0, 0.01122679215254],
53 52 [ 0, -0.01122679215254]
54 53 ]);
55 54
56 #print a2.shape
57 55
58 56 af = np.array([ [a1,a2] ], dtype=object)
59 57
60 58 s1 = a1[::-1]
61 59 s2 = a2[::-1]
62 60
63 61 sf = np.array([ [s1,s2] ], dtype=object)
64 62
65 63 return af, sf
Show file before | Show file after | Hide comments
@@ -1,468 +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 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 LINE BELOW (?)
317 Psi = deb4_basis(Nr); ###### REPLACED BY LINEs BELOW (?)
318
319 print 'FINALLY!'
320 print Psi.shape
318 321
319 322 # REMOVE THIS?--------------------------------
320 323 #wavelet1 = pywt.Wavelet('db4')
321 324 #Phi, Psi, x = wavelet1.wavefun(level=3)
322 325 # --------------------------------------------
323 326
324 327 # add "sum to 1" constraint
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 H2 = H2.T.conj();
329
330 print 'H2 shape', H2.shape
331 print 'Psi shape', Psi.shape
332
333 s = irls_dn2(H2*Psi,g2,0.5,G);
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);
334 337 # f_cs = Psi*s;
335 338 #
336 339 # # plot
337 340 # plot(thetar,f_cs,'r.-');
338 341 # hold on;
339 342 # plot(thetat,fact,'k-');
340 343 # hold off;
341 344
342 345
343 346 # # # Scaling and shifting
344 347 # # # Only necessary for capon solution
345 348
346 349
347 350 f_capon = f_capon/np.max(f_capon)*np.max(fact);
348 351
349 352
350 353 ### analyze stuff ######################
351 354 # calculate MSE
352 355 rmse_fourier = np.sqrt(np.mean((f_fourier - factr)**2));
353 356 rmse_capon = np.sqrt(np.mean((f_capon - factr)**2));
354 357 rmse_maxent = np.sqrt(np.mean((f_maxent - factr)**2));
355 358 #rmse_cs = np.sqrt(np.mean((f_cs - factr).^2));
356 359
357 360
358 361 relrmse_fourier = rmse_fourier / np.linalg.norm(fact);
359 362 relrmse_capon = rmse_capon / np.linalg.norm(fact);
360 363 relrmse_maxent = rmse_maxent / np.linalg.norm(fact);
361 364 #relrmse_cs = rmse_cs / np.norm(fact);
362 365
363 366 # To be able to perform dot product (align matrices) done below within the dot calculations
364 367
365 368
366 369 #f_fourier = f_fourier.T.conj()
367 370 #f_capon = f_capon.T.conj()
368 371 #f_maxent = f_maxent.T.conj()
369 372
370 373 #factr = factr.T.conj()
371 374
372 375 # calculate correlation
373 376
374 377 corr_fourier = np.dot(f_fourier.T.conj(),factr) / (np.linalg.norm(f_fourier)*np.linalg.norm(factr));
375 378 corr_capon = np.dot(f_capon.T.conj(),factr) / (np.linalg.norm(f_capon)*np.linalg.norm(factr));
376 379 corr_maxent = np.dot(f_maxent.T.conj(),factr) / (np.linalg.norm(f_maxent)*np.linalg.norm(factr));
377 380 #corr_cs = np.dot(f_cs,factr) / (norm(f_cs)*norm(factr));
378 381
379 382
380 383 # calculate centered correlation
381 384 f0 = factr - np.mean(factr);
382 385 f1 = f_fourier - np.mean(f_fourier);
383 386
384 387 corrc_fourier = np.dot(f0.T.conj(),f1) / (np.linalg.norm(f0)*np.linalg.norm(f1));
385 388 f1 = f_capon - np.mean(f_capon);
386 389 corrc_capon = np.dot(f0.T.conj(),f1) / (np.linalg.norm(f0)*np.linalg.norm(f1));
387 390 f1 = f_maxent - np.mean(f_maxent);
388 391 corrc_maxent = np.dot(f0.T.conj(),f1) / (np.linalg.norm(f0)*np.linalg.norm(f1));
389 392 #f1 = f_cs - mean(f_cs);
390 393 #corrc_cs = dot(f0,f1) / (norm(f0)*norm(f1));
391 394
392 395
393 396
394 397 # # # plot stuff #########################
395 398
396 399 #---- Capon----
397 400 plt.figure(4)
398 401 plt.subplot(2, 1, 1)
399 402 plt.plot(180/math.pi*thetar, f_capon, 'r', label='Capon')
400 403 plt.plot(180/math.pi*thetat,fact, 'k--', label='Truth')
401 404 plt.ylabel('Power (arbitrary units)')
402 405 plt.legend(loc='upper right')
403 406
404 407 # formatting y-axis
405 408 locs,labels = plt.yticks()
406 409 plt.yticks(locs, map(lambda x: "%.1f" % x, locs*1e4))
407 410 plt.text(0.0, 1.01, '1e-4', fontsize=10, transform = plt.gca().transAxes)
408 411
409 412
410 413 #---- MaxEnt----
411 414 plt.subplot(2, 1, 2)
412 415 plt.plot(180/math.pi*thetar, f_maxent, 'r', label='MaxEnt')
413 416 plt.plot(180/math.pi*thetat,fact, 'k--', label='Truth')
414 417 plt.ylabel('Power (arbitrary units)')
415 418 plt.legend(loc='upper right')
416 419
417 420 # formatting y-axis
418 421 locs,labels = plt.yticks()
419 422 plt.yticks(locs, map(lambda x: "%.1f" % x, locs*1e4))
420 423 plt.text(0.0, 1.01, '1e-4', fontsize=10, transform = plt.gca().transAxes)
421 424
422 425 plt.show()
423 426
424 427
425 428 # # PLOT PARA COMPRESSED SENSING
426 429 # #
427 430 # # subplot(3,1,3);
428 431 # # plot(180/pi*thetar,f_cs,'r-');
429 432 # # hold on;
430 433 # # plot(180/pi*thetat,fact,'k--');
431 434 # # hold off;
432 435 # # ylim([min(f_cs) 1.1*max(fact)]);
433 436 # # # title(sprintf('rel. RMSE: #.2e\tCorr: #.3f Corrc: #.3f', relrmse_cs, corr_cs, corrc_cs));
434 437 # # # title 'Compressed Sensing - Debauchies Wavelets'
435 438 # # xlabel 'Degrees'
436 439 # # ylabel({'Power';'(arbitrary units)'})
437 440 # # legend('Comp. Sens.','Truth');
438 441 # #
439 442 # # # set(gcf,'Position',[749 143 528 881]); # CSL
440 443 # # # set(gcf,'Position',[885 -21 528 673]); # macbook
441 444 # # pause(0.01);
442 445
443 446
444 447 # # Store Results
445 448 corr[0, snri, Ni] = corr_fourier;
446 449 corr[1, snri, Ni] = corr_capon;
447 450 corr[2, snri, Ni] = corr_maxent;
448 451 #corr[3, snri, Ni] = corr_cs;
449 452
450 453 rmse[0,snri,Ni] = relrmse_fourier;
451 454 rmse[1,snri,Ni] = relrmse_capon;
452 455 rmse[2,snri,Ni] = relrmse_maxent;
453 456 #rmse[3,snri,Ni] = relrmse_cs;
454 457
455 458 corrc[0,snri,Ni] = corrc_fourier;
456 459 corrc[1,snri,Ni] = corrc_capon;
457 460 corrc[2,snri,Ni] = corrc_maxent;
458 461 #corrc[3,snri,Ni] = corrc_cs;
459 462
460 463
461 464 print 'Capon:\t', elapsed_time_capon, 'sec';
462 465 print 'Maxent:\t',elapsed_time_maxent, 'sec';
463 466 #print 'CS:\t%3.3f sec\n',elapsed_time_cs;
464 467
465 468 print (NN*(snri+1) + Ni), '/', (SNRdBvec.size*NN);
466 469
467 470 print corr
468 471 No newline at end of file
Show file before | Show file after | Hide comments
@@ -1,59 +1,72
1 1 '''
2 2 Created on May 29, 2014
3 3
4 4 @author: Yolian Amaro
5 5 '''
6 6
7 7 #from sp import multirate
8 import cshift
8 from cshift import *
9 9 from multirate import upfirdn
10 10
11 11 def afb(x, af):
12 12
13 13 # Analysis filter bank
14 14 #
15 15 # USAGE:
16 16 # [lo, hi] = afb(x, af)
17 17 # INPUT:
18 18 # x - N-point vector, where
19 19 # 1) N is even
20 20 # 2) N >= length(af)
21 21 # af - analysis filters
22 22 # af(:, 1) - lowpass filter (even length)
23 23 # af(:, 2) - highpass filter (even length)
24 24 # OUTPUT:
25 25 # lo - Low frequecy output
26 26 # hi - High frequency output
27 27 # EXAMPLE:
28 28 # [af, sf] = farras;
29 29 # x = rand(1,64);
30 30 # [lo, hi] = afb(x, af);
31 31 # y = sfb(lo, hi, sf);
32 32 # err = x - y;
33 33 # max(abs(err))
34 34 #
35 35 # WAVELET SOFTWARE AT POLYTECHNIC UNIVERSITY, BROOKLYN, NY
36 36 # http://taco.poly.edu/WaveletSoftware/
37 37
38 38 N = x.size;
39 L = (af).size/2;
40 x = cshift(x,-L);
39 L = (af).size/4; #L should be = 5
40 #print af
41 #print 'L', L
42 x = cshift(x,-(L-1));
43
44 # print 'afb x', x.shape
45 # print 'af[:,0]',af[:,0].shape
46 # print 'af[:,1]',af[:,1].shape
47 # print '-----------------------'
41 48
42 49 # lowpass filter
43 50 lo = upfirdn(x, af[:,0], 1, 2);
44 51
52
45 53 # VERIFY THIS!!!!!!!!!!!!
46 54 for i in range(0, L):
47 lo[i] = lo[N/2+[i]] + lo[i];
55 lo[i] = lo[N/2+i] + lo[i];
48 56
49 lo = lo[1:N/2];
57 lo = lo[0:N/2];
58
50 59
51 60 # highpass filter
52 hi = upfirdn(x, af[:,2], 1, 2);
61 hi = upfirdn(x, af[:,1], 1, 2);
53 62
54 63 for j in range(0, L):
55 hi[j] = hi(N/2+[j]) + hi[j];
64 hi[j] = hi[N/2+j] + hi[j];
56 65
57 hi = hi[1:N/2];
66 hi = hi[0:N/2];
67
68 # Reshape from 1D to 2D
69 lo = lo.reshape(1, lo.size)
70 hi = hi.reshape(1, hi.size)
58 71
59 72 return lo, hi
Show file before | Show file after | Hide comments
@@ -1,29 +1,34
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 # m - amount of shift
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 N = x.size;
25 n = np.arange(N-1);
26 n = np.arange(N);
26 27 n = np.mod(n-m, N);
27 y = x[n];
28
29 print x.shape
30
31 y = x[0,n];
32
28 33
29 34 return y No newline at end of file
Show file before | Show file after | Hide comments
@@ -1,38 +1,41
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 import FSfarras
9 import dualfilt1
8 from FSfarras import *
9 from dualfilt1 import *
10 from dualtree import *
11 from idualtree import *
10 12
11 13 def deb4_basis(N):
12 14
13 15 Psi = np.zeros(shape=(N,2*N+1));
14 16 idx = 1;
15 17
16 18 J = 4;
17 [Faf, Fsf] = FSfarras;
18 [af, sf] = dualfilt1;
19 [Faf, Fsf] = FSfarras();
20 [af, sf] = dualfilt1();
19 21
20 22 # compute transform of zero vector
21 23 x = np.zeros(shape=(1,N));
22 #w = dualtree(x, J, Faf, af);
23 # #
24 # # # Uses both real and imaginary wavelets
25 # # for i in range (1, J+1):
26 # # for j in range (1, 2):
27 # # for k in range (1, (w[i][j]).size):
28 # # w[i][j](k) = 1;
29 # # y = idualtree(w, J, Fsf, sf);
30 # # w[i][j](k) = 0;
31 # # # store it
32 # # Psi(:,idx) = y.T.conj();
33 # # idx = idx + 1;
34 # #
35 # # # Add uniform vector (seems to be useful if there's a background
36 # # Psi(:,2*N+1) = 1/np.sqrt(N);
37 #
38 # return Psi No newline at end of file
24 w = dualtree(x, J, Faf, af);
25
26
27 # Uses both real and imaginary wavelets
28 for i in range (0, J):
29 for j in range (0, 1):
30 for k in range (0, (w[i][j]).size):
31 w[i][j][0,k] = 1;
32 y = idualtree(w, J, Fsf, sf);
33 w[i][j][0,k] = 0;
34 # store it
35 Psi[:,idx] = y.T.conj();
36 idx = idx + 1;
37
38 # Add uniform vector (seems to be useful if there's a background
39 Psi[:,2*N+1] = 1/np.sqrt(N);
40
41 return Psi No newline at end of file
Show file before | Show file after | Hide comments
@@ -1,63 +1,76
1 1 '''
2 2 Created on May 29, 2014
3 3
4 4 @author: Yolian Amaro
5 5 '''
6 6
7 7 import numpy as np
8 import afb
8 from afb import *
9 9
10 10 def dualtree(x, J, Faf, af):
11 11
12 12 # Dual-tree Complex Discrete Wavelet Transform
13 13 #
14 14 # USAGE:
15 15 # w = dualtree(x, J, Faf, af)
16 16 # INPUT:
17 17 # x - N-point vector
18 18 # 1) N is divisible by 2^J
19 19 # 2) N >= 2^(J-1)*length(af)
20 20 # J - number of stages
21 21 # Faf - filters for the first stage
22 22 # af - filters for the remaining stages
23 23 # OUTPUT:
24 24 # w - DWT coefficients
25 25 # w{j}{1}, j = 1..J - real part
26 26 # w{j}{2}, j = 1..J - imaginary part
27 27 # w{J+1}{d} - lowpass coefficients, d = 1,2
28 28 # EXAMPLE:
29 29 # x = rand(1, 512);
30 30 # J = 4;
31 31 # [Faf, Fsf] = FSfarras;
32 32 # [af, sf] = dualfilt1;
33 33 # w = dualtree(x, J, Faf, af);
34 34 # y = idualtree(w, J, Fsf, sf);
35 35 # err = x - y;
36 36 # max(abs(err))
37 37 #
38 38 # WAVELET SOFTWARE AT POLYTECHNIC UNIVERSITY, BROOKLYN, NY
39 39 # http://taco.poly.edu/WaveletSoftware/
40 40
41 # ---------Trees Structure---------------#
42 # w [ 0 1 2 .... J ] #
43 # | | | | #
44 # [0 1] [0 1] [0 1] [0 1] #
45 #----------------------------------------#
46
41 47 # normalization
42 48 x = x/np.sqrt(2);
49
43 50
44 w = np.zeros(shape=(J,2)) ### VERIFY THIS DEFINITION
51 w = np.zeros(shape=(J+1), dtype=object)
45 52
53 for j in range (0, w.size):
54 w[j] = np.zeros(shape=(J+1), dtype=object)
55
46 56 # Tree 1
47 [x1,w[1,0]] = afb(x, Faf[0,1]); # w{1}{1}
57 [x1, w[0][0]] = afb(x, Faf[0,0]); # w{1}{1}
48 58
49 for j in range (2,J):
50 [x1,w[j,0]] = afb(x1, af[0,1]); #check this
51
52 w[J+1,1] = x1;
53 59
60 for j in range (1,J):
61 [x1,w[j][0]] = afb(x1, af[0,0]); ### or 0,1????
62
63
64
65 w[J][0] = x1;
66
54 67 # Tree 2
55 [x2,w[1,2]] = afb(x, Faf[0,1]);
56
57 for j in range (2,J):
58 [x2,w[j,2]] = afb(x2, af[0,1]);
59
60 w[J+1,2] = x2;
68 [x2,w[0][1]] = afb(x, Faf[0,1]);
69
70 for j in range (1,J):
71 [x2,w[j][1]] = afb(x2, af[0,1]);
72
73 w[J][1] = x2;
61 74
62 75 return w
63 76
General Comments 0
You need to be logged in to leave comments. Login now