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