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